#!/usr/bin/env python3
"""Zero/pole construction checks: exact Fraction certificates + 70-digit diagnostics.

Requires mpmath. The proof pages justify infinite limits; finite high-precision
checks are not interval enclosures or replacements for those proofs.
Use --output to write JSON; without it, print only.
"""
from fractions import Fraction as F
from collections import Counter
from pathlib import Path
from math import comb
import argparse,json
import mpmath as mp
mp.mp.dps=70

def elementary(p,w):
 return (1-w)*mp.exp(sum(w**k/k for k in range(1,p+1)))
def product_sine(N,z):return mp.pi*z*mp.fprod(1-z*z/(n*n) for n in range(1,N+1))
def partial_cot(N,z):return 1/z+sum(2*z/(z*z-n*n) for n in range(1,N+1))
def growing_log(N,z):return sum(n*(mp.log1p(-z/(n*n))+z/(n*n)) for n in range(1,N+1))
def double_pole_term(n,z):return 1/(z-n)**2+z*z/(n*(n-z))
def double_pole_sum(N,z):return sum(double_pole_term(n,z) for n in range(1,N+1))
def jensen_function(w):return w*w*mp.exp(w)*(1-2*w)*(1-w/2)
def angular_average(R,M):
 return sum(mp.log(abs(jensen_function(R*mp.exp(2j*mp.pi*j/M)))) for j in range(M))/M
def circle_integral(fn,center,radius,M=192):
 # Returns the trapezoidal diagnostic for integral/(2*pi*i), not an enclosure.
 return sum(radius*mp.exp(2j*mp.pi*j/M)*fn(center+radius*mp.exp(2j*mp.pi*j/M)) for j in range(M))/M
def quad_segment(fn,a,b):
 # Break long sides into pieces of length <= 1/2. A single 11-unit segment
 # underresolved the rapid tanh transition near its midpoint at 70 digits.
 steps=max(1,int(mp.ceil(2*abs(b-a))))
 return mp.quad(lambda t:fn(a+t*(b-a))*(b-a),[mp.mpf(i)/steps for i in range(steps+1)])
def s(x):return mp.nstr(x,48)

