Re: Does the order matter to add a sequence of floating numbers

On Oct 30, 9:29 am, Richard Heathfield <r...@xxxxxxxxxxxxxxx> wrote:
JosephLee said:

Let's say float arr[100], we want to add them up, does the order
matter? what is the different between arr[0] +..+ arr[99], and
arr[99]+...+arr[0]?

You wouldn't think it would make a difference, would you? But it can,
nevertheless, especially when some or all of the values are near the
limits of what can be represented.

I tried fairly hard to write a program to demonstrate this, and failed. :-)

Take 2 and add DBL_EPSILON to it 1 million times.
Then, take DBL_EPSILON, add it up 1 million times, and add 2 to it.

In general, the most accurate way to add floating point numbers is
from smallest to largest.

You can get nearly the same accuracy by using Kahan's compensated
summation (Kahan's adder).
For example:

//
// Kahan Summation (also known as compensated summation) template
// Implementation by Dann Corbit 8-11-1998
// Thanks also to Tom Lippincott & Nathan Myers who gave valuable
advice.
// But the crappy stuff in here is my fault, not theirs! ;-)
//
// This technique is due to Dr. William Kahan,
// a Professor at University of California, Berkeley
//
// Personal note (DRC): Prof. Kahan is a very nice guy and a
gentleman.
//
// Ref:
// "What Every Computer Scientist Should Know About Floating Point
Arithmetic,"
// by David Goldberg, published in the March, 1991 issue of Computing
Surveys.
// Copyright 1991, ACM
//
// Theorem 8 (Kahan Summation Formula) on page 203 shows:
// error in iterative summation can be reduced to 2 * epsilon,
// whereas naive summation can produce n * epsilon.
//
template <class Etype> class KahanAdder
{

private:
Etype sum_; //The current working sum.
Etype carry_; //The carry from the previous operation

// You may be tempted to move these into the body of member
add() as
// local variables. DON'T DO IT! For some extended precision
types
// the constructor time may dominate your application, or even
cause
// it to fail due to memory fragmentation.
Etype temp_; //A temporary used to capture residual bits of
precision
Etype y_; //A temporary used to capture residual bits of
precision

// This is the heart of the template. After each add operation,
// the carry will contain the additional bits that could be left
out
// from the previous operation.
void add(const Etype & datum) {
y_ = datum - carry_;
temp_ = sum_ + y_;
carry_ = (temp_ - sum_) - y_;
sum_ = temp_;
}

public:
// Constructors. Both temp_ and y_ are initialized
// on first use. Note that these FORCE a constructor for zero.
// This was done to make _asolutely sure_ that the components
// are initialized to zero and not arbitrary values.
// I'm a bit torn, however, because I may inherit a type for
// which an integral constructor makes no sense.
// I would like to hear more debate on this.
KahanAdder() : sum_(0), carry_(0) {}
KahanAdder(const Etype& e): sum_(e), carry_(0) {}
KahanAdder(const Etype& e, const Etype& c): sum_(e), carry_(c)
{}

// Default destructor. Nothing special happens.
// ~KahanAdder() {} // add stuff here if you want to build one.

// Accessor method to return the working sum.
Etype get_sum() const { return sum_; }

// Accessor method to return the current carry.
Etype get_carry() const { return carry_; }

// I have powerful pro's and con's over this one.
// Opinions?
operator const Etype&() const { return sum_; }

// Syntactic sugar so that we can use the adder in a transparent
way.
// Here, we add an element to the sum with +=, which really just
calls
// the add() method.
KahanAdder & operator += (const Etype & g) { add(g); return
*this; }

// Syntactic sugar so that we can use the adder in a transparent
way.
// Here, we subtract an element to the sum with -=, which really
just
// calls the add() method with the arg negated.
KahanAdder & operator -= (const Etype & g) { add(-g); return
*this; }

// Add two KahanAdders together.
template< class Etype >
KahanAdder<Etype> operator+(const KahanAdder<Etype>& right)
{
*this += right.get_carry();
return *this += right.get_sum();
}

// Subtract one KahanAdder from another.
template< class Etype >
KahanAdder<Etype> operator-(const KahanAdder<Etype>& right)
{
*this -= right.get_carry();
return *this -= right.get_sum();
}
};

.