"""Check the actual NumPy cell-table implementation against Arb at every
J_m(w) used by the two frozen grids and JS. Record direct-convolution source.
"""
import os
for k in ['OPENBLAS_NUM_THREADS','OMP_NUM_THREADS','MKL_NUM_THREADS','NUMEXPR_NUM_THREADS']:
    os.environ[k]='1'
import inspect, math, json
import numpy as np
from flint import arb, ctx
import audit_v6 as A
ctx.prec=256; ctx.threads=1
source=inspect.getsource(np.convolve)
open('numpy_convolve_source.txt','w').write(source)
assert 'multiarray.correlate' in source
A.cells()
H=200; nmax=59*H
c=(1.0/(H+np.arange(nmax+1,dtype=np.float64)))*(1+1e-15)
assert all(arb(float(v))>=arb(1)/(H+i) for i,v in enumerate(c))
p=c.copy(); comparisons=0; min_excess=None
for m in range(1,A.WMAX):
    if m>1: p=np.convolve(p,c)[:nmax+1]*(1+1e-10)
    for w in range(m,A.WMAX+1):
        q=arb(float(np.sum(p[:(w-m)*H+1]))*(1+1e-10))
        exact=A.J[m,w]
        assert q>=exact
        excess=q/exact-1
        if min_excess is None or excess<min_excess: min_excess=excess
        comparisons+=1
u=arb(1)/2**53
print('weights verified',len(c),'J values verified',comparisons,'least relative excess',min_excess)
for n in [7601,11801,15201]:
    gamma=2*n*u/(1-2*n*u)
    print('length',n,'gamma_2n',gamma,'inflation margin',(arb(float(1+1e-10))*(1-u)*(1-gamma)-1))
print('WMAX60 normality: min product >= 12000^-59 =',arb(12000)**-59)
print('WMAX60 magnitude: coefficient,sum < 60^59 =',arb(60)**59)