I'm reading the book Prime and Programming and I'm stuck on one of the computer exercises.
I'm checking for Fermat Pseudoprimes and I've written a program that works for reasonably small numbers, e.g. $\le 10^6$. It just brute-force checks that $b^{n-1} \equiv 1\bmod n$.
Now I'm asked to find the smallest prime or base 2 pseudoprime greater that $10^{15}$.
Python is unable to calculate $2^{10^{15}-1} \bmod 10^{15}$, never mind check case after case.
I must be missing a piece of Maths that will simplify the calculation massively.
Any suggestions?
To calculate powers modulo n you can use Exponentiation by squaring together with Modular arithmetic; the Python built-in function pow does that (not math.pow, which accepts only 2 parameters).
pow (x, y, n) calculates (x**y) % n.
The paper New Implementations of Tabulating Strong Pseudoprimes and Liars, Wuyang Liu, Illinois Wesleyan University describes an algorithm that calculates (strong) Fermat Pseudoprimes up to a given bound. But for a range of big numbers as $n > 10^{15}$ that still requires much time.
The first 16-digit strong pseudoprime base 2 is 1000003918690669 (OEIS-A062568).
If you look for some, but not necessarily the next, pseudoprime you can calculate candidates ($q$) with the following formula in nested loops, e. g. k running from 1 to 99, m from k+1 to 100 and x such that $q > 10^{15}$ and then test if $q$ is a strong pseudoprime. Those number are more often strong pseudoprime than odd nunbers in general.
$\ q = (2\cdot k \cdot x + 1) \cdot (2 \cdot m \cdot x + 1)\ $ with
$\quad x \in \mathbb N$,
$\quad \gcd(k, m) = 1\ $ and
$\quad 0 < k \ll q,\ 0 < m \ll q\ $