def run():
 counts=Counter()
 def ck(name,condition):
  if not condition:raise ArithmeticError((name,counts[name]))
  counts[name]+=1
 tolerance=mp.mpf('1e-55')
 # Explicit logarithmic correction and quantitative bound, for complex w.
 for p in range(7):
  for r in [mp.mpf(1)/8,mp.mpf(1)/4,mp.mpf(1)/2]:
   for j in range(16):
    w=r*mp.exp(2j*mp.pi*j/16)
    L=mp.log1p(-w)+sum(w**k/k for k in range(1,p+1))
    B=r**(p+1)/((p+1)*(1-r))
    ck('elementary_log_bound',abs(L)<=B+tolerance)
    ck('elementary_exponent_identity',abs(mp.exp(L)-elementary(p,w))<tolerance)
    ck('paired_elementary_factor',abs(elementary(1,w)*elementary(1,-w)-(1-w*w))<tolerance)
 points=[mp.mpf(1)/4,mp.mpf(1)/2,mp.mpf(3)/4,mp.mpc(-.5,.25),mp.mpc(1.5,.5),mp.mpc(-1,.5)]
 for N in [5,10,40,100]:
  for z in points:
   PN=product_sine(N,z);r=abs(z)
   ck('sine_normal_product_bound',abs(mp.sin(mp.pi*z)-PN)<=abs(PN)*mp.expm1(r*r/N)+tolerance)
   B=2*r/(N*(1-(r/(N+1))**2))
   ck('cot_uniform_tail_bound',abs(mp.pi*mp.cot(mp.pi*z)-partial_cot(N,z))<=B+tolerance)
 # A reference on |z|<1 obtained by changing two absolutely convergent sums.
 # zeta(2k-1)<2 bounds the omitted k-tail explicitly, independently of n-cutoff.
 K=220;zetas=[None,None]+[mp.zeta(2*k-1) for k in range(2,K+1)]
 for z in [mp.mpf(1)/4,mp.mpf(1)/2,mp.mpf(-3)/4,mp.mpc(.5,.25),mp.mpc(-.5,.5)]:
  ref=-sum(zetas[k]*z**k/k for k in range(2,K+1));r=abs(z)
  delta=2*r**(K+1)/((K+1)*(1-r))
  for N in [3,10,30,100]:
   tail=ref-growing_log(N,z);B=r*r/(2*N*N)
   ck('growing_multiplicity_log_tail',abs(tail)<=B+delta+tolerance)
   ck('growing_multiplicity_relative_tail',abs(mp.expm1(tail))<=mp.expm1(B+delta)+tolerance)
 for n in range(1,101):
  ck('uncorrected_growth_diverges',n*mp.log1p(mp.mpf(1)/(n*n))>=mp.mpf(1)/(2*n))
  # Residue of z/[n(z-n^2)] is n, checked exactly before taking a limit.
  ck('growing_zero_order_residue',F(n*n,n)==n)
 ck('capstone_truncation_budget',F(2,1415**2-2)<F(1,10**6))
 # Translated center and center multiplicity are incorporated in this G(w+1).
 averages=[]
 for R in [mp.mpf(1)/4,mp.mpf(3)/4,mp.mpf(1),mp.mpf(3)/2,mp.mpf(5)/2,mp.mpf(3)]:
  exact=2*mp.log(R)+sum(mp.log(R/a) for a in [mp.mpf(1)/2,mp.mpf(2)] if a<R)
  qs=[min(R/a,a/R) for a in [mp.mpf(1)/2,mp.mpf(2)]]
  direct=mp.quad(lambda theta:mp.log(abs(jensen_function(R*mp.exp(1j*theta)))),[0,mp.pi,2*mp.pi])/(2*mp.pi)
  ck('jensen_continuous_circle_average',abs(direct-exact)<mp.mpf('1e-50'))
  for M in [8,16,32,64]:
   actual=angular_average(R,M)
   predicted=exact+sum(mp.log1p(-q**M) for q in qs)/M
   B=sum(q**M/(1-q**M) for q in qs)/M
   ck('jensen_exact_discrete_average',abs(actual-predicted)<tolerance)
   ck('jensen_discrete_error_bound',-tolerance<=exact-actual<=B+tolerance)
   if M==16 and R in [1,3]:averages.append({'radius':s(R),'continuous_average':s(exact),'sample_average':s(actual),'error_bound':s(B)})
 # Explicit Taylor choices for genuinely complex poles and higher-order data.
 for index,(a,coeffs) in enumerate([(mp.mpc(2,1),[1]),(mp.mpc(3,-2),[1,-2]),(mp.mpc(-4,1),[1,2,3j])],1):
  radius=abs(a)/2;Mbound=sum(abs(c)/(abs(a)/4)**k for k,c in enumerate(coeffs,1));degree=0
  while 3*Mbound*mp.mpf(2)**(degree+1)/mp.mpf(3)**(degree+1)>mp.mpf(2)**(-index):degree+=1
  def principal(z):return sum(c/(z-a)**k for k,c in enumerate(coeffs,1))
  coeff=[sum(c*(-a)**(-k)*comb(k+j-1,j)*a**(-j) for k,c in enumerate(coeffs,1)) for j in range(degree+1)]
  B=3*Mbound*(mp.mpf(2)/3)**(degree+1)
  for rho in [radius/3,2*radius/3,radius]:
   for j in range(24):
    z=rho*mp.exp(2j*mp.pi*j/24);q=sum(c*z**j for j,c in enumerate(coeff))
    ck('mittag_taylor_constructive_bound',abs(principal(z)-q)<=B+tolerance)
  ck('mittag_taylor_choice_terminates',B<=mp.mpf(2)**(-index))
 # Contour identity uses all four oriented sides, not a reused partial sum.
 for N in [2,5]:
  R=mp.mpf(N)+mp.mpf(1)/2;vertices=[-R-1j*R,R-1j*R,R+1j*R,-R+1j*R]
  for z in [mp.mpc(.25,.5),mp.mpc(1.25,-.25)]:
   fn=lambda w:mp.pi*mp.cot(mp.pi*w)/(w*w-z*z)
   I=sum(quad_segment(fn,a,b) for a,b in zip(vertices,vertices[1:]+vertices[:1]))
   error=z*I/(2j*mp.pi)
   ck('cot_finite_oriented_contour_identity',abs(mp.pi*mp.cot(mp.pi*z)-partial_cot(N,z)-error)<mp.mpf('1e-50'))
   ck('cot_square_contour_bound',abs(error)<=8*abs(z)*R/(R*R-abs(z)**2))
 # Double poles and varying residues: inspect the two Laurent coefficients.
 pole_records=[]
 for n in [1,5,12]:
  residue=circle_integral(lambda z:double_pole_sum(24,z),n,mp.mpf(1)/4)
  leading=circle_integral(lambda z:(z-n)*double_pole_sum(24,z),n,mp.mpf(1)/4)
  ck('double_pole_residue',abs(residue+n)<mp.mpf('1e-50'))
  ck('double_pole_leading_coefficient',abs(leading-1)<mp.mpf('1e-50'))
  pole_records.append({'pole':n,'residue':s(residue),'leading_coefficient':s(leading)})
 for N in [5,20,50]:
  for z in points:
   r=abs(z);tail=sum(double_pole_term(n,z) for n in range(N+1,401))
   ck('double_pole_finite_tail_bound',abs(tail)<=(4+2*r*r)/N+tolerance)
 # Exact rational enclosure of cos(pi/4), independent of high precision.
 C10=F(1)
 for n in range(1,11):C10*=1-F(1,4*(2*n-1)**2)
 low=F(151,152)*C10;high=C10
 ck('cosine_rational_enclosure_lower',2*low*low<1)
 ck('cosine_rational_enclosure_upper',2*high*high>1)
 for N in [3,10,30,100]:
  for z in points:
   PN=mp.fprod(1-4*z*z/(2*n-1)**2 for n in range(1,N+1))
   ck('cosine_normal_product_tail',abs(mp.cos(mp.pi*z)-PN)<=abs(PN)*mp.expm1(2*abs(z)**2/(2*N-1))+tolerance)
 cap={'growing_zero_budget_N':1415,'growing_zero_relative_bound':str(F(2,1415**2-2)),
      'translated_jensen_at_16_samples':averages,'jensen_R1_exact_rational_bound':str(F(1,524280)),
      'jensen_R3_exact_rational_bound':str((F(1,6**16-1)+F(2**16,3**16-2**16))/16),
      'double_pole_laurent_diagnostics':pole_records,'cos_pi_quarter_rational_lower':str(low),'cos_pi_quarter_rational_upper':str(high)}
 return {'status':'PASS','assertions':sum(counts.values()),'groups':dict(counts),'precision_digits':70,'capstone':cap,
         'limits':['Fraction inequalities are exact; complex numerical quadrature is a 70-digit diagnostic, not an interval certificate.','Infinite convergence, zero orders, and meromorphic existence are established by the written proofs.','Finite tail checks do not claim to enumerate an infinite series.']}

if __name__=='__main__':
 parser=argparse.ArgumentParser();parser.add_argument('--output',type=Path);args=parser.parse_args();out=json.dumps(run(),ensure_ascii=False,indent=2)+'\n'
 if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(out)
 print(out,end='')