Commit cfcf94ae authored by Simon Marlow's avatar Simon Marlow

update to version from Seq no More paper

parent ebb8287d
{-
$Id: CRA.hs,v 1.1 1996/01/08 20:07:34 partain Exp $
This is revision: $Revision: 1.1 $
Module realizing the Chinese Remainder Algorithm (on 2 inputs as well as
on 2 input lists).
Changelog:
$Log: CRA.hs,v $
Revision 1.1 1996/01/08 20:07:34 partain
Initial revision
--# Revision 1.2 1994/12/13 19:53:03 hwloidl
--# tree_IMCRA0 is a tree-based version of CRA over lists.
--# Nomally it works over infinite lists (with -DNO_FILTER over finite lists;
--# NOTE: Depending on NO_FILTER different tree_IMCRA0 functions are compiled
--# with *different typings*. So, any module using tree_IMCRA0 must be
--# recompiled, too, when changing this flag).
--# s2IMCRA is the same as s1IMCRA (i.e. fold-based) but works on infinite lists.
--# tree_IMCRA0 handles the det-check now (therefore list of dets as input).
--# It tries to push the det-check as much backwards as possible by using the
--# handle fails fct (possible problem: only fails is immediately needed when calling
--# handle_fails)
--#
--# Revision 1.1 1994/12/08 19:24:21 hwloidl
--# Initial revision
--#
---------------------------------------------------------------------------- -}
module CRA (binCRA,s1IMCRA,s2IMCRA,tree_IMCRA0,
isPrime) where
import ModArithm (modHom, modDif, modProd, modInv)
-- Time-stamp: <2010-11-03 10:31:32 simonmar>
--
-- Chinese Remainder Algorithm.
-- Works over lists of Integers (2 input lists).
----------------------------------------------------------------------------
import Random (randomInts) -- Just for testing
import ParForce
-- @node Chinese Remainder Algorithm, , ,
-- @chapter Chinese Remainder Algorithm
#ifdef SEQ
module CRA (binCRA, seq_list_CRA0, seq_list_CRA,
-- par_binCRA, tree_IMCRA0,
isPrime) where
par :: a -> b -> b
par x y = y
-- @menu
-- * Imports::
-- * Auxiliary functions::
-- * Basic binary CRA operation::
-- * CRA over lists::
-- @end menu
seq :: a -> b -> b
seq x y = y
-- @node Imports, Auxiliary functions, Chinese Remainder Algorithm, Chinese Remainder Algorithm
-- @section Imports
_parGlobal_ :: Int# -> a -> b -> b
_parGlobal_ n x y = y
import ModArithm (modHom, modDif, modProd, modInv)
_parLocal_ :: Int# -> a -> b -> b
_parLocal_ n x y = y
#if defined(STRATEGIES)
import Control.Parallel
import Control.Parallel.Strategies
#endif
#ifdef GUM
-- GUM has only par
import {-fool mkdependHS; ToDo: rm-}
Parallel
_parGlobal_ :: Int# -> a -> b -> b
_parGlobal_ n x y = par x y
_parLocal_ :: Int# -> a -> b -> b
_parLocal_ n x y = par x y
#endif
-- @node Auxiliary functions, Basic binary CRA operation, Imports, Chinese Remainder Algorithm
-- @section Auxiliary functions
isPrime :: Integer -> Bool
--isPrime :: Integer -> Bool
isPrime 2 = True
isPrime n
| even n = False
| otherwise = isPrime' n 3
isPrime' :: Integer -> Integer -> Bool
--isPrime' :: Integer -> Integer -> Bool
isPrime' n l1 | n < l1*l1 = True
| n `mod` l1 == 0 = False
| otherwise = isPrime' n (l1+2)
-- @node Basic binary CRA operation, CRA over lists, Auxiliary functions, Chinese Remainder Algorithm
-- @section Basic binary CRA operation
#if 0 && defined(STRATEGIES)
par_binCRA :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer
par_binCRA m1 m2 inv a1 a2 =
(if d == 0 then a1
else a)
`sparking` rnf a
where ab = modHom m2 a1
d = _parGlobal_ 0# a
modDif m2 a2 a1
d = modDif m2 a2 a1
b = modProd m2 d inv
b0 = if (b+b>m2) then b-m2
else b
a = m1*b0+a1 {- not HWL version -}
a = m1*b0+a1
#endif
binCRA :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer
......@@ -102,25 +71,34 @@ binCRA m1 m2 inv a1 a2 = (if d == 0 then a1
else b
a = m1*b0+a1 {- not HWL version -}
s1IMCRA :: [Integer] -> [Integer] -> (Integer,Integer)
-- @node CRA over lists, , Basic binary CRA operation, Chinese Remainder Algorithm
-- @section CRA over lists
s1IMCRA (m:ms) (r:rs) = foldl
(\ (m0,r0) (m1,r1) -> (m0*m1,binCRA m0 m1 (modInv m1 m0) r0 r1) )
(m,r) (zip ms rs)
-- @node Sequential, , CRA over lists, CRA over lists
-- @subsection Sequential
seq_list_CRA0 :: [Integer] -> [Integer] -> (Integer,Integer)
seq_list_CRA0 (m:ms) (r:rs) =
foldl (\ (m0,r0) (m1,r1) -> (m0*m1, binCRA m0 m1 (modInv m1 m0) r0 r1) )
(m,r)
(zip ms rs)
-----------------------------------------------------------------------------
-- Same as s1IMCRA but applicable for infinite lists
-----------------------------------------------------------------------------
s2IMCRA :: Integer -> [Integer] -> [Integer] -> [Integer] -> (Integer,Integer)
s2IMCRA mBound (m:ms) (r:rs) (d:ds) = s2IMCRA' mBound ms rs ds m r
-- Same as above but applicable for infinite lists
s2IMCRA' :: Integer -> [Integer] -> [Integer] -> [Integer] -> Integer -> Integer ->
seq_list_CRA :: Integer -> [Integer] -> [Integer] -> [Integer] -> (Integer,Integer)
seq_list_CRA mBound (m:ms) (r:rs) (d:ds) = seq_list_CRA' mBound ms rs ds m r
seq_list_CRA' :: Integer -> [Integer] -> [Integer] -> [Integer] -> Integer -> Integer ->
(Integer,Integer)
s2IMCRA' mBound (m:ms) (r:rs) (d:ds) accM accR
seq_list_CRA' mBound (m:ms) (r:rs) (d:ds) accM accR
| accM > mBound = (accM, accR)
| d == 0 = s2IMCRA' mBound ms rs ds accM accR
| otherwise = s2IMCRA' mBound ms rs ds (m*accM) (binCRA accM m (modInv m accM) accR r)
| d == 0 = seq_list_CRA' mBound ms rs ds accM accR
| otherwise = -- trace ("seq_list_CRA': (m,accM) = " ++ (show (m,accM))) $
seq_list_CRA' mBound ms rs ds (m*accM)
(my_cra accM m (modInv m accM) accR r)
where my_cra = binCRA {- or par_binCRA -}
-----------------------------------------------------------------------------
-- Note:
......@@ -135,102 +113,55 @@ s2IMCRA' mBound (m:ms) (r:rs) (d:ds) accM accR
-- Note: ms as ds are infinte lists
tree_IMCRA0 :: Integer -> [Integer] -> [Integer] -> [Integer] ->
(Integer,Integer)
(Integer,Integer)
tree_IMCRA0 n ms as ds =
let
res@(m, a, fails) = {-_parGlobal_ 11# (forcelist ms')
(_parGlobal_ 12# (forcelist as')
(_parGlobal_ 13# (forcelist ds') -}
tree_IMCRA0' ms' as' ds' -- ))
where ms' = take n ms
as' = take n as
ds' = take n ds
handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
-- _parGlobal_ 46# m $
-- _parGlobal_ 47# a $
if n == 0 then (m, a)
else if d1 == 0 then (handle_fails n m a ms as ds)
else handle_fails (n-1) m' a' ms as ds
where
m' = m * m1
a' = par_binCRA m m1 inv a a1
inv = modInv m1 m
in
handle_fails fails m a ms as ds
res@(m, a, fails) = {-parGlobal 11 11 (forcelist ms')
(parGlobal 12 12 (forcelist as')
(parGlobal 13 13 (forcelist ds') -}
tree_IMCRA0' ms' as' ds' -- ))
where ms' = take (fromInteger n) ms
as' = take (fromInteger n) as
ds' = take (fromInteger n) ds
handle_fails :: Integer -> Integer -> Integer ->
[Integer] -> [Integer] -> [Integer] -> (Integer, Integer)
handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
-- parGlobal 46 46 1 0 m $
-- parGlobal 47 47 1 0 a $
(if n == 0
then (m, a)
else if d1 == 0
then (handle_fails n m a ms as ds)
else handle_fails (n-1) m' a' ms as ds
)
#if defined(STRATEGIES)
`using` \ x -> do par m (pseq a (return x)) --{ m' <- rpar m; a' <- rseq a; return x }
#endif
where
m' = m * m1
a' = {-par_binCRA-} binCRA m m1 inv a a1
inv = modInv m1 m
in
handle_fails fails m a ms as ds
tree_IMCRA0' [m] [a] [0] = (1, 1, 1) -- FAIL due to unlucky prime
tree_IMCRA0' [m] [a] [_] = (m, a, 0) -- normal case
tree_IMCRA0' ms as ds =
let
n = length ms
(left_ms, right_ms) = splitAt (n `div` 2) ms
(left_as, right_as) = splitAt (n `div` 2) as
(left_ds, right_ds) = splitAt (n `div` 2) ds
left@(left_M, left_CRA, left_fails) =
tree_IMCRA0' left_ms left_as left_ds
right@(right_M, right_CRA, right_fails) =
tree_IMCRA0' right_ms right_as right_ds
inv = modInv right_M left_M
cra = par_binCRA left_M right_M inv left_CRA right_CRA
in
#if defined(TU_CRA)
_parGlobal_ 5# left
(_parGlobal_ 6# right
(_parLocal_ 9# inv
--seq cra
(left_M * right_M,
cra,
left_fails + right_fails ) -- ))
#elif defined(COMP_CRA)
_parGlobal_ 5# left_CRA
(_parGlobal_ 6# left_M
(_parGlobal_ 7# right_CRA
(_parGlobal_ 8# right_M
(_parGlobal_ 9# inv
(left_M * right_M,
cra,
left_fails + right_fails)))))
#elif defined(ASYM_CRA)
_parGlobal_ 5# left
(_parGlobal_ 7# right
(_parGlobal_ 9# inv
(seq cra
(left_M * right_M,
cra,
left_fails + right_fails)) ) )
#endif
-- +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-- Experimental pat starts here
-- +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-- This is a brain damaged version of doing a binary-tree CRA over lists.
tree_IMCRA :: [Integer] -> [Integer] -> (Integer,Integer)
tree_IMCRA ms as = _parGlobal_ 0# inv
(m1*m2, par_binCRA m1 m2 inv a1 a2)
where --( [(m1,a1),(m2,a2)] : _ ) =
--(m1,a1) = (head foo) !! 0
--(m2,a2) = (head foo) !! 1
-- foo :: Integer -- (Integer,Integer) -> (Integer,Integer) -> Boolean
[(m1, a1), (m2, a2)] =
head (dropWhile (\ l -> (length l) /= 2)
(iterate tree_IMCRA' (zip ms as)))
inv = modInv m2 m1
-- cra = par_binCRA m1 m2 inv a1 a2
tree_IMCRA' :: [(Integer,Integer)] -> [(Integer,Integer)]
tree_IMCRA' [] = []
tree_IMCRA' [(m1,a1)] = [(m1,a1)]
tree_IMCRA' ((m1,a1) : (m2,a2) : rest) =
_parGlobal_ 11# cra ({-_parGlobal_ 12# (forcelist rest')-} (cra : rest') )
where cra = (m1*m2, binCRA m1 m2 (modInv m2 m1) a1 a2)
rest' = tree_IMCRA' rest
let
n = length ms
(left_ms, right_ms) = splitAt (n `div` 2) ms
(left_as, right_as) = splitAt (n `div` 2) as
(left_ds, right_ds) = splitAt (n `div` 2) ds
left@(left_M, left_CRA, left_fails) =
tree_IMCRA0' left_ms left_as left_ds
right@(right_M, right_CRA, right_fails) =
tree_IMCRA0' right_ms right_as right_ds
inv = modInv right_M left_M
cra = {-par_binCRA-} binCRA left_M right_M inv left_CRA right_CRA
in
(left_M * right_M, cra, left_fails + right_fails )
-- IMCRA with list of products of the modules as additional argument i.e.
-- in s2IMCRA m M r the following holds M_i = m_1 * ... * m_{i-1}
......@@ -251,72 +182,3 @@ s3IMCRA (m:ms) pp (r:rs) = foldl
(\ (m0,r0) (m1,r1,p) -> (p,binCRA m0 m1 (modInv m1 m0) r0 r1) )
(m,r) (zip3 ms rs pp)
-----------------------------------------------------------------------------
-- Old cheating version
-----------------------------------------------------------------------------
#if 0 /* defined(NO_FILTER) */
-- Note: ms as are finte lists
-- no ds needed as no test (== 0) is performed
tree_IMCRA0 :: [Integer] -> [Integer] ->
(Integer,Integer)
-- tree_IMCRA0 [] [] = (1, 1)
tree_IMCRA0 [m] [a] = (m, a)
tree_IMCRA0 ms as =
let
n = length ms
(left_ms, right_ms) = splitAt (n `div` 2) ms
(left_as, right_as) = splitAt (n `div` 2) as
left@(left_M, left_CRA) = tree_IMCRA0 left_ms left_as
right@(right_M, right_CRA) = tree_IMCRA0 right_ms right_as
inv = modInv right_M left_M
in
#if defined(TU_CRA)
_parGlobal_ 5# left
(_parGlobal_ 6# right
(_parLocal_ 9# inv
(left_M * right_M,
par_binCRA left_M right_M inv left_CRA right_CRA)))
#elif defined(COMP_CRA)
_parGlobal_ 5# left_CRA
(_parGlobal_ 6# left_M
(_parGlobal_ 7# right_CRA
(_parGlobal_ 8# right_M
(_parGlobal_ 9# inv
(left_M * right_M,
par_binCRA left_M right_M inv left_CRA right_CRA)))))
#elif defined(ASYM_CRA)
_parGlobal_ 5# left
--(_parLocal_ 6# right
(_parGlobal_ 9# inv
(left_M * right_M,
par_binCRA left_M right_M inv left_CRA right_CRA)) --)
#endif
{-
--(left_CRA + right_CRA,
--binCRA left_M right_M inv left_CRA right_CRA) )
-}
{- Test check for equal length lists.
| length ms /= length as = error "tree_IMCRA0: different lengths of input lists"
-}
#endif
-- Old fail handler
{-
handle_fails 0 m a _ _ _ = (m, a)
handle_fails n m a (m1:ms) (a1:as) (0:ds) =
handle_fails n m a ms as ds
handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
handle_fails (n-1) m' a' ms as ds
where
m' = m * m1
a' = par_binCRA m m1 inv a a1
inv = modInv m1 m
-}
{-
$Id: LinSolv.hs,v 1.1 1997/07/31 22:15:08 sof Exp $
This is revision: $Revision: 1.1 $
Module for computing the solution of a system of linear system of equations
in several variables.
The rhs of the system of equations is the first parameter (given as a
square matrix). The lhs of the system is the second parameter (given as a
vector).
ChangeLog:
$Log: LinSolv.hs,v $
Revision 1.1 1997/07/31 22:15:08 sof
Renamed LinSolv-par to LinSolv
Revision 1.1 1996/01/08 20:07:34 partain
Initial revision
--# Revision 1.4 1994/12/13 19:49:38 hwloidl
--# Version based on infinite lists.
--# Has seq steps in all variants (even without det check)
--# -DNOFILTER uses finite part of inifinte xList and avoids any det check
--# -DSEQ_CRA uses a fold-based CRA instead of the tree-based one (both over
--# infinte lists)
--#
--# Revision 1.3 1994/11/23 01:04:43 hwloidl
--# Test version of fwd mapping and hom sol (no CRA in this version)
--# New: prime list computed in par and packet forced
--# modDet computed in par to hom sol
--# No modDet == 0 test in gen_xList to avoid seq in that step ->
--# count and filter are used to detect number of necessary hom ims
--# xList is now again an infinite list (it is not forced in this version)
--# Note: The result of linsolv just triggers computation up to hom sols but
--# not until CRAs. The value is not important, just the tiggering
--# NO_FILTER flag: if set any modDet == 0 test is avoided (i.e. no filtering
--# either; this doesn't always yield a correct solution but
--# shows how this test seq the computation
--#
--# Revision 1.2 1994/11/19 21:49:38 hwloidl
--# Compute list of primes in parallel to fwd mapping and sol in hom ims
--#
--# Revision 1.1 1994/11/19 02:00:05 hwloidl
--# Initial revision
--#
---------------------------------------------------------------------------- -}
module LinSolv(linSolv) where
import ModArithm {- (modHom, Hom(hom))-}
import Matrix (SqMatrix, Vector, vector, matBounds,
determinant, transp, replaceColumn, size,
maxElem, maxElemVec, scalarMult, vecScalarQuot,
matGcd, vecGcd, matHom, vecHom)
import CRA (s1IMCRA,s2IMCRA,tree_IMCRA0)
#ifndef SEQ
import ParForce
import {-fool mkdependHS; ToDo: rm-}
Parallel
#else
parmap = map
par :: a -> b -> b
par x y = y
seq :: a -> b -> b
seq x y = y
_parGlobal_ :: Int# -> a -> b -> b
_parGlobal_ n x y = y
_parLocal_ :: Int# -> a -> b -> b
_parLocal_ n x y = y
#endif
#ifdef GUM
-- GUM has only par
_parGlobal_ :: Int# -> a -> b -> b
_parGlobal_ n x y = par x y
_parLocal_ :: Int# -> a -> b -> b
_parLocal_ n x y = par x y
#endif
-- ---------------------------------------------------------------------------
-- Auxiliary functions
-- ---------------------------------------------------------------------------
primes :: [Integer]
primes = 2 : oddPrimesFrom 3
where oddPrimesFrom n | isPrime n primes = n : oddPrimesFrom (n+2)
| otherwise = oddPrimesFrom (n+2)
isPrime :: Integer -> [Integer] -> Bool
isPrime n (l1:ls) | n < l1*l1 = True
| n `mod` l1 == 0 = False
| otherwise = isPrime n ls
projection :: Integer -> [[Integer]] -> [Integer]
projection i [] = []
projection i (ls1:lss) = ls1!!(i-1) : projection i lss
fact :: Integer -> Integer
fact 0 = 1
fact (n+1) = (n+1) * fact n
-- ---------------------------------------------------------------------------
-- Function for solving a system of linear equations
-- ---------------------------------------------------------------------------
-- linSolv :: (Integral a) => SqMatrix a -> Vector a -> (a,a,Vector a)
linSolv a b =
let
{- Step1: -}
n = toInteger (size a)
s = max (maxElem a) (maxElemVec b)
pBound = 2 * s^n * fact n
{- Step2: -}
xList :: [[Integer]]
xList = _parGlobal_ 19# (forcelist1 100 primes)
_parGlobal_ 18# noOfPrimes
(gen_xList primes)
gen_xList :: [Integer] -> [[Integer]]
gen_xList (p:ps)
= let
b0 = vecHom p b
a0 = matHom p a
modDet = modHom p (determinant a0)
{- Parallelism inside a homomorphic solution -}
homSol = _parGlobal_ 8# modDet
(p : modDet : pmx)
pmx = parmap ( \ j ->
let
a1 = replaceColumn j a0 b0
in
modHom p (determinant a1) )
[jLo..jHi]
((iLo,jLo),(iHi,jHi)) = matBounds a
restList = gen_xList ps
in
_parGlobal_ 3# homSol
(homSol : restList)
{- Step3: -}
{- s1IMCRA variant -}
count (l1:ls) i accum = if accum>pBound
then i
else count ls (i+1) (l1*accum)
-- noOfPrimes = count (projection 1 xList) 0 1
-- Just for TESTING:
-- noOfPrimes_0 = count primes 0 1
-- xList_0 = take noOfPrimes xList
#if defined(NO_FILTER)
{-
noOfPrimes = count (projection 1 xList) 0 1
xList'' = take noOfPrimes xList
luckyPrimes :: [Integer]
luckyPrimes = projection 1 xList''
-}
-- 1. The lists passed to CRA are finite
-- 2. No test for det == 0 is ever performed
noOfPrimes = count (projection 1 xList) 0 1
primeList = take noOfPrimes (projection 1 xList)
detList = take noOfPrimes (projection 2 xList)
-- Force the lists?
det = snd (tree_IMCRA0 primeList detList)
x_i i = snd (tree_IMCRA0 primeList x_i_List)
where
x_i_List = take noOfPrimes (projection (i+2) xList)
x = vector (parmap x_i [1..n])
#else
-- No more cheating ...
{- The naive, inefficient version
xList' = filter (\ (x1:x2:xs) -> x2 /= 0) xList
noOfPrimes = count (projection 1 xList') 0 1
xList'' = take noOfPrimes xList'
-}
noOfPrimes = count (projection 1 xList) 0 1
primeList = projection 1 xList
detList = projection 2 xList -- list of dets in hom ims if det==0
-- we will throw the result away in
-- tree_IMCRA0 => needed in there
#if defined(SEQ_CRA)
det = a
where (_, a) = s2IMCRA pBound {-noOfPrimes-} primes detList detList
x_i i = a
where (_, a) = s2IMCRA pBound {-noOfPrimes-} primes x_i_List detList