-- verify 2: given a possible Ramanujan number, generate all its components;
-- try to check random failures of computers, on the theory that even if one
-- or two sums-of-components get corrupted during a run, there will still be
-- something; thus, if the output of this test is the same as the input,
-- then there was (probably) no such corruption.

module Main where
import Data.Ratio
import Text.ParserCombinators.Parsec

-- data describing a Ramanujan sum: the Ramanumber, then individual
-- components which when cubed and added combine to Ramanumber

data Rama = Rama Integer [(Integer,Integer)] deriving Show

-- parser for a file containing a bunch of Ramanujan sums:
-- this generates a list of Rama objects

int_num =
  do s <- option '+' (oneOf "+-")
     d <- many1 digit
     let sgn = if s == '+' then 1 else -1
         val = read d
     return (sgn*val)

pair_comp =
  do oneOf [' ', '\t']
     val1 <- int_num
     oneOf [' ']
     val2 <- int_num
     return (val1, val2)

rama_line =
  do sum <- many1 digit
     pairs <- many1 pair_comp
     newline
     return (Rama (read sum) pairs)

rama_lines =
  do lines <- many1 rama_line
     eof
     return lines

get_rsum (Rama i (xs)) = i

-- functions to regenerate the list of components given the Ramanumber

-- various things dealing with square and cube roots

sqrti x = floor (0.5 + (sqrt (fromInteger x)))
cbrti x = 10 + ceiling (exp ((1.0/3.0) * log (fromInteger (4*x))))
is_square x = x == ((sqrti x)^2)
dif r k =
  let n = 4*r - k^3
      d = 3*k
      (qq,qr) = quotRem n d
  in (r, k, qq, qr)
isc (i,k,l,r) = (r == 0) && (is_square l)
mog1 (r,k,d,_) = (r, k, (sqrti d))
fint (r,k,l) = (even k) == (even l)
mog2 (r,k,l) = (r, (k - l) `div` 2, (k + l) `div` 2)
srt (i,j,k) = (j < k)

{- this is where all the heavy lifting happens: from right to left

	generate list of k values:
	    [1 .. cbrti x]
	transform list entries into 4-tuples of (r, k, maybe-l, remainder):
	    map (dif x) ...
	filter out those that can't be right:
	    remainder non-zero, or maybe-l isn't a perfect square:
	    filter isc ...
	transform 4-tuples into 3-tuples
	    map mog1 ...
	filter out those that don't generate integral r
	    filter fint
	final transformation to turn (r,k,l) into (r,i,j):
	    map mog2
	done! so simple, really! :-}

rama_gen x =
  map mog2 (filter fint (map mog1 (filter isc (map (dif x) [1 .. cbrti x]))))

-- A bunch of stuff to format the output exactly the way that ari* does it.

tailOrNot [] = "meh"
tailOrNot lst = tail lst

gl m a b = a ++ m ++ b
gls = gl " "
gln = gl "\n"
comp1 (i,j,k) = gls (show j) (show k)
comp2 l = tailOrNot (foldl gls "" (map comp1 l))

glom [] [] = []
glom (r:rs) (c:cs) = (gl "\t" (show r) c) : glom rs cs

main = do input <- getContents
          let val = case parse rama_lines "stdin" input of
                         Left err -> error $ "Input:\n" ++ show input ++ 
                                             "\nError:\n" ++ show err
                         Right result -> result
              ramas = map get_rsum val
              internal = map comp2 (map rama_gen ramas)
              display = glom ramas internal
          putStrLn (foldl gln "Haskell sez:" display)