# Primality test in Python

Here is a working Python implementation of primality test. Is there something that I could change in code to achieve a better running time? For more information about algorithm see this post.

from sympy import *
from sympy.abc import x

n=int(input("Enter a number : "))

def xmat(r,n):
return Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)

def smallestr(n):
if n==1 or n%2==0:
return 0
else:
r=3
while r<1000000:
u=n%r
if u==0 and r<n:
return 0
elif not(u==0) and not(u==1) and not(u==r-1):
return r
else:
r=nextprime(r)

def myisprime(n):
r=smallestr(n)
if r==0:
return n==2
else:
xp=(xmat(r,n)**n)*Matrix([[x],[1]])
return trunc(xp[1],n)==(rem(x*(1%n),x**r-1))**n

if myisprime(n):
print("prime")
else:
print("composite")


You can run this code here.

EDIT

To clarify things, this Python code represents my attempt to translate the following PARI/GP code which is very fast.

Ooooh boy, a prime number finder with the question, "how do I make this faster." It's like asking a bartender what their favorite drink is: it really depends. However, before we get to performance, let's tackle some of the stylistic considerations in this code.

1. The algorithm you're implementing is fairly complicated. Add some brief comments and docstrings to explain what's going on.

3. Name your constants. Why r < 1000000? This would be clearer if we had something like MAX_PRIME_VALUE = 1000000, r < MAX_PRIME_VALUE.

Now in terms of performance. Are there ways to make this faster? Yes, of course. If you had a dictionary of all prime numbers <= 1000000, then you could do a simple lookup and be done. The ultimate trade-off between runtime and storage space is to simply pre-compute all the answers to your question. Also, it depends on what hardware/software you have available (pure no-library python? Only sympy? Numpy? Numba? Cython? CUDA?). So in order to answer this question, you really need to ask what your use case is and what trade-offs you're able to make.

For the sake of argument, let's use your code above as our basic requirements. Let's say we have access to sympy and basic Python, but nothing else. As a baseline, let's consider a trivial prime number finder that use nothing but pure Python. We'll just count up from 3 up to sqrt(n): it's naive, it's dead simple to write, and it's actually reasonably fast just because Python is a decent language (and it even has an O(sqrt(N)) runtime, which is not a bad place to start):

def naive_is_prime(n):
if n == 2:
return True
if n < 2 or n % 2 == 0:
return False

# Consolidated code thanks to GZ0's comment
return all(n % i != 0 for i in range(3, int(n ** 0.5 + 1), 2))


And a couple benchmarks:

%timeit myisprime(11)
1.23 ms ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit naive_is_prime(11)
1.11 µs ± 47.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)


Uh oh, that's not a good place to begin, with the naive code giving us a 1000x speedup. Well, maybe that's just because the numbers are so small; maybe the overhead of using sympy will shine when we start considering larger primes. Unfortunately, it slows down rapidly, and waiting for it to compute if tiny numbers are prime takes an unbearable amount of time:

%time myisprime(61)
Wall time: 28.3 s


whereas the naive implementation above continues to be almost instantaneous up until fairly large primes:

%timeit naive_is_prime(2**31 - 1)  # 2147483647
2.7 ms ± 160 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


I think the moral of the story is that "premature optimization is the root of all evil". I fall into this all the time, and even did so while answering this question (I had a more complicated baseline to show off that performed much worse than that trivial code), and it's very tempting to use a fancy package like sympy to accomplish what's not a very fancy problem. Start simple, then look into some common approaches to the problem (like the Sieve approach) and figure out if they actually help you or not.

• The overall review is great however there is a minor bug in the implemention: the upper bound should not be int(math.ceil(math.sqrt(n))), it should be int(n**0.5) + 1 (this does not require math import), or alternatively math.floor(math.sqrt(n)) + 1. Otherwise it fails on some square numbers such as 9. Also, the loop can be written in a one-liner using the built-in all function. – GZ0 Sep 11 '19 at 18:10
• For performance simply writing all(n % i for i ...) is slightly faster despite that it hurts readablity a bit. – GZ0 Sep 11 '19 at 18:41
• Please, see EDIT – Peđa Terzić Sep 12 '19 at 7:57
• @PeđaTerzić It's likely that the PARI/GP has a faster implementation than the nextprime you seem to be importing from sympy. It's possible that you could speed up your code if you create a generator for your primes there. But so far, this answer is very valid to your question, which was "how do I speed up this primality check" – Gloweye Sep 12 '19 at 8:37
• @Gloweye, nextprime doesn't even register in the profile output. – Peter Taylor Sep 12 '19 at 9:37

