How to use randomized algorithm within the implementation of pure data structures?

Hi cafe, I'm currently implementing the data structure representing algebraic numbers. Implementing algebraic number arithmetic requires polynomial factorisation. Pure algorithm for factoring (Berlekamp's algorithm) is more expensive than the randomized one (Cantor-Zassenhaus algorithm) , so I want to use the latter algorithm. Since C-Z algorithm is a randomized algorithm, we have to have an access for random number generator when calculating algebraic number arithmetics e.g. writing Num instance for algebraic numbers. Here is the problem: how to pass-around random number generator throughout pure computaion? I think one immediate solution is to create global state with `newStdGen` and `unsafePerformIO` like below: ``` randSrc :: IORef StdGen randSrc = unsafePerformIO $ newIORef =<< newStdGen {-# NOINLINE randSrc #-} getRand :: IO Int getRand = atomicModifyIORef' randSrc (swap . next) -- We can probably use the following function to avoid calling -- `unsafePerformIO` whenever calling 'getRand`, -- assuming some optimization flag enabled: getRandUnsafe :: Int getRandUnsafe = unsafePerformIO getRand {-# NOINLINE getRandUnsafe #-} {-# RULES "getRandUnsafe" getRandUnsafe = unsafePerformIO getRand #-} ``` But this hack seems rather dirty and unsafe. Is there any workaround to achieve the same thing? -- Hiromi ISHII konn.jinro@gmail.com

On 2014-11-01 10:25, Hiromi ISHII wrote:
Hi cafe,
I'm currently implementing the data structure representing algebraic numbers.
Implementing algebraic number arithmetic requires polynomial factorisation. Pure algorithm for factoring (Berlekamp's algorithm) is more expensive than the randomized one (Cantor-Zassenhaus algorithm) , so I want to use the latter algorithm.
Since C-Z algorithm is a randomized algorithm, we have to have an access for random number generator when calculating algebraic number arithmetics e.g. writing Num instance for algebraic numbers.
Here is the problem: how to pass-around random number generator throughout pure computaion?
I think one immediate solution is to create global state with `newStdGen` and `unsafePerformIO` like below:
You can just use the State monad to thread the StdGen around and "update" it when you need to. You can get a pure interface by hiding away the runState behind a function: -- The is the function you export from your module. myAlgorithm :: StdGen -> ... -> ... myAlgorithm g .... = fst . runState (myAlgorithm' ...) g myAlgorithm' :: ... -> State StdGen Result myAlgorithm' ... = do ... x <- rand ... return $ ... rand :: State StdGen a rand = do x <- get (a, g') <- random -- Here "random" is from System.Random put $! g' return a (The above probably contains typos, and can probably also be prettified, but hopefully you get the gist.) Regards,

Hi Bardur,
You can just use the State monad to thread the StdGen around and "update" it when you need to. You can get a pure interface by hiding away the runState behind a function:
Thank you for your rapid response! Unfortunately, I didn't describe my problem accurately. This approach (or using MonadRandom) to pass around random generator with Monad, works fine when it's just enough to feed generator to the algorithm. But my situation is slightly different: random generator has to be passed around to implement the instance method for `Num`, so it can't take random generator as its argument. So I need some way to hide random generator from function type signatures. Fortunately, your response suggested me the alternative approach: converting the data-type into continuation-passing style. This should work fine when we just do some operations on data-type, but we have to feed the generator when we want to inspect its value, so it's not sufficient, though... -- Hiromi ISHII konn.jinro@gmail.com

Just an idea here, but would implicit-params work? It only gives you
Reader-monad capabilities, but you can always split random generators.
There might be repercussions for the quality of the generated numbers,
though, for which I have no idea.
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
import Data.Implicit
import System.Random
import Text.Printf
import Data.Default.Class
instance Default StdGen where
def = mkStdGen 10
newtype AlgebraicNumber = AlgebraicNumber {- your data here -} String
deriving Show
instance Implicit_ StdGen => Num AlgebraicNumber where
AlgebraicNumber x + AlgebraicNumber y =
-- You could tidy this up with a State monad.
let g = param_ :: StdGen
(g1, g2) = split g
gx = fst $ next g1 -- compute on x using left generator
gy = fst $ next g2 -- compute on y using right generator
in
AlgebraicNumber (printf ("%s computed with rand = %d,"
++ "%s computed with rand = %d")
x gx y gy)
On Sat, Nov 1, 2014 at 11:28 AM, Hiromi ISHII
Hi Bardur,
You can just use the State monad to thread the StdGen around and "update" it when you need to. You can get a pure interface by hiding away the runState behind a function:
Thank you for your rapid response! Unfortunately, I didn't describe my problem accurately.
This approach (or using MonadRandom) to pass around random generator with Monad, works fine when it's just enough to feed generator to the algorithm.
But my situation is slightly different: random generator has to be passed around to implement the instance method for `Num`, so it can't take random generator as its argument. So I need some way to hide random generator from function type signatures.
Fortunately, your response suggested me the alternative approach: converting the data-type into continuation-passing style. This should work fine when we just do some operations on data-type, but we have to feed the generator when we want to inspect its value, so it's not sufficient, though...
-- Hiromi ISHII konn.jinro@gmail.com
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
-- Jun Inoue

Hi Ishii-san, お久しぶりです。 On 2014年11月01日 18:25, Hiromi ISHII wrote:
Since C-Z algorithm is a randomized algorithm, we have to have an access for random number generator when calculating algebraic number arithmetics e.g. writing Num instance for algebraic numbers.
Here is the problem: how to pass-around random number generator throughout pure computaion?
I think one immediate solution is to create global state with `newStdGen` and `unsafePerformIO` like below: <SNIP> But this hack seems rather dirty and unsafe.
Is there any workaround to achieve the same thing?
If it works with the algorithm, you could use a pseudo-random number generator with a fixed seed. For example, here is a program to estimate the value of π (purely) using a Monte Carlo simulation: {-# LANGUAGE BangPatterns #-} module Main where import System.Random (mkStdGen, randomRs) -- | Estamate pi via monte-carlo simulation mcpi :: Int -- ^ number of iterations -> Double -- ^ estimated value of pi mcpi count = step (randomRs (0.0, 1.0) (mkStdGen 1331)) 0 count where step :: [Double] -> Int -> Int -> Double step (x:y:rs) !qrt !i | i < 1 = 4.0 * fromIntegral qrt / fromIntegral count | hit x y = step rs (qrt + 1) (i - 1) | otherwise = step rs qrt (i - 1) step _ _ _ = error "impossible" hit :: Double -> Double -> Bool hit x y = x ^ (2 :: Int) + y ^ (2 :: Int) <= 1.0 main :: IO () main = putStrLn $ "pi ~= " ++ show (mcpi 1000000) Cheers, Travis
participants (4)
-
Bardur Arantsson
-
Hiromi ISHII
-
Jun Inoue
-
Travis Cardwell