4n+1型の素数を平方数の和に分解する。

e.g. 5 = 1+4=1^2+4^2

1249 = 15^2+32^2

17669 = 70^2+113^2

225221 = 410^2+239^2 

import Number
import Data.List
pRoot p = g.foldl1' f $ [1..div (p-1) 2]
where f a b = mod (a*b) p
g x = min x $ p-x
sqSum p = let r = pRoot p
in fst.until((==1).snd) next$ ((1,r),div (r*r+1) p)
next ((x,y),k) =
let u = g$ x `mod` k
v = g$ y `mod` k
x' = abs$ (x*v-y*u) `div` k
y' = abs$ (x*u+y*v) `div` k
k' = (u*u+v*v) `div` k
in ((x',y'),k')
where g q = if q < div k 2 then q else q-k

X^2=-1 (mod p) を得るために pRoot を実行してるが、この関数は計算時間がよろしくない、約p/2回の乗算が必要。

別のいい方法はないか？

つまり、log(p)の多項式程度の計算回数で X^2 = -1 (mod p) を求めたい。

[追記]

pの平方非剰余が1つわかっていればできるらしい。

http://books.google.com/books?id=WZIrv5NDzE0C&pg=PA69&lpg=PA69&dq=%E5%B9%B3%E6%96%B9%E5%89%B0%E4%BD%99+%E5%B9%B3%E6%96%B9%E6%A0%B9+%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0&source=bl&ots=DORHN4Ln3s&sig=HfvbJ4nW__oRjLtqL93Wc_5x8_8&hl=ja&sa=X&oi=book_result&resnum=7&ct=result

[追記]

アルゴリズムにすると shanks-tonelli というらしい。

http://en.wikipedia.org/wiki/Shanks-Tonelli_algorithm

やはり、平方非剰余が必要。

実装してみた。 

import Data.List
gcd' m n | m < n = let (a,b) = gcd'  n m in (b,a)
| r ==  = (1,1-q)
| otherwise = let (a,b) = gcd' n r
in (b,a-b*q)
where (q,r) = divMod m n
powMod a n p | n ==  = 1
| even n = powMod (mod (a*a) p) (div n 2) p
| otherwise = (a * powMod a (n-1) p)`mod` p
jacobi a n | a <  && mod n 4 == 1 = jacobi (-a) n
| a <  && mod n 4 == 3 = - jacobi (-a) n
| a < 2 = a
| even a && (mod n 8 == 1 || mod n 8 == 7) = jacobi (div a 2) n
| even a && (mod n 8 == 3 || mod n 8 == 5) = - jacobi (div a 2) n
| mod a 4 == 3 && mod n 4 == 3 = - jacobi (mod n a) a
| otherwise = jacobi (mod n a) a
sqrtMod n p =
let m = fst.gcd' n$ p
w = head.filter((==(-1)).flip jacobi p)$[2..]
(s,q) = until(odd.snd) div2$ (,p-1)
v = powMod w q p
step u = case findIndex (==1).take s.iterate sqMod $ mod (u^2*m) p of
Just  -> u
Just i -> step $ (u*(powMod v (2^(s-i-1)) p)) `mod` p
Nothing -> error$ "There is no such x that x^2 ="++show n++" (mod "++show p++")"
in step $ powMod n (div (q+1) 2) p
where div2 (s,q) = (s+1,div q 2)
sqMod x = powMod x 2 p
sqSum p = let r = sqrtMod (p-1) p
in fst.until((==1).snd) next$ ((1,r),div (r*r+1) p)
where next ((x,y),k) =
let u = g$ x `mod` k
v = g$ y `mod` k
x' = abs$ (x*v-y*u) `div` k
y' = abs$ (x*u+y*v) `div` k
k' = (u*u+v*v) `div` k
in ((x',y'),k')
where g q = if q < div k 2 then q else q-k

平方剰余の平方根を求めることがメインになってしまった。