Erf.hs

-- GSL.SpecFunc.Erf - translated from specfunc/erfc.c in GSL sources.
-- Original header from specfunc/erfc.c follows.

-- specfunc/erfc.c
-- 
-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003 Gerard Jungman
-- 
-- This program is free software; you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation; either version 3 of the License, or (at
-- your option) any later version.
-- 
-- This program is distributed in the hope that it will be useful, but
-- WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-- General Public License for more details.
-- 
-- You should have received a copy of the GNU General Public License
-- along with this program; if not, write to the Free Software
-- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
module GSL.SpecFunc.Erf (erf, erfc, erff, erfcf) where

import qualified Data.Vector as V
import NR.Ch5.S8

-- I'm lazy: these should be implemented as native Float operations
-- with truncated series, etc., but they aren't.
erff :: Float -> Float
erff = realToFrac . erf . realToFrac
erfcf :: Float -> Float
erfcf = realToFrac . erfc . realToFrac


erf :: Double -> Double
erf x
    | abs x < 1 = erfseries x
    | otherwise = 1 - erfc x

erfseries :: Double -> Double
erfseries x = go 1 x x
    where
        go k coef e
            | k < 30    = 
                let
                    k'      = k + 1
                    coef'   = coef * negate (x*x/k)
                    e'      =  e + coef' / (2 * k + 1)
                in go k' coef' e'
            | otherwise = 2 / sqrt pi * e
 
erfc :: Double -> Double
erfc x 
    | x >= 0    = estimate
    | otherwise = 2 - estimate
    where
        ax = abs x
        estimate
            | ax <= 1   = chebEval erfc_xlt1_cs (2 * ax - 1)
            | ax <= 5   = 
                let t = 0.5 * (ax - 3)
                in exp (negate x*x) * chebEval erfc_x15_cs t
            | ax < 10   = 
                let t = (2 * ax - 15) / 5
                in (exp (negate x*x) / ax) * chebEval erfc_x510_cs t
            | otherwise = erfc8 ax

erfc8 :: Double -> Double
erfc8 x = erfc8_sum x * exp (negate x * x)

-- |estimates erfc(x) valid for 8 < x < 100
-- This is based on index 5725 in Hart et al
erfc8_sum :: Double -> Double
erfc8_sum x = num / den
    where
        p = [ 2.97886562639399288862
            , 7.409740605964741794425
            , 6.1602098531096305440906
            , 5.019049726784267463450058
            , 1.275366644729965952479585264
            , 0.5641895835477550741253201704
            ]
        q = [ 3.3690752069827527677
            , 9.608965327192787870698
            , 17.08144074746600431571095
            , 12.0489519278551290360340491
            , 9.396034016235054150430579648
            , 2.260528520767326969591866945
            , 1.0
            ]
        num = poly p x
        den = poly q x

poly p x = go p
    where
        go []       = 0
        go (c:cs)   = c + x * go cs 

erfc_xlt1_cs = Chebyshev
    { chebCoeffs    = V.fromList
        [   1.06073416421769980345174155056
        , -0.42582445804381043569204735291
        ,  0.04955262679620434040357683080
        ,  0.00449293488768382749558001242
        , -0.00129194104658496953494224761
        , -0.00001836389292149396270416979
        ,  0.00002211114704099526291538556
        , -5.23337485234257134673693179020e-7
        , -2.78184788833537885382530989578e-7
        ,  1.41158092748813114560316684249e-8
        ,  2.72571296330561699984539141865e-9
        , -2.06343904872070629406401492476e-10
        , -2.14273991996785367924201401812e-11
        ,  2.22990255539358204580285098119e-12
        ,  1.36250074650698280575807934155e-13
        , -1.95144010922293091898995913038e-14
        , -6.85627169231704599442806370690e-16
        ,  1.44506492869699938239521607493e-16
        ,  2.45935306460536488037576200030e-18
        , -9.29599561220523396007359328540e-19
        ]
    , chebOrder     = 19
    , chebA         = -1
    , chebB         = 1
    }
erfc_x15_cs = Chebyshev
    { chebCoeffs    = V.fromList
        [  0.44045832024338111077637466616
        , -0.143958836762168335790826895326
        ,  0.044786499817939267247056666937
        , -0.013343124200271211203618353102
        ,  0.003824682739750469767692372556
        , -0.001058699227195126547306482530
        ,  0.000283859419210073742736310108
        , -0.000073906170662206760483959432
        ,  0.000018725312521489179015872934
        , -4.62530981164919445131297264430e-6
        ,  1.11558657244432857487884006422e-6
        , -2.63098662650834130067808832725e-7
        ,  6.07462122724551777372119408710e-8
        , -1.37460865539865444777251011793e-8
        ,  3.05157051905475145520096717210e-9
        , -6.65174789720310713757307724790e-10
        ,  1.42483346273207784489792999706e-10
        , -3.00141127395323902092018744545e-11
        ,  6.22171792645348091472914001250e-12
        , -1.26994639225668496876152836555e-12
        ,  2.55385883033257575402681845385e-13
        , -5.06258237507038698392265499770e-14
        ,  9.89705409478327321641264227110e-15
        , -1.90685978789192181051961024995e-15
        ,  3.50826648032737849245113757340e-16
        ]
    , chebOrder     = 24
    , chebA         = -1
    , chebB         = 1
    }
erfc_x510_cs = Chebyshev
    { chebCoeffs    = V.fromList
        [  1.11684990123545698684297865808
        ,  0.003736240359381998520654927536
        , -0.000916623948045470238763619870
        ,  0.000199094325044940833965078819
        , -0.000040276384918650072591781859
        ,  7.76515264697061049477127605790e-6
        , -1.44464794206689070402099225301e-6
        ,  2.61311930343463958393485241947e-7
        , -4.61833026634844152345304095560e-8
        ,  8.00253111512943601598732144340e-9
        , -1.36291114862793031395712122089e-9
        ,  2.28570483090160869607683087722e-10
        , -3.78022521563251805044056974560e-11
        ,  6.17253683874528285729910462130e-12
        , -9.96019290955316888445830597430e-13
        ,  1.58953143706980770269506726000e-13
        , -2.51045971047162509999527428316e-14
        ,  3.92607828989125810013581287560e-15
        , -6.07970619384160374392535453420e-16
        ,  9.12600607264794717315507477670e-17
  ]
    , chebOrder     = 19
    , chebA         = -1
    , chebB         = 1
    }