9

I am attempting to form a double precision floating point number (64-bit) by taking the ratio of one product of integers divided by another product of integers. I wish to do so in a way that reduces the rounding error.

I am familiar with Kahan summation for addition and subtraction. What techniques work for division?

The numerator is a product of many long values (tens of thousands), likewise the denominator. I wish to prevent overflow and underflow, too. (One application is estimating infinite products by stopping after a sufficient number of terms.)

One thing I have tried is to factor the easily factorable numbers (using trial division by known primes up to a million) and cancel common factors, which helps, but not enough. My errors are approximately 1.0E-13.

I am working in C#, but any code that works with IEEE standard floating point numbers is welcome.

RESEARCH:

I came across a good paper that discusses EFT (Error Free Transformations) for + - x /, Horner's Rule (polynomials), and square root. The title is "4ccurate 4lgorithms in Floating Point 4rithmetic" by Philippe Langlois. See http://www.mathematik.hu-berlin.de/~gaggle/S09/AUTODIFF/projects/papers/langlois_4ccurate_4lgorithms_in_floating_point_4rithmetic.pdf

The above pointed me to Karp and Markstein (for division): https://cr.yp.to/bib/1997/karp.pdf

1
  • Indeed, Python has built in support for high precision arithmetic. To do the same, I would need to use an extended precision library in .Net, which I prefer to avoid. May 16, 2016 at 14:52

4 Answers 4

3

What techniques work for division?

For the division a/b, you can evaluate the residue (remainder):

a = b*q + r

This remainder r is easily accessible if you have fused-multiply-add

q = a/b ;
r = fma(b,q,-a) ;

The same fma trick can be applied on multiplication:

y = a*b ;
r = fma(a,b,-y) ; // the result is y+r

Then if you end up with two approximate operands after the products (a0+ra) / (b0+rb), you are interested in (a0+ra) = q*(b0+rb) + r.
You can first evaluate:

q0 = a0/b0 ;
r0 = fma(b0,q0,-a0);

Then approximate the remainder as:

r = fma(q0,rb,r0-ra);

Then correct the quotient as:

q = q0 + r/b0;

EDIT: What if fma is not available?

We can emulate the fma using an exact product à la Dekker, which is decomposed into the exact sum of 2 floating point, then a Boldo-Melquiond roundToOdd trick to be sure to have the sum of 3 floating point exactly rounded.

But it's going to be overkill. We use fma only for evaluating the residual error, so we generally have c very close to -ab. In this case, ab+c is exact, and we only have 2 floating points to sum, not 3.

Anyway, we only roughly estimate the residual error of a bunch of operations, so the last bit of this residue would not have been that important.

So the fma can be written like this:

/* extract the high 26 bits of significand */
double upperHalf( double x ) {
    double secator = 134217729.0; /* 1<<27+1 */
    double p = x * secator; /* simplified... normally we should check if overflow and scale down */
    return p + (x - p);
}

/* emulate a fused multiply add: roundToNearestFloat(a*b+c)
   Beware: use only when -c is an approximation of a*b
   otherwise there is NO guaranty of correct rounding */
double emulated_fma(a,b,c) {
    double aup = upperHalf(a);
    double alo = a-aup;
    double bup = upperHalf(b);
    double blo = b-bup;

    /* compute exact product of a and b
       which is the exact sum of ab and a residual error resab */
    double high = aup*bup; 
    double mid  = aup*blo + alo*bup;
    double low  = alo*blo;
    double ab = high + mid;
    double resab = (high - ab) + mid + low;

    double fma = ab + c; /* expected to be exact, so don't bother with residual error */
    return resab + fma;
}

Well, a bit less overkill than the general emulated fma, but it might be more clever to use a language that provides native fma for this part of the job...

3
  • 1
    The fma() call in r = fma(q0*rb,r0-ra); has only two instead of the required three operands, I would suggest double-checking that line.
    – njuffa
    May 18, 2016 at 10:49
  • From another question, I learned that FMA is not available in C#. May 26, 2016 at 19:22
  • @PaulChernoch in this case we can still emulate it, I completed the answer
    – aka.nice
    May 26, 2016 at 23:10
