{-# LANGUAGE TypeFamilies #-}

-- |
-- Module     : Simulation.Aivika.RealTime.Generator
-- Copyright  : Copyright (c) 2016-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 8.0.1
--
-- Here is defined a random number generator, where
-- the 'RT' monad can be an instance of 'MonadGenerator'.
--
module Simulation.Aivika.RealTime.Generator () where

import Control.Monad
import Control.Monad.Trans

import System.Random
import qualified System.Random.MWC as MWC

import Data.IORef
import Data.Vector

import Simulation.Aivika.Trans.Generator
import Simulation.Aivika.Trans.Generator.Primitive

import Simulation.Aivika.RealTime.Internal.RT

instance (Functor m, Monad m, MonadIO m) => MonadGenerator (RT m) where

  {-# SPECIALISE instance MonadGenerator (RT IO) #-}

  data Generator (RT m) =
    Generator { Generator (RT m) -> RT m Double
generator01 :: RT m Double,
                -- ^ the generator of uniform numbers from 0 to 1
                Generator (RT m) -> RT m Double
generatorNormal01 :: RT m Double,
                -- ^ the generator of normal numbers with mean 0 and variance 1
                Generator (RT m) -> RT m Int
generatorSequenceNo :: RT m Int
                -- ^ the generator of sequence numbers
              }

  {-# INLINE generateUniform #-}
  generateUniform :: Generator (RT m) -> Double -> Double -> RT m Double
generateUniform = RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateUniform01 (RT m Double -> Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateUniformInt #-}
  generateUniformInt :: Generator (RT m) -> Int -> Int -> RT m Int
generateUniformInt = RT m Double -> Int -> Int -> RT m Int
forall (m :: * -> *). Monad m => m Double -> Int -> Int -> m Int
generateUniformInt01 (RT m Double -> Int -> Int -> RT m Int)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Int
-> Int
-> RT m Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateTriangular #-}
  generateTriangular :: Generator (RT m) -> Double -> Double -> Double -> RT m Double
generateTriangular = RT m Double -> Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> Double -> m Double
generateTriangular01 (RT m Double -> Double -> Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Double
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateNormal #-}
  generateNormal :: Generator (RT m) -> Double -> Double -> RT m Double
generateNormal = RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateNormal01 (RT m Double -> Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generatorNormal01

  {-# INLINE generateLogNormal #-}
  generateLogNormal :: Generator (RT m) -> Double -> Double -> RT m Double
generateLogNormal = RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateLogNormal01 (RT m Double -> Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generatorNormal01

  {-# INLINE generateExponential #-}
  generateExponential :: Generator (RT m) -> Double -> RT m Double
generateExponential = RT m Double -> Double -> RT m Double
forall (m :: * -> *). Monad m => m Double -> Double -> m Double
generateExponential01 (RT m Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateErlang #-}
  generateErlang :: Generator (RT m) -> Double -> Int -> RT m Double
generateErlang = RT m Double -> Double -> Int -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Int -> m Double
generateErlang01 (RT m Double -> Double -> Int -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Int
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generatePoisson #-}
  generatePoisson :: Generator (RT m) -> Double -> RT m Int
generatePoisson = RT m Double -> Double -> RT m Int
forall (m :: * -> *). Monad m => m Double -> Double -> m Int
generatePoisson01 (RT m Double -> Double -> RT m Int)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> RT m Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateBinomial #-}
  generateBinomial :: Generator (RT m) -> Double -> Int -> RT m Int
generateBinomial = RT m Double -> Double -> Int -> RT m Int
forall (m :: * -> *). Monad m => m Double -> Double -> Int -> m Int
generateBinomial01 (RT m Double -> Double -> Int -> RT m Int)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Int
-> RT m Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateGamma #-}
  generateGamma :: Generator (RT m) -> Double -> Double -> RT m Double
generateGamma Generator (RT m)
g = RT m Double -> RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateGamma01 (Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generatorNormal01 Generator (RT m)
g) (Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01 Generator (RT m)
g)

  {-# INLINE generateBeta #-}
  generateBeta :: Generator (RT m) -> Double -> Double -> RT m Double
generateBeta Generator (RT m)
g = RT m Double -> RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateBeta01 (Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generatorNormal01 Generator (RT m)
g) (Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01 Generator (RT m)
g)

  {-# INLINE generateWeibull #-}
  generateWeibull :: Generator (RT m) -> Double -> Double -> RT m Double
generateWeibull = RT m Double -> Double -> Double -> RT m Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateWeibull01 (RT m Double -> Double -> Double -> RT m Double)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> Double
-> Double
-> RT m Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateDiscrete #-}
  generateDiscrete :: Generator (RT m) -> DiscretePDF a -> RT m a
generateDiscrete = RT m Double -> DiscretePDF a -> RT m a
forall (m :: * -> *) a. Monad m => m Double -> DiscretePDF a -> m a
generateDiscrete01 (RT m Double -> DiscretePDF a -> RT m a)
-> (Generator (RT m) -> RT m Double)
-> Generator (RT m)
-> DiscretePDF a
-> RT m a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator (RT m) -> RT m Double
forall (m :: * -> *). Generator (RT m) -> RT m Double
generator01

  {-# INLINE generateSequenceNo #-}
  generateSequenceNo :: Generator (RT m) -> RT m Int
generateSequenceNo = Generator (RT m) -> RT m Int
forall (m :: * -> *). Generator (RT m) -> RT m Int
generatorSequenceNo

  {-# INLINABLE newGenerator #-}
  newGenerator :: GeneratorType (RT m) -> RT m (Generator (RT m))
newGenerator GeneratorType (RT m)
tp =
    case GeneratorType (RT m)
tp of
      GeneratorType (RT m)
SimpleGenerator ->
        do let g :: IO (IO Double)
g = Gen (PrimState IO) -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen (PrimState IO) -> IO Double)
-> IO (Gen (PrimState IO)) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   IO (Gen (PrimState IO))
MWC.createSystemRandom
           IO Double
g' <- IO (IO Double) -> RT m (IO Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO (IO Double)
g
           RT m Double -> RT m (Generator (RT m))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 (IO Double -> RT m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO Double
g')
      SimpleGeneratorWithSeed Word32
x ->
        do let g :: IO (IO Double)
g = Gen RealWorld -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen RealWorld -> IO Double)
-> IO (Gen RealWorld) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   Vector Word32 -> IO (Gen (PrimState IO))
forall (m :: * -> *) (v :: * -> *).
(PrimMonad m, Vector v Word32) =>
v Word32 -> m (Gen (PrimState m))
MWC.initialize (Word32 -> Vector Word32
forall a. a -> Vector a
singleton Word32
x)
           IO Double
g' <- IO (IO Double) -> RT m (IO Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO (IO Double)
g
           RT m Double -> RT m (Generator (RT m))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 (IO Double -> RT m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO Double
g')
      CustomGenerator RT m (Generator (RT m))
g ->
        RT m (Generator (RT m))
g
      CustomGenerator01 RT m Double
g ->
        RT m Double -> RT m (Generator (RT m))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 RT m Double
g

  {-# INLINABLE newRandomGenerator #-}
  newRandomGenerator :: g -> RT m (Generator (RT m))
newRandomGenerator g
g = 
    do IORef g
r <- IO (IORef g) -> RT m (IORef g)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef g) -> RT m (IORef g)) -> IO (IORef g) -> RT m (IORef g)
forall a b. (a -> b) -> a -> b
$ g -> IO (IORef g)
forall a. a -> IO (IORef a)
newIORef g
g
       let g01 :: RT m Double
g01 = do g
g <- IO g -> RT m g
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO g -> RT m g) -> IO g -> RT m g
forall a b. (a -> b) -> a -> b
$ IORef g -> IO g
forall a. IORef a -> IO a
readIORef IORef g
r
                    let (Double
x, g
g') = g -> (Double, g)
forall a g. (Random a, RandomGen g) => g -> (a, g)
random g
g
                    IO () -> RT m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> RT m ()) -> IO () -> RT m ()
forall a b. (a -> b) -> a -> b
$ IORef g -> g -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef g
r g
g'
                    Double -> RT m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
x
       RT m Double -> RT m (Generator (RT m))
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 RT m Double
g01

  {-# INLINABLE newRandomGenerator01 #-}
  newRandomGenerator01 :: RT m Double -> RT m (Generator (RT m))
newRandomGenerator01 RT m Double
g01 =
    do RT m Double
gNormal01 <- RT m Double -> RT m (RT m Double)
forall (m :: * -> *). MonadIO m => m Double -> m (m Double)
newNormalGenerator01 RT m Double
g01
       IORef Int
gSeqNoRef <- IO (IORef Int) -> RT m (IORef Int)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Int) -> RT m (IORef Int))
-> IO (IORef Int) -> RT m (IORef Int)
forall a b. (a -> b) -> a -> b
$ Int -> IO (IORef Int)
forall a. a -> IO (IORef a)
newIORef Int
0
       let gSeqNo :: RT m Int
gSeqNo =
             do Int
x <- IO Int -> RT m Int
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Int -> RT m Int) -> IO Int -> RT m Int
forall a b. (a -> b) -> a -> b
$ IORef Int -> IO Int
forall a. IORef a -> IO a
readIORef IORef Int
gSeqNoRef
                IO () -> RT m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> RT m ()) -> IO () -> RT m ()
forall a b. (a -> b) -> a -> b
$ IORef Int -> (Int -> Int) -> IO ()
forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
gSeqNoRef (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                Int -> RT m Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
x
       Generator (RT m) -> RT m (Generator (RT m))
forall (m :: * -> *) a. Monad m => a -> m a
return Generator :: forall (m :: * -> *).
RT m Double -> RT m Double -> RT m Int -> Generator (RT m)
Generator { generator01 :: RT m Double
generator01 = RT m Double
g01,
                          generatorNormal01 :: RT m Double
generatorNormal01 = RT m Double
gNormal01,
                          generatorSequenceNo :: RT m Int
generatorSequenceNo = RT m Int
gSeqNo }

-- | Create a normal random number generator with mean 0 and variance 1
-- by the specified generator of uniform random numbers from 0 to 1.
newNormalGenerator01 :: MonadIO m
                        => m Double
                        -- ^ the generator
                        -> m (m Double)
{-# INLINABLE newNormalGenerator01 #-}
newNormalGenerator01 :: m Double -> m (m Double)
newNormalGenerator01 m Double
g =
  do IORef Double
nextRef <- IO (IORef Double) -> m (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> m (IORef Double))
-> IO (IORef Double) -> m (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Bool
flagRef <- IO (IORef Bool) -> m (IORef Bool)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Bool) -> m (IORef Bool))
-> IO (IORef Bool) -> m (IORef Bool)
forall a b. (a -> b) -> a -> b
$ Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
xi1Ref  <- IO (IORef Double) -> m (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> m (IORef Double))
-> IO (IORef Double) -> m (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
xi2Ref  <- IO (IORef Double) -> m (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> m (IORef Double))
-> IO (IORef Double) -> m (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
psiRef  <- IO (IORef Double) -> m (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> m (IORef Double))
-> IO (IORef Double) -> m (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     let loop :: m ()
loop =
           do Double
psi <- IO Double -> m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> m Double) -> IO Double -> m Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
              if (Double
psi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1.0) Bool -> Bool -> Bool
|| (Double
psi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0)
                then do Double
g1 <- m Double
g
                        Double
g2 <- m Double
g
                        let xi1 :: Double
xi1 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            xi2 :: Double
xi2 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            psi :: Double
psi = Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi2
                        IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
xi1
                        IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
xi2
                        IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
psi
                        m ()
loop
                else IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (- Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
psi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
psi)
     m Double -> m (m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (m Double -> m (m Double)) -> m Double -> m (m Double)
forall a b. (a -> b) -> a -> b
$
       do Bool
flag <- IO Bool -> m Bool
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Bool -> m Bool) -> IO Bool -> m Bool
forall a b. (a -> b) -> a -> b
$ IORef Bool -> IO Bool
forall a. IORef a -> IO a
readIORef IORef Bool
flagRef
          if Bool
flag
            then do IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
False
                    IO Double -> m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> m Double) -> IO Double -> m Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
nextRef
            else do IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
0.0
                    IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
0.0
                    IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
0.0
                    m ()
loop
                    Double
xi1 <- IO Double -> m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> m Double) -> IO Double -> m Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi1Ref
                    Double
xi2 <- IO Double -> m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> m Double) -> IO Double -> m Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi2Ref
                    Double
psi <- IO Double -> m Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> m Double) -> IO Double -> m Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
                    IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
True
                    IO () -> m ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> m ()) -> IO () -> m ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
nextRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi
                    Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi