-- t5.hs: nhc98 Segmentation fault
-- 2003-Dec-03 / TN

module Main where

main
  = do
      putStr "t5: 2003-Dec-03 23.28\n"
      putStr ( "t5: " ++ show ( ec 2 ((2^127-1)^2) 50 1 ) ++ "\n" )

data ECPoint
  = ECUndefined { ecd :: Integer }
  | ECZero
  | ECNonzero { ecx, ecy :: Integer }
  deriving ( Show, Eq )

ec a n l1 d = osA0 (\(xs,a)->(1,show a)) (ecadd a n) ECZero pone ecext n l1 d

ecadd a n p1@( ECUndefined { ecd = d } ) p2 = p1
ecadd a n p1 p2@( ECUndefined { ecd = d } ) = p2

ecadd a n ECZero p = p
ecadd a n p ECZero = p

ecadd a n p1@( ECNonzero { ecx = x1, ecy = y1 } ) p2@( ECNonzero { ecx = x2, ecy = y2 } )
  = let
      (a1,b1,g1) = chinrem n (x1 - x2)
    in
      if g1 == 1 then
        eclambda n ((y1 - y2)*b1) p1 p2
      else
        if g1 == n then
          let
            (a2,b2,g2) = chinrem n (y1 + y2)
          in
            if g2 == 1 then
              eclambda n ((3*x1*x1 + a)*b2) p1 p2
            else
              if g2 == n then
                ECZero
              else
                ECUndefined { ecd = g2 }
        else
          ECUndefined { ecd = g1 }

ecext (ECUndefined { ecd=d }) = d
ecext _ = 1

eclambda n lambda ( ECNonzero { ecx = x1, ecy = y1 } ) ( ECNonzero { ecx = x2, ecy = y2 } )
  = let
      x3 = mod (lambda*lambda - x1 - x2) n
    in
      ECNonzero { ecx = x3, ecy = mod (lambda*(x1 - x3) - y1) n }

pone = ECNonzero 2 3

dikke a b = b `mod` a /= 0

pp l = [ last (takeWhile (<= max l p) (iterate (*p) p)) |
         p<-2:3:5:dropWhile (<=1) [ i*30+j | i<-[0..], j<-filter (dikke 5) (filter (dikke 3) (filter (dikke 2) [0..29] ) ) ]
       ]

osA0 cont op pzero pone ext n l
  = split
      cont
      n
      ( \(x:xs,_)-> x>l )
      ( \(x:xs,a)->(xs,repeter op pzero x a) )
      (\(_,p)->ext p)
      ( pp l, pone )

chinrem n 0 = (signum n,0,abs n)
chinrem n m
  = let
      (q,r) = divMod n m
      (a,b,g) = chinrem m r
    in
      (b,a-b*q,g)

repeter op p0 0 p = p0
repeter op p0 k p
  = let
      kpd2 = repeter op p0 (div k 2) p
      kp2 = op kpd2 kpd2
    in
      if mod k 2 == 0 then
        kp2
      else
        op p kp2

maxx ended j (x:xs) = if j <= 0 || ended x then x else maxx ended (j-1) xs

split cont n ended next ext s0 d
  = if ended s0 then
      cont s0
    else
      let
        si = maxx ended d (iterate next s0)
        f = gcd n (ext si)
      in
        if f /= 1 then
          if d > 1 then
            split cont n ended next ext s0 (d `div` 2)
          else
            (f,take 50 (show si))
        else
          split cont n ended next ext si d