summaryrefslogtreecommitdiff
path: root/testsuite/tests/dph/nbody/Randomish.hs
blob: 7aeefa33d84e35c6e5a8ad6b2e567cad15b849bb (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
{-# LANGUAGE BangPatterns #-}

module Randomish
        ( randomishInts
        , randomishDoubles)
where
import Data.Word
import Data.Vector.Unboxed                      (Vector)
import qualified Data.Vector.Unboxed.Mutable    as MV
import qualified Data.Vector.Unboxed            as V
import qualified Data.Vector.Generic            as G            


-- | Use the "minimal standard" Lehmer generator to quickly generate some random
--   numbers with reasonable statistical properties. By "reasonable" we mean good
--   enough for games and test data, but not cryptography or anything where the
--   quality of the randomness really matters.
-- 
--   From "Random Number Generators: Good ones are hard to find"
--   Stephen K. Park and Keith W. Miller.
--   Communications of the ACM, Oct 1988, Volume 31, Number 10.
--
randomishInts 
        :: Int                  -- Length of vector.
        -> Int                  -- Minumum value in output.
        -> Int                  -- Maximum value in output.
        -> Int                  -- Random seed. 
        -> Vector Int           -- Vector of random numbers.

randomishInts !len !valMin' !valMax' !seed'
        
 = let  -- a magic number (don't change it)
        multiplier :: Word64
        multiplier = 16807

        -- a merzenne prime (don't change it)
        modulus :: Word64
        modulus = 2^(31 :: Integer) - 1

        -- if the seed is 0 all the numbers in the sequence are also 0.
        seed    
         | seed' == 0   = 1
         | otherwise    = seed'

        !valMin = fromIntegral valMin'
        !valMax = fromIntegral valMax' + 1
        !range  = valMax - valMin

        {-# INLINE f #-}
        f x             = multiplier * x `mod` modulus
 in G.create 
     $ do       
        vec             <- MV.new len

        let go !ix !x 
                | ix == len     = return ()
                | otherwise
                = do    let x'  = f x
                        MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin
                        go (ix + 1) x'

        go 0 (f $ f $ f $ fromIntegral seed)
        return vec


-- | Generate some randomish doubles with terrible statistical properties.
--   This is good enough for test data, but not much else.
randomishDoubles 
        :: Int                  -- Length of vector
        -> Double               -- Minimum value in output
        -> Double               -- Maximum value in output
        -> Int                  -- Random seed.
        -> Vector Double        -- Vector of randomish doubles.

randomishDoubles !len !valMin !valMax !seed
 = let  range   = valMax - valMin

        mx      = 2^(30 :: Integer) - 1
        mxf     = fromIntegral mx
        ints    = randomishInts len 0 mx seed
        
   in   V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints