読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

(15) Lenstraの楕円曲線法による素因数分解アルゴリズム

6番目のFermat数の因数を楕円曲線法で探してみましょう。

/* 4.3.13m */

/* 因数発見を示すグローバルスイッチ */
notFoundYet : true$

/* Maximaの組み込み関数gcdexは、最大公約数の符号が不定なので、正となるように規格化する関数を用意。 */
extendedGCD(m, n) := block([t],
  t : gcdex(m, n),
  return(t * signum(t[3])))$

ellipticAdd(t, u, n) := block([v : [0, 0], slope, iList],
  iList : extendedGCD(u[1] - t[1], n),
  if iList[3] # 1 then (
    notFoundYet : false,
    print(iList[3]),
    return([0, 0])),
  slope : mod((u[2] - t[2]) * iList[1], n),
  v[1] : mod(slope ^ 2 - u[1] - t[1], n),
  v[2] : mod(-t[2] + slope * (t[1] - v[1]), n),
  return(v))$

ellipticDouble(t, a, n) := block([v : [0, 0], slope, iList],
  iList : extendedGCD(2 * t[2], n),
  if iList[3] # 1 then (
    notFoundYet : false,
    print(iList[3]),
    return([0, 0])),
  slope : mod((3 * t[1] ^ 2 + a) * iList[1], n),
  v[1] : mod(slope ^ 2 - 2 * t[1], n),
  v[2] : mod(-t[2] + slope * (t[1] - v[1]), n),
  return(v))$

ellipticMul(m, z, a, n) := block([k, t, u],
  k : m,
  u : z,
  while mod(k, 2) = 0 do (
    u : ellipticDouble(u, a, n),
    k : k / 2),
  k : floor(k / 2),
  t : u,
  while k > 0 do (
    t : ellipticDouble(t, a, n),
    if mod(k, 2) = 1 then u : ellipticAdd(t, u, n),
    k : floor(k / 2)),
  return(u))$

n : 2 ^ (2 ^ 6) + 1$
a : 23$
blim : 100$
logc : log(sqrt(n))$
t : [1, 1]$
for p : 2 next next_prime(p) while notFoundYet and p < blim do (
  d : floor(logc / log(p)),
  print(p, " ", d),
  t : ellipticMul(p ^ d, t, a, n))$

使用した機能

274177がの因数であることがわかりました。