"""Rigorous lower bound for the variational quotient (Section 5.2 of the paper) of F_psi(t) = psi(t_1+...+t_k) prod_i g(t_i), g(t) = 1_{[0,T]}(t)/(c+bt), c = theta/log k, T = beta/log k, b = k-1 (theta is gamma in the paper; h is psi). psi is piecewise constant: psi = h_j on [E_j D, E_{j+1} D), j = 0..nc-1, with lattice spacing D (= 2^-18 or 2^-19, passed in from the trial file; the module default is not used), E_0 = 0, E_nc D = 1 (eps_cells = 0 in the paper), h_0 = 0, h_j in [0,1] a multiple of 2^-40. With X ~ g^2/m2 on [0,T], m2 = int g^2, S_n = X_1+...+X_n (Lemma 4.1): Q = (k/m2) E[H(S_{k-1})^2] / E[h(S_k)^2], H(s) = int_0^T h(s+t) g(t) dt = h_{nc-1} Gt(E_nc D - s) + sum_{j=1}^{nc-1} (h_{j-1}-h_j) Gt(E_j D - s), Gt(y) = G(min(T, max(0,y))), G(y) = log(1+by/c)/b. Numerator: Gt is nondecreasing, so on each lattice cell [i D, (i+1) D] one has H >= l_i (each term evaluated at the appropriate end, using rigorous enclosures of Gt on the lattice), and E[H^2] >= sum_i l_i^2 P(S_{k-1} in cell_i) (Abel summation with the two-sided CDF bounds of fcore.py). Denominator: Abel summation over the h-cells.""" import sys, json, math, time import numpy as np from flint import arb, ctx import os; sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from fcore import CDF DELTA = 2.0 ** -17 U = 2.0 ** -53 def gtable(k, theta, beta, nmax, DELTA=DELTA): """Lower/upper binary64 bounds of Gt(n*DELTA), n = 0..nmax (Gt = G(T) beyond T).""" ctx.prec = 128 kk = arb(k); L = kk.log(); c = arb(theta) / L; T = arb(beta) / L; b = kk - 1 GT = (1 + b * T / c).log() / b lo = np.empty(nmax + 1); hi = np.empty(nmax + 1) for n in range(nmax + 1): y = arb(n) * arb(DELTA) # exact lattice point if y < T: # Arb comparisons are rigorous ("certainly less") v = (1 + b * y / c).log() / b elif y > T: v = GT else: v = GT.union((1 + b * y / c).log() / b) # cannot decide: enclose both lo[n] = float(v.lower()); hi[n] = float(v.upper()) # float() of arb.lower()/upper() rounds to nearest: widen by one ulp to be safe lo = np.nextafter(lo, -np.inf); hi = np.nextafter(hi, np.inf) lo = np.maximum(lo, 0.0) m2 = (1 / b) * (1 / c - 1 / (c + b * T)) return lo, hi, m2, GT def hcert(k, theta, beta, E, h, eps_cells=0, verbose=True, cdf=None, ilo=None, DELTA=DELTA): t0 = time.time() E = np.asarray(E, dtype=np.int64); h = np.asarray(h, dtype=float) assert E[0] == 0 and np.all(np.diff(E) > 0) and len(h) == len(E) - 1 assert np.all(h >= 0) and np.all(h <= 1) assert np.all(h * 2.0 ** 40 == np.round(h * 2.0 ** 40)) # differences below are exact ONE = int(round(1 / DELTA)) assert E[-1] == ONE + eps_cells nc = len(h) if cdf is None: cdf = CDF(k, theta, beta, ns=(k - 1, k), verbose=verbose) Tf = cdf.Tf nmax = int(math.ceil(Tf / DELTA)) + 2 glo, ghi, m2, GT = gtable(k, theta, beta, nmax, DELTA) # ---------- denominator: E[h(S_k)^2] = h_{nc-1}^2 C(E_nc) + sum_{j=1}^{nc-1} (h_{j-1}^2 - h_j^2) C(E_j) ye = E[1:] * DELTA Lk, Uk = cdf.bounds(k, ye) den = arb(h[-1]) ** 2 * arb(Uk[-1]) denl = arb(h[-1]) ** 2 * arb(Lk[-1]) for j in range(1, nc): coef = arb(h[j - 1]) ** 2 - arb(h[j]) ** 2 if float(coef.mid()) >= 0: den += coef * arb(Uk[j - 1]); denl += coef * arb(Lk[j - 1]) else: den += coef * arb(Lk[j - 1]); denl += coef * arb(Uk[j - 1]) den_up = float(np.nextafter(float(den.upper()), np.inf)); den_lo = float(np.nextafter(float(denl.lower()), -np.inf)) t_den = time.time() - t0 # ---------- H lower bound on lattice cells [iD, (i+1)D], i = ilo .. itop-1 itop = ONE - eps_cells if ilo is None: ilo = int(E[1]) if nc > 1 else int(0.3 / DELTA) # mass of S_{k-1} below this is dropped (lower bound) I = np.arange(ilo, itop) def tab(arr, idx): out = np.empty(idx.shape) neg = idx < 0; big = idx > nmax mid = ~(neg | big) out[neg] = 0.0 out[big] = arr[nmax] out[mid] = arr[idx[mid]] return out # term for the top edge: h_{nc-1} Gt(E_nc D - s), lower bound at s = (i+1)D acc = h[-1] * tab(glo, E[-1] - (I + 1)) absacc = np.abs(acc) for j in range(1, nc): d = h[j - 1] - h[j] # exact if d == 0: continue if d > 0: v = d * tab(glo, E[j] - (I + 1)) else: v = d * tab(ghi, E[j] - I) acc = acc + v absacc = absacc + np.abs(v) l = acc - 2.0 * (nc + 2) * U * absacc * (1 + 1e-10) l = np.maximum(l, 0.0) # sanity: H <= max(h) G(T) # ---------- numerator: sum_i l_i^2 (C(s_{i+1}) - C(s_i)), s_i = i D pts = np.arange(ilo, itop + 1) * DELTA L1, U1 = cdf.bounds(k - 1, pts) t_cdf = time.time() - t0 - t_den # Abel: sum_{i} l_i^2 (C_{i+1} - C_i) = l_last^2 C_top - l_first^2 C_lo + sum_{i interior} (l_{i-1}^2 - l_i^2) C_i ctx.prec = 128 l2 = [arb(x) ** 2 for x in l] num = l2[-1] * arb(L1[-1]) - l2[0] * arb(U1[0]) for q in range(1, len(l)): coef = l2[q - 1] - l2[q] Cb = L1[q] if float(coef.mid()) >= 0 else U1[q] num += coef * arb(Cb) num_lo = float(np.nextafter(float(num.lower()), -np.inf)) # float(arb) rounds to nearest: step outward kk = arb(k) Q = kk / m2 * arb(num_lo) / arb(den_up) out = dict(k=k, theta=theta, beta=beta, eps_cells=eps_cells, nc=nc, Xi1=cdf.Xi1, cells=len(l), delta=DELTA, num_lower=num_lo, den_upper=den_up, den_lower=den_lo, Pk_lower=float(Lk[-1]), Q_lower=float(np.nextafter(float(Q.lower()), -np.inf)), max_cdf_width=float(max(np.max(U1 - L1), np.max(Uk - Lk))), secs=dict(den=round(t_den, 1), cdf=round(t_cdf, 1), total=round(time.time() - t0, 1))) if verbose: print(json.dumps(out), flush=True) return out if __name__ == "__main__": import pickle k, th, be = int(sys.argv[1]), float(sys.argv[2]), float(sys.argv[3]) import os; sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from hopt_lat import optimise r = optimise(k, th, be, ncell=int(sys.argv[4]) if len(sys.argv) > 4 else 200) print(json.dumps(dict(Q_fft=r["Q_fft"])), flush=True) hcert(k, th, be, r["E"], r["h"]) # control: h == 1 (plain product trial) hcert(k, th, be, np.array([0, int(round(1 / DELTA))]), np.array([1.0]))