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つわかっていればできるらしい。

##### [追記]

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

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

```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
(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
```