Opened 10 years ago
Last modified 10 years ago
#13516 closed defect
prime_powers doesn't work with start very well — at Version 26
Reported by: | Karl-Dieter Crisman | Owned by: | William Stein |
---|---|---|---|
Priority: | major | Milestone: | sage-5.6 |
Component: | number theory | Keywords: | beginner |
Cc: | William Stein | Merged in: | |
Authors: | Kevin Halasz | Reviewers: | Punarbasu Purkayastha |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
See this sage-support thread.
In Sage 5.3, the function prime_powers behaves a little strange: sage: prime_powers(4,10) [4, 5, 7, 8, 9] # As expected sage: prime_powers(5,10) [7, 8, 9] # 5 isn't a prime power anymore??? # And now things become even worse: sage: prime_powers(7,10) --------------------------------------------------------------------------- IndexError Traceback (most recent call last) /home/mueller/<ipython console> in <module>() /home/mueller/local/sage-5.3/local/lib/python2.7/site-packages/sage/rings/arith.pyc in prime_powers(start, stop) 743 i = bisect(v, start) 744 if start > 2: --> 745 if v[i] == start: 746 i -= 1 747 w = list(v[i:]) IndexError: list index out of range
Yeah, this seems problematic. The code in question is old, too, so perhaps there is a more efficient way to do it in the meantime...
Apply to devel/sage
: 13516_reviewer.2.patch.
Change History (26)
comment:1 Changed 10 years ago by
Keywords: | beginner added |
---|
comment:2 Changed 10 years ago by
comment:3 Changed 10 years ago by
Description: | modified (diff) |
---|---|
Status: | new → needs_review |
comment:4 follow-up: 5 Changed 10 years ago by
I found the code to be riddled with errors, so I decided to completely rework it. Also, I have qualms about calling 1 a prime power, but did so because the old function did. If you think its fine to drop this, let me know.
comment:5 Changed 10 years ago by
Thanks for your work - hopefully someone will review it soon. You can put your real name in the author area.
Also, I have qualms about calling 1 a prime power, but did so because the old function did. If you think its fine to drop this, let me know.
Well, John Horton Conway does call -1 a prime, in which case every nonzero integer (not just positive) is a unique product of prime powers - not a unique product of primes, note, nor of the exponents, but of the prime powers themselves (I can't find a link for this right now, my apologies) in which case positives get the power 1 and and negatives -1. I think that's right... anyway, maybe they were thinking this?
comment:6 Changed 10 years ago by
I would prefer to leave 1 as a prime power because it is listed in Sloane's tables as a prime power: http://oeis.org/A000961
There he says "Since 1 = p^{0 does not have a well defined prime base p, it is sometimes not regarded as a prime power.", which might be where your misgivings come from. }
If by "prime power" one thinks of "power of a prime", the only question is in what set are we considering the prime powers. If we take the natural numbers, then the number 1 is definitely a power of a prime.
If by "prime power" one thinks "power of a specific canonical prime", then 1 is not such a thing.
In this case, the best thing to do is stick with what is there (to avoid introducing bugs in other people's code!) and clearly document/define what a prime power is in Sage.
comment:7 Changed 10 years ago by
Authors: | → Kevin Halasz |
---|
comment:8 Changed 10 years ago by
I just updated the docstring to make the fact that 1 is a prime power explicit.
comment:9 Changed 10 years ago by
Could you speed this up slightly by making s = stop.sqrt()
or something so that it's not computed for each prime. In fact, even that is a more expensive comparison each time because stop.sqrt()
is likely a symbolic element, so maybe even stop.sqrt().n()
would be appropriate... Also, once p >stop.sqrt()
, presumably all remaining p
are beyond it as well, so maybe there could be some speedup there too. Just some ideas.
comment:10 Changed 10 years ago by
I've changed the patch so that s=stop.sqrt() is calculated outside of the for loop. After some tests, I saw that this was faster than setting s=stop.sqrt().n().
Also, note that when p>s, the content of that if loop is a break command, meaning that the entire for loop ends. Thus, once a single p>s, no more p values are tried.
comment:11 Changed 10 years ago by
Status: | needs_review → needs_work |
---|
The comment on line 708 in sage/rings/arith.py
needs to be fixed, too - it talks about primes rather than prime powers.
I also think that the following:
sage: prime_powers(10,7) 761 Traceback (most recent call last): 762 ... 763 ValueError: the first input must be less than the second input, however, 10 > 7
i.e. the corresponding implementation logic is not right, in the sense that it should just return empty lists rather than throwing exceptions. And negative start
should be allowed too (cf. the semantics of range()
).
Then, in the following fragment
783 output = prime_range(start,stop) 784 if start == 1: 785 output.append(1) 786 787 s = stop.sqrt() 788 for p in prime_range(stop):
prime_range()
, which is not cheap, is basically called two times instead of one.
One can do with one call to prime_range(stop)
just fine.
comment:12 Changed 10 years ago by
Status: | needs_work → needs_review |
---|
Dimpase, I've addressed all of your suggestions. Let me know what you think of these changes/if you have other suggestions.
comment:13 Changed 10 years ago by
Description: | modified (diff) |
---|
comment:14 Changed 10 years ago by
Description: | modified (diff) |
---|
comment:15 Changed 10 years ago by
Status: | needs_review → needs_work |
---|
You always should coerce stop
into Integer. Indeed, currently one gets:
sage: prime_powers(1,int(9)) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /usr/local/src/sage/sage-5.4.rc0/devel/sage-main/<ipython console> in <module>() /usr/local/src/sage/sage-5.4.rc0/local/lib/python2.7/site-packages/sage/rings/arith.pyc in prime_powers(start, stop) 761 from sage.rings.integer import is_Integer 762 if not (is_Integer(start) and (is_Integer(stop) or stop == None)): --> 763 raise TypeError, "both inputs must be integers, but your inputs were %s and %s"%(start,stop) 764 765 # deal with the case in which only one input is given TypeError: both inputs must be integers, but your inputs were 1 and 9
This is because taking sqrt(int(9))
does not work too well in Sage...
Second issue: the comment # check that all inputs are positive integers
on line 760 is misleading!
comment:16 Changed 10 years ago by
In the documentation please write `start`, `stop`
with double backticks. The single backticks will make them be formatted as latex, which is not what is desired.
comment:17 follow-up: 18 Changed 10 years ago by
Dimpase,
I'm not sure if I understand what you're suggesting I do. Should I replace the check that the inputs are Integers with a coercion of the inputs into Integers? Once I coerce the elements the check is redundant, as either it raises an error in and of itself (say, if it was passed a string it will raise a TypeError?), or will fix the problem.
comment:18 Changed 10 years ago by
Replying to khalasz:
Dimpase,
I'm not sure if I understand what you're suggesting I do. Should I replace the check that the inputs are Integers with a coercion of the inputs into Integers? Once I coerce the elements the check is redundant, as either it raises an error in and of itself (say, if it was passed a string it will raise a TypeError?), or will fix the problem.
I think I would prefer the input of type int
to be coerced into Integer
, and throw an error if it's neither int
nor Integer
. The reason is that one could potentially try to apply this function to different from ZZ
rings, with strange results, if one just blindly coerces stuff into Integer
.
comment:19 Changed 10 years ago by
I think you can write it like this
from sage.rings.integer import Integer if not isinstance(start, (int, Integer)): raise TypeError("start must be an integer") if stop is not None and not isinstance(stop, (int, Integer)): raise TypeError("stop must be an integer")
comment:20 follow-up: 21 Changed 10 years ago by
Status: | needs_work → needs_review |
---|
I've made the changes. Let me know what you think!
comment:21 Changed 10 years ago by
Status: | needs_review → needs_work |
---|
Replying to khalasz:
I've made the changes. Let me know what you think!
I think there is a bug in the code, coming from the fact that
sage: Integer(None)==None False
You also should have test cases (doctests) for all the different combinations of start/stop, and test them, too. You know that you can run Sage so that it tests doctests in a particular file, right?
E.g. the bug above would have been caught by the proper doctest.
Also, for some reason you removed the test sage: v = prime_powers(10)
, but it was there for good reason.
comment:22 Changed 10 years ago by
I forgot to rebuild sage before doctesting before posting this patch. Sorry for putting it up with such a stupid mistake. I've fixed it, and added back the test sage: v = prime_powers(1)
.
I think I've covered all the possible basic scenarios in the doctests, in both the EXAMPLES and the TESTS, do you disagree?
comment:23 Changed 10 years ago by
Status: | needs_work → needs_review |
---|
comment:24 Changed 10 years ago by
Some more nitpicks. :)
- Don't use
==
or!=
when comparing againstNone
. See PEP 8. - When you describe the arguments
start
andstop
, then describe the default value too.- ``start`` - an integer. If two inputs are given, .... - ``stop`` - an integer (default: ``None``). An upper bound for ....
- Please don't remove the trac numbers from the examples. They were put there after someone fixed some bug.
sage: type(v[0]) # trac #922
- Can you write the
TypeError
in python 3 style? Every small bit will help in the migration to python 3 later.raise TypeError("start must be an integer, %s is not an integer"%start) raise TypeError("stop must be an integer, %s is not an integer"%stop)
comment:26 Changed 10 years ago by
Description: | modified (diff) |
---|---|
Reviewers: | → Punarbasu Purkayastha |
Thanks a lot for addressing my concerns. I have made some changes to your patch.
- Fixed trailing whitespaces.
- Made sure
prime_powers(-1, positive integer)
works. - Fixed
TypeError
.
The changes can be seen in 13516_reviewer.patch. All these changes have been merged with your patch and the new patch is now 13516_primepowers.2.patch.
Aside from the above corrections, the changes introduced by your patch has positive review from my side. If you think my changes are ok, feel free to change the ticket to positive review.
as
prime_powers(m)
works, an easy workaround is to re-implementprime_powers(m,n)
as the difference ofprime_powers(n)
andprime_powers(m-1)
.