Gaussian sampler improvements

[osdnk]

1. Recalibrate Ramez coeffs

HI = RealField(200)
DBL = RealField(53)     
ln2 = log(HI(2))
L = ln2   # target: x in [-ln2, 0]

# The reconstruction 1 - (x*c/(c-2) - x) equals exp(x) exactly iff
# c = x*coth(x/2) - 2.  The code stores c = x - x^2*P(x^2), so P must match:
def g(x):
    x = HI(x)
    if x == 0:
        return HI(1)/6
    return (x/tanh(x/2) - 2)/x^2

def fit(d):                  # degree-d poly in z=x^2, interpolated at Chebyshev nodes
    nodes = [-L/2 + L/2*cos(HI(pi)*(2*k+1)/(2*(d+1))) for k in range(d+1)]
    A = matrix(HI, d+1, d+1, lambda i, j: nodes[i]^(2*j))
    b = vector(HI, [g(x) for x in nodes])
    return A.solve_right(b)

def exp_small(x, C, F): 
    x = F(x); z = x*x
    P = F(0)
    for c in reversed(list(C)):
        P = P*z + F(c)
    cc = x - z*P
    return 1 - ((x*cc)/(cc - 2) - x)

def max_rel_err(C, F, N=40000):
    m = HI(0)
    for k in range(N+1):
        x = -L*k/N
        m = max(m, abs(HI(exp_small(x, C, F)) - exp(HI(x))) / exp(HI(x)))
    return float(log(m)/ln2)

for d in [4, 5, 6]:
    C = fit(d)
    print("  degree %d (%d coeffs):  approx err 2^%.1f   double-precision err 2^%.1f"
          % (d, d+1, max_rel_err(C, HI), max_rel_err([DBL(c) for c in C], DBL)))

d = 5
for j, c in enumerate(fit(d)):
    print("#define C%d (%s)" % (j+1, DBL(c).hex()))

Re-calibrated to match the precision of double

2. Constant time exp_small

I start from a guess y = -0.4193, which is in the middle of x \in [−ln2, 0]. The worst-case initial error is then 0.16. After 5 steps, we square the error to 2^-84, below the double precision.

3. Constant-time sweep over the CDF table using __builtin_sub_overflow -- standard trick for ct comparison.