gauge 0.2.0 → 0.2.1

raw patch · 18 files changed

+1162/−105 lines, 18 filesdep −math-functionsdep −mwc-randomdep ~basedep ~deepseqPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: math-functions, mwc-random

Dependency ranges changed: base, deepseq

API changes (from Hackage documentation)

Files

• Gauge/Analysis.hs +14/−16
• Gauge/Benchmark.hs +17/−21
• Gauge/Format.hs +119/−0
• Gauge/Main.hs +32/−30
• Gauge/Optional.hs +6/−3
• README.markdown +12/−2
• cbits/cycles.h +1/−0
• cbits/time-posix.c +18/−3
• changelog.md +12/−0
• gauge.cabal +11/−13
• math-functions/Numeric/MathFunctions/Comparison.hs +57/−0
• math-functions/Numeric/MathFunctions/Constants.hs +115/−0
• math-functions/Numeric/SpecFunctions.hs +3/−0
• math-functions/Numeric/SpecFunctions/Internal.hs +116/−0
• math-functions/Numeric/Sum.hs +198/−0
• mwc-random/System/Random/MWC.hs +428/−0
• statistics/Statistics/Regression.hs +1/−2
• statistics/Statistics/Resampling.hs +2/−15

@@ -37,20 +37,18 @@     , benchmarkWith'     ) where --- Temporary: to support pre-AMP GHC 7.8.4:-import Data.Monoid- import Control.Applicative import Control.Arrow (second) import Control.DeepSeq (NFData(rnf)) import Control.Monad (forM_, when)-import Gauge.Benchmark (Benchmark (..), Benchmarkable, runBenchmark)+import Gauge.Benchmark (Benchmark (..), BenchmarkAnalysis(..), Benchmarkable, runBenchmark) import Gauge.IO.Printf (note, printError, prolix, rewindClearLine) import Gauge.Main.Options (defaultConfig, Config(..), Verbosity (..),                            DisplayMode (..)) import Gauge.Measurement (Measured(measTime), secs, rescale, measureKeys,                           measureAccessors_, validateAccessors, renderNames) import Gauge.Monad (Gauge, askConfig, gaugeIO, Crit(..), askCrit, withConfig)+import Gauge.Format import qualified Gauge.CSV as CSV import Data.Data (Data, Typeable) import Data.Int (Int64)@@ -106,9 +104,8 @@  instance NFData OutlierEffect -instance Monoid Outliers where-    mempty  = Outliers 0 0 0 0 0-    mappend = addOutliers+outliersEmpty :: Outliers+outliersEmpty = Outliers 0 0 0 0 0  addOutliers :: Outliers -> Outliers -> Outliers addOutliers (Outliers s a b c d) (Outliers t w x y z) =@@ -192,9 +189,10 @@       rnf reportName `seq` rnf reportKeys `seq`       rnf reportMeasured `seq` rnf reportAnalysis `seq` rnf reportOutliers `seq`       rnf reportKDEs+ -- | Classify outliers in a data set, using the boxplot technique. classifyOutliers :: Sample -> Outliers-classifyOutliers sa = U.foldl' ((. outlier) . mappend) mempty ssa+classifyOutliers sa = U.foldl' ((. outlier) . addOutliers) outliersEmpty ssa     where outlier e = Outliers {                         samplesSeen = 1                       , lowSevere = if e <= loS && e < hiM then 1 else 0@@ -412,7 +410,7 @@          case displayMode of             StatsTable -> do-              _ <- note "%sbenchmarked %s\n" rewindClearLine desc+              _ <- note "%sbenchmarked %s%s%s\n" rewindClearLine green desc reset               let r2 n = printf "%.3f R\178" n               forM_ builtin $ \Regression{..} ->                 case Map.lookup "iters" regCoeffs of@@ -436,7 +434,7 @@               _ <- note "\n"               pure ()             Condensed -> do-              _ <- note "%s%-40s " rewindClearLine desc+              _ <- note "%s%s%-40s%s " rewindClearLine green desc reset               bsSmall secs "mean" anMean               bsSmall secs "( +-" anStdDev               _ <- note ")\n"@@ -445,8 +443,9 @@         return rpt       where bs :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge ()             bs f metric e@Estimate{..} =-              note "%-20s %-10s (%s .. %s%s)\n" metric-                   (f estPoint) (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)+              note "%s%-20s%s %s%-10s%s (%s .. %s%s)\n" yellow metric reset+                   red (f estPoint) reset+                   (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)                    (let cl = confIntCL estError                         str | cl == cl95 = ""                             | otherwise  = printf ", ci %.3f" (confidenceLevel cl)@@ -454,7 +453,7 @@                    )             bsSmall :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge ()             bsSmall f metric Estimate{..} =-              note "%s %-10s" metric (f estPoint)+              note "%s %s%-10s%s" metric red (f estPoint) reset              reportStat :: Verbosity                        -> (Measured -> Maybe Double)@@ -487,9 +486,8 @@  -- | Run a benchmark interactively and analyse its performance. benchmarkWith' :: Config -> Benchmarkable -> IO ()-benchmarkWith' cfg bm =-  withConfig cfg $-    runBenchmark (const True) (Benchmark "function" bm) analyseBenchmark+benchmarkWith' cfg bm = withConfig cfg $+    runBenchmark (const True) (Benchmark "function" bm) (BenchmarkNormal analyseBenchmark)  -- | Run a benchmark interactively and analyse its performanc. benchmark' :: Benchmarkable -> IO ()

@@ -59,7 +59,7 @@      -- * Running Benchmarks     , runBenchmark-    , runBenchmarkIters+    , BenchmarkAnalysis(..)     ) where  import Control.Applicative@@ -652,29 +652,25 @@ -- Running benchmarks ------------------------------------------------------------------------------- +-- | The function to run after measurement+data BenchmarkAnalysis =+      forall a . BenchmarkNormal (String -> V.Vector Measured -> Gauge a)+    | BenchmarkIters Int64+ -- | Run benchmarkables, selected by a given selector function, under a given -- benchmark and analyse the output using the given analysis function.-runBenchmark-  :: (String -> Bool) -- ^ Select benchmarks by name.-  -> Benchmark-  -> (String -> V.Vector Measured -> Gauge a) -- ^ Analysis function.-  -> Gauge ()+runBenchmark :: (String -> Bool) -- ^ Select benchmarks by name.+             -> Benchmark+             -> BenchmarkAnalysis -- ^ Analysis function+             -> Gauge () runBenchmark selector bs analyse =-  for selector bs $ \_idx desc bm ->-      runBenchmarkable desc bm >>= analyse desc >>= \_ -> return ()---- XXX For consistency, this should also use a separate process when--- --measure-with is specified.--- | Run a benchmark without analysing its performance.-runBenchmarkIters-  :: (String -> Bool) -- ^ Select benchmarks by name.-  -> Benchmark-  -> Int64            -- ^ Number of iterations to run.-  -> Gauge ()-runBenchmarkIters selector bs iters =-  for selector bs $ \_idx desc bm -> do-    _ <- note "benchmarking %s\r" desc-    gaugeIO $ iterateBenchmarkable_ bm iters+    for selector bs $ \_idx desc bm ->+        case analyse of+            BenchmarkNormal step ->+                runBenchmarkable desc bm >>= step desc >>= \_ -> return ()+            BenchmarkIters iters -> do+                _ <- note "benchmarking %s\r" desc+                gaugeIO $ iterateBenchmarkable_ bm iters  -- | Iterate over benchmarks. for :: (String -> Bool)

