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
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.