-- verify 3: independent generation of Ramanujan numbers

module Main where
import Prelude
import System
import Data.List
import Text.ParserCombinators.Parsec

cbrtf x = floor (exp ((1.0/3.0) * log (fromInteger x)))
cbrtc x = ceiling (exp ((1.0/3.0) * log (fromInteger x)))
cbrth x = ceiling (exp ((1.0/3.0) * log (0.5*(fromInteger x))))

cbrtfb x = if x < 1 then 1 else cbrtf x
cbrtcb x = if x < 1 then 1 else cbrtc x

-- generate a list of triples (s, i, j), where
-- s = i^3 + j^3, i < j, and lower <= s <= upper

gen_j l u i = [cbrtfb (l - i^3) .. cbrtcb (u - i^3)]
gen_triple i j = (i^3 + j^3,i,j)
is_good l u (s,i,j) = (i < j) && (l <= s) && (s <= u)
gen_trips l u i = filter (is_good l u) (map (gen_triple i) (gen_j l u i))

-- sort the triples in order of increasing s

cmp1 (s1,i1,j1) (s2,i2,j2) =
  if (s1 == s2) then
        compare i1 i2
     else
        compare s1 s2
gen_sums l u = sortBy cmp1 (concat (map (gen_trips l u) [1 .. cbrth u]))

-- group the triples into groups with the same s and with the gcd of
-- all the (i,j) equal to 1; ie, primitive Ramanujan numbers

cmp2 (s1,i1,j1) (s2,i2,j2) = (s1 == s2)
gcd1 (s,i,j) = gcd i j
ggcd [] = 0
ggcd (x:xs) = gcd (gcd1 x) (ggcd xs)
bylen x = ((length x) > 1 && (ggcd x) == 1)

-- reshuffle the numbers so that there is just one s for each group

sval (s,i,j) = s
ijval (s,i,j) = (i,j)
get_comps x = (sval (head x), map ijval x)

-- format everything nicely

gl a b = a ++ " " ++ b
ppr (i,j) = (show i) ++ " " ++ (show j)
print_num (s,sl) =
  (show s) ++ "\t" ++ (tail (foldl gl "" (map ppr sl))) ++ "\n"
do_interval l u =
  map print_num ( map get_comps (filter bylen (groupBy cmp2 (gen_sums l u))))

-- main iterator over intervals

iloop s i e =
  if (s > e) then
        do putStrLn "All done!"
     else
        do putStrLn ("# start " ++ show s ++ " end " ++ show (s + i))
           putStr (foldl (++) "" (do_interval s (s + i)))
           iloop (s + i) i e

-- parser for a bigint: \d+(e\d+)?
-- this could conceivably be done more easily with a regexp?

expo =
  do oneOf ['e', 'E']
     eee <- many1 digit
     return eee

bignum =
  do mant <- many1 digit
     eee <- option "0" expo
     eof
     return (mant, eee)

convert x =
  let (m,e) = case parse bignum "" x of
                   Left err -> error $ "Input: " ++ x ++ "\n" ++ show err
                   Right result -> result
  in (read m) * (10^(read e))

main =
  do args <- getArgs
     let len = (length args)
         start = convert (args !! 0)
         end = convert (args !! 1)
         interval = convert (args !! 2)
     if (len < 3) then
           putStrLn "usage: verify3 start end interval"
        else
           iloop start interval end