{-# LANGUAGE NoMonomorphismRestriction #-}
{-
   The Taylor series for the exponential function e^x:
    e^z = lim n|->inf \Sum_1^inf z^n / n!
   See: http://en.wikipedia.org/wiki/Power_series

   See the corresponding C+GMP version for comparison

   Compile:  ghc -O2 -threaded -rtsopts -o power_circ_par power_circ_par.hs
   Run:      ./power_circ_par 4000 10 1 +RTS -N7 -sstderr
   Measure:  time ./power_circ_par 4000 10 1 +RTS -N7 -sstderr 1>/dev/null

   Parallel version obtained by replacing sum function with parfoldr.
-}

import System.Environment
import Control.Monad
import Data.Ratio
import Control.Parallel
import Control.Parallel.Strategies
import Control.DeepSeq
import Data.List (foldl')
import GHC.Conc (numCapabilities)

fact 0 = 1
fact n = n*(fact (n-1))

-- compute e^z up to d digits 
power_e :: Integer -> Integer -> Integer -> Rational
power_e z d v = case v of
                     1 -> parfoldr (+) 0  (takeWhile (>(1 % (10^d))) taylor)
                     2 -> parfoldl (+) 0  (takeWhile (>(1 % (10^d))) taylor)
                     3 -> parfoldl' (+) 0  (takeWhile (>(1 % (10^d))) taylor)
                     _ -> sum  (takeWhile (>(1 % (10^d))) taylor)
              where -- infinite list of entire Taylor series
            	    taylor =  [ (pow_ints!!n % 1) / (factorials!!n % 1) | n<-[0,1..] ] 
	            -- 2 circular lists, to memoise earlier computations 
	            factorials = 1:1:[ (toInteger i)*(factorials!!(i-1)) | i<-[2..] ] 
                    pow_ints = 1:[ z*(pow_ints!!(i-1)) | i<-[1..] ]

-- TODO: change based on numCapabilities
chunksize xs = 100 --l `div` sparksPerCapability
          where
                l = length xs
                sparksPerCapability = 2 * numCapabilities

chunk n [] = [] 
chunk n xs = ys : chunk n zs 
  where (ys,zs) = splitAt n xs 


parfold fold f z [] = z 
parfold fold f z xs = res 
  where 
        parts = chunk (chunksize xs) xs 
        partsRs = map (fold f z) parts `using` parList rdeepseq 
        res = fold f z partsRs 

parfoldr = parfold foldr
parfoldl = parfold foldl
parfoldl' = parfold foldl'


main = do
         args <- getArgs
         let (z:d:v:_) = map read args :: [Integer]
         let res = power_e z d v
         putStrLn $ "e^"++(show z)++" with a precision of "++(show d)++" digits is "++(res `pseq` " omitted")
         {-
         putStrLn $ "e^"++(show z)++" with a precision of "++(show d)++" digits is "++(show ((fromRational res)::Double))  
         putStrLn $ " or as Rational " ++(show res)
         putStrLn $ "Expected result: " ++(show (exp (fromInteger z)))
         -}