-- Generate ramanujan numbers according to theorem 412 of Hardy & Wright

module Main where
import Prelude
import IO
import Data.Ratio
import System

-- the basic step from one rational point along the elliptic curve to another:
-- generate a list of r rational points along the same curve

rat_step :: Integer -> (Rational, Rational) -> [(Rational,Rational)]
rat_step r (x,y) | r <= 0     = []
                 | r == 1     = [(x,y)]
                 | otherwise  = (x,y) : (rat_step (r - 1) (gnew x y))
  where gnew x1 y1 = let x13 = x1^3
                         y13 = y1^3
                         d1 = x13 - y13
                         xi = x1*(x13 + 2*y13)/d1
                         yi = y1*(2*x13 + y13)/d1
                         xi3 = xi^3
                         yi3 = yi^3
                         d2 = xi3 + yi3
                         x2 = xi*(xi3 - 2*yi3)/d2
                         y2 = yi*(2*xi3 - yi3)/d2
                     in (x2,y2)

-- extract the denominators of all the pairs of rational numbers

denoms :: [(Rational, Rational)] -> [Integer]
denoms [] = []
denoms ((x,y):rs) = (denominator x):(denominator y):(denoms rs)

-- find the GCD and LCM of a list of integers, not just a pair of them

lgcd :: [Integer] -> Integer
lgcd l = foldl1 gcd l

llcm :: [Integer] -> Integer
llcm l = div (product l) (lgcd l)

-- rescale the rational numbers to integers: the argument s is constructed
-- such that the overall scaling will produce an integer with no remainder
-- in the division

rscale :: Integer -> [(Rational, Rational)] -> [(Integer, Integer)]
rscale s [] = []
rscale s ((x,y):rs) =
  ((rsc s x),(rsc s y)):(rscale s rs)
  where rsc s r = div ((numerator r)*s) (denominator r)

flatten :: [(Integer, Integer)] -> [Integer]
flatten [] = []
flatten ((i,j):lst) = i:j:(flatten lst)

iscale :: Integer -> [(Integer, Integer)] -> [(Integer, Integer)]
iscale s [] = []
iscale s ((x,y):rs) = ((div x s),(div y s)):(iscale s rs)

-- sanity-check the final result: compute sums of cubes, see that they all
-- add up to the same value

scube :: [(Integer, Integer)] -> [Integer]
scube [] = []
scube ((x,y):rs) = (x^3 + y^3):(scube rs)

checkit lst =
  if (ae (head lst) (tail lst)) then "good" else "bad!"
  where ae s [] = True
        ae s (x:xs) = if s /= x then False else ae s xs

-- do the whole shebang!

doit :: Integer -> Rational -> Rational -> String
doit r x y =
  let rlist = rat_step r (x,y)
      ilist = rscale (llcm (denoms rlist)) rlist
      jlist = iscale (lgcd (flatten ilist)) ilist
  in (show jlist) ++ "\n" ++ (checkit (scube jlist))

-- extract the Nth argument as an integer
rda alst n = read (alst !! n)

main =
  do args <- getArgs
     if (length args) >= 5
        then putStrLn (doit (rda args 0)
                           ((rda args 1) % (rda args 2))
                           ((rda args 3) % (rda args 4)))
        else putStrLn "usage: prog r xnum xden ynum yden"