カーマイケル数（Carmichael number）
あまりなじみの無い名前ですがカーマイケル数（Wikipedia）によると、
\begin{align}
&カーマイケル数とは、nと互いに素である任意の自然数 a に対し\\
&\large a^{n-1} \equiv 1 \pmod n \\
&となるnのこと
\end{align}
この式を見るとフェルマーの小定理の合成数版と考えると分かりやすいかもしれません。
今回はこのカーマイケル数をできるだけ多く探していこうと思います。
カーマイケル数（Carmichael number）の判定方法
\begin{align}
&カーマイケル数の必要十分条件は\\
&1. 素数ではない\\
&2. 無平方数（\mathrm{squarefree\ number}）である\\
&3. すべての素因子pに対して p-1 | n-1 となる\\
&（ p-1 | n-1 はp-1がn-1の約数だということをを表す）
\end{align}
この判定方法で$n \le 2000$のカーマイケル数を求めます。
import sympy as sp
N = 2*10**3
total = 0
for n in range(3, N+1, 2):
f = sp.factorint(n)
if len(f)>1 and max(f.values()) == 1 and all((n-1)%(p-1)==0 for p in f):
print("#", n, f)
total += n
print(f"# N = {N}, total={total}")
# 561 {3: 1, 11: 1, 17: 1}
# 1105 {5: 1, 13: 1, 17: 1}
# 1729 {7: 1, 13: 1, 19: 1}
# N = 2000, total=3395
アルゴリズムの高速化（無平方数を生成）
さすがにこの方法だとすべての奇数について因数分解をしているので、もっと良い方法がありそうです。
判定方法の２．に 無平方数（squarefree number）というのがあるので素数のリストからすべてのN以下の無平方数を生成する方法が考えられます。
ただ問題はその最大素因数pですが、よく考えると$p \le \sqrt{N}$であることが証明できます（下に証明を添付しました。）
この考え方を使ったコードが以下になります。なんとか$N=10^8$ぐらいまでは求めることが出来ました。
from sympy import nextprime
def Carmichael(n, plist):
ret, p1 = 0, 2 if plist == [] else plist[-1]
while p1 <= int(N**0.5): # prime factor is less than sqrt(N)
p1 = nextprime(p1)
n1, plist1 = n * p1, plist+[p1]
if n1 > N: break
if len(plist1) > 1 and all((n1-1)%(q-1) == 0 for q in plist1):
ret += n1
ret += Carmichael(n1, plist1)
return ret
N = 10**8
total = Carmichael(1, [])
print(f"N = {N}, total = {total}")
# N = 100000000, total = 7145594043
カーマイケル数の最大素因数pがNの平方根以下になる証明
\begin{align}
&カーマイケル数をn、最大素因数pとすると\\
&n = pr = (p-1)r + r\\
&n-1 = (p-1)r + r-1 \\
&p-1 | n-1 なので p-1 | r-1 \\
&p(p-1) | p(r-1) = n -p \\
&p(p-1) | n -p \\
&p(p-1) \le n - p \\
&p^2 \le n\\
\end{align}
（開発環境：Google Colab）
この考え方はProject Euler, Problem 536: Modulo Power Identityを解くのに