[boost] [rfc] Inexact quantities arithmetics

From: Marat Khalili (0x8.0p15_at_[hidden])
Date: 2007-07-08 14:07:29

Hello to everyone,

I'm developing a library for working with inexact quantities, like
(0.3162 +- 0.0079), in other words, for auto-estimating the propagation
of errors in calculations. I plan to release it under Boost license when
ready. Any comments (including sensibility of the whole project) are
welcome.

I. Some basic principles and limitations I established so far

1. Inexact quantity is a pair of 'value' and 'error' components.

2. Error component is always non-negative, but can be infinite or NaN.

3. Value and error are of the same type (at least for the user).
For example, value can not be 'double' while error is 'float', nor value
can be 'complex<double>' while error is 'double'.

4. Value result of calculation does not depend on error and coincides
with result of the same calculation for bare value. Correspondingly it
is biased from math expectation for any non-linear function.
For example, cos(0 +- 0.1).value() gives 1, not exp(-0.05).

5. Error calculation is limited to the first order. This means f(mx +-
dx) gives (f(mx) +- abs(dx * f'(mx))), f'' and subsequent terms
proportional to greater powers of dx are ignored.
For example:
     cos(0 +- 0.1) gives exactly 1;
     sqrt(0 +- 0.1) gives (0 +- inf).

Later, template parameter for altering behavior described in the
previous two items can be added.

6. Correlations are not considered; instead, worst case is always assumed.
For example:
     (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx +
my) +- sqrt(dx*dx + dy*dy);
     x *= x (where x is of inexact type) produces correct result, while
     x -= x will not produce 0 but, incorrectly, (0 +- 2*x.error()).

In the distant future additional operations with correlations (and
distributions) explicitly specified by user can probably be added.

NOTE: Items 4-6 make specific distribution (normal, lognormal or
whatever) irrelevant.

II. Proposed IO format

1. Currently IO format is specified by template parameter. I really need
your comments here. I want IO format to be both tweakable and
replaceable by user with minimum headache. Ideally I'd like to see my
specific manipulators similar to std::scientific and
std::precision(int), but I'm not sure it's possible.

2. Already, the following format for inexact<T> is implemented:

a) When T is simple floating-point type, the format is:
'(' value dummy error ')' [('e' | 'E') common_exponent]

Here dummy is any whitespace-terminated string. For wchar_t streams
sensible value is L" \xB1 ". For char streams default value is " +- ",
but one might decide to use " \pm ", or " &plusmn; ", or UTF-8
representation of L" \xB1 ", or probably something else as long as it
ends on space. The point is to read parenthesized part with:
     in >> value >> ws >> dummy >> ws >> error;

In output, presence of common exponent and specific format of value and
error (fixed or scientific, precision) depend on stream flags and values
of value and error in a complex way I will not elaborate on here.
Results look like this:
Default: (-0.0010676 +- 0.0000010)  (-1.16767e+06 +- 0.0011)
Fixed:    (-0.001068 +- 0.000001)  (-1167672.1657 +- 0.0011)
Scientific: (-1.0676 +- 0.0010)e-3 (-1.167672e+06 +- 1.1e-03)

b) Alternatively, just the bare value (without parentheses) can be used
in place of (value +- 0).

c) When T is collection C<E> of elements E, the format is by definition
the same as for C<inexact<E> >. Corresponding converters are provided.

For example, output of both
inexact<complex<double> >(
        complex<double>(1, 2),
        complex<double>(3, 4))
and
complex<inexact<double> >(
        inexact<double>(1, 3),
        inexact<double>(2, 4))
is "((1 +- 2),(3 +- 4))".

3. Format where error is stored in and retrieved from the number of
value digits can be useful as well.

III. Some notes about implementation

1. In implementation the type of error is wrapped in order to guarantee
non-negativeness. I'm working on this right now.

a) When value type T is simple floating-point type, error will be of
type absolute<T>, that trivially converts to T but calls abs() when
converting from T.

b) When T is collection C<E> of elements E, error will be of type C<A>,
where A is error type of E.

For example, in inexact<vector<complex<double> > > the value type will
be vector<complex<double> >, and error type will be
vector<complex<absolute<double> > >.

Hopefully all this stuff will get optimized out. I tried to do without
this wrapping first, but code is much more transparent with.

2. I'll need derivatives for all cmath functions at least. It is
probable that some template-based derivation-finder will be created as a
by-product. Or, does one exist already?

With Best Regards,
Marat

P.S. Pardon for possible double-post. I'm really lost how boost list + 
gmane combination works.