<commit_before>
module Main where

import Numeric.AERN.DoubleBasis.RealIntervalApprox
import Numeric.AERN.DoubleBasis.MRealIntervalApprox

import Control.Monad.ST (ST, runST)

type R = RealIntervalApprox
type RI = RealIntervalApprox

main =
  do
  putStrLn $ "riemann 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemann 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erf e x = 2/(sqrt pi) * riemann e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erf 0.001 1 = " ++ show (erf 0.001 1) 
  putStrLn $ "riemannInPlace 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemannInPlace 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erfInPlace e x = 2/(sqrt pi) * riemannInPlace e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erfInPlace 0.001 1 = " ++ show (erfInPlace 0.001 1) 
  putStrLn $ "riemannInPlace 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemannInPlace 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erfInPlace e x = 2/(sqrt pi) * riemannInPlace e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erfInPlace 0.001 1 = " ++ show (erf 0.001 1) 

-- compute the error function to within e accuracy
erf :: R -> R -> R
erf e x =
  2/(sqrt pi) * riemann e (\t -> exp (-t^2)) (0 </\> x)

-- compute the integral of f over d to within e accuracy
riemann :: R -> (R -> R) -> RI -> R
riemann e f d =
  riemann' e f (width d) [d] 0

riemann' _ _ _ []     result = result 
riemann' e f initWidth (d:ds) result
  | reachedSufficientAccuracy =
      riemann' e f initWidth ds (dArea + result)
  | otherwise = 
      riemann' e f initWidth (dl:dr:ds) result
  where
  reachedSufficientAccuracy =
    case dError <? e * dWidth / initWidth of
      Just True -> True
      _ -> False
  dError = dWidth * (width fd)
  dArea = dWidth * fd
  fd = f d
  dWidth = width d
  (dl, dr) = bisect Nothing d

-- erf using riemannInPlace
erfInPlace :: R -> R -> R
erfInPlace e x =
  2/(sqrt pi) * riemannInPlace e (\t -> exp (-t^2)) (0 </\> x)

-- riemann using in-place accumulator
riemannInPlace :: R -> (R -> R) -> RI -> R
riemannInPlace e f d =
  runST $
    do
    resM <- makeMutable 0
    riemannInPlace' resM e f (width d) [d]
    result <- unsafeReadMutable resM
    return result

riemannInPlace' resM _ _ _ [] =
  return () 
riemannInPlace' resM e f initWidth (d:ds)
  | reachedSufficientAccuracy =
      do 
      dAreaM <- unsafeMakeMutable dArea
      resM <+>= dAreaM
      riemannInPlace' resM e f initWidth ds 
  | otherwise = 
      riemannInPlace' resM e f initWidth (dl:dr:ds)
  where
  reachedSufficientAccuracy =
    case dError <? e * dWidth / initWidth of
      Just True -> True
      _ -> False
  dError = dWidth * (width fd)
  dArea = dWidth * fd
  fd = f d
  dWidth = width d<commit_msg>Fixed old Riemann demo so that it compiles with recent versions.<commit_after>
module Main where

import Numeric.AERN.DoubleBasis.RealIntervalApprox
import Numeric.AERN.DoubleBasis.MRealIntervalApprox

import Control.Monad.ST (ST, runST)

type R = RealIntervalApprox
type RI = RealIntervalApprox

main =
  do
  putStrLn $ "riemann 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemann 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erf e x = 2/(sqrt pi) * riemann e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erf 0.001 1 = " ++ show (erf 0.001 1) 
  putStrLn $ "riemannInPlace 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemannInPlace 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erfInPlace e x = 2/(sqrt pi) * riemannInPlace e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erfInPlace 0.001 1 = " ++ show (erfInPlace 0.001 1) 
  putStrLn $ "riemannInPlace 0.1 (\\x -> x^2) (Interval 0 1) = " ++
    show (riemannInPlace 0.1 (\x -> x^2) (0 </\> 1))
  putStrLn "erfInPlace e x = 2/(sqrt pi) * riemannInPlace e (\\t -> exp (-t^2)) (0 </\\> x)"
  putStrLn $ "erfInPlace 0.001 1 = " ++ show (erf 0.001 1) 

-- compute the error function to within e accuracy
erf :: R -> R -> R
erf e x =
  2/(sqrt piOut) * riemann e (\t -> exp (-t^2)) (0 </\> x)

-- compute the integral of f over d to within e accuracy
riemann :: R -> (R -> R) -> RI -> R
riemann e f d =
  riemann' e f (width d) [d] 0

riemann' _ _ _ []     result = result 
riemann' e f initWidth (d:ds) result
  | reachedSufficientAccuracy =
      riemann' e f initWidth ds (dArea + result)
  | otherwise = 
      riemann' e f initWidth (dl:dr:ds) result
  where
  reachedSufficientAccuracy =
    case dError <? e * dWidth / initWidth of
      Just True -> True
      _ -> False
  dError = dWidth * (width fd)
  dArea = dWidth * fd
  fd = f d
  dWidth = width d
  (dl, dr) = bisect Nothing d

-- erf using riemannInPlace
erfInPlace :: R -> R -> R
erfInPlace e x =
  2/(sqrt piOut) * riemannInPlace e (\t -> exp (-t^2)) (0 </\> x)

-- riemann using in-place accumulator
riemannInPlace :: R -> (R -> R) -> RI -> R
riemannInPlace e f d =
  runST $
    do
    resM <- makeMutable 0
    riemannInPlace' resM e f (width d) [d]
    result <- unsafeReadMutable resM
    return result

riemannInPlace' resM _ _ _ [] =
  return () 
riemannInPlace' resM e f initWidth (d:ds)
  | reachedSufficientAccuracy =
      do 
      dAreaM <- unsafeMakeMutable dArea
      resM <+>= dAreaM
      riemannInPlace' resM e f initWidth ds 
  | otherwise = 
      riemannInPlace' resM e f initWidth (dl:dr:ds)
  where
  reachedSufficientAccuracy =
    case dError <? e * dWidth / initWidth of
      Just True -> True
      _ -> False
  dError = dWidth * (width fd)
  dArea = dWidth * fd
  fd = f d
  dWidth = width d
{'fraction_non_alphanumeric': 0.14733542319749215, 'fraction_numerical': 0.028397565922920892, 'mean_word_length': 2.5084087968952136, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4199537', 'n_tokens_mistral': 2270, 'n_tokens_neox': 2079, 'n_words': 685}