### Implementation

def xmat(r,n):
return Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)


Am I correct in thinking that n will always be a positive integer greater than 2? If so, *rem(1,x**r-1)*(1%n) can be optimised away entirely.

def myisprime(n):
r=smallestr(n)
if r==0:
return n==2
else:
xp=(xmat(r,n)**n)*Matrix([[x],[1]])
...


How to parse smallestr? Is it small_e_str? PEP8 promotes using underscores to separate words in names.

I inserted here a debug line:

        print(r, xp)


When n = 3 we get

5 Matrix([[-4*x**2 + x*(2*x*(4*x**2 - 1) - 2*x) + 1], [x*(4*x**2 - 1) - 2*x]])


It seems that sympy is manipulating expression trees rather than polynomials, and not simplifying at all ever. If you can't find a way to make it auto-simplify, it would probably be faster to roll your own polynomial class.

        return trunc(xp[1],n)==(rem(x*(1%n),x**r-1))**n


Simplifying:

        return trunc(xp[1],n)==x**n


Now, trunc implies that you should be working over $$\\mathbb{Z} / n\mathbb{Z}\$$ rather than $$\\mathbb{Z}\$$. xmat(r,n)**n is going to generate coefficients which are exponentially larger than n: working modulo n all the way through would be far far faster.

(I implemented a naïve version out of curiosity, and it tests 61 in about 2.4 milliseconds, vs 42 seconds for your code).

### Deep dive into the algorithm

By diagonalisation, if I haven't messed up the algebra,

$$\begin{pmatrix}2x & -1 \\ 1 & 0\end{pmatrix}^n \begin{pmatrix}x \\ 1\end{pmatrix} %= \left(S J S^{-1} \right)^n \begin{pmatrix}x \\ 1\end{pmatrix} %= S \left(J \right)^n S^{-1} \begin{pmatrix}x \\ 1\end{pmatrix} %= \begin{pmatrix} x-\sqrt{x^2-1} & x+\sqrt{x^2-1} \\ 1 & 1 \end{pmatrix} \begin{pmatrix} x-\sqrt{x^2-1} & 0 \\ 0 & x+\sqrt{x^2-1} \end{pmatrix}^n \begin{pmatrix} -\frac{1}{2\sqrt{x^2-1}} & \frac{x}{2\sqrt{x^2-1}}+\frac12 \\ \frac{1}{2\sqrt{x^2-1}} & \frac12 - \frac{x}{2\sqrt{x^2-1}} \end{pmatrix} \begin{pmatrix}x \\ 1\end{pmatrix} \\ %= \begin{pmatrix} x-\sqrt{x^2-1} & x+\sqrt{x^2-1} \\ 1 & 1 \end{pmatrix} \begin{pmatrix} (x-\sqrt{x^2-1})^n & 0 \\ 0 & (x+\sqrt{x^2-1})^n \end{pmatrix} \begin{pmatrix}\frac12 \\ \frac12\end{pmatrix} \\ %= \begin{pmatrix} x-\sqrt{x^2-1} & x+\sqrt{x^2-1} \\ 1 & 1 \end{pmatrix} \begin{pmatrix}\frac12 (x-\sqrt{x^2-1})^n \\ \frac12 (x+\sqrt{x^2-1})^n \end{pmatrix} \\ = \frac12 \begin{pmatrix} (x-\sqrt{x^2-1})^{n+1} + (x+\sqrt{x^2-1})^{n+1} \\ (x-\sqrt{x^2-1})^n + (x+\sqrt{x^2-1})^n \end{pmatrix}$$

So what this test boils down to is $$(x-\sqrt{x^2-1})^n + (x+\sqrt{x^2-1})^n \equiv 2x^n \pmod n$$

Expanding the LHS, $$(x-\sqrt{x^2-1})^n + (x+\sqrt{x^2-1})^n = \sum_{i=0}^n \binom{n}{i} x^{n-i} \left( (-\sqrt{x^2-1})^{i} + (\sqrt{x^2-1})^{i}\right) \\ = \sum_{i=0}^n \binom{n}{i} x^{n-i} (\sqrt{x^2-1})^{i} \left( (-1)^i + 1^i\right) \\ = \sum_{i=0}^{n/2} \binom{n}{2i} x^{n-2i} (\sqrt{x^2-1})^{2i} \\ = \sum_{i=0}^{n/2} \binom{n}{2i} x^{n-2i} (x^2-1)^{i} \\ = \sum_{i=0}^{n/2} \binom{n}{2i} x^{n-2i} \sum_{j=0}^i \binom{i}{j} (x^2)^{i-j} (-1)^j \\ = \sum_{j=0}^{n/2} (-1)^j x^{n-2j} \sum_{i=j}^{n/2} \binom{n}{2i} \binom{i}{j} \\$$

