Uh oh!
There was an error while loading.Please reload this page.
- Notifications
You must be signed in to change notification settings - Fork33.3k
Description
Proposal:
Replace the1 + erf(s) computation witherfc(-s) as suggested in thisStackOverflow discussion and hinted in thisJohn Cook blog post. The core idea is to exploit the identity,1 + erf(x) == erfc(-x), to eliminate the addition step thus avoiding loss of precision.
Empirical analysis
Here is some code for empircal analysis. For 10,000 random values in each range bin, it shows the maximum difference in between the two approaches as measured in ULPs. It also shows how often each technique beats or agrees with the other as compared to a reference implementation usingmpmath:
importmathfrommpmathimportmpfromcollectionsimportCounterfromrandomimportuniformfromitertoolsimportpairwisemp.dps=50defbest(x):'Which is better? erf, erfc, or are they the same?'ref=1+mp.erf(x)e=1+math.erf(x)c=math.erfc(-x)de=abs(e-ref)dc=abs(c-ref)return'erf'ifde<dcelse'erfc'ifdc<deelse'='deferr(x):'Diffence in ulp'e=1+math.erf(x)c=math.erfc(-x)returnabs(round((e-c)/math.ulp(c)))deffrange(lo,hi,steps=10):step= (hi-lo)/stepsreturn [lo+i*stepforiinrange(steps+1)]forlo,hiinpairwise(frange(-5,5,steps=20)):xarr= [uniform(lo,hi)foriinrange(10_000)]winners=Counter(map(best,xarr)).most_common()max_err=max(map(err,xarr))print(f'{lo:-5.2f} to{hi:-5.2f}:{max_err:12d} ulp',winners)
On macOS with clang-1600.0.26.6, this outputs:
-5.00 to -4.50: 411273532899 ulp [('erfc', 10000)]-4.50 to -4.00: 3211997254 ulp [('erfc', 10000)]-4.00 to -3.50: 25164209 ulp [('erfc', 10000)]-3.50 to -3.00: 784219 ulp [('erfc', 9999), ('=', 1)]-3.00 to -2.50: 24563 ulp [('erfc', 9996), ('erf', 3), ('=', 1)]-2.50 to -2.00: 1535 ulp [('erfc', 9921), ('erf', 45), ('=', 34)]-2.00 to -1.50: 95 ulp [('erfc', 9490), ('erf', 298), ('=', 212)]-1.50 to -1.00: 11 ulp [('erfc', 7540), ('=', 1246), ('erf', 1214)]-1.00 to -0.50: 2 ulp [('erfc', 4211), ('=', 4039), ('erf', 1750)]-0.50 to 0.00: 1 ulp [('=', 9764), ('erfc', 155), ('erf', 81)] 0.00 to 0.50: 1 ulp [('=', 9228), ('erf', 698), ('erfc', 74)] 0.50 to 1.00: 1 ulp [('=', 9553), ('erfc', 364), ('erf', 83)] 1.00 to 1.50: 1 ulp [('=', 8708), ('erfc', 1156), ('erf', 136)] 1.50 to 2.00: 1 ulp [('=', 8792), ('erfc', 1167), ('erf', 41)] 2.00 to 2.50: 1 ulp [('=', 8722), ('erfc', 1271), ('erf', 7)] 2.50 to 3.00: 1 ulp [('=', 8780), ('erfc', 1220)] 3.00 to 3.50: 1 ulp [('=', 8728), ('erfc', 1272)] 3.50 to 4.00: 1 ulp [('=', 8730), ('erfc', 1270)] 4.00 to 4.50: 1 ulp [('=', 8742), ('erfc', 1258)] 4.50 to 5.00: 1 ulp [('=', 8756), ('erfc', 1244)]The results show massive improvement for negative inputs. For positive inputs, the difference is no more than 1 ulp anderfc wins in every bucket except for0.00 to 0.50.
@tim-one ran this on a Windows build (which uses a different math library) and found that "on Windows too there was no bin in which erf won more often than erfc".