1import math, json, random
 2import numpy as np
 3import torch
 4
 5SEED=417
 6random.seed(SEED); np.random.seed(SEED); torch.manual_seed(SEED)
 7
 8def f(r): return 1.0/((r+2.0)*math.log(r+2.0)**2)
 9
10def math_check():
11    # Add the analytic integral remainder after the finite summation cutoff.
12    M=2_000_000
13    rs=np.arange(1,M+1,dtype=np.float64)
14    vals=1.0/((rs+2.0)*np.log(rs+2.0)**2)
15    tails=[]
16    for R in [8,16,32,64,128,256,1024]:
17        finite=float(vals[R:].sum())
18        remainder=1.0/math.log(M+2.0)
19        corrected=finite+remainder
20        estimate=1.0/math.log(R+2.0)
21        tails.append({'R':R,'tail_with_remainder':corrected,'integral_estimate':estimate,'ratio':corrected/estimate})
22    bands=[]
23    for k in range(1,9):
24        lo,hi=2**k,2**(k+1)
25        exact=sum(f(r) for r in range(lo,hi))
26        integral=1.0/math.log(lo+2.0)-1.0/math.log(hi+2.0)
27        asymptotic=1.0/(k*(k+1)*math.log(2.0))
28        bands.append({'k':k,'exact':exact,'shifted_integral':integral,'asymptotic_dyadic_mass':asymptotic,'exact_over_integral':exact/integral})
29    return {'tails':tails,'bands':bands}
30
31def exact_kernel(L):
32    # symmetric normalized off-diagonal kernel, including only distances available in length L
33    a=np.array([0.0]+[f(r) for r in range(1,L)],dtype=np.float64)
34    return a/(2*a[1:].sum())
35
36def approx_dyadic(L):
37    # Effective coefficient for each distance induced by band averaging and fixed band masses.
38    q=np.zeros(L)
39    for k in range(int(math.log2(L))):
40        lo,hi=2**k,min(2**(k+1),L)
41        mass=1.0/((k+1)*(k+2)*math.log(2.0))
42        q[lo:hi]+=mass/(hi-lo)
43    q=q/(2*q[1:].sum())
44    return q
45
46def local(L,R=8):
47    q=np.zeros(L); q[1:R+1]=1.0/(2*R); return q
48
49def compare():
50    rows=[]
51    for L in [64,256,1024,4096]:
52        p=exact_kernel(L)
53        d=approx_dyadic(L); c=local(L)
54        # relative L1 discrepancy and mass beyond the local radius
55        rows.append({'L':L,'dyadic_l1_error':float(np.abs(p-d).sum()),'local_l1_error':float(np.abs(p-c).sum()),'true_mass_r_gt_8':float(p[9:].sum()),'dyadic_mass_r_gt_8':float(d[9:].sum())})
56    return rows
57
58print(json.dumps({'math':math_check(),'kernel_approximation':compare()},indent=2))