Opened 10 years ago

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

Reported by: Owned by: Karl-Dieter Crisman William Stein major sage-5.6 number theory beginner William Stein Kevin Halasz N/A

```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. http://trac.sagemath.org/sage_trac/raw-attachment/ticket/13516/13516_primepowers.patch

### 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) new → needs_review

### comment:4 follow-up:  5 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_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 Kevin Halasz

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

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