cash-0.1.0.0: src/Math/ComputerAlgebra/Cash.hs
{-# 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