diff --git a/erf-native.cabal b/erf-native.cabal --- a/erf-native.cabal +++ b/erf-native.cabal @@ -1,8 +1,8 @@ name: erf-native -version: 1.0.0.0 -stability: experimental +version: 1.0.0.1 +stability: provisional -cabal-version: >= 1.2 +cabal-version: >= 1.6 build-type: Simple author: James Cook @@ -18,9 +18,12 @@ so I'm throwing this one out there. It incorporates code translated from GSL's C source, and so is licensed under the GPL. +source-repository head + type: git + location: git://github.com/mokus0/erf-native.git + Library hs-source-dirs: src exposed-modules: Data.Number.Erf other-modules: GSL.SpecFunc.Erf - NR.Ch5.S8 - build-depends: base >= 3 && <5, vector + build-depends: base >= 3 && <5, polynomial diff --git a/src/GSL/SpecFunc/Erf.hs b/src/GSL/SpecFunc/Erf.hs --- a/src/GSL/SpecFunc/Erf.hs +++ b/src/GSL/SpecFunc/Erf.hs @@ -20,8 +20,7 @@ -- 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 +import Math.Polynomial.Chebyshev -- I'm lazy: these should be implemented as native Float operations -- with truncated series, etc., but they aren't. @@ -55,13 +54,13 @@ where ax = abs x estimate - | ax <= 1 = chebEval erfc_xlt1_cs (2 * ax - 1) + | ax <= 1 = evalChebyshevSeries erfc_xlt1_cs (2 * ax - 1) | ax <= 5 = let t = 0.5 * (ax - 3) - in exp (negate x*x) * chebEval erfc_x15_cs t + in exp (negate x*x) * evalChebyshevSeries erfc_x15_cs t | ax < 10 = let t = (2 * ax - 15) / 5 - in (exp (negate x*x) / ax) * chebEval erfc_x510_cs t + in (exp (negate x*x) / ax) * evalChebyshevSeries erfc_x510_cs t | otherwise = erfc8 ax erfc8 :: Double -> Double @@ -95,89 +94,74 @@ 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 - } +erfc_xlt1_cs = + [ 1.06073416421769980345174155056 / 2 + , -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 + ] +erfc_x15_cs = + [ 0.44045832024338111077637466616 / 2 + , -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 + ] +erfc_x510_cs = + [ 1.11684990123545698684297865808 / 2 + , 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 + ] diff --git a/src/NR/Ch5/S8.hs b/src/NR/Ch5/S8.hs deleted file mode 100644 --- a/src/NR/Ch5/S8.hs +++ /dev/null @@ -1,34 +0,0 @@ -{-# LANGUAGE RecordWildCards, BangPatterns, FlexibleContexts #-} --- |Evaluating Chebyshev polynomials, translated from code given --- in Numerical Recipes, 3d Ed., Section 5.8. -module NR.Ch5.S8 where - -import Data.Vector.Generic (Vector, generate, (!)) - -data Chebyshev v a b = Chebyshev - { chebOrder :: Int - -- ^ Order of the (truncated) approximation, which may be different - -- from the order of the fitted approximation, which is - -- @length chebCoeffs - 1@ - , chebCoeffs :: v b - , chebA :: a - , chebB :: a - } deriving (Eq, Show) - -chebEval :: (Real a, Fractional b, Vector v b) => Chebyshev v a b -> a -> b -chebEval c = chebEvalM (chebOrder c) c - -chebEvalM :: (Real a, Fractional b, Vector v b) => Int -> Chebyshev v a b -> a -> b -chebEvalM order Chebyshev{..} x = go order 0 0 - where - go !j !dd !d - | j > 0 = go (j-1) d (y2*d - dd + (chebCoeffs!j)) - | otherwise = y * d - dd + 0.5 * (chebCoeffs!0) - - y = (2 * realToFrac x - a - b) / (b - a) - y2 = 2 * y - - a = realToFrac chebA - b = realToFrac chebB - -