Miller-Rabin 素数判定法

PyCryptodomeのisPrime はRandom BasisのMiller-Rabin

miller rabinをかいくぐる合成数の作りかた

• baseを固定した時: baseについてのstrong pseudoprime。 Arnault's method

• baseがランダムな時:

• ランダムな base についての 1/4 の誤判定率(理論値)をほぼ達成する数は lambda(x) をいい感じにするように適当に生成すれば出てきた。factordb に誤判定したリストがあったので眺めて規則性を見つけたが、 なんでこれがいい感じなのかは分からない……

• いい感じ: xが大きな素数数個の積で表されており、かつ lambda(x) が大きな素因数を一つもつ

• 誤判定リスト: http://www.factordb.com/prooffailed.php

def mr_test(n, d, r, a):
  a = Zmod(n)(a)
  y = a ^ d
  if y == 1 or y == -1: return True
  for _ in range(r - 1):
    y ^= 2
    if y == -1: break
  else:
    return False
  return True

def liar_prob(n):
  d = n - 1
  r = 0
  while d % 2 == 0:
    d /= 2
    r += 1
  if n <= 1000:
    l = list(range(1, n))
  elif n <= 10000:
    l = list(range(1, n))
    random.shuffle(l)
    l = l[:1000]
  else:
    s = set()
    R = Zmod(n)
    while len(s) < 1000:
      s.add(random.randrange(1, n))
    l = list(s)
  res = float(len([i for i in l if mr_test(n, d, r, i)]) / len(l))
  return res

def find_spsp(B):
    while True:
        p = getPrime(B)
        a = p * 2 * 3 + 1
        b = p * 4 * 3 + 1
        if isPrime(a): print(a, b)
        if isPrime(a) and isPrime(b):
            pr = liar_prob(a * b)
            print(pr)
            if 0.2 <= pr:
                print(a*b, pr)
                return a * b
               
find_spsp(128)
# 4237060777988884017800671780056666555374693122698792516704419635488944279425403 0.249