|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-03-28 07:27:39
Author: johnmaddock
Date: 2008-03-28 07:27:38 EDT (Fri, 28 Mar 2008)
New Revision: 43908
URL: http://svn.boost.org/trac/boost/changeset/43908
Log:
Updated minimax docs, and program.
Changed erf/erfc approximations to more efficient versions.
Text files modified:
sandbox/math_toolkit/boost/math/special_functions/erf.hpp | 965 ++++++++++++++++++++++++---------------
sandbox/math_toolkit/libs/math/doc/sf_and_dist/minimax.qbk | 12
sandbox/math_toolkit/libs/math/minimax/f.cpp | 24
3 files changed, 623 insertions(+), 378 deletions(-)
Modified: sandbox/math_toolkit/boost/math/special_functions/erf.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/erf.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/erf.hpp 2008-03-28 07:27:38 EDT (Fri, 28 Mar 2008)
@@ -207,9 +207,27 @@
}
else
{
- static const T n[7] = { 0.00337916709551257778174L, -0.000147024115786688745475L, -0.37463022236812520164L, 0.0163061594494816999803L, -0.0534354147807331748737L, 0.00161898096813581982844L, -0.0059528010489182840404L };
- static const T d[7] = { 1, -0.0435089806536379531594L, 0.442761965043509204727L, -0.017375974533016704678L, 0.0772756490303260060769L, -0.00210552465858669941879L, 0.00544772980263244037286L };
- result = static_cast<T>(z * 1.125f + z * tools::evaluate_polynomial(n, z) / tools::evaluate_polynomial(d, z));
+ // Maximum Deviation Found: 1.561e-17
+ // Expected Error Term: 1.561e-17
+ // Maximum Relative Change in Control Points: 1.155e-04
+ // Max Error found at double precision = 2.961182e-17
+
+ static const T Y = 1.044948577880859375f;
+ static const T P[] = {
+ 0.0834305892146531832907L,
+ -0.338165134459360935041L,
+ -0.0509990735146777432841L,
+ -0.00772758345802133288487L,
+ -0.000322780120964605683831L,
+ };
+ static const T Q[] = {
+ 1L,
+ 0.455004033050794024546L,
+ 0.0875222600142252549554L,
+ 0.00858571925074406212772L,
+ 0.000370900071787748000569L,
+ };
+ result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
}
}
else if((z < 14) || ((z < 28) && invert))
@@ -218,90 +236,113 @@
// We'll be calculating erfc:
//
invert = !invert;
- T r, b;
- if(z < 0.75)
+ if(z < 1.5f)
{
- // Worst case absolute error found: 8.554649561e-018
- static const T n[5] = { -0.0361790390718262468222L, 0.301888464724047222196L, 0.201731143672633894981L, 0.0659353268087389983319L, 0.00721876720062364930761L };
- static const T d[6] = { 1, 1.58814245739127341535L, 0.99354580430196422336L, 0.291753007176902027213L, 0.033994791234913855515L, -0.000104234653166533504303L };
- static const float f0 = 0.3440242112F;
- r = tools::evaluate_polynomial(n, z - 0.5) / tools::evaluate_polynomial(d, z - 0.5);
- b = f0;
- }
- else if(z < 1.25)
- {
- // Worst case absolute error found: 6.50251514e-018
- static const T n[6] = { -0.039787689261113685983L, 0.160309168830518003303L, 0.163049978514596540313L, 0.0710685660158400750009L, 0.01497188097404877543L, 0.00130080628375002584279L };
- static const T d[6] = { 1, 1.77564880074171280407L, 1.31438791181040008779L, 0.509359151038517059748L, 0.103958527905812829559L, 0.00901292460643094469406L };
- static const float f0 = 0.419990927F;
- r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75);
- b = f0;
- }
- else if(z < 2.25)
- {
- // Worst case absolute error found: 1.132743504e-017
- static const T n[6] = { -0.0300838560557949724172L, 0.0592886319615167248092L, 0.0622294724048409148736L, 0.0248575228109427909578L, 0.00463781847004901844581L, 0.000347305179334822548368L };
- static const T d[7] = { 1, 1.57915060645728571344L, 1.03342495188878679417L, 0.35158678814344218974L, 0.062469256580984456783L, 0.00466640448020624599948L, 0.290106403940303572448e-6L };
- static const float f0 = 0.4898625016F;
- r = tools::evaluate_polynomial(n, z - 1.25) / tools::evaluate_polynomial(d, z - 1.25);
- b = f0;
- }
- else if(z < 3.5)
- {
- // Worst case absolute error found: 3.446364609e-018
- static const T n[6] = { -0.0117907570137227857015L, 0.0162667227692515660221L, 0.0175329212378413544794L, 0.00620897681269247137578L, 0.000986614895094589251706L, 0.601354618401624353425e-4L };
- static const T d[6] = { 1, 1.33374851361555383557L, 0.73227756904205983415L, 0.207410266363727673685L, 0.0304034048466731110163L, 0.00185296959991832048613L };
- static const float f0 = 0.5317370892F;
- r = tools::evaluate_polynomial(n, z - 2.25) / tools::evaluate_polynomial(d, z - 2.25);
- b = f0;
- }
- else if(z < 5.5)
- {
- // Worst case absolute error found: 1.579588208e-018
- static const T n[6] = { -0.00588219091116732271979L, 0.00434428684527812140098L, 0.00466899990542371512895L, 0.00139937567253199794533L, 0.000179205902444982389766L, 0.845033527560949509345e-5L };
- static const T d[6] = { 1, 1.07389345953392962127L, 0.470965611895885060643L, 0.105594730223366124873L, 0.0121252833787344059719L, 0.000571755036133730341579L };
- static const float f0 = 0.5494099855F;
- r = tools::evaluate_polynomial(n, z - 3.5) / tools::evaluate_polynomial(d, z - 3.5);
- b = f0;
- }
- else if(z < 9)
- {
- // Worst case absolute error found: 1.410768708e-017
- static const T n[5] = { -0.00273864253749621265032L, 0.0013089921066773026803L, 0.000775841526778089659703L, 0.000110909476102006410909L, 0.472577590124068298534e-5L };
- static const T d[6] = { 1, 0.650694792327863647878L, 0.161126734432670927888L, 0.0180081468446110640846L, 0.000767341359508884026192L, -0.287636719206664167616e-9L };
- static const float f0 = 0.5580308437F;
- r = tools::evaluate_polynomial(n, z - 5.5) / tools::evaluate_polynomial(d, z - 5.5);
- b = f0;
- }
- else if(z < 14)
- {
- // Worst case absolute error found: 1.458310511e-018
- static const T n[5] = { -0.000995856413171151859346L, 0.000320252910249376187643L, 0.000129085624923151780987L, 0.121577881306587454509e-4L, 0.33293110334156470348e-6L };
- static const T d[5] = { 1, 0.428034987547594828954L, 0.0692297359775940896439L, 0.00501515176145997560701L, 0.00013733589151338416322L };
- static const float f0 = 0.5617653728F;
- r = tools::evaluate_polynomial(n, z - 9) / tools::evaluate_polynomial(d, z - 9);
- b = f0;
- }
- else if(z < 21)
- {
- // Worst case absolute error found: 1.08182873e-019
- static const T n[5] = { -0.000395463268432048215535L, 0.91155953112698182321e-4L, 0.237451641259281193813e-4L, 0.145759953022524466816e-5L, 0.259395907606548998142e-7L };
- static const T d[5] = { 1, 0.281604524251560309285L, 0.0298468482900092392397L, 0.00141114575715338885136L, 0.251128951158576064819e-4L };
- static const float f0 = 0.5631566644F;
- r = tools::evaluate_polynomial(n, z - 14) / tools::evaluate_polynomial(d, z - 14);
- b = f0;
+ // Maximum Deviation Found: 3.702e-17
+ // Expected Error Term: 3.702e-17
+ // Maximum Relative Change in Control Points: 2.845e-04
+ // Max Error found at double precision = 4.841816e-17
+ static const T Y = 0.405935764312744140625f;
+ static const T P[] = {
+ -0.098090592216281240205L,
+ 0.178114665841120341155L,
+ 0.191003695796775433986L,
+ 0.0888900368967884466578L,
+ 0.0195049001251218801359L,
+ 0.00180424538297014223957L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.84759070983002217845L,
+ 1.42628004845511324508L,
+ 0.578052804889902404909L,
+ 0.12385097467900864233L,
+ 0.0113385233577001411017L,
+ 0.337511472483094676155e-5L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 0.5) / tools::evaluate_polynomial(Q, z - 0.5);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 2.5f)
+ {
+ // Max Error found at double precision = 6.599585e-18
+ // Maximum Deviation Found: 3.909e-18
+ // Expected Error Term: 3.909e-18
+ // Maximum Relative Change in Control Points: 9.886e-05
+ static const T Y = 0.50672817230224609375f;
+ static const T P[] = {
+ -0.0243500476207698441272L,
+ 0.0386540375035707201728L,
+ 0.04394818964209516296L,
+ 0.0175679436311802092299L,
+ 0.00323962406290842133584L,
+ 0.000235839115596880717416L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.53991494948552447182L,
+ 0.982403709157920235114L,
+ 0.325732924782444448493L,
+ 0.0563921837420478160373L,
+ 0.00410369723978904575884L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 1.5) / tools::evaluate_polynomial(Q, z - 1.5);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 4.5f)
+ {
+ // Maximum Deviation Found: 1.512e-17
+ // Expected Error Term: 1.512e-17
+ // Maximum Relative Change in Control Points: 2.222e-04
+ // Max Error found at double precision = 2.062515e-17
+ static const T Y = 0.5405750274658203125f;
+ static const T P[] = {
+ 0.00295276716530971662634L,
+ 0.0137384425896355332126L,
+ 0.00840807615555585383007L,
+ 0.00212825620914618649141L,
+ 0.000250269961544794627958L,
+ 0.113212406648847561139e-4L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.04217814166938418171L,
+ 0.442597659481563127003L,
+ 0.0958492726301061423444L,
+ 0.0105982906484876531489L,
+ 0.000479411269521714493907L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 3.5) / tools::evaluate_polynomial(Q, z - 3.5);
+ result *= exp(-z * z) / z;
}
else
{
- // Worst case absolute error found: 7.010370259e-018
- static const T n[4] = { -0.000139182098873874523526L, 0.395254617101737287826e-4L, 0.376801239136290345387e-5L, 0.629017242098850415839e-7L };
- static const T d[4] = { 1, 0.15077096006891495258L, 0.00756136203065884121997L, 0.000126226197336507576933L };
- static const float f0 = 0.5636912584F;
- r = tools::evaluate_polynomial(n, z - 21) / tools::evaluate_polynomial(d, z - 21);
- b = f0;
+ // Max Error found at double precision = 2.997958e-17
+ // Maximum Deviation Found: 2.860e-17
+ // Expected Error Term: 2.859e-17
+ // Maximum Relative Change in Control Points: 1.357e-05
+ static const T Y = 0.5579090118408203125f;
+ static const T P[] = {
+ 0.00628057170626964891937L,
+ 0.0175389834052493308818L,
+ -0.212652252872804219852L,
+ -0.687717681153649930619L,
+ -2.5518551727311523996L,
+ -3.22729451764143718517L,
+ -2.8175401114513378771L,
+ };
+ static const T Q[] = {
+ 1L,
+ 2.79257750980575282228L,
+ 11.0567237927800161565L,
+ 15.930646027911794143L,
+ 22.9367376522880577224L,
+ 13.5064170191802889145L,
+ 5.48409182238641741584L,
+ };
+ result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+ result *= exp(-z * z) / z;
}
- T g = exp(-z * z) / z;
- result = g * b + g * r;
}
else
{
@@ -360,10 +401,28 @@
}
else
{
- // Worst case absolute error found: 6.688618532e-21
- static const T n[8] = { 0.00337916709551257388990745L, -0.00073695653048167948530905L, -0.374732337392919607868241L, 0.0817442448733587196071743L, -0.0421089319936548595203468L, 0.0070165709512095756344528L, -0.00495091255982435110337458L, 0.000871646599037922480317225L };
- static const T d[8] = { 1L, -0.218088218087924645390535L, 0.412542972725442099083918L, -0.0841891147873106755410271L, 0.0655338856400241519690695L, -0.0120019604454941768171266L, 0.00408165558926174048329689L, -0.000615900721557769691924509L };
- result = z * 1.125 + z * tools::evaluate_polynomial(n, z) / tools::evaluate_polynomial(d, z);
+ // Max Error found at long double precision = 1.623299e-20
+ // Maximum Deviation Found: 4.326e-22
+ // Expected Error Term: -4.326e-22
+ // Maximum Relative Change in Control Points: 1.474e-04
+ static const T Y = 1.044948577880859375f;
+ static const T P[] = {
+ 0.0834305892146531988966L,
+ -0.338097283075565413695L,
+ -0.0509602734406067204596L,
+ -0.00904906346158537794396L,
+ -0.000489468651464798669181L,
+ -0.200305626366151877759e-4L,
+ };
+ static const T Q[] = {
+ 1L,
+ 0.455817300515875172439L,
+ 0.0916537354356241792007L,
+ 0.0102722652675910031202L,
+ 0.000650511752687851548735L,
+ 0.189532519105655496778e-4L,
+ };
+ result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
}
}
else if((z < 110) || ((z < 110) && invert)) // TODO FIXME!!!
@@ -372,134 +431,128 @@
// We'll be calculating erfc:
//
invert = !invert;
- T r, b;
- if(z < 0.75)
- {
- // Worst case absolute error found: 5.582813374e-21
- static const T n[6] = { -0.0361790390718262471360258L, 0.292251883444882683221149L, 0.281447041797604512774415L, 0.125610208862766947294894L, 0.0274135028268930549240776L, 0.00250839672168065762786937L };
- static const T d[6] = { 1L, 1.8545005897903486499845L, 1.43575803037831418074962L, 0.582827658753036572454135L, 0.124810476932949746447682L, 0.0113724176546353285778481L };
- static const float f0 = 0.3440242112F;
- r = tools::evaluate_polynomial(n, z - 0.5) / tools::evaluate_polynomial(d, z - 0.5);
- b = f0;
- }
- else if(z < 1.25)
- {
- // Worst case absolute error found: 4.01854729e-21
- static const T n[7] = { -0.0397876892611136856954425L, 0.153165212467878293257683L, 0.191260295600936245503129L, 0.10276327061989304213645L, 0.029637090615738836726027L, 0.0046093486780275489468812L, 0.000307607820348680180548455L };
- static const T d[7] = { 1L, 1.95520072987627704987886L, 1.64762317199384860109595L, 0.768238607022126250082483L, 0.209793185936509782784315L, 0.0319569316899913392596356L, 0.00213363160895785378615014L };
- static const float f0 = 0.419990927F;
- r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75);
- b = f0;
- }
- else if(z < 2.25)
- {
- // Worst case absolute error found: 2.866005373e-21
- static const T n[7] = { -0.0300838560557949717328341L, 0.0538578829844454508530552L, 0.0726211541651914182692959L, 0.0367628469888049348429018L, 0.00964629015572527529605267L, 0.00133453480075291076745275L, 0.778087599782504251917881e-4L };
- static const T d[8] = { 1L, 1.75967098147167528287343L, 1.32883571437961120556307L, 0.552528596508757581287907L, 0.133793056941332861912279L, 0.0179509645176280768640766L, 0.00104712440019937356634038L, -0.106640381820357337177643e-7L };
- static const float f0 = 0.4898625016F;
- r = tools::evaluate_polynomial(n, z - 1.25) / tools::evaluate_polynomial(d, z - 1.25);
- b = f0;
- }
- else if(z < 3.5)
- {
- // Worst case absolute error found: 1.045355789e-21
- static const T n[7] = { -0.0117907570137227847827732L, 0.014262132090538809896674L, 0.0202234435902960820020765L, 0.00930668299990432009042239L, 0.00213357802422065994322516L, 0.00025022987386460102395382L, 0.120534912219588189822126e-4L };
- static const T d[7] = { 1L, 1.50376225203620482047419L, 0.965397786204462896346934L, 0.339265230476796681555511L, 0.0689740649541569716897427L, 0.00771060262491768307365526L, 0.000371421101531069302990367L };
- static const float f0 = 0.5317370892F;
- r = tools::evaluate_polynomial(n, z - 2.25) / tools::evaluate_polynomial(d, z - 2.25);
- b = f0;
- }
- else if(z < 5.25)
- {
- // Worst case absolute error found: 8.300028706e-22
- static const T n[7] = { -0.00546954795538729307482955L, 0.00404190278731707110245394L, 0.0054963369553161170521356L, 0.00212616472603945399437862L, 0.000394984014495083900689956L, 0.365565477064442377259271e-4L, 0.135485897109932323253786e-5L };
- static const T d[8] = { 1L, 1.21019697773630784832251L, 0.620914668221143886601045L, 0.173038430661142762569515L, 0.0276550813773432047594539L, 0.00240625974424309709745382L, 0.891811817251336577241006e-4L, -0.465528836283382684461025e-11L };
- static const float f0 = 0.5489973426F;
- r = tools::evaluate_polynomial(n, z - 3.5) / tools::evaluate_polynomial(d, z - 3.5);
- b = f0;
- }
- else if(z < 8)
- {
- // Worst case absolute error found: 1.700157534e-21
- static const T n[6] = { -0.00270722535905778347999196L, 0.0013187563425029400461378L, 0.00119925933261002333923989L, 0.00027849619811344664248235L, 0.267822988218331849989363e-4L, 0.923043672315028197865066e-6L };
- static const T d[7] = { 1L, 0.814632808543141591118279L, 0.268901665856299542168425L, 0.0449877216103041118694989L, 0.00381759663320248459168994L, 0.000131571897888596914350697L, 0.404815359675764138445257e-11L };
- static const float f0 = 0.5571740866F;
- r = tools::evaluate_polynomial(n, z - 5.25) / tools::evaluate_polynomial(d, z - 5.25);
- b = f0;
- }
- else if(z < 11.5)
+ if(z < 1.5)
{
- //Worst case absolute error found: 3.002278011e-22
- static const T n[6] = { -0.00109946720691742196814323L, 0.000406425442750422675169153L, 0.000274499489416900707787024L, 0.465293770646659383436343e-4L, 0.320955425395767463401993e-5L, 0.778286018145020892261936e-7L };
- static const T d[6] = { 1L, 0.588173710611846046373373L, 0.139363331289409746077541L, 0.0166329340417083678763028L, 0.00100023921310234908642639L, 0.24254837521587225125068e-4L };
- static const float f0 = 0.5609807968F;
- r = tools::evaluate_polynomial(n, z - 8) / tools::evaluate_polynomial(d, z - 8);
- b = f0;
- }
- else if(z < 17)
- {
- //Worst case absolute error found: 6.741114695e-21
- static const T n[5] = { -0.00056907993601094962855594L, 0.000169498540373762264416984L, 0.518472354581100890120501e-4L, 0.382819312231928859704678e-5L, 0.824989931281894431781794e-7L };
- static const T d[6] = { 1L, 0.339637250051139347430323L, 0.043472647870310663055044L, 0.00248549335224637114641629L, 0.535633305337152900549536e-4L, -0.117490944405459578783846e-12L };
- static const float f0 = 0.5626493692F;
- r = tools::evaluate_polynomial(n, z - 11.5) / tools::evaluate_polynomial(d, z - 11.5);
- b = f0;
- }
- else if(z < 24)
- {
- // Worst case absolute error found: 7.802346984e-22
- static const T n[5] = { -0.000241313599483991337479091L, 0.574224975202501512365975e-4L, 0.115998962927383778460557e-4L, 0.581762134402593739370875e-6L, 0.853971555085673614607418e-8L };
- static const T d[5] = { 1L, 0.233044138299687841018015L, 0.0204186940546440312625597L, 0.000797185647564398289151125L, 0.117019281670172327758019e-4L };
- static const float f0 = 0.5634598136F;
- r = tools::evaluate_polynomial(n, z - 17) / tools::evaluate_polynomial(d, z - 17);
- b = f0;
- }
- else if(z < 38)
- {
- // Worst case absolute error found: 2.414228989e-22
- static const T n[5] = { -0.000146674699277760365803642L, 0.162666552112280519955647e-4L, 0.269116248509165239294897e-5L, 0.979584479468091935086972e-7L, 0.101994647625723465722285e-8L };
- static const T d[5] = { 1L, 0.165907812944847226546036L, 0.0103361716191505884359634L, 0.000286593026373868366935721L, 0.298401570840900340874568e-5L };
- static const float f0 = 0.5638477802F;
- r = tools::evaluate_polynomial(n, z - 24) / tools::evaluate_polynomial(d, z - 24);
- b = f0;
- }
- else if(z < 60)
- {
- // Worst case absolute error found: 5.896543869e-24
- static const T n[5] = { -0.583905797629771786720406e-4L, 0.412510325105496173512992e-5L, 0.431790922420250949096906e-6L, 0.993365155590013193345569e-8L, 0.653480510020104699270084e-10L };
- static const T d[5] = { 1L, 0.105077086072039915406159L, 0.00414278428675475620830226L, 0.726338754644523769144108e-4L, 0.477818471047398785369849e-6L };
- static const float f0 = 0.5640528202F;
- r = tools::evaluate_polynomial(n, z - 38) / tools::evaluate_polynomial(d, z - 38);
- b = f0;
- }
- else if(z < 85)
- {
- // Worst case absolute error found: 3.080612264e-21
- static const T n[4] = { -0.196457797609229579459841e-4L, 0.157243887666800692441195e-5L, 0.543902511192700878690335e-7L, 0.317472492369117710852685e-9L };
- static const T d[5] = { 1L, 0.052803989240957632204885L, 0.000926876069151753290378112L, 0.541011723226630257077328e-5L, 0.535093845803642394908747e-15L };
- static const float f0 = 0.5641309023F;
- r = tools::evaluate_polynomial(n, z - 60) / tools::evaluate_polynomial(d, z - 60);
- b = f0;
+ // Max Error found at long double precision = 3.239590e-20
+ // Maximum Deviation Found: 2.241e-20
+ // Expected Error Term: -2.241e-20
+ // Maximum Relative Change in Control Points: 5.110e-03
+ static const T Y = 0.405935764312744140625f;
+ static const T P[] = {
+ -0.0980905922162812031672L,
+ 0.159989089922969141329L,
+ 0.222359821619935712378L,
+ 0.127303921703577362312L,
+ 0.0384057530342762400273L,
+ 0.00628431160851156719325L,
+ 0.000441266654514391746428L,
+ 0.266689068336295642561e-7L,
+ };
+ static const T Q[] = {
+ 1L,
+ 2.03237474985469469291L,
+ 1.78355454954969405222L,
+ 0.867940326293760578231L,
+ 0.248025606990021698392L,
+ 0.0396649631833002269861L,
+ 0.00279220237309449026796L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 2.5)
+ {
+ // Max Error found at long double precision = 3.686211e-21
+ // Maximum Deviation Found: 1.495e-21
+ // Expected Error Term: -1.494e-21
+ // Maximum Relative Change in Control Points: 1.793e-04
+ static const T Y = 0.50672817230224609375f;
+ static const T P[] = {
+ -0.024350047620769840217L,
+ 0.0343522687935671451309L,
+ 0.0505420824305544949541L,
+ 0.0257479325917757388209L,
+ 0.00669349844190354356118L,
+ 0.00090807914416099524444L,
+ 0.515917266698050027934e-4L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.71657861671930336344L,
+ 1.26409634824280366218L,
+ 0.512371437838969015941L,
+ 0.120902623051120950935L,
+ 0.0158027197831887485261L,
+ 0.000897871370778031611439L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 4.5)
+ {
+ // Maximum Deviation Found: 1.107e-20
+ // Expected Error Term: -1.106e-20
+ // Maximum Relative Change in Control Points: 1.709e-04
+ // Max Error found at long double precision = 1.446908e-20
+ static const T Y = 0.5405750274658203125f;
+ static const T P[] = {
+ 0.0029527671653097284033L,
+ 0.0141853245895495604051L,
+ 0.0104959584626432293901L,
+ 0.00343963795976100077626L,
+ 0.00059065441194877637899L,
+ 0.523435380636174008685e-4L,
+ 0.189896043050331257262e-5L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.19352160185285642574L,
+ 0.603256964363454392857L,
+ 0.165411142458540585835L,
+ 0.0259729870946203166468L,
+ 0.00221657568292893699158L,
+ 0.804149464190309799804e-4L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 3.5f) / tools::evaluate_polynomial(Q, z - 3.5f);
+ result *= exp(-z * z) / z;
}
else
{
- // Worst case absolute error found: 8.094633491e-22
- static const T n[4] = { -0.789224703978722689089794e-5L, 0.622088451660986955124162e-6L, 0.145728445676882396797184e-7L, 0.603715505542715364529243e-10L };
- static const T d[4] = { 1L, 0.0375328846356293715248719L, 0.000467919535974625308126054L, 0.193847039275845656900547e-5L };
- static const float f0 = 0.5641584396F;
- r = tools::evaluate_polynomial(n, z - 85) / tools::evaluate_polynomial(d, z - 85);
- b = f0;
- }
- T g = exp(-z * z) / z;
- result = g * b + g * r;
- BOOST_MATH_INSTRUMENT_CODE("r = " << r);
- BOOST_MATH_INSTRUMENT_CODE("b = " << b);
- BOOST_MATH_INSTRUMENT_CODE("g = " << g);
+ // Max Error found at long double precision = 7.961166e-21
+ // Maximum Deviation Found: 6.677e-21
+ // Expected Error Term: 6.676e-21
+ // Maximum Relative Change in Control Points: 2.319e-05
+ static const T Y = 0.55825519561767578125f;
+ static const T P[] = {
+ 0.00593438793008050214106L,
+ 0.0280666231009089713937L,
+ -0.141597835204583050043L,
+ -0.978088201154300548842L,
+ -5.47351527796012049443L,
+ -13.8677304660245326627L,
+ -27.1274948720539821722L,
+ -29.2545152747009461519L,
+ -16.8865774499799676937L,
+ };
+ static const T Q[] = {
+ 1L,
+ 4.72948911186645394541L,
+ 23.6750543147695749212L,
+ 60.0021517335693186785L,
+ 131.766251645149522868L,
+ 178.167924971283482513L,
+ 182.499390505915222699L,
+ 104.365251479578577989L,
+ 30.8365511891224291717L,
+ };
+ result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+ result *= exp(-z * z) / z;
+ }
}
else
{
//
- // Any value of z larger than 28 will underflow to zero:
+ // Any value of z larger than 110 will underflow to zero:
//
result = 0;
invert = !invert;
@@ -553,184 +606,352 @@
}
else
{
- // Worst case absolute error found: 1.928180863e-35
- static const T n[13] = { 0.0033791670955125738961589031215451706772L, -0.000356604747854533671135323429762519216044L, -0.374476838669183581687167228866769133591L, 0.0395338132469809122364498388174446488042L, -0.070405473508786506375820161461872523315L, 0.00575264725772369303419496752516485264994L, -0.0122324470706306942925087773122510971344L, 0.000982833333252586078523570049842642796291L, -0.000937806155615159592441487275938040285833L, 0.485407838108763091860415874932955355755e-4L, -0.50171236926234625577876479444632561922e-4L, 0.19406068817888598455243350289053451571e-5L, -0.119351103792049576459000102632508734863e-5L };
- static const T d[13] = { 1L, -0.105530368216503232473476334993759958083L, 0.488152943026846232046726653294817930988L, -0.0470361716364117780901924633553851211874L, 0.107663671943702835026199580597519084906L, -0.00919493879447389180633447493128337242362L, 0.0138231121717229362691899919242806829805L, -0.000994048559663865788847688218108232247441L, 0.00109769834527023265969224251892094019735L, -0.600458401801636062015615549258555311545e-4L, 0.51530723974502946291624848874654212384e-4L, -0.164121264470361558910636548509486296153e-5L, 0.112643498977070218963888579607359294396e-5L };
-
- result = z * 1.125 + z * tools::evaluate_rational(n, d, z);
+ // Max Error found at long double precision = 2.342380e-35
+ // Maximum Deviation Found: 6.124e-36
+ // Expected Error Term: -6.124e-36
+ // Maximum Relative Change in Control Points: 3.492e-10
+ static const T Y = 1.0841522216796875f;
+ static const T P[] = {
+ 0.0442269454158250738961589031215451778L,
+ -0.35549265736002144875335323556961233L,
+ -0.0582179564566667896225454670863270393L,
+ -0.0112694696904802304229950538453123925L,
+ -0.000805730648981801146251825329609079099L,
+ -0.566304966591936566229702842075966273e-4L,
+ -0.169655010425186987820201021510002265e-5L,
+ -0.344448249920445916714548295433198544e-7L,
+ };
+ static const T Q[] = {
+ 1L,
+ 0.466542092785657604666906909196052522L,
+ 0.100005087012526447295176964142107611L,
+ 0.0128341535890117646540050072234142603L,
+ 0.00107150448466867929159660677016658186L,
+ 0.586168368028999183607733369248338474e-4L,
+ 0.196230608502104324965623171516808796e-5L,
+ 0.313388521582925207734229967907890146e-7L,
+ };
+ result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
}
}
- else if((z < 9) || ((z < 110) && invert)) // TODO FIXME!!
+ else if((z < 9) || ((z < 110) && invert))
{
//
// We'll be calculating erfc:
//
invert = !invert;
- T r, b;
- if(z < 0.75)
+ if(z < 1)
+ {
+ // Max Error found at long double precision = 3.246278e-35
+ // Maximum Deviation Found: 1.388e-35
+ // Expected Error Term: 1.387e-35
+ // Maximum Relative Change in Control Points: 6.127e-05
+ static const T Y = 0.371877193450927734375f;
+ static const T P[] = {
+ -0.0640320213544647969396032886581290455L,
+ 0.200769874440155895637857443946706731L,
+ 0.378447199873537170666487408805779826L,
+ 0.30521399466465939450398642044975127L,
+ 0.146890026406815277906781824723458196L,
+ 0.0464837937749539978247589252732769567L,
+ 0.00987895759019540115099100165904822903L,
+ 0.00137507575429025512038051025154301132L,
+ 0.0001144764551085935580772512359680516L,
+ 0.436544865032836914773944382339900079e-5L,
+ };
+ static const T Q[] = {
+ 1L,
+ 2.47651182872457465043733800302427977L,
+ 2.78706486002517996428836400245547955L,
+ 1.87295924621659627926365005293130693L,
+ 0.829375825174365625428280908787261065L,
+ 0.251334771307848291593780143950311514L,
+ 0.0522110268876176186719436765734722473L,
+ 0.00718332151250963182233267040106902368L,
+ 0.000595279058621482041084986219276392459L,
+ 0.226988669466501655990637599399326874e-4L,
+ 0.270666232259029102353426738909226413e-10L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 1.5)
+ {
+ // Max Error found at long double precision = 2.215785e-35
+ // Maximum Deviation Found: 1.539e-35
+ // Expected Error Term: 1.538e-35
+ // Maximum Relative Change in Control Points: 6.104e-05
+ static const T Y = 0.45658016204833984375f;
+ static const T P[] = {
+ -0.0289965858925328393392496555094848345L,
+ 0.0868181194868601184627743162571779226L,
+ 0.169373435121178901746317404936356745L,
+ 0.13350446515949251201104889028133486L,
+ 0.0617447837290183627136837688446313313L,
+ 0.0185618495228251406703152962489700468L,
+ 0.00371949406491883508764162050169531013L,
+ 0.000485121708792921297742105775823900772L,
+ 0.376494706741453489892108068231400061e-4L,
+ 0.133166058052466262415271732172490045e-5L,
+ };
+ static const T Q[] = {
+ 1L,
+ 2.32970330146503867261275580968135126L,
+ 2.46325715420422771961250513514928746L,
+ 1.55307882560757679068505047390857842L,
+ 0.644274289865972449441174485441409076L,
+ 0.182609091063258208068606847453955649L,
+ 0.0354171651271241474946129665801606795L,
+ 0.00454060370165285246451879969534083997L,
+ 0.000349871943711566546821198612518656486L,
+ 0.123749319840299552925421880481085392e-4L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 1.0f) / tools::evaluate_polynomial(Q, z - 1.0f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 2.25)
+ {
+ // Maximum Deviation Found: 1.418e-35
+ // Expected Error Term: 1.418e-35
+ // Maximum Relative Change in Control Points: 1.316e-04
+ // Max Error found at long double precision = 1.998462e-35
+ static const T Y = 0.50250148773193359375f;
+ static const T P[] = {
+ -0.0201233630504573402185161184151016606L,
+ 0.0331864357574860196516686996302305002L,
+ 0.0716562720864787193337475444413405461L,
+ 0.0545835322082103985114927569724880658L,
+ 0.0236692635189696678976549720784989593L,
+ 0.00656970902163248872837262539337601845L,
+ 0.00120282643299089441390490459256235021L,
+ 0.000142123229065182650020762792081622986L,
+ 0.991531438367015135346716277792989347e-5L,
+ 0.312857043762117596999398067153076051e-6L,
+ };
+ static const T Q[] = {
+ 1L,
+ 2.13506082409097783827103424943508554L,
+ 2.06399257267556230937723190496806215L,
+ 1.18678481279932541314830499880691109L,
+ 0.447733186643051752513538142316799562L,
+ 0.11505680005657879437196953047542148L,
+ 0.020163993632192726170219663831914034L,
+ 0.00232708971840141388847728782209730585L,
+ 0.000160733201627963528519726484608224112L,
+ 0.507158721790721802724402992033269266e-5L,
+ 0.18647774409821470950544212696270639e-12L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if (z < 3)
+ {
+ // Maximum Deviation Found: 3.575e-36
+ // Expected Error Term: 3.575e-36
+ // Maximum Relative Change in Control Points: 7.103e-05
+ // Max Error found at long double precision = 5.794737e-36
+ static const T Y = 0.52896785736083984375f;
+ static const T P[] = {
+ -0.00902152521745813634562524098263360074L,
+ 0.0145207142776691539346923710537580927L,
+ 0.0301681239582193983824211995978678571L,
+ 0.0215548540823305814379020678660434461L,
+ 0.00864683476267958365678294164340749949L,
+ 0.00219693096885585491739823283511049902L,
+ 0.000364961639163319762492184502159894371L,
+ 0.388174251026723752769264051548703059e-4L,
+ 0.241918026931789436000532513553594321e-5L,
+ 0.676586625472423508158937481943649258e-7L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.93669171363907292305550231764920001L,
+ 1.69468476144051356810672506101377494L,
+ 0.880023580986436640372794392579985511L,
+ 0.299099106711315090710836273697708402L,
+ 0.0690593962363545715997445583603382337L,
+ 0.0108427016361318921960863149875360222L,
+ 0.00111747247208044534520499324234317695L,
+ 0.686843205749767250666787987163701209e-4L,
+ 0.192093541425429248675532015101904262e-5L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 2.25f) / tools::evaluate_polynomial(Q, z - 2.25f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 3.5)
+ {
+ // Maximum Deviation Found: 8.126e-37
+ // Expected Error Term: -8.126e-37
+ // Maximum Relative Change in Control Points: 1.363e-04
+ // Max Error found at long double precision = 1.747062e-36
+ static const T Y = 0.54037380218505859375f;
+ static const T P[] = {
+ -0.0033703486408887424921155540591370375L,
+ 0.0104948043110005245215286678898115811L,
+ 0.0148530118504000311502310457390417795L,
+ 0.00816693029245443090102738825536188916L,
+ 0.00249716579989140882491939681805594585L,
+ 0.0004655591010047353023978045800916647L,
+ 0.531129557920045295895085236636025323e-4L,
+ 0.343526765122727069515775194111741049e-5L,
+ 0.971120407556888763695313774578711839e-7L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.59911256167540354915906501335919317L,
+ 1.136006830764025173864831382946934L,
+ 0.468565867990030871678574840738423023L,
+ 0.122821824954470343413956476900662236L,
+ 0.0209670914950115943338996513330141633L,
+ 0.00227845718243186165620199012883547257L,
+ 0.000144243326443913171313947613547085553L,
+ 0.407763415954267700941230249989140046e-5L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 3.0f) / tools::evaluate_polynomial(Q, z - 3.0f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 5.5)
+ {
+ // Maximum Deviation Found: 5.804e-36
+ // Expected Error Term: -5.803e-36
+ // Maximum Relative Change in Control Points: 2.475e-05
+ // Max Error found at long double precision = 1.349545e-35
+ static const T Y = 0.55000019073486328125f;
+ static const T P[] = {
+ 0.00118142849742309772151454518093813615L,
+ 0.0072201822885703318172366893469382745L,
+ 0.0078782276276860110721875733778481505L,
+ 0.00418229166204362376187593976656261146L,
+ 0.00134198400587769200074194304298642705L,
+ 0.000283210387078004063264777611497435572L,
+ 0.405687064094911866569295610914844928e-4L,
+ 0.39348283801568113807887364414008292e-5L,
+ 0.248798540917787001526976889284624449e-6L,
+ 0.929502490223452372919607105387474751e-8L,
+ 0.156161469668275442569286723236274457e-9L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.52955245103668419479878456656709381L,
+ 1.06263944820093830054635017117417064L,
+ 0.441684612681607364321013134378316463L,
+ 0.121665258426166960049773715928906382L,
+ 0.0232134512374747691424978642874321434L,
+ 0.00310778180686296328582860464875562636L,
+ 0.000288361770756174705123674838640161693L,
+ 0.177529187194133944622193191942300132e-4L,
+ 0.655068544833064069223029299070876623e-6L,
+ 0.11005507545746069573608988651927452e-7L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 4.5f) / tools::evaluate_polynomial(Q, z - 4.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 7.5)
+ {
+ // Maximum Deviation Found: 1.007e-36
+ // Expected Error Term: 1.007e-36
+ // Maximum Relative Change in Control Points: 1.027e-03
+ // Max Error found at long double precision = 2.646420e-36
+ static const T Y = 0.5574436187744140625f;
+ static const T P[] = {
+ 0.000293236907400849056269309713064107674L,
+ 0.00225110719535060642692275221961480162L,
+ 0.00190984458121502831421717207849429799L,
+ 0.000747757733460111743833929141001680706L,
+ 0.000170663175280949889583158597373928096L,
+ 0.246441188958013822253071608197514058e-4L,
+ 0.229818000860544644974205957895688106e-5L,
+ 0.134886977703388748488480980637704864e-6L,
+ 0.454764611880548962757125070106650958e-8L,
+ 0.673002744115866600294723141176820155e-10L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.12843690320861239631195353379313367L,
+ 0.569900657061622955362493442186537259L,
+ 0.169094404206844928112348730277514273L,
+ 0.0324887449084220415058158657252147063L,
+ 0.00419252877436825753042680842608219552L,
+ 0.00036344133176118603523976748563178578L,
+ 0.204123895931375107397698245752850347e-4L,
+ 0.674128352521481412232785122943508729e-6L,
+ 0.997637501418963696542159244436245077e-8L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z - 6.5f) / tools::evaluate_polynomial(Q, z - 6.5f);
+ result *= exp(-z * z) / z;
+ }
+ else if(z < 11.5)
{
- // Worst case absolute error found: 9.46579566e-36
- static const T n[10] = { -0.0361790390718262471349157886581290316118L, 0.268773785250238404882137450640472787307L, 0.46350995084084251624649426251701042395L, 0.368375435727102373204587584306335625665L, 0.177618123820303858190236222513516291818L, 0.0566304173556669007529719743050764079095L, 0.0121631149481817424284077180037019529004L, 0.00171397353209314111395429418066990259845L, 0.000144662387395699594624184141956722488753L, 0.559870522050008635715382724858714587198e-5L };
- static const T d[10] = { 1L, 2.50344259590701770420935329380375393716L, 2.84905597172139276093882199286535521011L, 1.93691730181297099541395314232750876411L, 0.868059574796050528229446630538462280596L, 0.266360035323208212078527036132085926692L, 0.0560555526482963925944703505114360693216L, 0.0078174400311465420803366235814673576269L, 0.000657067309046405057499687417839930873806L, 0.254293850077789079098316521097979388983e-4L };
- static const float f0 = 0.3440242112F;
- r = tools::evaluate_rational(n, d, z - 0.5);
- b = f0;
- }
- else if(z < 1.25)
- {
- // Worst case absolute error found: 1.222145602e-35
- static const T n[10] = { -0.03978768926111368569548863384587917014L, 0.136218360681765349252731304877153919181L, 0.252782160406474440925641829129129001834L, 0.198264231106182362320012632943145619752L, 0.0923045825293507328801206570363391760624L, 0.0281157216148097885766639832985410722743L, 0.00573041663561645197870019701493117161792L, 0.000762341440133027349203518836487137709687L, 0.60471020134417423449877859375492618899e-4L, 0.219005333943510376644902615714724932217e-5L };
- static const T d[11] = { 1L, 2.38113277319993574121349184069891082204L, 2.57380422881476860215664207822277590181L, 1.65937045609044738941173490190122101824L, 0.704055811320312044285417250966993014161L, 0.20414913933328592198279939394283925451L, 0.0405162285360227740710964820549709038107L, 0.00531638867177288975915820230980317499728L, 0.000419364368135139398723983192742319455284L, 0.151874665979234971229096136924566078234e-4L, 0.807869459506748684117962248796937508011e-11L };
- static const float f0 = 0.419990927F;
- r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75);
- b = f0;
- }
- else if(z < 2)
- {
- // Worst case absolute error found: 5.893842955e-36
- static const T n[11] = { -0.0255063683486569102096736247449691465143L, 0.045782379672906795594927072060091308408L, 0.113248439610400562258072020811195716817L, 0.0996016254422112410086711272219455446695L, 0.0508749250027894453228337309651895478017L, 0.0171081937013828309576540212196644542209L, 0.00395354196550210630440706013523069756354L, 0.000629022203390154585475081628606234279007L, 0.664903286194855400689101617763591823345e-4L, 0.423935693893425355108071655059640137208e-5L, 0.124304036910852727351487636048151737214e-6L };
- static const T d[11] = { 1L, 2.39207679390801118396945702674440915308L, 2.62237869972649377524874287442154430843L, 1.73645189911172997548091140085423843725L, 0.769812706091926741262865732006953282036L, 0.238986814887891859065369830215615790694L, 0.0526759147274379214313767032352419949829L, 0.00814993801398361741777997755108018659382L, 0.00084829993036396244429607826663252633817L, 0.537276435448416921594616417908288527881e-4L, 0.157537193656690184073389824392669625417e-5L };
- static const float f0 = 0.4852850139F;
- r = tools::evaluate_rational(n, d, z - 1.25);
- b = f0;
- }
- else if(z < 2.75)
- {
- // Worst case absolute error found: 4.024770853e-36
- static const T n[10] = { -0.0108897177067473013323228381829144739013L, 0.0202210475357865979950082670101965963435L, 0.0403242149078991892599316678797913295452L, 0.0288492313188655084113941326565482276741L, 0.0116982647742533555237861890442866083526L, 0.00301908913020213436098518520436147263177L, 0.000511140165864993121203730804407968689429L, 0.555507897975436549741754647662158917103e-4L, 0.354571088276496768574495922471690102061e-5L, 0.101789333060641482357520298518780163915e-6L };
- static const T d[11] = { 1L, 1.98184247277299581801610266575921553552L, 1.77518826058984218945272444617044495028L, 0.943934346956188464279312722940302202684L, 0.328630002685235061519039528479761588138L, 0.0777535542351039388345270792222146705189L, 0.0125143974181120800829065248546370953609L, 0.00132270605931460450441108147393979771563L, 0.834118048375722123506409257130329786209e-4L, 0.239456257167492104073765911366304033453e-5L, 0.197067742893785800814802969598122120027e-13L };
- static const float f0 = 0.5216810703F;
- r = tools::evaluate_polynomial(n, z - 2) / tools::evaluate_polynomial(d, z - 2);
- b = f0;
- }
- else if(z < 3.75)
- {
- // Worst case absolute error found: 2.119317982e-36
- static const T n[10] = { -0.00669534469911386821762042893742722704704L, 0.00779239067529714323524154862288379264056L, 0.0167670669587476509267036865033136655094L, 0.0113887916348251443051357686146040093464L, 0.00426976750247047946700539147728477144579L, 0.00100469100574064832665606356894416652764L, 0.000153533145320881108157829902752192859922L, 0.149337551064665413462766906201269176262e-4L, 0.846377837919513024118176704010972579138e-6L, 0.214045773545256889299689737489755489478e-7L };
- static const T d[11] = { 1L, 1.78724215851193637544287795626580411105L, 1.44052576962222794702178612920219772782L, 0.687639905608366315245841860669607532265L, 0.214374367225822611754443822822738563207L, 0.0452948320968245754796139856381504201504L, 0.00649108394178118005887118777181540680812L, 0.000608904720665003139414993591868256489088L, 0.33959064390570911588709483563995284603e-4L, 0.858811916085393051026834431509997486704e-6L, 0.618878592093890510233502654703683447468e-15L };
- static const float f0 = 0.5392661095F;
- r = tools::evaluate_polynomial(n, z - 2.75) / tools::evaluate_polynomial(d, z - 2.75);
- b = f0;
- }
- else if(z < 5)
- {
- // Worst case absolute error found: 3.131502824e-36
- static const T n[10] = { -0.00378088626017041998499190989910098718437L, 0.0029008905292996011997575492874095588328L, 0.00662431938391549599757244232386689480515L, 0.00417809740399798845564363621020984935218L, 0.00142019749135652688012034919213168974543L, 0.000299107980170253223293530710056814995102L, 0.405161538841561396150507786831930770839e-4L, 0.346344371670880861875666253626975779945e-5L, 0.171091054330494778613793054233437928605e-6L, 0.373924433717749484258186454458704819755e-8L };
- static const T d[10] = { 1L, 1.5810750672399887547849540367499389454L, 1.12479852885403050655678225945856872694L, 0.47277272679268851560291322980574597267L, 0.129444913616967588584693095240544707208L, 0.0239544490709674941887988416318107990646L, 0.00299775294496053944060700963645084591246L, 0.00024478412843088575835960648397300177201L, 0.118424712755145205302405346348931402917e-4L, 0.2588206250858483868392167535345681119e-6L };
- static const float f0 = 0.549742341F;
- r = tools::evaluate_rational(n, d, z - 3.75);
- b = f0;
- }
- else if(z < 6.5)
- {
- // Worst case absolute error found: 3.352877573e-35
- static const T n[9] = { -0.00210683958249012010206456754123471415706L, 0.00146329021314062287670019911742786780092L, 0.00242029064025351202243048169807220157512L, 0.0011321990764681390160708960047630195582L, 0.000277123780759366982673218537550876769487L, 0.401236501288775561636453586216146028714e-4L, 0.347799938323835778216424009916867086167e-5L, 0.167678812729975154456097184107934455429e-6L, 0.346722083660429051057284107535869165775e-8L };
- static const T d[10] = { 1L, 1.22334833521124956366395053444841951468L, 0.661433457507589455018784737495591428263L, 0.206503622658778732280997770028712044451L, 0.0407323027711252752353388616742333806362L, 0.00519969889874877079704615143005539754407L, 0.000419679230772141031030427156828631265963L, 0.195896640254171597965013007459411704085e-4L, 0.405070207572551760879797790899826058473e-6L, -0.949400883467250846930389103621356900319e-17L };
- static const float f0 = 0.5556300282F;
- r = tools::evaluate_polynomial(n, z - 5) / tools::evaluate_polynomial(d, z - 5);
- b = f0;
- }
- else if(z < 8)
- {
- // Worst case absolute error found: 2.10254551e-36
- static const T n[9] = { -0.00107224589975857477185569028693588970638L, 0.00081159959093417892787006088639848404179L, 0.00105587689576932891666032146026668833287L, 0.000416243954540394829165805666548948771809L, 0.861189599093384016322579078144012057531e-4L, 0.105064862265647286966467927732505059558e-4L, 0.763674245263385902692134637353517251296e-6L, 0.306964079269190247526141442183490066292e-7L, 0.525762928524110511201313708396204710874e-9L };
- static const T d[9] = { 1L, 1.03391233015873996503551085347368889767L, 0.471295765635712684194436077437130977978L, 0.123736066749315618886080242926593910851L, 0.0204690897886919138685460664198600282119L, 0.00218521816743913946855947853274936296576L, 0.000147057386621380823003258590658747813774L, 0.570514093706434168568509838021466564264e-5L, 0.977166974697066620826028345712327325748e-7L };
- static const float f0 = 0.5588091016F;
- r = tools::evaluate_rational(n, d, z - 6.5);
- b = f0;
- }
- else if(z < 10)
- {
- // Worst case absolute error found: 8.006848023e-37
- static const T n[9] = { -0.000764310289345400483607004463638641680792L, 0.000375959928342415987920641866513058424701L, 0.000477297615927227258961911005347799033473L, 0.000165573027183931250708334564452626717999L, 0.296617493157889183515359094859055088531e-4L, 0.310848843632513098624143817615592253324e-5L, 0.192932185180269434271810693046412848027e-6L, 0.658630702075882625552897504189436319318e-8L, 0.952880249971934343233104137698620618851e-10L };
- static const T d[9] = { 1L, 0.885953297549629585611202995885653291163L, 0.345435500008639080844920938390800739845L, 0.0774289910213823638414558561872084410517L, 0.0109142290775434200609901181618181718948L, 0.000990815749814504617347317658063197107511L, 0.565801964505889288566864514277149126893e-4L, 0.185846832474795168475989869562411740416e-5L, 0.268875677412705938842038805073576012915e-7L };
- static const float f0 = 0.5606456399F;
- r = tools::evaluate_rational(n, d, z - 8);
- b = f0;
- }
- else if(z < 12.5)
- {
- // Worst case absolute error found: 1.920487881e-37
- static const T n[9] = { -0.00049563543942917072170112779531688430711L, 0.000181627479476470045833245371263435213396L, 0.000205326116055098869964168775605689851547L, 0.602131211119027087508128340443602490578e-4L, 0.904046610725767796892834734953040318357e-5L, 0.790026495730279360618915285828083295597e-6L, 0.407154634490633661408148126521656550974e-7L, 0.114947371191075676274146563385045083492e-8L, 0.136959786076422777905476122283384578138e-10L };
- static const T d[9] = { 1L, 0.738936412939629363226408445024616124644L, 0.239911614287295587917806937612822134282L, 0.0447038494121568822243673246874110377585L, 0.00522920850086320874490830611785659031521L, 0.000393238307986133717620560371515858231357L, 0.18566831172454022627187937328433090172e-4L, 0.503267648562793696253090536560738864711e-6L, 0.599643373938283798258195761814169705225e-8L };
- static const float f0 = 0.5619055629F;
- r = tools::evaluate_rational(n, d, z - 10);
- b = f0;
- }
- else if(z < 15.5)
- {
- // Worst case absolute error found: 1.257535398e-36
- static const T n[8] = { -0.000310716603972278501158850534578560332617L, 0.00011678083970731888844953949114132723885L, 0.800190190430061077253903477346761222901e-4L, 0.156297703728913451581100601056534652076e-4L, 0.151432948873760306237990776980248434808e-5L, 0.798843820137664551611401320879346402013e-7L, 0.219398460602279142004550137383605410162e-8L, 0.244774638611231694720102050826438123042e-10L };
- static const T d[9] = { 1L, 0.536396467095662522242258924420035847546L, 0.12368534422025248132153213926057643819L, 0.0158935327926481775020421628983323726346L, 0.00122923842249710594941031294559763761829L, 0.572250268558063378795115723535491980973e-4L, 0.148480028895780193980507551658637328542e-5L, 0.165653602646537759781637790321962585489e-7L, 0.159564067829807076335030582566349996674e-21L };
- static const float f0 = 0.5627119541F;
- r = tools::evaluate_polynomial(n, z - 12.5) / tools::evaluate_polynomial(d, z - 12.5);
- b = f0;
- }
- else if(z < 20)
- {
- // Worst case absolute error found: 8.757869781e-37
- static const T n[8] = { -0.000232165799411403604511899928682939571679L, 0.473280617901953093430938763139716601257e-4L, 0.322197983787538821545512424866901113822e-4L, 0.53341606003892328294004958671866936958e-5L, 0.429392533762301420884052643595872389588e-6L, 0.186707830002841677949013638707964861935e-7L, 0.420531779988891521855765212230340499168e-9L, 0.38311153882270641561622347190918577449e-11L };
- static const T d[9] = { 1L, 0.440698415779467873589164469370853517393L, 0.0834079018864546179293964148377602724235L, 0.00878845776227345123101968908701592510214L, 0.000556792236817609981676018545344816931862L, 0.212109697877283363015270621439889468055e-4L, 0.449886743173619367191275104721917344569e-6L, 0.409854405005280582884236774230760402868e-8L, 0.185071069100878210011727114952196630136e-23L };
- static const float f0 = 0.5632548332F;
- r = tools::evaluate_polynomial(n, z - 15.5) / tools::evaluate_polynomial(d, z - 15.5);
- b = f0;
- }
- else if(z < 26)
- {
- // Worst case absolute error found: 6.571456853e-37
- static const T n[8] = { -0.000143129243915019877016909310584833592972L, 0.203555078616578752774553439209122798188e-4L, 0.116664173863253360297276766667046754002e-4L, 0.153607643549360281089355497044679566329e-5L, 0.976486412462035616905934994399021437658e-7L, 0.33416295362298759817564985788402188186e-8L, 0.590812660432887787644458409396105030781e-10L, 0.421471133628743008309458424282421874073e-12L };
- static const T d[8] = { 1L, 0.346848503043261151874606241871432055165L, 0.0516236057301875770334953837585483539519L, 0.0042740199483843978391191804633398315544L, 0.000212586135154357046245694533825075007902L, 0.635258188334431428038271842518161892901e-5L, 0.105600415847309067601860633252823051505e-6L, 0.753327238218310630683194279382721367469e-9L };
- static const float f0 = 0.5636301041F;
- r = tools::evaluate_polynomial(n, z - 20) / tools::evaluate_polynomial(d, z - 20);
- b = f0;
- }
- else if(z < 34)
- {
- // Worst case absolute error found: 2.529847851e-38
- static const T n[8] = { -0.863162280463127451272161303767299107489e-4L, 0.876391266908317792353253474461536821127e-5L, 0.407624004719762912183133597708988715137e-5L, 0.418371930808379615690956857824418800194e-6L, 0.206376792034344913360458422974245946873e-7L, 0.546878311460852031076829190722479684e-9L, 0.74760305098509923829341087347518769626e-11L, 0.411832046806658129073165530819472782663e-13L };
- static const T d[9] = { 1L, 0.268714425336129161136892060316889824437L, 0.0309686025544536788982104017649851516554L, 0.00198426609900573235086828549632290797514L, 0.763402107420631006526499294753645914808e-4L, 0.176354119917411263184292422389304506735e-5L, 0.226504369338582253171306523992412653547e-7L, 0.124774448034213281307712525982862926228e-9L, 0.541581693417048102342931921476367282462e-28L };
- static const float f0 = 0.5638595223F;
- r = tools::evaluate_polynomial(n, z - 26) / tools::evaluate_polynomial(d, z - 26);
- b = f0;
- }
- else if(z < 46)
- {
- // Worst case absolute error found: 4.880939969e-36
- static const T n[7] = { -0.552701065525823960321509114250962372959e-4L, 0.459798238451342341380837226471283129117e-5L, 0.117462487430397988459985818542278619172e-5L, 0.704965053290383355071079647747711714696e-7L, 0.191250482739845326510859812663255140286e-8L, 0.245913419605775857221974833985059356312e-10L, 0.120466123381598216554945834019040289508e-12L };
- static const T d[8] = { 1L, 0.175852634004021068883919193145643406305L, 0.0128923775281953424344037205817061020944L, 0.000504385604406829419856756492198778141939L, 0.111061123996047697713589874603922096536e-4L, 0.130499191758882778721577274330215931713e-6L, 0.639279131688964342759306361922751772829e-9L, -0.10503012469906917938721140272475203112e-26L };
- static const float f0 = 0.564001143F;
- r = tools::evaluate_polynomial(n, z - 34) / tools::evaluate_polynomial(d, z - 34);
- b = f0;
- }
- else if(z < 62)
- {
- // Worst case absolute error found: 8.634961697e-37
- static const T n[7] = { -0.299551937061912926289705712581858190494e-4L, 0.188783643501597286680709990049243153264e-5L, 0.353403900815908094401075506391935032602e-6L, 0.156779149815977342177830075875441013335e-7L, 0.31456307256618424444443489905810774564e-9L, 0.29912544832845265905351204765862702307e-11L, 0.108360038290329929702659864116372774364e-13L };
- static const T d[7] = { 1L, 0.13020345740128026085404079010013007005L, 0.00706598549835633149505923515407700416484L, 0.000204578844720147762058725212299415091588L, 0.333280847765611305843836201217690887394e-5L, 0.289666303356425428524772241712503072453e-7L, 0.104933404691983708511908027657434756019e-9L };
- static const float f0 = 0.564086318F;
- r = tools::evaluate_polynomial(n, z - 46) / tools::evaluate_polynomial(d, z - 46);
- b = f0;
- }
- else if(z < 80)
- {
- // Worst case absolute error found: 6.10700166e-39
- static const T n[7] = { -0.146162619512779884168960178459825270831e-4L, 0.952303834226435420147786300516273108344e-6L, 0.114559890033396819882468040960469980168e-6L, 0.368994353517438072494331750992706039868e-8L, 0.545158660461829568567388818070830588533e-10L, 0.383569328729331398089160922704933035281e-12L, 0.103104113324271568678309317039606225294e-14L };
- static const T d[7] = { 1L, 0.0966822058944670599111360985490948186413L, 0.00389546596914826027568119425228340291579L, 0.837237328321088890541798035513749762375e-4L, 0.101236677684940943809478588316884539423e-5L, 0.652985528810044575089151077574382356898e-8L, 0.175523663960756825512727785416137345625e-10L };
- static const float f0 = 0.5641308427F;
- r = tools::evaluate_polynomial(n, z - 62) / tools::evaluate_polynomial(d, z - 62);
- b = f0;
+ // Maximum Deviation Found: 8.380e-36
+ // Expected Error Term: 8.380e-36
+ // Maximum Relative Change in Control Points: 2.632e-06
+ // Max Error found at long double precision = 9.849522e-36
+ static const T Y = 0.56083202362060546875f;
+ static const T P[] = {
+ 0.000282420728751494363613829834891390121L,
+ 0.00175387065018002823433704079355125161L,
+ 0.0021344978564889819420775336322920375L,
+ 0.00124151356560137532655039683963075661L,
+ 0.000423600733566948018555157026862139644L,
+ 0.914030340865175237133613697319509698e-4L,
+ 0.126999927156823363353809747017945494e-4L,
+ 0.110610959842869849776179749369376402e-5L,
+ 0.55075079477173482096725348704634529e-7L,
+ 0.119735694018906705225870691331543806e-8L,
+ };
+ static const T Q[] = {
+ 1L,
+ 1.69889613396167354566098060039549882L,
+ 1.28824647372749624464956031163282674L,
+ 0.572297795434934493541628008224078717L,
+ 0.164157697425571712377043857240773164L,
+ 0.0315311145224594430281219516531649562L,
+ 0.00405588922155632380812945849777127458L,
+ 0.000336929033691445666232029762868642417L,
+ 0.164033049810404773469413526427932109e-4L,
+ 0.356615210500531410114914617294694857e-6L,
+ };
+ result = Y + tools::evaluate_polynomial(P, z / 2 - 4.75f) / tools::evaluate_polynomial(Q, z / 2 - 4.75f);
+ result *= exp(-z * z) / z;
}
else
{
- // Worst case absolute error found: 3.409761622e-39
- static const T n[7] = { -0.103600733768855845278685109220598569282e-4L, 0.324846376725138276091774803419773168357e-6L, 0.376580042454826796817160322889753111347e-7L, 0.971125540805193703472871419997820785081e-9L, 0.112499588573918670534994853431614458564e-10L, 0.6161310085325929027114924903443564594e-13L, 0.128358125353302055468778305481781957985e-15L };
- static const T d[7] = { 1L, 0.0749579802028981336679035635535776767532L, 0.00234137768051079846068630120744888560113L, 0.390095348825177571970348898222511084593e-4L, 0.365628111271063883320224395719085602867e-6L, 0.182790705425846241876560215158605843026e-8L, 0.380806548535012144489621218246876843627e-11L };
- static const float f0 = 0.5641558766F;
- r = tools::evaluate_polynomial(n, z - 80) / tools::evaluate_polynomial(d, z - 80);
- b = f0;
+ // Maximum Deviation Found: 1.132e-35
+ // Expected Error Term: -1.132e-35
+ // Maximum Relative Change in Control Points: 4.674e-04
+ // Max Error found at long double precision = 1.162590e-35
+ static const T Y = 0.5632686614990234375f;
+ static const T P[] = {
+ 0.000920922048732849448079451574171836943L,
+ 0.00321439044532288750501700028748922439L,
+ -0.250455263029390118657884864261823431L,
+ -0.906807635364090342031792404764598142L,
+ -8.92233572835991735876688745989985565L,
+ -21.7797433494422564811782116907878495L,
+ -91.1451915251976354349734589601171659L,
+ -144.1279109655993927069052125017673L,
+ -313.845076581796338665519022313775589L,
+ -273.11378811923343424081101235736475L,
+ -271.651566205951067025696102600443452L,
+ -60.0530577077238079968843307523245547L,
+ };
+ static const T Q[] = {
+ 1L,
+ 3.49040448075464744191022350947892036L,
+ 34.3563592467165971295915749548313227L,
+ 84.4993232033879023178285731843850461L,
+ 376.005865281206894120659401340373818L,
+ 629.95369438888946233003926191755125L,
+ 1568.35771983533158591604513304269098L,
+ 1646.02452040831961063640827116581021L,
+ 2299.96860633240298708910425594484895L,
+ 1222.73204392037452750381340219906374L,
+ 799.359797306084372350264298361110448L,
+ 72.7415265778588087243442792401576737L,
+ };
+ result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+ result *= exp(-z * z) / z;
}
- T g = exp(-z * z) / z;
- result = g * b + g * r;
}
else
{
Modified: sandbox/math_toolkit/libs/math/doc/sf_and_dist/minimax.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/sf_and_dist/minimax.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/sf_and_dist/minimax.qbk 2008-03-28 07:27:38 EDT (Fri, 28 Mar 2008)
@@ -94,12 +94,20 @@
try a higher value if an approximation won't converge,
or a lower value to get speedier convergence.]]
[[x-offset val][Sets the x-offset to /val/: the approximation will
- be generated for `f(x + X) + Y` where /X/ is the x-offset
+ be generated for `f(S * (x + X)) + Y` where /X/ is the
+ x-offset, /S/ is the x-scale
and /Y/ is the y-offset. Defaults to zero. To avoid
rounding errors, take care to specify a value that can
be exactly represented as a floating point number.]]
+[[x-scale val][Sets the x-scale to /val/: the approximation will
+ be generated for `f(S * (x + X)) + Y` where /S/ is the
+ x-scale, /X/ is the x-offset
+ and /Y/ is the y-offset. Defaults to one. To avoid
+ rounding errors, take care to specify a value that can
+ be exactly represented as a floating point number.]]
[[y-offset val][Sets the y-offset to /val/: the approximation will
- be generated for `f(x + X) + Y` where /X/ is the x-offset
+ be generated for `f(S * (x + X)) + Y` where /X/
+ is the x-offset, /S/ is the x-scale
and /Y/ is the y-offset. Defaults to zero. To avoid
rounding errors, take care to specify a value that can
be exactly represented as a floating point number.]]
Modified: sandbox/math_toolkit/libs/math/minimax/f.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/minimax/f.cpp (original)
+++ sandbox/math_toolkit/libs/math/minimax/f.cpp 2008-03-28 07:27:38 EDT (Fri, 28 Mar 2008)
@@ -21,13 +21,18 @@
switch(variant)
{
case 0:
- return boost::math::erfc(x) * x / exp(-x * x);
+ {
+ boost::math::ntl::RR x_ = sqrt(x == 0 ? 1e-80 : x);
+ return boost::math::erf(x_) / x_;
+ }
case 1:
- return boost::math::erf(x);
+ {
+ boost::math::ntl::RR x_ = 1 / x;
+ return boost::math::erfc(x_) * x_ / exp(-x_ * x_);
+ }
case 2:
{
- boost::math::ntl::RR x_ = x == 0 ? 1e-80 : x;
- return boost::math::erf(x_) / x_;
+ return boost::math::erfc(x) * x / exp(-x * x);
}
case 3:
{
@@ -294,6 +299,17 @@
{
return (boost::math::expint(x) - x) * x * exp(-x);
}
+ case 26:
+ {
+ //
+ // erf_inv in range [0, 0.5]
+ //
+ boost::math::ntl::RR y = x;
+ if(y == 0)
+ y = boost::math::tools::epsilon<boost::math::ntl::RR>() / 64;
+ y = sqrt(y);
+ return boost::math::erf_inv(y) / (y);
+ }
}
return 0;
}
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk