Opened 10 years ago

Last modified 10 years ago

#13516 closed defect

prime_powers doesn't work with start very well — at Version 14

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:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by Dima Pasechnik)

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

  1. patch

Change History (14)

comment:1 Changed 10 years ago by Karl-Dieter Crisman

Keywords: beginner added

comment:2 Changed 10 years ago by Dima Pasechnik

as prime_powers(m) works, an easy workaround is to re-implement prime_powers(m,n) as the difference of prime_powers(n) and prime_powers(m-1).

comment:3 Changed 10 years ago by Kevin Halasz

Description: modified (diff)
Status: newneeds_review

comment:4 Changed 10 years ago by Kevin Halasz

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 in reply to:  4 Changed 10 years ago by Karl-Dieter Crisman

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 William Stein

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 = p0 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 Kevin Halasz

Authors: Kevin Halasz

comment:8 Changed 10 years ago by Kevin Halasz

I just updated the docstring to make the fact that 1 is a prime power explicit.

comment:9 Changed 10 years ago by Karl-Dieter Crisman

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 Kevin Halasz

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 Dima Pasechnik

Status: needs_reviewneeds_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 Kevin Halasz

Status: needs_workneeds_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 Dima Pasechnik

Description: modified (diff)

comment:14 Changed 10 years ago by Dima Pasechnik

Description: modified (diff)
Note: See TracTickets for help on using tickets.