Re: Comments on Special functions, Quaternions, Octonions

From: Hubert HOLIN (Hubert.Holin_at_[hidden])
Date: 2001-05-25 17:49:08

Paris (U.E.), le 26/05/2001

--- In boost_at_y..., Peter Schmitteckert (boost) <boost_at_s...> wrote:
> Salut,
> 
> first of all I would see the library in boost, but I suggest that problems 

        Thanks!

> with accuracy/overflows  should be revised, or at least be mentioned.
> 
> Special Functions:
> ==============
> 
> There are calculations of the following form
> 
> if (abs(x) <= numeric_limits<T>::epsilon())
> {
>   return(static_cast<T>(1));
> }
> else
> {
>  return(sin(x)/x);
> }
> 
> I would like to suggest to use something similar to
> 
> const T  eps = sqrt(sqrt(numeric_limits<T>::epsilon()));
> if	(abs(x) <= eps) // or  "<" ? which is better for performance?
> {
> 	return(static_cast<T>(1) - x*x/6);
> } else
> {
> 	return(sin(x)/x);
> }
> 
> e.g. using t.C with g++ 2.9.3  one can see a small inaccuracy with long 
> double:
> #include <iostream.h>
> #include <math.h>
> 
> using namespace std;
> 
> int main()
> {
> 
>     cout << sizeof(double) << '\t' << sizeof(long double) << endl << endl;
> 
>     cout.precision(20);
> 
>     cout << "double: " << endl;
>     for( double x = 1e-7; x >1e-8; x *=0.8 )
>     {
>         cout << x << '\t' << sin(x) << '\t' << sin(x) / x << '\t' << 
> double(1.0) - x*x/6.0
>             << '\t' << double(1.0) - x*x/6.0 + x*x*x*x/120.0 - x*x*x*x*x*x / 
> 5040. <<endl;
>     }
> 
>     cout << endl << "long double: " << endl;
> 
>     for( long double x = 1e-7; x >1e-8; x *=0.8 )
>     {
>         cout << x << '\t' << sin(x) << '\t' << sin(x) / x << '\t' << 1.0 - 
> x*x/6.0
>             << '\t' << 1.0 - x*x/6.0 + x*x*x*x/120.0 - x*x*x*x*x*x / 5040. 
> <<endl;
>    }
> 
>   return 0;
> }
> 
> Same applies to  sinhc, similar to atanh().

        Two remarks here:

1)	could you compare the discrepancy to numeric_limits<T>::epsilon() 
for small values of the argument.

2)	this is, as noted in the docs, a "stop-gap" library, we will have 
to fine-tune it for various hardware/software combinations... so this 
situation will probably require some extensive list of ifdefs.

> By the way, why is no-one using the correct name:
> arctan == arcus tangens; artanh = area tangens hyperbolicus

        Histerical raisins ;-)

> Inaddition, what is the deeper meening of sinc_pi, it behaves
> like sinc. It confused me quite some time.

        "sinc" is sometimes the name of a one-parameter family of 
functions, of which "sinc_pi" is the most famous (this is the approach 
we have taken here,  see docs). At other times (f.i. MATLAB), "sinc" 
just designates what I called here "sinc_pi". Just one more unfortunate 
ambiguity of mathematics...

> Quaternions:
> ===========

[SNIP]

        (for this part see another post, will be done).

        Looks like I forgot to turn my escapes into ordinary white spaces 
for this update, sorry...

                Hubert Holin
                Hubert.Holin_at_[hidden]