Hello,
I noticed that Todd left a hatch open for increasing the precision of
numeric summations through the use of NumericTypeTraits<>. I'm still
investigating how to take advantage of this within the context of the Blitz
source code.
My first step was to sketch out an adaptor class which provides
high-precision floating point summation support. Within this class,
intermediate sums are represented by a pair of floating point values. An
ordinary addition (or subtraction) is replaced by _four_ addition steps.
Todd's design should allow for compiling a Blitz program with high precision
support enabled for all the summation loops. This leaves open a very quick
way to test if one suspects numerical instabilities.
Here is the code which shows how the high precision addition is implemented.
#include <iostream>
#include <iomanip>
using namespace std;
template <class T=double>
struct StableSum {
typedef StableSum<T> Self;
StableSum (T v = 0) : x(v), e(0) {}
operator T() const { return x+e; }
Self& operator+= (const Self& v) {
e += v.x;
T sum = x + e;
T gain = sum - x;
e -= gain;
x = sum;
return *this;
}
Self& operator-= (const Self& v) {
x += -v.x;
return *this;
}
Self operator+ (const Self& v) const {
return Self(x) += v;
}
Self operator- (const Self& v) const {
return Self(x) -= v;
}
T x;
T e; // excess term
};
void stable_sum_test (void)
{
cout << setprecision (20);
double base = 1;
for (int k = 0; k < 15; ++k, base *= 10) {
double tot = 1;
double sum = base;
StableSum<double> stable = base;
for (int i = 1; i != 10; ++i, tot /= 10)
for (int j = 0; j != 1000*i; ++j) {
sum += tot; stable += tot;
}
sum -= base; stable -= base;
cout << "double= " << setw(20) << sum
<< " stable= " << setw(10) << stable << endl;
}
}
int main ()
{
stable_sum_test ();
return 0;
}
Program Output:
double= 1234.5678899999064 stable= 1234.56789
double= 1234.5678899998961 stable= 1234.56789
double= 1234.5678899998802 stable= 1234.56789
double= 1234.5678900005623 stable= 1234.56789
double= 1234.5678900128405 stable= 1234.56789
double= 1234.5678899873747 stable= 1234.56789
double= 1234.5678898273036 stable= 1234.56789
double= 1234.5678918063641 stable= 1234.56789
double= 1234.5679700374603 stable= 1234.56789
double= 1234.5679998397827 stable= 1234.56789
double= 1234.5657348632812 stable= 1234.56789
double= 1234.649658203125 stable= 1234.56789
double= 1234.4970703125 stable= 1234.56789
double= 1236.328125 stable= 1234.56789
double= 1234.375 stable= 1234.56789
You can see that the stable sums have maintained full precision long after
the ordinary double is shaving off digits.
Allan
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to majordomo@oonumerics.org, with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 04:30:08 EST