2

The multiplication equivalent of Kahan summation you are looking for is “double-double multiplication”. Here, if your integers are representable as double values, the function Mul122 from crlibm is mostly enough.

#define Mul122(resh,resl,a,bh,bl)                 \
{                                                 \
    double _t1, _t2, _t3, _t4;                    \
                                                  \
    Mul12(&_t1,&_t2,(a),(bh));                    \
    _t3 = (a) * (bl);                             \
    _t4 = _t2 + _t3;                              \
    Add12((*(resh)),(*(resl)),_t1,_t4);           \
}

bh and bl are the running product stored with additional precision as the sum of two double values. a is the next integer (we are assuming it is exactly converted to a double). resh and resl receive the next running product, in which the factor a has been taken into account.

In order to avoid underflow and overflow, you can externalize the exponent to an integer of the width you wish. This is done by periodically applying the frexp function to the high part of the running product, and then normalizing the running product by dividing both components by the same power of two (tracking the total power of two by which the running product has been divided can be done on the side with an integer variable of the desired width).

How often to apply frexp depends on the bound you have on the integers you are multiplying. If the integers are below 253, which would help with them being exactly representable as double values, you can do about 19 multiplications before having to normalize the running product, because the double-precision exponent goes up to 1023.

Once you have computed the products corresponding to the numerator and to the denominator, throw away the low components, and divide the high components. This will only introduce an error of about 1ULP. You were not aiming for an error of less than a double-precision ULP, were you?

Do not forget the powers of two that you left on the side for both numerator and denominator! Subtract them and apply the difference to the quotient with the ldexp function.

2
  • @PaulChernoch If what you call “ratio of one product by another” is in fact a product of ratios that tend towards one, like the “infinite product” example, the powers of two extracted from numerator and denominator can be added to and subtracted from the same variable, for which a simple int should then suffice.compute the numerator product and the denominator product at the same rhythm in order to do this. May 16, 2016 at 20:58
  • I am using infinite products as a unit test, and those terms tend towards one. But my intended usage involves products of terms that may not always tend towards one. May 16, 2016 at 21:31
2

Division doesn't suffer from the same catastrophic cancellation effects as addition and subtraction, and using IEEE floats is correctly rounded, and so should have have relative error of about 1/2 ulps (~2e-16). Any errors larger than that are most likely a result of intermediate products, so care needs to be taken with these.

Dekker (1971) has some algorithms for extending the precision of elementary mathematical operations: as pointed out by the other answer, these can be simplified if you have access to an fma operation.

1
  • I do not know how C# deals with FMA. Perhaps the jitter can perform the optimization if I arrange my operations in an order that it can recognize. May 18, 2016 at 12:51
0

The other answers are good if you have access to FMA (fused multiply-add), but C# does not make use of it. I continue to look for a fast solution, but I have found an accurate one.

Step 1: Collect the numerators and denominators separately.

Step 2: Strip off the sign and count how many multipliers were negative to learn the sign of the answer.

Step 3: Loop over all numbers, computing the natural log of each.

Step 4: Accumulate separate compensated sums for the logs of the numerators and denominators. (Use Kahan summation.)

Step 5: Take the difference between the two sums and compute the exponential.

Step 6: Restore the sign.

I tested this against 100,000 random integers in the numerator and the same numbers in the denominator, but with both sets shuffled in a different random order. If I use the naive approach of regular multiplication and division, my cumulative error is about 2x10^-15. Using my compensated log approach, the error is zero. (I got lucky maybe?) I will do more testing of harder cases. Nevertheless, by compensating the sum of the logs, I get almost twice the precision before the final rounding.

I am surprised that it worked so well. Obviously performing 200,000 logarithms is not ideal.

Theory note:

Cumulative rounding error is like a random walk. After N computations, you can expect an error of sqrt(N)*ULP/2. If ULP/2 is 5.0E-18 and N is 200,000, then you get 2.2E-15, which is close to what I got for the naive approach.

1
  • You don't need a FMA to compute Mul122, it is only more efficient if you have one. But the inefficiency of FMA-less Mul122 pales in comparison of that of log. May 27, 2016 at 6:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.