Opened 9 years ago

Closed 9 years ago

# prime_powers doesn't work with start very well

Reported by: Owned by: kcrisman was major sage-5.6 number theory beginner was sage-5.6.beta1 Kevin Halasz Dmitrii Pasechnik, Punarbasu Purkayastha, Karl-Dieter Crisman 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_primepowers.2.patch.

### comment:2 Changed 9 years ago by dimpase

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 9 years ago by khalasz

• Description modified (diff)
• Status changed from new to needs_review

### comment:4 follow-up: ↓ 5 Changed 9 years ago by khalasz

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 9 years ago by kcrisman

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 9 years ago by was

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 9 years ago by khalasz

• Authors set to Kevin Halasz

### comment:8 Changed 9 years ago by khalasz

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

### comment:9 Changed 9 years ago by kcrisman

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 9 years ago by khalasz

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 9 years ago by dimpase

• Status changed from needs_review to 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 9 years ago by khalasz

• Status changed from needs_work to 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 9 years ago by dimpase

• Description modified (diff)

### comment:14 Changed 9 years ago by dimpase

• Description modified (diff)

### comment:15 Changed 9 years ago by dimpase

• Status changed from needs_review to 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 9 years ago by ppurka

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 9 years ago by 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.

### comment:18 in reply to: ↑ 17 Changed 9 years ago by dimpase

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 9 years ago by ppurka

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 9 years ago by khalasz

• Status changed from needs_work to needs_review

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

### comment:21 in reply to: ↑ 20 Changed 9 years ago by dimpase

• Status changed from needs_review to 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 9 years ago by khalasz

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 9 years ago by khalasz

• Status changed from needs_work to needs_review

### comment:24 Changed 9 years ago by ppurka

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 9 years ago by ppurka

• Description modified (diff)
• Reviewers set to Punarbasu Purkayastha

1. Fixed trailing whitespaces.
2. Made sure `prime_powers(-1, positive integer)` works.
3. 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_reviewer.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.

Patchbot: apply 13516_reviewer.2.patch

Version 0, edited 9 years ago by ppurka (next)

### comment:27 Changed 9 years ago by ppurka

• Reviewers changed from Punarbasu Purkayastha to Dmitrii Pasechnik, Punarbasu Purkayastha

### comment:28 Changed 9 years ago by ppurka

• Description modified (diff)

Argggh! Uploaded the wrong patch earlier X(

Patchbot: apply 13516_primepowers.2.patch

### comment:29 Changed 9 years ago by khalasz

• Status changed from needs_review to positive_review

Looks good to me. Thanks for your help/review!

### comment:30 Changed 9 years ago by kcrisman

Nice work. Below, things almost certainly not worth addressing now that ppurka and khalasz have worked so hard to get this in, but I'll still point them out if ppurka is really bored and wants to fix them:

• Typo
```# inserted to prevent an error from occuring
```
• Might be worth using Python 3 string formatting in the error messages
• Conceivable to replace the reference to #922 with a nice new `:trac:` thingie

### comment:31 Changed 9 years ago by ppurka

@kcrisman - fixed. Thanks. :)

### comment:32 Changed 9 years ago by kcrisman

```sage: type(v[0])      # :trac:`922`
```

and have it formatted, because it's in a literal block, it would have to be

```TESTS:

Check :trac:`922`is fixed::

sage: v = prime_powers(10)
sage: type(v[0])
<type 'sage.rings.integer.Integer'>
```

or something like that. But I don't know if that's worth it...

### Changed 9 years ago by ppurka

apply only this to devel/sage

### comment:33 follow-up: ↓ 34 Changed 9 years ago by ppurka

Yes, you are right. I have reverted the last change.

### comment:34 in reply to: ↑ 33 Changed 9 years ago by kcrisman

Yes, you are right. I have reverted the last change.

Sorry to be so picky, anyway nice work by both of you!

### comment:35 Changed 9 years ago by ppurka

• Reviewers changed from Dmitrii Pasechnik, Punarbasu Purkayastha to Dmitrii Pasechnik, Punarbasu Purkayastha, Karl-Dieter Crisman

Sorry to be so picky, anyway nice work by both of you!

No problem. But you do deserve to be in the reviewers list! ;)

### comment:36 Changed 9 years ago by jdemeyer

• Merged in set to sage-5.6.beta1
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.