1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
{-# LANGUAGE CPP #-}
{-# LANGUAGE ForeignFunctionInterface #-}
-- |
-- Functions which have different implementations on different platforms
module Numeric.SpecFunctions.Compat (
erf
, erfc
, log1p
, expm1
) where
#if !defined(USE_SYSTEM_ERF) || !defined(USE_SYSTEM_EXPM1)
import qualified Data.Vector.Unboxed as U
#endif
#if !defined(USE_SYSTEM_ERF)
import Numeric.Polynomial.Chebyshev (chebyshev)
import Numeric.Polynomial (evaluateOddPolynomial)
#endif
#if !defined(USE_SYSTEM_EXPM1)
import Control.Applicative (liftA2)
import Numeric.Polynomial.Chebyshev (chebyshevBroucke)
import Numeric.Series (scanSequence,sumSeries,enumSequenceFrom)
import Numeric.MathFunctions.Constants
#endif
#if defined(USE_SYSTEM_EXPM1)
import GHC.Float (log1p,expm1)
#endif
----------------------------------------------------------------
-- erf & erfc
--
-- We provide pure haskell implementation for GHCJS and accessible on
-- GHC via flag
----------------------------------------------------------------
#if defined(USE_SYSTEM_ERF)
erf :: Double -> Double
erf = c_erf
{-# INLINE erf #-}
erfc :: Double -> Double
erfc = c_erfc
{-# INLINE erfc #-}
foreign import ccall unsafe "erf" c_erf :: Double -> Double
foreign import ccall unsafe "erfc" c_erfc :: Double -> Double
#else
erf :: Double -> Double
erf x
-- Computing erf as 1-erfc loses precision near 0 so we switch to
-- Taylor expansion here
| abs x < 0.1 = 0.56418958354775629
* evaluateOddPolynomial x erfTaylorSeries
| x < 0 = (-1) + erfcCheb (-x)
| otherwise = 1 - erfcCheb x
erfTaylorSeries :: U.Vector Double
{-# NOINLINE erfTaylorSeries #-}
erfTaylorSeries = U.fromList
[ 2
, -2/3
, 1/5
, -1/21
, 1/108
, -1/660
, 1/4680
]
erfc :: Double -> Double
erfc x | x < 0 = 2 - erfcCheb (-x)
| otherwise = erfcCheb x
-- Adapted from Numerical Recipes §6.2.2
erfcCheb :: Double -> Double
erfcCheb z
= t * exp( -z * z + chebyshev ty erfcCoef )
where
-- We're using approximation:
--
-- erfc(z) ≈ t·exp(-z² + P(t))
-- t = 2 / (2 + z)
t = 2 / (2 + z)
ty = 2 * t - 1
erfcCoef :: U.Vector Double
{-# NOINLINE erfcCoef #-}
erfcCoef = U.fromList
[ -0.6513268598908546 , 6.4196979235649026e-1 , 1.9476473204185836e-2
, -9.561514786808631e-3 , -9.46595344482036e-4 , 3.66839497852761e-4
, 4.2523324806907e-5 , -2.0278578112534e-5 , -1.624290004647e-6
, 1.303655835580e-6 , 1.5626441722e-8 , -8.5238095915e-8
, 6.529054439e-9 , 5.059343495e-9 , -9.91364156e-10
, -2.27365122e-10 , 9.6467911e-11 , 2.394038e-12
, -6.886027e-12 , 8.94487e-13 , 3.13092e-13
, -1.12708e-13 , 3.81e-16 , 7.106e-15
, -1.523e-15 , -9.4e-17 , 1.21e-16
, -2.8e-17
]
#endif
----------------------------------------------------------------
-- expm1 & log1p
--
-- We use one provided by base of for GHCJS use hand-coded one
----------------------------------------------------------------
#if !defined(USE_SYSTEM_EXPM1)
-- | Compute @exp x - 1@ without loss of accuracy for x near zero.
expm1 :: Double -> Double
-- NOTE: this is simplest implementation and not terribly efficient.
expm1 x
| x < (-37.42994775023705) = -1
| x > m_max_log = m_pos_inf
| abs x > 0.5 = exp x - 1
| otherwise = sumSeries $ liftA2 (*) (scanSequence (*) x (pure x))
(1 / scanSequence (*) 1 (enumSequenceFrom 2))
-- | Compute the natural logarithm of 1 + @x@. This is accurate even
-- for values of @x@ near zero, where use of @log(1+x)@ would lose
-- precision.
log1p :: Double -> Double
log1p x
| x == 0 = 0
| x == -1 = m_neg_inf
| x < -1 = m_NaN
| x' < m_epsilon * 0.5 = x
| (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0)
= x * (1 - x * 0.5)
| x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs)
| otherwise = log (1 + x)
where
x' = abs x
coeffs = U.fromList [
0.10378693562743769800686267719098e+1,
-0.13364301504908918098766041553133e+0,
0.19408249135520563357926199374750e-1,
-0.30107551127535777690376537776592e-2,
0.48694614797154850090456366509137e-3,
-0.81054881893175356066809943008622e-4,
0.13778847799559524782938251496059e-4,
-0.23802210894358970251369992914935e-5,
0.41640416213865183476391859901989e-6,
-0.73595828378075994984266837031998e-7,
0.13117611876241674949152294345011e-7,
-0.23546709317742425136696092330175e-8,
0.42522773276034997775638052962567e-9,
-0.77190894134840796826108107493300e-10,
0.14075746481359069909215356472191e-10,
-0.25769072058024680627537078627584e-11,
0.47342406666294421849154395005938e-12,
-0.87249012674742641745301263292675e-13,
0.16124614902740551465739833119115e-13,
-0.29875652015665773006710792416815e-14,
0.55480701209082887983041321697279e-15,
-0.10324619158271569595141333961932e-15
]
#endif
|
