Opened 10 years ago

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

Reported by: Owned by: Karl-Dieter Crisman William Stein major sage-5.6 number theory beginner William Stein Kevin Halasz Punarbasu Purkayastha 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 to `devel/sage`: 13516_reviewer.2.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)

### comment:14 Changed 10 years ago by Dima Pasechnik

Description: modified (diff)

### comment:15 Changed 10 years ago by Dima Pasechnik

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 Punarbasu Purkayastha

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

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 in reply to:  17 Changed 10 years ago by Dima Pasechnik

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 Punarbasu Purkayastha

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

Status: needs_work → needs_review

I've made the changes. Let me know what you think!

### comment:21 in reply to:  20 Changed 10 years ago by Dima Pasechnik

Status: needs_review → needs_work

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

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

Status: needs_work → needs_review

### comment:24 Changed 10 years ago by Punarbasu Purkayastha

Some more nitpicks. :)

1. Don't use `==` or `!=` when comparing against `None`. See PEP 8.
2. When you describe the arguments `start` and `stop`, then describe the default value too.
```    - ``start`` - an integer. If two inputs are given, ....
- ``stop`` - an integer (default: ``None``). An upper bound for ....
```
3. 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
```
4. 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 Punarbasu Purkayastha

Description: modified (diff) → Punarbasu Purkayastha

2. Made sure `prime_powers(-1, positive integer)` works.
3. Fixed `TypeError`.