4n+1素数の平方和分解

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 == 0 = (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 == 0 = 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 < 0 && mod n 4 == 1 = jacobi (-a) n
           | a < 0 && 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$ (0,p-1)
        v = powMod w q p
        step u = case findIndex (==1).take s.iterate sqMod $ mod (u^2*m) p of
                   Just 0 -> 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

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