c++ - How to calculate floating-point precision after round-off errors in +, -, *, and /? -

for purpose of verification, able calculate reasonably tight upper bound on accumulated error due rounding representable values during specific arithmetic computation.

assume have function foo() claims perform specific arithmetic computation. assume there implied guarantee maximum error (due rounding) stemming involved type being either float or double , implied (or stated) way in foo() caries out computation.

i able verify result foo() specific set of input values carrying out computation in fashion keeps tracks accumulated worst-case error, , check if 2 results close final worst-case error demands.

i imagine possible introducing new arithmetic class track_prec<t> adds precision tracking 1 of fundamental floating-point types, , let implementation of arithmetic operators of class calculate worst-case errors of each subexpression. problem don't know how compute these worst-case errors in these general cases:

// t = float or double template<class t> class track_prec { public:     t value;     t ulp; // http://en.wikipedia.org/wiki/unit_in_the_last_place      track_prec& operator+=(const track_prec& v)     {         value += v.value;         ulp = ???; // how this? , -=, *=, , /=?     }      friend bool operator==(t, const track_prec&)     {         // how should comparison done?     } }; 

assume, example, foo() simple summation on sequence of numbers. use track_prec<t> follows:

std::vector<t> values { 0.4, -1.78, 1.3e4, -9.29e3, ... }; check_equal(std::accumulate(values.begin(), values.end(), track_prec<t>()),             foo(values.begin(), values.end())); 

of course, kind of welcome, pointers free , working code nice.

i found these links on subject, not seem provide direct answers questions.

• http://www.cl.cam.ac.uk/~jrh13/papers/fparith.pdf
• http://docs.oracle.com/cd/e19957-01/806-3568/ncg_goldberg.html#1125
• http://www.holoborodko.com/pavel/mpfr/#features

the simplest way track precision of floating-point computations called interval arithmetic. not require ieee 754 arithmetic, computations can switched round upwards or downwards, bounds of computed intervals @ each step contain reals have resulted same computations had been done real arithmetic.

you should able find many existing implementations of interval arithmetic.

the accuracy of computation @ given step width of interval computed @ step.

beware of constants: if wish approximate πx, need multiply floating-point interval contains π interval x. if multiply interval x double denoted 3.1415926535897932, error multiplication of x double (which equals neither π nor 3.1415926535897932). in code in question, constant 0.4 means “the double nearest 0.4”. if need rational number 4/10, use interval bounds 2 doubles denoted respectively 3.9999999999999997e-01 , 4.0000000000000002e-01.