Comments on Special functions, Quaternions, Octonions

From: boost (boost_at_[hidden])
Date: 2001-05-25 07:39:33

Salut,

first of all I would see the library in boost, but I suggest that problems 
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().
By the way, why is no-one using the correct name:
arctan == arcus tangens; artanh = area tangens hyperbolicus
Inaddition, what is the deeper meening of sinc_pi, it behaves
like sinc. It confused me quite some time.

Quaternions:
===========

There seem to be problem with unneccessary overflows,
but I didn't have the time to look into detail:

template<typename X>											\
quaternion<type> &		operator /= (quaternion<X> const & rhs)	\
{																\
        type	ar = static_cast<type>(rhs.R_component_1());		\
        type	br = static_cast<type>(rhs.R_component_2());		\
        type	cr = static_cast<type>(rhs.R_component_3());		\
        type	dr = static_cast<type>(rhs.R_component_4());		\
                                                                                                                                \
        type	denominator = ar*ar+br*br+cr*cr+dr*dr;				\
                                                                                                                                \
        type	at = (+a*ar+b*br+c*cr+d*dr)/denominator;			\
        type	bt = (-a*br+b*ar-c*dr+d*cr)/denominator;			\
        type	ct = (-a*cr+b*dr+c*ar-d*br)/denominator;			\
        type	dt = (-a*dr-b*cr+c*br+d*ar)/denominator;			\
                                                                                                                                                \
        a = at;														\
        b = bt;														\
        c = ct;														\
        d = dt;														\
                                                                                                                                \
        return(*this);												\
}

Can this and similar expressions be reorganized to avoid overflows if 
at,bt,ct,dt ? If this is not done for speed reasons, it should be noted,
that there is a danger of overflow in the denominator, allthough the result
itself won't overflow.

Best wishes,
Peter