{-# OPTIONS_GHC -cpp -XParallelListComp -XScopedTypeVariables #-}
-- Time-stamp: <Sat Jun 26 2010 03:18:15 Stardate: Stardate: [-28]3280.27 hwloidl>
--
-- cash: the computer algebra shell
-- This is the start of a ghci-based shell, using SCSCP calls for computer algebra
-- Compiling sequentially picks the (older) SCSCP_API.hs module.
-- Compiling in parallel picks the multi-threaded client ParSCSCP.hs module.
-----------------------------------------------------------------------------

-- are we talking directly to a GAP server or a Haskell-side Coordination server

module Math.ComputerAlgebra.Cash (server, initServer, call0, call1, call2, mkRand, mult, resultant, sumEuler_seq, euler, hcf, relprime, eulerWithSCSCP, sumEulerWithSCSCP
                                 ,sumEulerListSCSCP, sumEulerFromToWithSCSCP, doClient, toOMMatrix
                                 ,toOMList, fromOMList, a_sparse_9, b_sparse_9, polyFromString, bagInter, myGcd)
where

import Math.ComputerAlgebra.Cash.Date

#undef GAP_SERVER

#undef  RUN_BUILTINS
#undef  RUN_POLYS 
#undef  RUN_SUMEULER 
#undef  RUN_SUMEULER_PAR 

import List (delete)
import Network
import System
import System.Exit
import System.IO
import System.IO.Unsafe
-- import Time
import Control.Monad
import Control.Concurrent
import qualified Control.Exception as C

import Math.ComputerAlgebra.Cash.SGPTypes
import Math.ComputerAlgebra.Cash.SCSCP_API
import Math.ComputerAlgebra.Cash.HS_SCSCP
import Math.ComputerAlgebra.Cash.HS2SCSCP

-- examples how to use it
-- import SCSCP_Ex -- cut down version of SCSCP_Examples

#ifdef HAVE_SKELETONS
-- abstractions over patterns of parallel coordination
import Math.ComputerAlgebra.Cash.CA_Skeletons
#endif

-- These are services, known to the client from the start
import Math.ComputerAlgebra.Cash.BaseServices

#ifdef __PARALLEL_HASKELL__ 
-- import EdenHelpers -- helper functions
-- import FoldDM
import Eden
#endif

-- import Poly
-- import TestPolys
-- import Karatsuba 

-- boilerplate setup, including pre-shared info
#if 0
main :: IO()
main = do
        args <- getArgs
        let portNum = if null args 
                         then 12321
                         else fromInteger (read (head args))
        doClient portNum
#endif

ngoq = initServer (server "localhost" (Just 12321))

-------------------------------------------------------
-- Util fcts (from ParSCSCP.hs)

call0 :: CAName -> OMObj 
call0 name = unsafePerformIO (putStrLn $ "call0 of "++(show name)) `seq`   
             (callSCSCP name [])

call1 :: (OMData a, OMData b) =>
         CAName -> a -> b
call1 name x  = fromOM (callSCSCP name [toOM x])

call2 :: (OMData a, OMData b, OMData c) =>
         CAName -> a -> b -> c
call2 name x y = fromOM (callSCSCP name [toOM x,toOM y])

mkRand n = callSCSCP scscp_CS_RandomPolynomialAsString  ( map toOM [n] )

mult p1 p2 = callSCSCP scscp_CS_KaratsubaStr_x [p1, p2]

resultant p1OM_str p2OM_str = callSCSCP scscp_CS_Resultant ( map toOM [p1OM_str, p2OM_str] )

-------------------------------------------------------
-- purely Haskell side code

-- 
fact :: Integer -> Integer
fact 0 = 1
fact n = n*(fact (n-1))

factAcc :: Integer -> Integer -> Integer
factAcc 0 acc = acc
factAcc n acc = factAcc (n-1) (n*acc)

sumEuler_seq :: Int -> Int
sumEuler_seq = sum . map euler . enumFromTo (1::Int)

---------------------------------------------------------------------------
-- main fct

euler :: Int -> Int
euler n = length (filter (relprime n) [1..(n-1)])

---------------------------------------------------------------------------
-- aux fcts
hcf     :: Int -> Int -> Int
hcf x 0 = x
hcf x y = hcf y (rem x y)

relprime     :: Int -> Int -> Bool
relprime x y = hcf x y == 1

-----------------------------------------------------------------------------

eulerWithSCSCP :: Int -> IO Int
eulerWithSCSCP n = do
  -- do the computation ...
  let fName = scscp_WS_Phi
  let args  = map toOM [n]
  let resOM = callSCSCP fName args
  -- print the result
  let x = (fromOM resOM) :: Int
  return x

sumEulerWithSCSCP :: Int -> IO Int
sumEulerWithSCSCP n = do
   xs <- mapM eulerWithSCSCP [1..n]
   return (sum xs)

sumEulerListSCSCP :: [Int] -> IO Int
sumEulerListSCSCP ns = do
   xs <- mapM eulerWithSCSCP ns
   return (sum xs)
                      
sumEulerFromToWithSCSCP :: Int -> Int -> IO Int
sumEulerFromToWithSCSCP m n = do
   xs <- mapM eulerWithSCSCP [m..n]
   return (sum xs)

-----------------------------------------------------------------------------

-- client for calling Euler totient function 
-- currently only works with the purely Haskell dummyServer
-- GAP server works fine up to Phi function, but needs name scscp_WS_Phi 
-- polynomials only supported for the Haskell server
doClient :: PortNumber -> IO ()
doClient portNum = do
                    --  init ...
	            putStrLn ("starting up client, opening port " ++ show portNum)
                    initServer (server "localhost" (Just portNum))
#ifdef RUN_BUILTINS
                    -------------------------------------------------------
                    -- ask for available services
	            putStrLn $ "Request for GetServiceDescr ..."
	            let fName = Right GetServiceDescr
                    let args = []
                    let resOM = callSCSCP fName args
                    putStrLn $ "Reply: "++(show resOM)
                    -------------------------------------------------------
# ifdef GAP_SERVER
                    -- ask for available services
	            putStrLn $ "Request for GetTransientCD ..."
	            let fName = Right GetTransientCD
                    let args = map toOM ["scscp_transient_1"]
                    let resOM = callSCSCP fName args
                    putStrLn $ "Reply: "++(show resOM)
# endif
                    -------------------------------------------------------
                    -- ask for available services
	            putStrLn $ "Request for GetAllowedHeads ..."
	            let fName = Right GetAllowedHeads
                    let args = []
                    let resOM = callSCSCP fName args
                    putStrLn $ "Reply: "++(show resOM)
                    -------------------------------------------------------
# ifdef GAP_SERVER
                    -- ask for available services
	            putStrLn $ "Request for GetSignature ..."
	            let fName = Right GetSignature 
                    let args = map toOM ["scscp_transient_1", "WS_Phi"] -- HWL: BROKEN: wrong encoding of args; should be in OMS
                    let resOM = callSCSCP fName args
                    putStrLn $ "Reply: "++(show resOM)
# endif
#endif
                    -------------------------------------------------------
                    -- do the computation ...
	            putStrLn $ "Running phi 12 ..."
	            let fName = 
#  ifdef GAP_SERVER
                                scscp_WS_Phi 
#  else
                                scscp_CS_Phi
#  endif
                    let args = map toOM [12::Int]
                    let resOM = callSCSCP fName args
                    -- print the result
                    let x = (fromOM resOM) :: Int
	            putStrLn $ "Result: "++(show x)
#if 0
	            putStrLn $ "Running factorial 5 ..."
	            let fName = 
#  ifdef GAP_SERVER
                                scscp_WS_Factorial
#  else
                                scscpFact
#  endif
                    let args = map toOM [5::Int]
                    let resOM = callSCSCP fName args
                    -- print the result
                    let x = (fromOM resOM) :: Int
	            putStrLn $ "Result: "++(show x)
#endif
                    -- do the computation ...
#ifdef RUN_SUMEULER
                    let n = 87
	            putStrLn $ "Running sumEulerWithSCSCP "++(show n)++" ..."
                    x <- sumEulerWithSCSCP n
	            putStrLn $ "Result: "++(show x)
	            putStrLn $ "Running sumEuler_seq "++(show n)++" ..."
	            let y = sumEuler_seq n
	            putStrLn $ "Result: "++(show y)
                    putStrLn $ "Are the two results the same: "++(show (x==y))
#endif
                    -----------------------------------------------------------------------------         
                    -- do the computation ...
                    -- shutdown
                    releaseServer

#ifdef __PARALLEL_HASKELL__
linearSolver :: [ Integer ]        
            -> [ [Integer] ]      
            -> [ Integer ]
            -> Arith          
linearSolver ps ms vs
  = multHomImg (call2 scscp_WS_Mod) (call2 scscp_WS_Sol) (call2 scscp_WS_CRA) ps (toOMList ((toMatrix ms):[toOMList (map toNum vs)]))


toNum :: Integer -> Arith
toNum x = Num x

toVector :: [Integer] -> Arith
toVector xs = MatrixRow (map (\z -> Num z) xs)

toMatrix :: [ [Integer] ] -> Arith
toMatrix (x:xs) = Matrix (map (\z -> MatrixRow (map toNum z)) (x:xs)) 

instance NFData (Arith) 
instance Trans (Arith)


parZipWith f primes list2 -- list3
  = newTasks
    where
      workerProcs = [process (zip [n,n..] . (worker f)) | n <- [1..noPe] ]

      (newReqs, newTasks) = (unzip . concat) (zipWith ( # ) workerProcs (distributeLists (primes, list2) requests)) -- (zipWith (#) workerProcs (distributeLists (list1, list2) requests))

      requests = (concat (replicate 2 [1..noPe])) ++ newReqs

      -- worker f [] = []
      -- worker :: (Trans a, Trans b) => (a -> b -> c) -> [(a,b)] -> [c]
      worker f [] = []
      worker f ((p, t2) : ts) = (f p t2) : worker f ts

      -- distributeLists :: ([t], [t]) -> [Int] -> [[(t,t)]]
      distributeLists tasks reqs = [ taskList reqs tasks n | n <- [ 1 .. noPe ] ]
        where
           taskList (r:rs) ( p:ps, t2:ts2 ) pe | pe == r   = (p, t2) : (taskList rs (ps, ts2) pe)
                                               | otherwise =      taskList rs (ps, ts2) pe
           taskList _      _                _  = []

          
multHomImg :: -- :: (Trans p, Trans c', Trans b') =>
              (Integer -> Arith -> Arith ) ->  -- map input to homomorphic images
              (Integer -> Arith -> Arith) -> -- solve the problem in the hom. imgs.
              ([Integer] -> Arith -> Arith) ->  -- combine the results to an overall result              
              [ Integer ] ->             -- hom. imgs. to use. NOTE - will be supplied in a vector format!
              Arith ->               -- input (matrix, vector)
              Arith              -- result-}
multHomImg h f g ps x =
 res
 where
    xList    = zipWith h ps (repeat x)
    resL     = parZipWith f ps xList
    res      = g ps (toOMMatrix resL)    

#endif

toOMMatrix xs = Matrix xs
toOMList xs = Math.ComputerAlgebra.Cash.SGPTypes.List xs

fromOMList (List xs) = xs
fromOMList (Matrix xs) = xs
fromOMList (MatrixRow xs) = xs

primes :: [Integer]
primes = sieve [2..]
   where
    sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0]  

a_sparse_9 = 
        [  [ 763,    0,    0,    0,  633,    0,    0,  45,    0],
           [  0,    1,    0,    0,    0,    0,    0,    0,    0],
           [  0,    0,    3,    0,    0,    0,    0,    0,   42],
           [  0,  125,    0,    1,    0,    0,  572,    0,    0],
           [ 27,    0,   35,    0,    1,    0,    0,    0,   11],
           [  0,    0,    0,    0,    0,    1,   36,    0,    0],
           [  0,    0,   26,    0,    0,    0,    1,    0,   46],
           [  0,   19,    0,    0,    0,    0,    7,    1,    0],
           [ 57,    0,    0,    0,    2,    0,    0,    0,   92 ]]
b_sparse_9 = [0, 0, 0, 7, 0, 0, 0, 0, 11]

newtype Polynomial a = P OMObj

instance Show a => Show (Polynomial a) where
  show (P pOM) = fromOM pOM

instance Eq a => Eq (Polynomial a) where
  (P p1) == (P p2) = p1 == p2

instance  Num (Arith {-Polynomial a-}) where
  (*) p1@(Polynomial _) p2@(Polynomial _) = call2 scscp_WS_ProdPoly p1 p2
  (+) p1@(Polynomial _) p2@(Polynomial _) = call2 scscp_WS_SumPoly p1 p2
  abs = error "abs on Polynomial not implemented"
  signum = error "signum on Polynomial not implemented"
  fromInteger n = fromOM $ toOM ("0*x_1+"++(show n))

polyFromString :: String -> Arith -- Polynomial Int
polyFromString = fromOM . toOM

instance OMData (Polynomial a) where
  toOM (P pOM) = pOM
  fromOM pOM = P pOM

class Factorisable a where
  factors :: a -> [a]

instance Factorisable Int where
  factors = call1 scscp_WS_FactorsInt

instance Factorisable Integer where
  factors = call1 scscp_WS_FactorsInt

instance Factorisable Arith where
  factors p1@(Polynomial _) = call1 scscp_WS_Factors p1



bagInter :: (Eq a) => [a] -> [a] -> [a]
bagInter [] _ = []
bagInter (x:xs) ys | elem x ys = x:(bagInter xs (delete x ys))
                   | otherwise = bagInter xs ys

myGcd :: (Num a, Factorisable a) => a -> a -> a
-- myGcd :: Arith -> Arith -> Arith
myGcd x y = let 
              xs = factors x
              ys = factors y
              zs = xs `bagInter` ys
            in
            if null zs then fromInteger 1 else  foldl1 (*) zs