"""Why parallel sums disagree: tracing FP32 reduction order
Deterministic, standard-library experiment. RESULTS contains the full data.
"""
import math
import struct

def f32(value):
    return struct.unpack('f', struct.pack('f', value))[0]

def serial(values):
    total = f32(0)
    for value in values:
        total = f32(total + value)
    return total

def tree(values):
    work = [f32(value) for value in values]
    while len(work) > 1:
        if len(work) % 2:
            work.append(0.0)
        work = [f32(work[i] + work[i + 1])
                for i in range(0, len(work), 2)]
    return work[0] if work else 0.0

for values in ([1e8, 1, -1e8, 1], [1e8, -1e8, 1, 1]):
    reference = math.fsum(values)
    print(f'serial={serial(values):.1f} tree={tree(values):.1f} reference={reference:.1f}')
assert serial([1e8, 1, -1e8, 1]) == 1
assert tree([1e8, 1, -1e8, 1]) == 0
assert tree([]) == 0
print('Rounding-order demonstration, not a GPU reduction benchmark.')


from itertools import permutations
rows=[]
for values in sorted(set(permutations([1e8,1.,-1e8,1.]))):
    reference=math.fsum(values); a=serial(values); b=tree(values)
    rows.append(dict(inputs=list(values),serial=a,pairwise=b,reference=reference,serial_abs_error=abs(a-reference),pairwise_abs_error=abs(b-reference),sum_abs=sum(abs(v) for v in values)))
assert len(rows)==12
RESULTS=dict(permutations=rows,serial_results=sorted(set(r['serial'] for r in rows)),pairwise_results=sorted(set(r['pairwise'] for r in rows)),empty_tree=tree([]))
print('Distinct input permutations:',len(rows))
print('Serial outcomes:',RESULTS['serial_results'],'pairwise outcomes:',RESULTS['pairwise_results'])