Testing special functions

ROOT
1

In the presentation regarding MathCore special function accuracy, linked from the ROOT 5.19 release notes, I noticed that M. Slawinska is testing the accuracy using absolute and relative differences.

However, this is not the proper way to test mathematical functions, and one should test and plot the accuracy across the input range using ulp (unit in the last place) instead. See e.g. J.-M. Muller, Elementary Functions, 2nd Ed. (Birkhäuser, Boston, 2006), p. 14. The unit ulp is preferred among virtually the entire field of computer arithmetics, e.g. CPU datasheets and standard library implementations. Only by this form of normalization one can read off the number of error bits in the machine representation, and how far the implementation is away from the minimum error (with correct rounding) of 0.5 ulp.

2

Another remark regarding testing against Mathematica is to make sure that you actually evaluate in Mathematica with the sufficient internal accuracy. For performance reasons, Mathematica defaults to $MachinePrecision in all its algorithms, which means it does not guarantee even the last bit accuracy.

You can see the behavior e.g.:

[code]In[1]:= Erfc[7.741] // InputForm

Out[1]//InputForm= 6.8361796318803375*^-28

In[2]:= N[Erfc[7741/1000], 20] // InputForm

Out[2]//InputForm=
6.8361796318803098220690529129963840378727151852417`20.*^-28
[/code]

Note the difference is 2.842 * 10^(-42), while 1 ulp at this value is 4.968 * 10^(-44) (= 2^(-143), since log_2(6.836 * 10^(-28)) = -90.241).

3

In case you wonder in the remark to Mathematica whether the difference is due to the machine representation of the initial argument 7.741: It is not. And you can see that by comparing against the double precision, binary fraction representation of 7.741:

[code]In[3]:= Erfc[7.741`] - N[Erfc[544724448679297/70368744177664], 20]

               -43

Out[3]= -8.96831 10[/code]

I also attached two plots (normalized at 1 ulp) that shows the global behavior for Erfc[] and Gamma[], both running on Mathematica 6.0.1.0 for Linux x86-64.



4

Hi,
thank you very much for your useful remarks and observations.
I agree with you that is better using ulp for testing instead of relative errors.
Thanks also for the remarks against Mathematica. We used the normal machine precision, but it would have been preferable as you suggested to use a larger internal accuracy. Anyway the older implementation we tested of a function like TMath::Erfc had errors much larger than the double machine precision.
We will eventually produce new test plots of the special functions using ULP and comparing with internal values of Mathematica used with better internal accuracy.

Best Regards

Lorenzo