Cornacchia Algorithm

1.5.2
The Algorithm of
Cornacchia
今回の内容
 Algorithm1.5.2(Cornacchia)
 Algorithm1.5.2の実装
 Algorithm1.5.2のExample
 Algorithm1.5.3(Modified-Cornacchia)
Algorithm1.5.2(Cornacchia)
p : prim enum ber( p  1(mod4)) を与えたとき ,
x2  y2  p
をみたす x, y  Zが存在する
より、一般的に・・・・
d  Z 0 , p : odd primeを与えたとき ,
x 2  dy2  p をみたす x, y  Z は存在するか ?
存在するとしたらそれは何か ?
Algorithm1.5.2(Cornacchia)
 Input: d , p (d  Z , p : Primenumber,0  d  p)
 Output: x, y  Z ( x2  dy2  p の整数解 x, y が存在すれば)
 Step1: k    d  ( Algorithm1.4.12を用いる )
 p 


k  1 なら , 整数解を持たない .
 Step2:
x0  d (mod p)をみたす
, x0 を求める
.
2
(Algorit hm1.5.1を用い
る)
x0 は, p / 2  x0  p をみたすように
する
(必要ならば, x  p  x0 とする)
a  p, b  x0 , l 
 p  とする
Algorithm1.5.2(Cornacchia)
 Step3: b  l の間, r  a(modb), a  b, b  r を繰り返す
このbがxの候補
になる
Euclid Algorithmの適用
2
d
が
p

b
で割り切れない ,または ,
 Step4:
c  ( p  b ) / d が平方数(Algorithm1.7.3)
2
ではないならば , 整数解は存在しない
そうでないとき , ( x, y )  (b, c )
Algorithm1.5.2(Cornacchia)
 証明について
 F. Morain, J.-L. Nicolas による証明がある.
(URL)
http://web.math.hr/~duje/tbkript/tbksem.html

x 2  dy2  p
について
 Primes of the form x2+ny2 (Cox, David A. 著)
に多くのことが書かれている.
Pythonでの実装
import math
import kro1
平方剰余の計算
import shanks
import square_test
a=p
b = x0
l = int(math.floor(math.sqrt(p)))
while b > l:
r=a%b
def cornacchia(d,p):
a=b
k = kro1.kro_b(-d,p)
Algorithm1.5.1
b=r
if k == -1:
c = (p - (b ** 2)) / d
return "no solution"
t = squaretest.square_test(c)
x0= shanks.shanks(-d,p)
if ((d % (p - (b ** 2)) != 0) or
if x0 == 0:
square_test.square_test(c) > 1):
return "no solution"
return (b,t)
if x0 < p/2:
return "no solution"
x0 = p – x0
Algorithm1.7.3
Pythonでの実装(Algorithm1.5.1)
import math
import random
import kro1
def shanks(c,g):
temp = g -1
e=0
while temp % 2 == 0:
temp = temp // 2
e=e+1
q = (g-1) // 2**e
k=0
while k != -1:
n =int(math.floor
(1000*random.random()))
k = kro1.kro_b(n,g)
z = (n ** q) % g
y=z
r=e
x = (c ** ((q - 1) // 2)) % g
b = (c * (x ** 2)) % g
x = (c * x) % g
while b % g != 1:
m=1
while ((b ** (2 ** m)) % g) != 1:
m=m+1
if m == r:
print “[a] is not a quadratic residue mod p"
return 0
t = (y ** (2 ** (r - m - 1))) % g
y = (t ** 2) % g
r=m%g
x = (x * t) % g
b = (b * y) % g
return x
Pythonでの実装(Algorithm1.7.1&3)
import math
def square_test(c):
q11 = []
q63 = []
q64 = []
q65 = []
for k in range(0,11):
q11.append(0)
for k in range(0,6):
q11[(k ** 2) % 11] =1
for k in range(0,63):
q63.append(0)
for k in range(0,32):
q63[(k ** 2) % 63] = 1
for k in range(0,64):
q64.append(0)
for k in range(0,32):
q64[(k ** 2) % 64] = 1
for k in range(0,65):
q65.append(0)
for k in range(0,33):
q65[(k ** 2) % 65] = 1
t = c % 64
if q64[t] == 0:
return "sono1"
r = c % 45045
if q63[r % 63] == 0:
return "sono2"
if q65[r % 65] == 0:
return "sono3"
if q11[r % 11] == 0:
return "sono4"
x=c
y=0
y = math.floor
((x + math.floor(c/x))/2)
while y < x:
x=y
y = math.floor
((x + math.floor(c/x))/2)
if int(math.floor(x ** 2))
== int(math.floor(c)):
return 1
else:
return 0
Example
x 2  2 y 2  97


2
Step1: 
  1 (algorithm1.4.12)
 97 
Step2: x0  17 (algorithm1.5.1)
p/2<x0<p
となるように
x0  97 17  80と置き換える
a  97, b  80, l 


Step3:
Step4:
a
97
80
17
12
b
80
17
12
5
 97  9
r
5
5
5
Euclid Algorithm
を適用
p  b 2 ( 97  25  72) d (2)
( p  b ) / d ( 72 / 2  36  6 )
2
 ( x, y )  (5,6)
2
平方数
になっている
Algorithm1.5.3(Modified Cornacchia)
x 2  Dy 2  4 p ( D  0 or 1, p : odd prime)
をみたす x,y を求める
(1)
4| D
2| x
x'  dy  p
2
2
( x'  x / 2, d  D / 4)
(2)
D  1 (mod8)
x 2  y 2  4 (mod8)
①
( p  2k  1とおく )
①かつ, x, y : even
x'  y '  p
2
2
( x'  x / 2, y '  y / 2, d  D )
(3) D  5(mod8)
x, y の偶奇は決定できない
Algorithm1.5.3(Modified Cornacchia)
 (3)の場合,Algorithm1.5.2を適用することができな
い
D  5(mod8)
のとき
 1.5.2を修正したAlgorithm1.5.3を用いることで,この
問題を解決することができる
Algorithm1.5.3(Modified Cornacchia)
 Input:
p : primenumber
D  Z 0 s.t. D  0 or 1 (mod4) and D  4 p

 Step1: p  2 のとき ,
2
2
x
,
y

Z
(
x

D
y
 4 p の解 x, y が存在すれば)
Output:
1.5.2と
の違い
D  8 が整数の平方であるならば ,
( x, y )  ( D  8 ,1)
それ以外のとき , 解は存在しない
 Step2:
D
k    ( Algorithm1.4.12を用いる .)
 p
k  1 なら , 解を持たない
Algorithm1.5.3(Modified Cornacchia)
, x0 を求める
 Step3: x0  D (mod p)をみたす
(Algorithm1.5.1を用い
る)
2
1.5.2との
違い
x0が mod2 で 0 と合同でないとき ,
x  p  x0 とする


a  2 p, b  x0 , l  2 p とする
 Step4: b  l の間, r  a(modb), a  b, b  r を繰り返す
Algorithm1.5.3(Modified Cornacchia)
 Step5: D が 4 p  b 2 で割り切れない ,または ,
c  ( p  b ) / D が平方数(Algorithm1.7.3)
2
でないならば , 解は存在しない
そうでないとき , ( x, y )  (b, c )

Download Report