00001 # ifndef CPPAD_ERF_INCLUDED
00002 # define CPPAD_ERF_INCLUDED
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 # include <cmath>
00190 # include <cppad/cppad.hpp>
00191
00192
00193 namespace CppAD {
00194
00195 template <class Type>
00196 Type erf_template(const Type &x)
00197 { using CppAD::exp;
00198
00199 static const Type
00200 tiny = static_cast<Type>( 1e-300 ),
00201 one = static_cast<Type>( 1.00000000000000000000e+00 ),
00202
00203 erx = static_cast<Type>( 8.45062911510467529297e-01 ),
00204
00205
00206
00207 pp0 = static_cast<Type>( 1.28379167095512558561e-01 ),
00208 pp1 = static_cast<Type>( -3.25042107247001499370e-01 ),
00209 pp2 = static_cast<Type>( -2.84817495755985104766e-02 ),
00210 pp3 = static_cast<Type>( -5.77027029648944159157e-03 ),
00211 pp4 = static_cast<Type>( -2.37630166566501626084e-05 ),
00212 qq1 = static_cast<Type>( 3.97917223959155352819e-01 ),
00213 qq2 = static_cast<Type>( 6.50222499887672944485e-02 ),
00214 qq3 = static_cast<Type>( 5.08130628187576562776e-03 ),
00215 qq4 = static_cast<Type>( 1.32494738004321644526e-04 ),
00216 qq5 = static_cast<Type>( -3.96022827877536812320e-06 ),
00217
00218
00219
00220 pa0 = static_cast<Type>( -2.36211856075265944077e-03 ),
00221 pa1 = static_cast<Type>( 4.14856118683748331666e-01 ),
00222 pa2 = static_cast<Type>( -3.72207876035701323847e-01 ),
00223 pa3 = static_cast<Type>( 3.18346619901161753674e-01 ),
00224 pa4 = static_cast<Type>( -1.10894694282396677476e-01 ),
00225 pa5 = static_cast<Type>( 3.54783043256182359371e-02 ),
00226 pa6 = static_cast<Type>( -2.16637559486879084300e-03 ),
00227 qa1 = static_cast<Type>( 1.06420880400844228286e-01 ),
00228 qa2 = static_cast<Type>( 5.40397917702171048937e-01 ),
00229 qa3 = static_cast<Type>( 7.18286544141962662868e-02 ),
00230 qa4 = static_cast<Type>( 1.26171219808761642112e-01 ),
00231 qa5 = static_cast<Type>( 1.36370839120290507362e-02 ),
00232 qa6 = static_cast<Type>( 1.19844998467991074170e-02 ),
00233
00234
00235
00236 ra0 = static_cast<Type>( -9.86494403484714822705e-03 ),
00237 ra1 = static_cast<Type>( -6.93858572707181764372e-01 ),
00238 ra2 = static_cast<Type>( -1.05586262253232909814e+01 ),
00239 ra3 = static_cast<Type>( -6.23753324503260060396e+01 ),
00240 ra4 = static_cast<Type>( -1.62396669462573470355e+02 ),
00241 ra5 = static_cast<Type>( -1.84605092906711035994e+02 ),
00242 ra6 = static_cast<Type>( -8.12874355063065934246e+01 ),
00243 ra7 = static_cast<Type>( -9.81432934416914548592e+00 ),
00244 sa1 = static_cast<Type>( 1.96512716674392571292e+01 ),
00245 sa2 = static_cast<Type>( 1.37657754143519042600e+02 ),
00246 sa3 = static_cast<Type>( 4.34565877475229228821e+02 ),
00247 sa4 = static_cast<Type>( 6.45387271733267880336e+02 ),
00248 sa5 = static_cast<Type>( 4.29008140027567833386e+02 ),
00249 sa6 = static_cast<Type>( 1.08635005541779435134e+02 ),
00250 sa7 = static_cast<Type>( 6.57024977031928170135e+00 ),
00251 sa8 = static_cast<Type>( -6.04244152148580987438e-02 ),
00252
00253
00254
00255 rb0 = static_cast<Type>( -9.86494292470009928597e-03 ),
00256 rb1 = static_cast<Type>( -7.99283237680523006574e-01 ),
00257 rb2 = static_cast<Type>( -1.77579549177547519889e+01 ),
00258 rb3 = static_cast<Type>( -1.60636384855821916062e+02 ),
00259 rb4 = static_cast<Type>( -6.37566443368389627722e+02 ),
00260 rb5 = static_cast<Type>( -1.02509513161107724954e+03 ),
00261 rb6 = static_cast<Type>( -4.83519191608651397019e+02 ),
00262 sb1 = static_cast<Type>( 3.03380607434824582924e+01 ),
00263 sb2 = static_cast<Type>( 3.25792512996573918826e+02 ),
00264 sb3 = static_cast<Type>( 1.53672958608443695994e+03 ),
00265 sb4 = static_cast<Type>( 3.19985821950859553908e+03 ),
00266 sb5 = static_cast<Type>( 2.55305040643316442583e+03 ),
00267 sb6 = static_cast<Type>( 4.74528541206955367215e+02 ),
00268 sb7 = static_cast<Type>( -2.24409524465858183362e+01 );
00269
00270 Type R,S,P,Q,s,y,z,r;
00271
00272 Type ax, value;
00273
00274 Type zero(0);
00275
00276
00277 ax = CondExpGe(x, zero, x, -x);
00278
00279
00280 z = ax*ax;
00281 r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
00282 s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
00283 y = r/s;
00284 Type Le0Dot84375 = ax + ax*y;
00285
00286
00287 s = ax-one;
00288 P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
00289 Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
00290 value = erx + P/Q;
00291
00292
00293 Type Le1Dot25 = CondExpGt(ax, Type(0.84375), value, Le0Dot84375);
00294
00295
00296 s = one/(ax*ax);
00297 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
00298 ra5+s*(ra6+s*ra7))))));
00299 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
00300 sa5+s*(sa6+s*(sa7+s*sa8)))))));
00301 r = exp(Type(-0.5625)) * exp(-ax*ax+R/S);
00302 value = one-r/ax;
00303
00304
00305 Type Le1OverDot35 = CondExpGt(ax, Type(1.25), value, Le1Dot25);
00306
00307
00308 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
00309 rb5+s*rb6)))));
00310 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
00311 sb5+s*(sb6+s*sb7))))));
00312 r = exp(Type(-0.5625)) * exp(-ax*ax+R/S);
00313 value = one-r/ax;
00314
00315
00316 Type Le6 = CondExpGt(ax, Type(1/.35), value, Le1OverDot35);
00317
00318
00319 Type Gt6 = one - tiny;
00320
00321
00322 Type Gt0 = CondExpGt(ax, Type(6.), Gt6, Le6);
00323
00324
00325
00326 return CondExpGe(x, zero, Gt0, -Gt0);
00327 }
00328
00329 inline float erf(const float &x)
00330 { return erf_template(x); }
00331
00332 inline double erf(const double &x)
00333 { return erf_template(x); }
00334
00335 template <class Base>
00336 inline AD<Base> erf(const AD<Base> &x)
00337 { return erf_template(x); }
00338
00339 template <class Base>
00340 inline AD<Base> erf(const VecAD_reference<Base> &x)
00341 { return erf_template( x.ADBase() ); }
00342
00343
00344 }
00345
00346 # endif