fixing a bug in perfect_power

[fchapoton]

Fix #40846

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.

[vincentmacri]

You beat me to it!

Looks good, just add a sentence before the new test saying Test small powers of even numbers that are not a power of 2 (see :issue:`40846`).

[vincentmacri]

You beat me to it!

Looks good, just add a sentence before the new test saying Test small powers of even numbers that are not a power of 2 (see :issue:`40846`).

Not too sure about how the docstring formatting works, this might require changing TESTS:: to TESTS: and adding a description to the previous tests which test powers of two.

[github-actions]

Documentation preview for this PR (built with commit 5649932; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[mantepse]

wouldn't it be better to test against len(_small_primes_table) (or store this in a constant)?

[vincentmacri]

wouldn't it be better to test against len(_small_primes_table) (or store this in a constant)?

_small_primes_table is kind of weird. It's not a list or table of primes themselves. Rather, _small_primes_table[n] is 1 if $2n + 1$ is prime, and 0 if $2n + 1$ is not prime (so _small_primes_table[0] = 0 because 1 is not prime, and _small_primes_table[1] = 1 because 3 is prime).

[whoami730]

This and similar logic is being employed in three different places in the code, ie, check if value < 1000 and is odd and if so, invoke _small_primes_table[value >> 1]. We should store this 1000 value in a constant.

[mantepse]

I think that 36.perfect_power() should also be a test. I don't see right now why the PR would fix it.

[vincentmacri]

I think that 36.perfect_power() should also be a test.

Sure.

I don't see right now why the PR would fix it.

With the PR 36 will fail all the checks for the faster version and reach the pari fallback (which works). The bug was caused by 36 reaching the case for small odd primes.

if not (n & 1): checks if n is even. If so, if mpz_popcount(self.value) == 1: checks if n is a power of 2.

Previously 36 passed the first check and failed the second, then it checked if it was less than 1000 and if so whether $2 \lfloor {\frac{n}{2}} \rfloor + 1$ prime. If $n$ is odd then this is expression is equal to n. The return statement operates under this assumption. But if you plug in $n = 36$ you get $37$ which is prime. It then thinks that n is prime and the way you write a prime as a power is just itself to the power of 1 hence the bug.

[whoami730]

Looks good, just add a sentence before the new test saying Test small powers of even numbers that are not a power of 2 (see :issue:`40846`).
+1.

[fchapoton]

comment added. Other changes left for other PRs.

[whoami730]

Other changes left for other PRs.

Let's add a TODO then to indicate what changes need to be made?

[fchapoton]

we should really quickly set this one to "positive review" ; the rest is not urgent

[mantepse]

ok for me

[vincentmacri]

Thank you for the quick fix!

[vincentmacri]

Not needed for this PR, but since we were talking about it, some food for thought on _small_primes_table if someone wants to work on cleaning it up:

• Why is _small_primes_table a list of int? Surely a tuple of bool would be more memory-efficient. (Or better yet bits, but I think working with bits is awkward unless there's some type for it that I'm not aware of.) And if it's more memory efficient then we could probably make it bigger (go up to 10000 instead of 1000 for example).
• Does Cython support a way to define a constant/macro that is computed at compile-time? _small_primes_table would be more readable if it were defined as _small_primes_table = tuple(is_prime(2 * n + 1) for n in range(500)), but this module is supposed to be efficient so we should only do that if it can be computed at compile-time.