Free42 possible accuracy flaw

[Albert Chan]

If you don't mind cost of fma, we could square with more precisions.

Code:

def normalize(x, y):
    z = x + y
    return complex(z, x-z+y)

def sqr(z):
    x, y = z.real, z.imag
    x2 = x*x
    return normalize(x2, fma(x,x,-x2) + 2*x*y)

Above use complex part to carry the error.

(x+ε)^2 ≈ x^2 + 2xε       , where x^2 = float(x^2) + fma(x, x, -float(x^2))

Lets do (1+123/2^26) ^ (2^26)

>>> from gmpy2 import *
>>> z = 1 + 123/2**26
>>> t1 = t2 = z
>>> for i in range(26): t1=t1*t1; t2=sqr(t2)
...
>>> t1
2.6192220598552021e+53
>>> t2
(2.6192220641937287e+53+1.3516279600598801e+36j)

Real part result is correct, with digits to spare.