import math
import mpmath
N = 100 # will test pi/2 +- N ulps
def ulp(x):
mant, exp = math.frexp(x)
return math.ldexp(0.5, exp - 52)
x = math.pi / 2
for _ in range(N):
x -= ulp(x)
mpmath.mp.prec = 200 # number of mantissa bits
print("trying tan(math.pi/2 +-%d" % N, "ulps)")
nfail = 0
maxdiff = 0.0
for i in range(-N, N+1):
want = mpmath.tan(mpmath.mpf(x))
if 0:
want = float(want) # OOPS!
else:
# float(want) doesn't round correctly - it chops:
# https://github.com/fredrik-johansson/mpmath/issues/324
# worm around that by making mpmath round to its own
# internal format with 53 bits first
want = float(mpmath.mpf(want, prec=53))
got = math.tan(x)
if want != got:
nfail += 1
print("at", i, "ulp", x, x.hex())
print(" want", want, want.hex())
print(" got ", got, got.hex())
diff = (got - want) / ulp(want)
print(" ", diff, "ulps from want to got")
maxdiff = max(maxdiff, abs(diff))
x += ulp(x)
print("failures", nfail, "max abs(diff)", maxdiff)