@@ -0,0 +1,119 @@+-- |+-- Module      : Gauge.Format+-- Copyright   : (c) 2017 Vincent Hanquez+-- +-- Formatting helpers+--+-- shame there's no leftPad package to use. /s+--+module Gauge.Format+    ( printNanoseconds+    , printSubNanoseconds+    , tableMarkdown+    , reset+    , green+    , red+    , yellow+    ) where++import Gauge.Time+import Data.List+import Data.Word+import Text.Printf+import qualified Basement.Terminal.ANSI as ANSI+import           Basement.Bounded (zn64)+import           GHC.Exts (toList)++-- | Print a NanoSeconds quantity with a human friendly format+-- that make it easy to compare different values+--+-- Given a separator Char of '_':+--+-- 0           -> "             0"+-- 1000        -> "         1_000"+-- 1234567     -> "     1_234_567"+-- 10200300400 -> "10_200_300_400"+--+-- Note that the seconds parameters is aligned considered+-- maximum of 2 characters (i.e. 99 seconds).+-- +printNanoseconds :: Maybe Char -> NanoSeconds -> String+printNanoseconds thousandSeparator (NanoSeconds absNs) =+    case divSub1000 0 absNs of+        [ns]         -> padLeft maxLength $ printSpace ns+        [ns,us]      -> padLeft maxLength $ addSeparators1000 [printSpace us,print3 ns]+        [ns,us,ms]   -> padLeft maxLength $ addSeparators1000 [printSpace ms,print3 us,print3 ns]+        [ns,us,ms,s] -> padLeft maxLength $ addSeparators1000 [printSpace s,print3 ms,print3 us,print3 ns]+        _            -> error "printNanoSeconds: internal error: invalid format"+  where+    maxLength = 3 + 3 + 3 + 2 + (sepLength * 3)+  +    (addSeparators1000, sepLength) =+        case thousandSeparator of+            Nothing -> (concat, 0)+            Just c  -> (intercalate [c], 1)++    printSpace :: Word64 -> String+    printSpace n = printf "%3d" n+    print3 :: Word64 -> String+    print3 n = printf "%03d" n++    divSub1000 :: Int -> Word64 -> [Word64]+    divSub1000 n i+        | n == 3    = [i]+        | otherwise =+            let (d,m) = i `divMod` 1000+             in if d == 0 then [m] else m : divSub1000 (n+1) d++printSubNanoseconds :: Maybe Char -> PicoSeconds100 -> String+printSubNanoseconds ts p =+    printNanoseconds ts ns ++ "." ++ show fragment+  where+    (ns, fragment) = picosecondsToNanoSeconds p+++-- | Produce a table in markdown+--+-- This is handy when wanting to copy paste to a markdown flavor destination.+tableMarkdown :: String     -- ^ top left corner label+              -> [String]   -- ^ columns labels+              -> [[String]] -- ^ a list of row labels followed by content rows+              -> String     -- ^ the resulting string+tableMarkdown name cols rows =+    let hdr = "| " ++ intercalate " | " (padList (name : cols)) ++ " |\n"+        sep = "|-" ++ intercalate "-|-" (map (map (const '-')) (padList (name : cols))) ++ "-|\n"+     in hdr ++ sep ++ concatMap printRow (map padList rows)+  where+    printRow :: [String] -> String+    printRow l = "| " ++ intercalate " | " l ++ " |\n"++    getColN n = map (flip (!!) n) rows++    sizeCols :: [Int]+    sizeCols = map (\(i, c) -> maximum $ map length (c : getColN i)) $ zip [0..] (name : cols)++    padList l = zipWith padCenter sizeCols l++padLeft :: Int -> String -> String+padLeft sz s+    | sz <= len = s+    | otherwise = replicate leftPad ' ' ++ s+  where+    len = length s+    leftPad = (sz - len)++padCenter :: Int -> String -> String+padCenter sz s+    | sz <= len = s+    | otherwise = replicate leftPad ' ' ++ s ++ replicate rightPad ' '+  where+    len = length s+    (leftPad, r) = (sz - len) `divMod` 2+    rightPad = leftPad + r++-- | reset, green, red, yellow ANSI escape+reset, green, red, yellow :: String+reset = toList ANSI.sgrReset+green = toList $ ANSI.sgrForeground (zn64 2) True+red = toList $ ANSI.sgrForeground (zn64 1) True+yellow = toList $ ANSI.sgrForeground (zn64 3) True

@@ -23,6 +23,7 @@     -- * Running Benchmarks Interactively     , benchmark     , benchmarkWith+    , module Gauge.Benchmark     ) where  import Control.Applicative@@ -109,7 +110,7 @@ benchmarkWith :: Config -> Benchmarkable -> IO () benchmarkWith cfg bm =   withConfig cfg $-    runBenchmark (const True) (Benchmark "function" bm) quickAnalyse+    runBenchmark (const True) (Benchmark "function" bm) (BenchmarkNormal quickAnalyse)  -- | Run a benchmark interactively with default config, and analyse its -- performance.@@ -160,35 +161,36 @@ -- one in your benchmark driver's command-line parser). runMode :: Mode -> Config -> [String] -> [Benchmark] -> IO () runMode wat cfg benches bs =-  -- TBD: This has become messy. We use mode as well as cfg options for the-  -- same purpose It is possible to specify multiple exclusive options.  We-  -- need to handle the exclusive options in a better way.-  case wat of-    List    -> mapM_ putStrLn . sort . concatMap benchNames $ bs-    Version -> putStrLn versionInfo-    Help    -> putStrLn describe-    DefaultMode ->-      case measureOnly cfg of-        Just outfile -> runWithConfig runBenchmark (\_ r ->-                          gaugeIO (writeFile outfile (show r)))-        Nothing ->-          case iters cfg of-          Just nbIters -> runWithConfig runBenchmarkIters nbIters-          Nothing ->-            case quickMode cfg of-              True  -> runWithConfig runBenchmark quickAnalyse-              False -> do+    -- TBD: This has become messy. We use mode as well as cfg options for the+    -- same purpose It is possible to specify multiple exclusive options.  We+    -- need to handle the exclusive options in a better way.+    case wat of+        List        -> mapM_ putStrLn . sort . concatMap benchNames $ bs+        Version     -> putStrLn versionInfo+        Help        -> putStrLn describe+        DefaultMode -> runDefault+  where+    runDefault = do+        CSV.write (csvRawFile cfg) $ CSV.Row $ map (CSV.string . fst) measureAccessors_+        CSV.write (csvFile cfg) $ CSV.Row $ map CSV.string+            ["Name", "Mean","MeanLB","MeanUB","Stddev","StddevLB","StddevUB"]++        hSetBuffering stdout NoBuffering+        selector <- selectBenches (match cfg) benches bsgroup++        -- if compiled without analysis step, then default to quickmode #ifdef HAVE_ANALYSIS-                  CSV.write (csvRawFile cfg) $ CSV.Row $ map CSV.string-                        (map fst measureAccessors_)-                  CSV.write (csvFile cfg) $ CSV.Row $ map CSV.string-                        ["Name", "Mean","MeanLB","MeanUB","Stddev","StddevLB","StddevUB"]-                  runWithConfig runBenchmark analyseBenchmark+        let compiledAnalyseStep = analyseBenchmark #else-                  runWithConfig runBenchmark quickAnalyse+        let compiledAnalyseStep = quickAnalyse #endif-  where bsgroup = BenchGroup "" bs-        runWithConfig f arg = do-          hSetBuffering stdout NoBuffering-          selector <- selectBenches (match cfg) benches bsgroup-          withConfig cfg $ f selector bsgroup arg++        let mode = case (measureOnly cfg, iters cfg, quickMode cfg) of+                (Just outfile, _           , _   )  -> BenchmarkNormal $ \_ r -> gaugeIO (writeFile outfile (show r))+                (Nothing     , Just nbIters, _   )  -> BenchmarkIters nbIters+                (Nothing     , Nothing     , True)  -> BenchmarkNormal quickAnalyse+                (Nothing     , Nothing     , False) -> BenchmarkNormal compiledAnalyseStep++        withConfig cfg $ runBenchmark selector bsgroup mode++    bsgroup = BenchGroup "" bs

