"""Empirical check of the accuracy of numpy's binary64 sin and cos, which fcore.bounds uses in its floating-point
sums.  The rounding-error budget there assumes errors of at most four units in the last place (ulp); this script
compares numpy against Arb (200-bit ball arithmetic) on 10^5 arguments covering the range used (xi*y with
xi <= Xi1 <= 1697 and y <= 1) and prints the largest error found, in ulps of the computed value."""
import math
import numpy as np
from flint import arb, ctx

ctx.prec = 200
rng = np.random.default_rng(1)
xs = np.concatenate([rng.uniform(0, 2000, 60000), rng.uniform(0, 10, 20000), rng.uniform(1e-3, 1, 20000)])
s, c = np.sin(xs), np.cos(xs)
worst = {"sin": 0.0, "cos": 0.0}
for x, sv, cv in zip(xs, s, c):
    a = arb(float(x))
    for name, val, ref in (("sin", sv, a.sin()), ("cos", cv, a.cos())):
        err = abs(arb(float(val)) - ref).upper()
        u = math.ulp(abs(float(val))) if val != 0 else 5e-324
        worst[name] = max(worst[name], float(err) / u)
print("numpy %s: max error over %d arguments: sin %.3f ulp, cos %.3f ulp" % (np.__version__, len(xs), worst["sin"], worst["cos"]))