[Python] Double-precision floating-point loss for triple sum

For prototyping I am using Python and Sympy but the eventual implementation will target C++. (I don't think there will be any difference between C and C++ for the purposes of this question.) In all cases the floating-point math is IEEE 754 double precision.

From Sympy, the exact Rational arbitrary-precision values for these three quantities are known to be

-21868
381874348283167701/8796093022208
81525135/65536


I need to calculate their sum in a C context for which arbitrary precision is not available and double is required. The quantities shown in C floating-point hex notation after division:

-0x1.55b0000000000p+14
+0x1.532c2c939020fp+15
+0x1.36fe63c000000p+10


The C implementation will only have access to those quantities and not the true rationals.

When the rationals are summed, the exact total is

200463504468116437/8796093022208


and when divided produces

0x1.64183f632041fp+14


but the naive sum of the individual doubles is

0x1.64183f632041ep+14


which has one bit of error. I have tried to eliminate the error using the following strategies:

I have attempted the sum over every permutation of the terms:

for permutation in itertools.permutations(terms):
print(sum(permutation))


This does not change the output.

I have attempted the Kahan sum over every permutation of the terms:

def kahan(terms: typing.Iterable[float]) -> float:
# https://en.wikipedia.org/wiki/Kahan_summation_algorithm#The_algorithm
total = 0
comp = 0
for term in terms:
y = term - comp
t = total + y
comp = (t - total) - y
total = t
return total


for permutation in itertools.permutations(terms):
print(kahan(permutation))


This does not change the output. Is there any way to modify or replace the Kahan sum to reduce this error?

Continue reading...