Fix is_euler_pseudoprime

[haru-44]

There are several things called Euler pseudoprime. However, the current implementation is different from the docstring description and wikipedia. So we modified the implementation to be in line with the docstring definition.

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

• ntheory
  • The behavior of is_euler_pseudoprime has changed. Composite numbers that used to return False may now also return True.

[sympy-bot]

✅

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• ntheory
  • The behavior of is_euler_pseudoprime has changed. Composite numbers that used to return False may now also return True. (#25629 by @haru-44 and @smichr)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13.

Click here to see the pull request description that was parsed.

There are several things called Euler pseudoprime. However, the current implementation is different from the docstring description and wikipedia. So we modified the implementation to be in line with the docstring definition.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* ntheory
  * The behavior of `is_euler_pseudoprime` has changed. Composite numbers that used to return `False` may now also return `True`.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

| Change   | Before [8059df73] <sympy-1.12^0>   | After [fe959796]    |   Ratio | Benchmark (Parameter)                                                |
|----------|------------------------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 102±1ms                            | 67.1±0.4ms          |    0.66 | integrate.TimeIntegrationRisch02.time_doit(10)                       |
| -        | 99.5±1ms                           | 65.9±0.4ms          |    0.66 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 26.0±0.09μs                        | 46.3±0.2μs          |    1.78 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 8.29±0.03ms                        | 4.49±0.01ms         |    0.54 | logic.LogicSuite.time_load_file                                      |
| -        | 2.46±0ms                           | 786±3μs             |    0.32 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 12.0±0.02ms                        | 2.34±0ms            |    0.2  | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 436±3μs                            | 97.1±0.5μs          |    0.22 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 5.58±0.01ms                        | 431±0.8μs           |    0.08 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 12.2±0.03ms                        | 1.30±0ms            |    0.11 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 7.34±0.02ms                        | 4.70±0.02ms         |    0.64 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 31.9±0.06ms                        | 14.3±0.04ms         |    0.45 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 8.29±0.03ms                        | 1.40±0.01ms         |    0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 19.2±0.07ms                        | 11.0±0.05ms         |    0.57 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 245±0.7ms                          | 84.7±0.2ms          |    0.35 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 7.38±0.03ms                        | 635±2μs             |    0.09 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 31.7±0.07ms                        | 1.02±0ms            |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 719±2μs                            | 234±2μs             |    0.33 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 7.49±0.02ms                        | 237±3μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 19.5±0.05ms                        | 244±0.9μs           |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 192±0.5μs                          | 107±0.5μs           |    0.56 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 365±1μs                            | 130±0.8μs           |    0.36 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 37.7±0.1ms                         | 16.2±0.05ms         |    0.43 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[smichr]

Any idea what version of definition the original code was giving?

[smichr]

Wouldn't it be better to exclude primes from return True? Although the docstring said otherwise, it seems the tests were structured in that way, showing the first Euler pseudoprime for given bases, not the first number that would pass the test.

In retrospect it might have been better to name this something like euler_test which will pass for primes or some composites relatively prime to a given base. You could allow the base to be 1 and in that case return n!=2 and not isprime(n).

[haru-44]

Any idea what version of definition the original code was giving?

Definition of Euler pseudoprime

1. $a^{\frac{n-1}{2}} \overset{?}{\equiv} \pm 1 \pmod{n}$ : wikipedia, oeis
2. $a^{\frac{n-1}{2}} \overset{?}{\equiv} \left( \frac{a}{n} \right) \pmod{n}$ : gmpy2 (Oh, the gmpy2 documentation is wrong!)

However, the current implementation does not fall into either of these cases. The current implementation runs the Miller-Rabin primality test twice, the first time calling mr and the second time implementing it on its own. At least the second one is unnecessary.

Wouldn't it be better to exclude primes from return True?

Is the question whether to include or exclude prime numbers in the definition of pseudoprime? As noted on wikipedia, both definitions are possible.

Some sources use the term pseudoprime to describe all probable primes, both composite numbers and actual primes.

When considering composite numbers that are mistakenly determined to be prime, it is rational not to include prime numbers in pseudoprime.
On the other hand, in programming for prime number determination, it is more natural to include primes in pseudoprime (e.g., gmpy2.is_**_prp).
So I think the name as it is now will not cause any misunderstanding, and the explanation in the docstring will make it clearer.

[smichr]

What do you think about handling the case

if a == 1:
    return n % 2 and not isprime(n)  # odd composite

[haru-44]

When considering is_euler_pseudoprime as a probabilistic algorithm for determining prime numbers, the use of isprime internally raises the classic dilemma of the chicken and the egg. Consequently, I do not advocate the utilization of isprime.

If a==1 should not be an error, then it is True when n is odd. So, I think this can be translated into the question "Should a==1 be treated as an error?" There is no reason to actively make it an error, so I will modify it so that it should not be an error.

[haru-44]

thanks