Re: Comments on Special functions, Quaternions, Octonions

From: lums_at_[hidden]
Date: 2001-05-25 12:53:52

Just a quick reply to the overflow question.  One way to handle this 
is to find the largest of ar, br, cr, dr, scale them all by that 
value, compute the scaled denominator, compute the scaled at, bt, ct, 
dt, and then re-scale the results.

If I count correctly, there are three comparisons, four divides, and 
one more multiply introduced.

Actually, to do this a little better, one would want to compute 
inverse quantities and do multiplies rather than divides throughout.  
Then there would be three comparisons, one divide, and six more 
multiplies introduced to deal with the overflow issue.  The original 
method should be modified to include one divide and four multiplies 
rather than four divides.  This is all assuming the basic type of the 
quaternion is something the compiler/cpu would deal with efficiently.

I would produce an example, but my interface to boost is via yahoo 
groups and it is quite a pain to put code in.  I hope the explanation 
conveys the right idea.
                                                                

--- 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 
> 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