So the test is that $$\forall 0 < j \le \tfrac n2: \sum_{i=j}^{n/2} \binom{n}{2i} \binom{i}{j} \equiv 0 \pmod n$$

But a direct calculation on that basis is worse than simply verifying that $$\forall 0 < j \le \tfrac n2: \binom{n}{2j} \equiv 0 \pmod n$$ which is a perfectly valid primality test for odd $$\n\$$. And that itself is clearly worse than simply testing $$\forall 1 < j \le \sqrt{n}: n \not \equiv 0 \pmod j$$

So we're relying on the matrix exponentiation being extremely fast. Roughly speaking, we do $$\O(1)\$$ polynomial multiplications for polynomials of order $$\1, 2, 4, \ldots, n\$$, so if a multiplication of polynomials of length $$\\ell\$$ takes $$\\Theta(\ell^\alpha)\$$ then the overall exponentiation also takes $$\\Theta(n^\alpha)\$$. There's no way that $$\\alpha < \tfrac12\$$, so no matter how well optimised the implementation is it won't be faster than a naïve trial division.

• The number r seems to be completely irrelevant in the algorithm. I think it is quite likely that there is something wrong in the original implementation and r should take effect in some way (e.g. bound the polynominal lengths). – GZ0 Sep 12 '19 at 11:56
• Playing around the GP version of the code, it seems that the actual intention of *rem(1,x**r-1)*(1%n) is to create a symbolic representation of % (x**r - 1) and % n for each matrix element, which takes effect during the computation of matrix power. – GZ0 Sep 12 '19 at 12:16
• For more information see this post: mathoverflow.net/q/286304/88804 – Peđa Terzić Sep 12 '19 at 12:29
• @PeđaTerzić The link should be put in the main post (ideally, from the very beginning) so that others do not need to spend time on inferring the algorithm from the code. And you might want to explain that your intention is to implement that specific algorithm rather than needing a primality test program that works. – GZ0 Sep 12 '19 at 12:42
• @GZ0, I could edit your observations into my answer, but I'd prefer it if you posted an answer pointing out how the Python code fails to port the Pari/GP code because you should get the rep for your efforts, not me. – Peter Taylor Sep 12 '19 at 13:32

The provided implementation in the post is flawed and fails to port the Pari/GP code. In this statement

Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)


the actual intentions of rem(1, x**r-1) and (1%n) are to create symbolic representations of
% (x**r - 1) and % n for each matrix element, which is going to take effect during the computation of matrix power.

However, under this implementation, rem(1, x**r-1) and (1%n) are both evaluated right away:

r = 3
n = 11
print(rem(1, x**r-1), (1 % n))


Output:

1 1


Therefore, these are not no longer symbolic representations and would be eliminated immediately after being multiplied to the matrix.

All in all I do not think the overall algorithm could be implemented easily using the sympy library. It would be more convenient to implement the matrix power-modulo operation yourself with the help of numpy to work with polynominal computations (see doc). By the way, I am not sure whether that helps but numpy actually have a chebyshev module.

• Hand-rolled it's still slower than a naïve approach, but notably faster than my previous rewrite of OP's code. – Peter Taylor Sep 12 '19 at 16:12
• @PeterTaylor Given all those complex computations, there is no way that the performance of a pure python implementation could beat the naive approach on small numbers. The advantage on time complexity can only outweigh the overhead on significantly large inputs. Using numpy would delegate all the polynomial manipulation jobs to C, which should lead to some speedup. – GZ0 Sep 12 '19 at 16:42
• The disclaimer was preemptive: someone would have complained that it was still slower than trial division. But I wasn't actually basing anything on pure Python: I tested with PyPy. There may be some size of n at which NumPy gets an advantage, but its far far slower than the code I linked above for small n. – Peter Taylor Sep 12 '19 at 19:05
• Yeah, the computation is still way too much for small ns and numpy definitely would not save that. – GZ0 Sep 12 '19 at 23:06
• My point was that a hand-rolled implementation of Chebyshev polynomials modulo (n, x^r-1) using PyPy was orders of magnitude faster than an implementation using numpy polynomials. I.e. using numpy makes things worse, not better. PyPy seems to be a more effective way of bringing the speedup of native code to the problem. – Peter Taylor Sep 13 '19 at 7:38