@@ -1,4 +1,6 @@ {-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ConstraintKinds #-} -- | -- Module      : Gauge.Optional -- Copyright   : (c) 2017-2018 Vincent Hanquez@@ -29,6 +31,7 @@ import Data.Word import Data.Data import GHC.Generics+import Basement.Compat.CallStack  -- | A type representing a sum-type free Maybe a -- where a specific tag represent Nothing@@ -50,7 +53,7 @@     isOptionalTag d = isInfinite d || isNaN d  -- | Create an optional value from a -toOptional :: OptionalTag a => a -> Optional a+toOptional :: (HasCallStack, OptionalTag a) => a -> Optional a toOptional v     | isOptionalTag v = error "Creating an optional valid value using the optional tag"     | otherwise       = Optional v@@ -70,7 +73,7 @@                      | otherwise       = Just v {-# INLINE toMaybe #-} -fromMaybe :: OptionalTag a => Maybe a -> Optional a+fromMaybe :: (HasCallStack, OptionalTag a) => Maybe a -> Optional a fromMaybe Nothing  = Optional optionalTag fromMaybe (Just v)     | isOptionalTag v = error "fromMaybe: creating an optional value using the optional tag"@@ -82,7 +85,7 @@                      | otherwise       = Optional (f v)  {-# INLINE map #-} -both :: OptionalTag a => (a -> a -> a) -> Optional a -> Optional a -> Optional a+both :: (HasCallStack, OptionalTag a) => (a -> a -> a) -> Optional a -> Optional a -> Optional a both f o1 o2     | isOmitted o1    = o2     | isOmitted o2    = o1

@@ -36,7 +36,7 @@  Number of total dependencies (direct & indirect): -* gauge: 18 dependencies+* gauge: 12 dependencies * criterion: 63 dependencies  Dependencies removed:@@ -56,30 +56,39 @@ * bytestring 0.10.8.1 * cassava 0.4.5.1 * cereal 0.5.4.0+* code-page 0.1.3+* containers 0.5.7.1 * criterion 1.2.2.0 * directory 1.3.0.0 * dlist 0.8.0.3 * erf 2.0.0.0 * exceptions 0.8.3 * filepath 1.4.1.1+* ghc-boot-th 8.0.2 * hashable 1.2.6.1+* integer-gmp 1.0.0.1 * integer-logarithms 1.0.2 * js-flot 0.8.3 * js-jquery 3.2.1+* math-functions 0.2.1.0 * microstache 1.0.1.1 * monad-par 0.3.4.8 * monad-par-extras 0.3.3 * mtl 2.2.1+* mwc-random 0.13.6.0 * optparse-applicative 0.13.2.0 * parallel 3.2.1.1 * parsec 3.1.11+* pretty 1.1.3.3 * process 1.4.3.0 * random 1.1 * scientific 0.3.5.2 * statistics 0.14.0.2 * stm 2.4.4.1 * tagged 0.8.5+* template-haskell 2.11.1.0 * text 1.2.2.2+* time 1.6.0.1 * time-locale-compat 0.1.1.3 * transformers-compat 0.5.1.4 * unix 2.7.2.1@@ -87,7 +96,8 @@ * uuid-types 1.0.3 * vector-algorithms 0.7.0.1 * vector-binary-instances 0.2.3.5-* code-page 0.1.3+* vector-th-unbox 0.2.1.6+  Criterion graph of dependencies: 

@@ -2,6 +2,7 @@ #define CYCLES_H  #include <stddef.h>+#include "Rts.h"  #if x86_64_HOST_ARCH || i386_HOST_ARCH 

@@ -2,19 +2,30 @@ #include <stdint.h>  #include <unistd.h>-#include <asm-generic/unistd.h>-#include <linux/perf_event.h>  #include "gauge-time.h" +// #define USE_PERF_EVENT++#ifdef USE_PERF_EVENT+#include <asm-generic/unistd.h>+#include <linux/perf_event.h>+#else+#include <cycles.h>+#endif++#ifdef USE_PERF_EVENT static int gauge_rdtsc_fddev = -1;+#endif  void gauge_inittime(void) {+#ifdef USE_PERF_EVENT     static struct perf_event_attr attr;     attr.type = PERF_TYPE_HARDWARE;     attr.config = PERF_COUNT_HW_CPU_CYCLES;     gauge_rdtsc_fddev = syscall (__NR_perf_event_open, &attr, 0, -1, -1, 0);+#endif }  #define timespec_to_uint64(x) (                      \@@ -25,14 +36,18 @@ void gauge_record(struct gauge_time *tr) {     struct timespec ts, ts2;-    uint64_t res;+    uint64_t res = 0;      clock_gettime(CLOCK_MONOTONIC, &ts);     clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts2);      tr->clock_nanosecs = timespec_to_uint64(ts);     tr->cpu_nanosecs = timespec_to_uint64(ts2);+#ifdef USE_PERF_EVENT     tr->rdtsc = (read (gauge_rdtsc_fddev, &res, sizeof(res)) < sizeof(res)) ? 0 : res;+#else+    tr->rdtsc = instruction_rdtsc();+#endif }  double gauge_gettime(void)

@@ -1,3 +1,15 @@+# 0.2.1++* Inline math-functions & mwc-random:+  * Remove most functions, instances and types, that are unnecessary for gauge+  * Remove unsafe seeding with partial seed (unused in gauge anyway)+  * Remove vector-th-unbox dependency (transitively template-haskell, pretty, ghc-boot-th)+  * Remove time dependency+* Re-add Gauge.Benchmark to Gauge.Main to keep the transition between criterion and gauge easy+* Fix cycles reporting on linux, osx and windows+* Add some extra callstack for reporting on partial function+* Fix compilation with Semigroup => Monoid (compilation on 8.4). still unsupported+* Add some color on terminal output  # 0.2.0 

@@ -1,5 +1,5 @@ name:           gauge-version:        0.2.0+version:        0.2.1 synopsis:       small framework for performance measurement and analysis license:        BSD3 license-file:   LICENSE@@ -19,7 +19,7 @@   GHC==7.8.4,   GHC==7.10.3,   GHC==8.0.2,-  GHC==8.2.1+  GHC==8.2.2  description:   This library provides a powerful but simple way to measure software@@ -46,6 +46,7 @@     Gauge.Time     Gauge.Optional     Gauge.CSV+    Gauge.Format      Gauge.Source.RUsage     Gauge.Source.GC@@ -75,9 +76,16 @@         Statistics.Transform         Statistics.Types         Statistics.Types.Internal+        System.Random.MWC+        Numeric.MathFunctions.Comparison+        Numeric.MathFunctions.Constants+        Numeric.SpecFunctions+        Numeric.SpecFunctions.Internal+        Numeric.Sum -  hs-source-dirs: . statistics +  hs-source-dirs: . statistics mwc-random math-functions+   include-Dirs: cbits   c-sources: cbits/cycles.c   if os(darwin)@@ -96,16 +104,10 @@     base >= 4.7 && < 5,     basement >= 0.0.4,     deepseq >= 1.1.0.0,-    mwc-random >= 0.8.0.3,     vector >= 0.7.1,     process,     directory -    -- formely statistics dependency that we need-  if flag(analysis)-      build-depends:-        math-functions >= 0.1.7-   default-language: Haskell2010   ghc-options: -O2 -Wall -funbox-strict-fields   if flag(analysis)@@ -123,7 +125,6 @@     base,     bytestring,     gauge,-    deepseq,     tasty,     tasty-hunit @@ -140,7 +141,6 @@     base,     bytestring,     gauge,-    deepseq,     tasty,     tasty-hunit @@ -157,7 +157,6 @@     base,     bytestring,     gauge,-    deepseq,     tasty,     tasty-hunit @@ -174,7 +173,6 @@     base,     bytestring,     gauge,-    deepseq,     tasty,     tasty-hunit 

@@ -0,0 +1,57 @@+-- |+-- Module    : Numeric.MathFunctions.Comparison+-- Copyright : (c) 2011 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Functions for approximate comparison of floating point numbers.+--+-- Approximate floating point comparison, based on Bruce Dawson's+-- \"Comparing floating point numbers\":+-- <http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm>+module Numeric.MathFunctions.Comparison+    ( within+    ) where++import Basement.Floating (doubleToWord)+import Data.Word (Word64)++-- |+-- Measure distance between two @Double@s in ULPs (units of least+-- precision). Note that it's different from @abs (ulpDelta a b)@+-- since it returns correct result even when 'ulpDelta' overflows.+ulpDistance :: Double+            -> Double+            -> Word64+ulpDistance a b =+  -- IEEE754 floats use most significant bit as sign bit (not+  -- 2-complement) and we need to rearrange representations of float+  -- number so that they could be compared lexicographically as+  -- Word64.+  let big     = 0x8000000000000000+      order i | i < big   = i + big+              | otherwise = maxBound - i+      ai = order ai0+      bi = order bi0+      d  | ai > bi   = ai - bi+         | otherwise = bi - ai+   in d+  where+    ai0 = doubleToWord a+    bi0 = doubleToWord b+++-- | Compare two 'Double' values for approximate equality, using+-- Dawson's method.+--+-- The required accuracy is specified in ULPs (units of least+-- precision).  If the two numbers differ by the given number of ULPs+-- or less, this function returns @True@.+within :: Int                   -- ^ Number of ULPs of accuracy desired.+       -> Double -> Double -> Bool+within ulps a b+  | ulps < 0  = False+  | otherwise = ulpDistance a b <= fromIntegral ulps

@@ -0,0 +1,115 @@+-- |+-- Module    : Numeric.MathFunctions.Constants+-- Copyright : (c) 2009, 2011 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Constant values common to much numeric code.++module Numeric.MathFunctions.Constants+    (+      -- * IEE754 constants+      m_epsilon+    , m_huge+    , m_tiny+    , m_max_exp+    , m_pos_inf+    , m_neg_inf+    , m_NaN+    , m_max_log+    , m_min_log+      -- * Mathematical constants+    , m_1_sqrt_2+    , m_2_sqrt_pi+    , m_ln_sqrt_2_pi+    , m_sqrt_2+    , m_sqrt_2_pi+    , m_eulerMascheroni+    ) where++----------------------------------------------------------------+-- IEE754 constants+----------------------------------------------------------------++-- | Largest representable finite value.+m_huge :: Double+m_huge = 1.7976931348623157e308+{-# INLINE m_huge #-}++-- | The smallest representable positive normalized value.+m_tiny :: Double+m_tiny = 2.2250738585072014e-308+{-# INLINE m_tiny #-}++-- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately+-- representable as a 'Double'.+m_max_exp :: Int+m_max_exp = 1024++-- | Positive infinity.+m_pos_inf :: Double+m_pos_inf = 1/0+{-# INLINE m_pos_inf #-}++-- | Negative infinity.+m_neg_inf :: Double+m_neg_inf = -1/0+{-# INLINE m_neg_inf #-}++-- | Not a number.+m_NaN :: Double+m_NaN = 0/0+{-# INLINE m_NaN #-}++-- | Maximum possible finite value of @log x@+m_max_log :: Double+m_max_log = 709.782712893384+{-# INLINE m_max_log #-}++-- | Logarithm of smallest normalized double ('m_tiny')+m_min_log :: Double+m_min_log = -708.3964185322641+{-# INLINE m_min_log #-}+++----------------------------------------------------------------+-- Mathematical constants+----------------------------------------------------------------++-- | @sqrt 2@+m_sqrt_2 :: Double+m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766+{-# INLINE m_sqrt_2 #-}++-- | @sqrt (2 * pi)@+m_sqrt_2_pi :: Double+m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299+{-# INLINE m_sqrt_2_pi #-}++-- | @2 / sqrt pi@+m_2_sqrt_pi :: Double+m_2_sqrt_pi = 1.1283791670955125738961589031215451716881012586579977136881+{-# INLINE m_2_sqrt_pi #-}++-- | @1 / sqrt 2@+m_1_sqrt_2 :: Double+m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883+{-# INLINE m_1_sqrt_2 #-}++-- | The smallest 'Double' ε such that 1 + ε ≠ 1.+m_epsilon :: Double+m_epsilon = encodeFloat (signif+1) expo - 1.0+    where (signif,expo) = decodeFloat (1.0::Double)++-- | @log(sqrt((2*pi))@+m_ln_sqrt_2_pi :: Double+m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171+{-# INLINE m_ln_sqrt_2_pi #-}++-- | Euler–Mascheroni constant (γ = 0.57721...)+m_eulerMascheroni :: Double+m_eulerMascheroni = 0.5772156649015328606065121+{-# INLINE m_eulerMascheroni #-}

@@ -0,0 +1,3 @@+module Numeric.SpecFunctions ( module X ) where++import Numeric.SpecFunctions.Internal as X

@@ -0,0 +1,116 @@+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-}+-- |+-- Module    : Numeric.SpecFunctions.Internal+-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Internal module with implementation of special functions.+module Numeric.SpecFunctions.Internal+    ( erf+    , erfc+    , invErf+    , invErfc+    , log2+    ) where++import Data.Bits       ((.&.), (.|.), shiftR)+import Data.Word       (Word64)+import qualified Data.Vector.Unboxed as U++import Numeric.MathFunctions.Constants++----------------------------------------------------------------+-- Error function+----------------------------------------------------------------++-- | Error function.+--+-- \[+-- \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt+-- \]+--+-- Function limits are:+--+-- \[+-- \begin{aligned}+--  &\operatorname{erf}(-\infty) &=& -1 \\+--  &\operatorname{erf}(0)       &=& \phantom{-}\,0 \\+--  &\operatorname{erf}(+\infty) &=& \phantom{-}\,1 \\+-- \end{aligned}+-- \]+erf :: Double -> Double+{-# INLINE erf #-}+erf = c_erf++-- | Complementary error function.+--+-- \[+-- \operatorname{erfc}(x) = 1 - \operatorname{erf}(x)+-- \]+--+-- Function limits are:+--+-- \[+-- \begin{aligned}+--  &\operatorname{erf}(-\infty) &=&\, 2 \\+--  &\operatorname{erf}(0)       &=&\, 1 \\+--  &\operatorname{erf}(+\infty) &=&\, 0 \\+-- \end{aligned}+-- \]+erfc :: Double -> Double+{-# INLINE erfc #-}+erfc = c_erfc++foreign import ccall "erf"  c_erf  :: Double -> Double+foreign import ccall "erfc" c_erfc :: Double -> Double+++-- | Inverse of 'erf'.+invErf :: Double -- ^ /p/ ∈ [-1,1]+       -> Double+invErf p = invErfc (1 - p)++-- | Inverse of 'erfc'.+invErfc :: Double -- ^ /p/ ∈ [0,2]+        -> Double+invErfc p+  | p == 2        = m_neg_inf+  | p == 0        = m_pos_inf+  | p >0 && p < 2 = if p <= 1 then r else -r+  | otherwise     = modErr $ "invErfc: p must be in [0,2] got " ++ show p+  where+    pp = if p <= 1 then p else 2 - p+    t  = sqrt $ -2 * log( 0.5 * pp)+    -- Initial guess+    x0 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t)+    r  = loop 0 x0+    --+    loop :: Int -> Double -> Double+    loop !j !x+      | j >= 2    = x+      | otherwise = let err = erfc x - pp+                        x'  = x + err / (1.12837916709551257 * exp(-x * x) - x * err) -- // Halley+                    in loop (j+1) x'++-- | /O(log n)/ Compute the logarithm in base 2 of the given value.+log2 :: Int -> Int+log2 v0+    | v0 <= 0   = modErr $ "log2: nonpositive input, got " ++ show v0+    | otherwise = go 5 0 v0+  where+    go !i !r !v | i == -1        = r+                | v .&. b i /= 0 = let si = U.unsafeIndex sv i+                                   in go (i-1) (r .|. si) (v `shiftR` si)+                | otherwise      = go (i-1) r v+    b = U.unsafeIndex bv+    !bv = U.fromList [ 0x02, 0x0c, 0xf0, 0xff00+                     , fromIntegral (0xffff0000 :: Word64)+                     , fromIntegral (0xffffffff00000000 :: Word64)]+    !sv = U.fromList [1,2,4,8,16,32]++modErr :: String -> a+modErr msg = error $ "Numeric.SpecFunctions." ++ msg

@@ -0,0 +1,198 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts,+    MultiParamTypeClasses, TypeFamilies #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+-- |+-- Module    : Numeric.Sum+-- Copyright : (c) 2014 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Functions for summing floating point numbers more accurately than+-- the naive 'Prelude.sum' function and its counterparts in the+-- @vector@ package and elsewhere.+--+-- When used with floating point numbers, in the worst case, the+-- 'Prelude.sum' function accumulates numeric error at a rate+-- proportional to the number of values being summed. The algorithms+-- in this module implement different methods of /compensated+-- summation/, which reduce the accumulation of numeric error so that+-- it either grows much more slowly than the number of inputs+-- (e.g. logarithmically), or remains constant.+module Numeric.Sum (+    -- * Summation type class+      Summation(..)+    , sumVector+    , kbn+    ) where++import Control.DeepSeq (NFData(..))+import Control.Monad+import Data.Data (Typeable, Data)+import Data.Vector.Generic (Vector(..), foldl')+import qualified Data.Vector.Generic.Mutable as M++import qualified Data.Foldable as F+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++-- | A class for summation of floating point numbers.+class Summation s where+    -- | The identity for summation.+    zero :: s++    -- | Add a value to a sum.+    add  :: s -> Double -> s++    -- | Sum a collection of values.+    --+    -- Example:+    -- @foo = 'sum' 'kbn' [1,2,3]@+    sum  :: (F.Foldable f) => (s -> Double) -> f Double -> Double+    sum  f = f . F.foldl' add zero+    {-# INLINE sum #-}++instance Summation Double where+    zero = 0+    add = (+)++-- | Kahan-Babuška-Neumaier summation. This is a little more+-- computationally costly than plain Kahan summation, but is /always/+-- at least as accurate.+data KBNSum = KBNSum {-# UNPACK #-} !Double {-# UNPACK #-} !Double+            deriving (Eq, Show, Typeable, Data)++newtype instance U.MVector s KBNSum = MV_KBNSum (U.MVector s (Double,Double))+newtype instance U.Vector    KBNSum = V_KBNSum  (U.Vector    (Double,Double))++instance M.MVector U.MVector KBNSum where+    {-# INLINE basicLength #-}+    {-# INLINE basicUnsafeSlice #-}+    {-# INLINE basicOverlaps #-}+    {-# INLINE basicUnsafeNew #-}+    {-# INLINE basicUnsafeReplicate #-}+    {-# INLINE basicUnsafeRead #-}+    {-# INLINE basicUnsafeWrite #-}+    {-# INLINE basicClear #-}+    {-# INLINE basicSet #-}+    {-# INLINE basicUnsafeCopy #-}+    {-# INLINE basicUnsafeGrow #-}+    basicLength (MV_KBNSum v) = M.basicLength v+    basicUnsafeSlice i n (MV_KBNSum v) = MV_KBNSum $ M.basicUnsafeSlice i n v+    basicOverlaps (MV_KBNSum v1) (MV_KBNSum v2) = M.basicOverlaps v1 v2+    basicUnsafeNew n = MV_KBNSum `liftM` M.basicUnsafeNew n+    basicUnsafeReplicate n (KBNSum a b) = MV_KBNSum `liftM` M.basicUnsafeReplicate n (a,b)+    basicUnsafeRead (MV_KBNSum v) i = uncurry KBNSum `liftM` M.basicUnsafeRead v i+    basicUnsafeWrite (MV_KBNSum v) i (KBNSum a b) = M.basicUnsafeWrite v i (a,b)+    basicClear (MV_KBNSum v) = M.basicClear v+    basicSet (MV_KBNSum v) (KBNSum a b) = M.basicSet v (a,b)+    basicUnsafeCopy (MV_KBNSum v1) (MV_KBNSum v2) = M.basicUnsafeCopy v1 v2+    basicUnsafeMove (MV_KBNSum v1) (MV_KBNSum v2) = M.basicUnsafeMove v1 v2+    basicUnsafeGrow (MV_KBNSum v) n = MV_KBNSum `liftM` M.basicUnsafeGrow v n+#if MIN_VERSION_vector(0,11,0)+    {-# INLINE basicInitialize #-}+    basicInitialize (MV_KBNSum v) = M.basicInitialize v+#endif++instance G.Vector U.Vector KBNSum where+    {-# INLINE basicUnsafeFreeze #-}+    {-# INLINE basicUnsafeThaw #-}+    {-# INLINE basicLength #-}+    {-# INLINE basicUnsafeSlice #-}+    {-# INLINE basicUnsafeIndexM #-}+    {-# INLINE elemseq #-}+    basicUnsafeFreeze (MV_KBNSum v) = V_KBNSum `liftM` G.basicUnsafeFreeze v+    basicUnsafeThaw (V_KBNSum v) = MV_KBNSum `liftM` G.basicUnsafeThaw v+    basicLength (V_KBNSum v) = G.basicLength v+    basicUnsafeSlice i n (V_KBNSum v) = V_KBNSum $ G.basicUnsafeSlice i n v+    basicUnsafeIndexM (V_KBNSum v) i = uncurry KBNSum `liftM` G.basicUnsafeIndexM v i+    basicUnsafeCopy (MV_KBNSum mv) (V_KBNSum v) = G.basicUnsafeCopy mv v+    elemseq _ = seq+++instance U.Unbox KBNSum++instance Summation KBNSum where+    zero = KBNSum 0 0+    add  = kbnAdd++instance NFData KBNSum where+    rnf !_ = ()++kbnAdd :: KBNSum -> Double -> KBNSum+kbnAdd (KBNSum sum c) x = KBNSum sum' c'+  where c' | abs sum >= abs x = c + ((sum - sum') + x)+           | otherwise        = c + ((x - sum') + sum)+        sum'                  = sum + x++-- | Return the result of a Kahan-Babuška-Neumaier sum.+kbn :: KBNSum -> Double+kbn (KBNSum sum c) = sum + c++-- | /O(n)/ Sum a vector of values.+sumVector :: (Vector v Double, Summation s) =>+             (s -> Double) -> v Double -> Double+sumVector f = f . foldl' add zero+{-# INLINE sumVector #-}++-- $usage+--+-- Most of these summation algorithms are intended to be used via the+-- 'Summation' typeclass interface. Explicit type annotations should+-- not be necessary, as the use of a function such as 'kbn' or 'kb2'+-- to extract the final sum out of a 'Summation' instance gives the+-- compiler enough information to determine the precise type of+-- summation algorithm to use.+--+-- As an example, here is a (somewhat silly) function that manually+-- computes the sum of elements in a list.+--+-- @+-- sillySumList :: [Double] -> Double+-- sillySumList = loop 'zero'+--   where loop s []     = 'kbn' s+--         loop s (x:xs) = 'seq' s' loop s' xs+--           where s'    = 'add' s x+-- @+--+-- In most instances, you can simply use the much more general 'sum'+-- function instead of writing a summation function by hand.+--+-- @+-- -- Avoid ambiguity around which sum function we are using.+-- import Prelude hiding (sum)+-- --+-- betterSumList :: [Double] -> Double+-- betterSumList xs = 'sum' 'kbn' xs+-- @++-- Note well the use of 'seq' in the example above to force the+-- evaluation of intermediate values.  If you must write a summation+-- function by hand, and you forget to evaluate the intermediate+-- values, you are likely to incur a space leak.+--+-- Here is an example of how to compute a prefix sum in which the+-- intermediate values are as accurate as possible.+--+-- @+-- prefixSum :: [Double] -> [Double]+-- prefixSum xs = map 'kbn' . 'scanl' 'add' 'zero' $ xs+-- @++-- $references+--+-- * Kahan, W. (1965), Further remarks on reducing truncation+--   errors. /Communications of the ACM/ 8(1):40.+--+-- * Neumaier, A. (1974), Rundungsfehleranalyse einiger Verfahren zur+--   Summation endlicher Summen.+--   /Zeitschrift für Angewandte Mathematik und Mechanik/ 54:39–51.+--+-- * Klein, A. (2006), A Generalized+--   Kahan-Babuška-Summation-Algorithm. /Computing/ 76(3):279-293.+--+-- * Higham, N.J. (1993), The accuracy of floating point+--   summation. /SIAM Journal on Scientific Computing/ 14(4):783–799.

@@ -0,0 +1,428 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts,+    MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies,+    ForeignFunctionInterface #-}+-- |+-- Module    : System.Random.MWC+-- Copyright : (c) 2009-2012 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Pseudo-random number generation.  This module contains code for+-- generating high quality random numbers that follow a uniform+-- distribution.+--+-- For non-uniform distributions, see the+-- 'System.Random.MWC.Distributions' module.+--+-- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222)+-- multiply-with-carry generator, which has a period of 2^8222 and+-- fares well in tests of randomness.  It is also extremely fast,+-- between 2 and 3 times faster than the Mersenne Twister.+--+-- The generator state is stored in the 'Gen' data type. It can be+-- created in several ways:+--+--   1. Using the 'withSystemRandom' call, which creates a random state.+--+--   2. Supply your own seed to 'initialize' function.+--+--   3. Finally, 'create' makes a generator from a fixed seed.+--      Generators created in this way aren't really random.+--+-- For repeatability, the state of the generator can be snapshotted+-- and replayed using the 'save' and 'restore' functions.+--+-- The simplest use is to generate a vector of uniformly distributed values:+--+-- @+--   vs \<- 'withSystemRandom' . 'asGenST' $ \\gen -> 'uniformVector' gen 100+-- @+--+-- These values can be of any type which is an instance of the class+-- 'Variate'.+--+-- To generate random values on demand, first 'create' a random number+-- generator.+--+-- @+--   gen <- 'create'+-- @+--+-- Hold onto this generator and use it wherever random values are+-- required (creating a new generator is expensive compared to+-- generating a random number, so you don't want to throw them+-- away). Get a random value using 'uniform' or 'uniformR':+--+-- @+--   v <- 'uniform' gen+-- @+--+-- @+--   v <- 'uniformR' (1, 52) gen+-- @+module System.Random.MWC+    (+    -- * Gen: Pseudo-Random Number Generators+      Gen+    , initialize+    , createSystemRandom+    , GenIO+    , splitGen++    -- * Variates: uniformly distributed values+    , Variate(..)+    , uniformVector++    ) where++#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__)+#include "MachDeps.h"+#endif++import Control.Monad           (liftM, replicateM)+import Data.Bits               ((.|.), shiftL, shiftR)+import Data.Int                (Int8, Int16, Int32, Int64)+import Data.Vector.Generic     (Vector)+import Data.Word               (Word8, Word16, Word32, Word64)+#if !MIN_VERSION_base(4,8,0)+import Data.Word               (Word)+#endif+import Foreign.Marshal.Alloc   (allocaBytes)+import Foreign.Marshal.Array   (peekArray)+import qualified Data.Vector.Generic         as G+import qualified Data.Vector.Unboxed         as I+import qualified Data.Vector.Unboxed.Mutable as M+#if defined(mingw32_HOST_OS)+import Foreign.Ptr+import Foreign.C.Types+#else+import System.IO        (IOMode(..), hGetBuf, withBinaryFile)+#endif++import Basement.Monad (PrimMonad(..))+++-- | The class of types for which we can generate uniformly+-- distributed random variates.+--+-- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222)+-- multiply-with-carry generator, which has a period of 2^8222 and+-- fares well in tests of randomness.  It is also extremely fast,+-- between 2 and 3 times faster than the Mersenne Twister.+--+-- /Note/: Marsaglia's PRNG is not known to be cryptographically+-- secure, so you should not use it for cryptographic operations.+class Variate a where+    -- | Generate a single uniformly distributed random variate.  The+    -- range of values produced varies by type:+    --+    -- * For fixed-width integral types, the type's entire range is+    --   used.+    --+    -- * For floating point numbers, the range (0,1] is used. Zero is+    --   explicitly excluded, to allow variates to be used in+    --   statistical calculations that require non-zero values+    --   (e.g. uses of the 'log' function).+    --+    -- To generate a 'Float' variate with a range of [0,1), subtract+    -- 2**(-33).  To do the same with 'Double' variates, subtract+    -- 2**(-53).+    uniform :: Gen -> IO a+    -- | Generate single uniformly distributed random variable in a+    -- given range.+    --+    -- * For integral types inclusive range is used.+    --+    -- * For floating point numbers range (a,b] is used if one ignores+    --   rounding errors.+    uniformR :: (a,a) -> Gen -> IO a++instance Variate Word32 where+    uniform  = uniform1 fromIntegral+    uniformR a b = uniformRange a b+    {-# INLINE uniform  #-}+    {-# INLINE uniformR #-}++instance Variate Word64 where+    uniform  = uniform2 wordsTo64Bit+    uniformR a b = uniformRange a b+    {-# INLINE uniform  #-}+    {-# INLINE uniformR #-}++instance Variate Int where+#if WORD_SIZE_IN_BITS < 64+    uniform = uniform1 fromIntegral+#else+    uniform = uniform2 wordsTo64Bit+#endif+    uniformR a b = uniformRange a b+    {-# INLINE uniform  #-}+    {-# INLINE uniformR #-}++instance Variate Word where+#if WORD_SIZE_IN_BITS < 64+    uniform = uniform1 fromIntegral+#else+    uniform = uniform2 wordsTo64Bit+#endif+    uniformR a b = uniformRange a b+    {-# INLINE uniform  #-}+    {-# INLINE uniformR #-}++wordsTo64Bit :: (Integral a) => Word32 -> Word32 -> a+wordsTo64Bit x y =+    fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64)+{-# INLINE wordsTo64Bit #-}++-- | State of the pseudo-random number generator. It uses mutable+-- state so same generator shouldn't be used from the different+-- threads simultaneously.+newtype Gen = Gen (M.MVector (PrimState IO) Word32)++-- | A shorter name for PRNG state in the 'IO' monad.+type GenIO = Gen++ioff, coff :: Int+ioff = 256+coff = 257++-- | Create a generator for variates using the given seed of 256 elements+--+-- @gen' <- 'initialize' . 'fromSeed' =<< 'save'@+initialize :: I.Vector Word32 -> IO Gen+initialize seed+    | fini /= 256 = error "mwc seed invalid size"+    | otherwise   = do+        q <- M.unsafeNew 258+        fill q+        M.unsafeWrite q ioff 255+        M.unsafeWrite q coff 362436+        return (Gen q)+  where+    fini = G.length seed+    fill q = go 0 where+          go i | i == 256  = return ()+               | otherwise = M.unsafeWrite q i (G.unsafeIndex seed i) >> go (i+1)+{-# INLINE initialize #-}++-- | Acquire seed from the system entropy source. On Unix machines,+-- this will attempt to use @/dev/urandom@. On Windows, it will internally+-- use @RtlGenRandom@.+acquireSeedSystem :: IO [Word32]+acquireSeedSystem = do+#if !defined(mingw32_HOST_OS)+  -- Read 256 random Word32s from /dev/urandom+  let nbytes = 1024+      random = "/dev/urandom"+  allocaBytes nbytes $ \buf -> do+    nread <- withBinaryFile random ReadMode $+               \h -> hGetBuf h buf nbytes+    peekArray (nread `div` 4) buf+#else+  let nbytes = 1024+  -- Generate 256 random Word32s from RtlGenRandom+  allocaBytes nbytes $ \buf -> do+    ok <- c_RtlGenRandom buf (fromIntegral nbytes)+    if ok then return () else fail "Couldn't use RtlGenRandom"+    peekArray (nbytes `div` 4) buf++-- Note: on 64-bit Windows, the 'stdcall' calling convention+-- isn't supported, so we use 'ccall' instead.+#if defined(i386_HOST_ARCH)+# define WINDOWS_CCONV stdcall+#elif defined(x86_64_HOST_ARCH)+# define WINDOWS_CCONV ccall+#else+# error Unknown mingw32 architecture!+#endif++-- Note: On Windows, the typical convention would be to use+-- the CryptoGenRandom API in order to generate random data.+-- However, here we use 'SystemFunction036', AKA RtlGenRandom.+--+-- This is a commonly used API for this purpose; one bonus is+-- that it avoids having to bring in the CryptoAPI library,+-- and completely sidesteps the initialization cost of CryptoAPI.+--+-- While this function is technically "subject to change" that is+-- extremely unlikely in practice: rand_s in the Microsoft CRT uses+-- this, and they can't change it easily without also breaking+-- backwards compatibility with e.g. statically linked applications.+--+-- The name 'SystemFunction036' is the actual link-time name; the+-- display name is just for giggles, I guess.+--+-- See also:+--   - http://blogs.msdn.com/b/michael_howard/archive/2005/01/14/353379.aspx+--   - https://bugzilla.mozilla.org/show_bug.cgi?id=504270+--+foreign import WINDOWS_CCONV unsafe "SystemFunction036"+  c_RtlGenRandom :: Ptr a -> CULong -> IO Bool+#endif++-- | Seed a PRNG with data from the system's fast source of pseudo-random+-- numbers.+createSystemRandom :: IO GenIO+createSystemRandom = do+  seed <- acquireSeedSystem+  initialize (I.fromList seed)++-- | Compute the next index into the state pool.  This is simply+-- addition modulo 256.+nextIndex :: Integral a => a -> Int+nextIndex i = fromIntegral j+    where j = fromIntegral (i+1) :: Word8+{-# INLINE nextIndex #-}++aa :: Word64+aa = 1540315826+{-# INLINE aa #-}++data DoubleWord32 = DoubleWord32 {-# UNPACK #-} !Word32 {-# UNPACK #-} !Word32++uniformWord32 :: Gen -> IO Word32+uniformWord32 (Gen q) = do+  i  <- nextIndex `liftM` M.unsafeRead q ioff+  c  <- fromIntegral `liftM` M.unsafeRead q coff+  qi <- fromIntegral `liftM` M.unsafeRead q i+  let t  = aa * qi + c+      c' = fromIntegral (t `shiftR` 32)+      x  = fromIntegral t + c'+      (DoubleWord32 x' c'')  | x < c'    = DoubleWord32 (x + 1) (c' + 1)+                             | otherwise = DoubleWord32 x c'+  M.unsafeWrite q i x'+  M.unsafeWrite q ioff (fromIntegral i)+  M.unsafeWrite q coff (fromIntegral c'')+  return x'+{-# INLINE uniformWord32 #-}++uniform1 :: (Word32 -> a) -> Gen -> IO a+uniform1 f gen = do+  i <- uniformWord32 gen+  return $! f i+{-# INLINE uniform1 #-}++uniform2 :: (Word32 -> Word32 -> a) -> Gen -> IO a+uniform2 f (Gen q) = do+  i  <- nextIndex `liftM` M.unsafeRead q ioff+  let j = nextIndex i+  c  <- fromIntegral `liftM` M.unsafeRead q coff+  qi <- fromIntegral `liftM` M.unsafeRead q i+  qj <- fromIntegral `liftM` M.unsafeRead q j+  let t   = aa * qi + c+      c'  = fromIntegral (t `shiftR` 32)+      x   = fromIntegral t + c'+      DoubleWord32 x' c'' | x < c'    = DoubleWord32 (x + 1) (c' + 1)+                          | otherwise = DoubleWord32 x c'+      u   = aa * qj + fromIntegral c''+      d'  = fromIntegral (u `shiftR` 32)+      y   = fromIntegral u + d'+      DoubleWord32 y' d'' | y < d'    = DoubleWord32 (y + 1) (d' + 1)+                          | otherwise = DoubleWord32 y d'+  M.unsafeWrite q i x'+  M.unsafeWrite q j y'+  M.unsafeWrite q ioff (fromIntegral j)+  M.unsafeWrite q coff (fromIntegral d'')+  return $! f x' y'+{-# INLINE uniform2 #-}++-- Type family for fixed size integrals. For signed data types it's+-- its unsigned couterpart with same size and for unsigned data types+-- it's same type+type family Unsigned a :: *++type instance Unsigned Int8  = Word8+type instance Unsigned Int16 = Word16+type instance Unsigned Int32 = Word32+type instance Unsigned Int64 = Word64++type instance Unsigned Word8  = Word8+type instance Unsigned Word16 = Word16+type instance Unsigned Word32 = Word32+type instance Unsigned Word64 = Word64++-- This is workaround for bug #25.+--+-- GHC-7.6 has a bug (#8072) which results in calculation of wrong+-- number of buckets in function `uniformRange'. Consequently uniformR+-- generates values in wrong range.+--+-- Bug only affects 32-bit systems and Int/Word data types. Word32+-- works just fine. So we set Word32 as unsigned counterpart for Int+-- and Word on 32-bit systems. It's done only for GHC-7.6 because+-- other versions are unaffected by the bug and we expect that GHC may+-- optimise code which uses Word better.+#if (WORD_SIZE_IN_BITS < 64) && (__GLASGOW_HASKELL__ == 706)+type instance Unsigned Int   = Word32+type instance Unsigned Word  = Word32+#else+type instance Unsigned Int   = Word+type instance Unsigned Word  = Word+#endif+++-- Subtract two numbers under assumption that x>=y and store result in+-- unsigned data type of same size+sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a+sub x y = fromIntegral x - fromIntegral y+{-# INLINE sub #-}++add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a+add m x = m + fromIntegral x+{-# INLINE add #-}++-- Generate uniformly distributed value in inclusive range.+--+-- NOTE: This function must be fully applied. Otherwise it won't be+--       inlined, which will cause a severe performance loss.+--+-- > uniformR     = uniformRange      -- won't be inlined+-- > uniformR a b = uniformRange a b  -- will be inlined+uniformRange :: ( Integral a, Bounded a, Variate a+                , Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a))+             => (a,a) -> Gen -> IO a+uniformRange (x1,x2) g+  | n == 0    = uniform g   -- Abuse overflow in unsigned types+  | otherwise = loop+  where+    -- Allow ranges where x2<x1+    (i, j) | x1 < x2   = (x1, x2)+           | otherwise = (x2, x1)+    n       = 1 + sub j i+    buckets = maxBound `div` n+    maxN    = buckets * n+    loop    = do x <- uniform g+                 if x < maxN then return $! add i (x `div` buckets)+                             else loop+{-# INLINE uniformRange #-}++-- | Generate a vector of pseudo-random variates.  This is not+-- necessarily faster than invoking 'uniform' repeatedly in a loop,+-- but it may be more convenient to use in some situations.+uniformVector :: (Variate a, Vector v a)+             => Gen -> Int -> IO (v a)+uniformVector gen n = G.replicateM n (uniform gen)+{-# INLINE uniformVector #-}++-- | Split a generator into several that can run independently.+splitGen :: Int -> GenIO -> IO [GenIO]+splitGen n gen+  | n <= 0    = return []+  | otherwise =+  fmap (gen:) . replicateM (n-1) $+  initialize =<< uniformVector gen 256++-- $references+--+-- * Marsaglia, G. (2003) Seeds for random number generators.+--   /Communications of the ACM/ 46(5):90&#8211;93.+--   <http://doi.acm.org/10.1145/769800.769827>+--+-- * Doornik, J.A. (2007) Conversion of high-period random numbers to+--   floating point.+--   /ACM Transactions on Modeling and Computer Simulation/ 17(1).+--   <http://www.doornik.com/research/randomdouble.pdf>

@@ -20,11 +20,10 @@ import Statistics.Function as F import Statistics.Matrix import Statistics.Matrix.Algorithms (qr)-import Statistics.Resampling (splitGen) import Statistics.Types      (Estimate(..),ConfInt,CL,estimateFromInterval,significanceLevel) import Statistics.Sample (mean) import Statistics.Sample.Internal (sum)-import System.Random.MWC (GenIO, uniformR)+import System.Random.MWC (GenIO, uniformR, splitGen) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U

@@ -23,16 +23,12 @@     , resample       -- * Jackknife     , jackknife-      -- * Helper functions-    , splitGen     ) where  import Control.Concurrent (forkIO, newChan, readChan, writeChan) import Control.Monad-import Basement.Monad (PrimMonad(..)) import Data.Data (Data, Typeable) import Data.Vector.Generic (unsafeFreeze)-import Data.Word (Word32) import qualified Data.Foldable as T import qualified Data.Traversable as T import qualified Data.Vector.Generic as G@@ -45,7 +41,7 @@ import Statistics.Function (indices, inplaceSortIO) import Statistics.Sample (mean, stdDev, variance, varianceUnbiased) import Statistics.Types (Sample)-import System.Random.MWC (Gen, GenIO, initialize, uniformR, uniformVector)+import System.Random.MWC (Gen, GenIO, uniformR, splitGen)   ----------------------------------------------------------------@@ -131,8 +127,7 @@   ests' = map estimate ests  -- | Create vector using resamples-resampleVector :: G.Vector v a-               => Gen (PrimState IO) -> v a -> IO (v a)+resampleVector :: G.Vector v a => Gen -> v a -> IO (v a) resampleVector gen v   = G.replicateM n $ do i <- uniformR (0,n-1) gen                         return $! G.unsafeIndex v i@@ -209,11 +204,3 @@ singletonErr :: String -> a singletonErr func = error $                     "Statistics.Resampling." ++ func ++ ": singleton input"---- | Split a generator into several that can run independently.-splitGen :: Int -> GenIO -> IO [GenIO]-splitGen n gen-  | n <= 0    = return []-  | otherwise =-  fmap (gen:) . replicateM (n-1) $-  initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))