Re: [boost] Digest Number 594

From: Kevin Lynch (krlynch_at_[hidden])
Date: 2001-05-25 12:15:11

boost_at_[hidden] wrote:
> 
> 
> Message: 13
>    Date: Fri, 25 May 2001 14:39:33 +0200
>    From: Peter Schmitteckert (boost) <boost_at_[hidden]>
> Subject: Comments on Special functions, Quaternions, Octonions
> 
> 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.
> 

I concur with this...I think this is a valuable library for boost, but I
would like to see the possibility of overflows addressed.  

[snip]

> 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) \
> {                                                                                                              

[code deleted]

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

The problem, in more detail, is that intermediate numerical results
(such as denominator in the above function) may overflow, ruining the
entire calculation, despite the fact that the result is perfectly
representable, and that there are alternative strategies for obtaining
the intermediate results such that no overflow will occur.  This may
crop up anywhere where you do things like taking the square root of a
square, or dividing a quantity that contains a variable by another
quantity that contains the square of that same variable.  Places to
watch for this are especially

operator/
operator/=
abs(x)
norm(x)
sqrt(x)
etc....

This is a problem in complex<T> as well... if you look at, say, the
source to the g++ library implementation, or section 5.4 of Numerical
Recipes in C, you'll see how to combat this problem.  Here is an example
replacement for one of your functions.

quaternion<T> &		operator /= (::std::complex<T> const & rhs)
{
        T	ar = rhs.real();
        T	br = rhs.imag();

        T at,bt,ct,dt,abdiv,denom;

        if( abs(ar)>abs(br) ) 
        {
                abdiv = br/ar;
                denom = ar + br*abdiv;

                at = (+a+b*abdiv)/denom;
                bt = (-a*abdiv+b)/denom;
                ct = (+c-d*abdiv)/denom;
                dt = (+c*abdiv+d)/denom;
        } else {
                abdiv = ar/br;
                denom = ar*abdiv + br;

                at = (+a*abdiv+b)/denom;
                bt = (-a+b*abdiv)/denom;
                ct = (+c*abdiv-d)/denom;
                dt = (+c+d*abdiv)/denom;
        }

                        
        a = at;
        b = bt;
        c = ct;
        d = dt;
                        
        return(*this);
}

Of course, I haven't checked this implementation, and it may not even
compile.  Look upon it as only an example!

Now, the problem:  you need to find some way to do this consistently for 
denom = a^2 + b^2 + c^2 + d^2

but I haven't spent the time yet to analyze it in detail.  And I haven't
figured out abs, sqrt, or any of the others that might be affected.  But
I'm sure that your current implementation of abs has the problem.  And I
didn't find a sqrt implementation (does sqrt even make sense for
quaternions?  my ignorance is shining through....)

> 
> Best wishes,
> Peter
> 

-- 
-------------------------------------------------------------------------------
Kevin Lynch				voice:	 [phone]
Physics Department			Fax: [phone]
Boston University			office:	 PRB-565
590 Commonwealth Ave.			e-mail:	 krlynch_at_[hidden]
Boston, MA 02215 USA			http://physics.bu.edu/~krlynch
-------------------------------------------------------------------------------