"""Independent integral checks for T563; run: sage -python checks.py. No independent method in this file uses the Selberg product or gamma values. For beta=2, Andreief's determinant expansion gives n! times the determinant of integral x^(i+j)w(x) dx. For beta=4, the confluent Vandermonde expansion gives n! times the Pfaffian of (j-i)*integral x^(i+j-1)w(x) dx. These moments are computed by integrating the polynomial (1-x)^(b-1). For beta=1, de Bruijn's ordered determinant identity gives n! times a Pfaffian of skew moments, augmented with single moments when n is odd. With x=sin(u)^2, the one-variable basis is 2 sin(u)^(2i+r) cos(u)^s, 0<=u<=pi/2. Its finite Fourier series is integrated term by term, twice for skew moments, with complex balls at 1024 bits. No quadrature, gamma calls, or rounded input arguments enter this method. The determinant identities follow by expanding determinants into permutations and integrating products, or products over the ordered chamber. Independent controls run first: unit interval, uniform two-particle Vandermonde moments, and the arcsine-weight cases pi and 4. A wrong n! normalization is rejected. The optional run(stored=records) verifies serialized values read back from the site. All checks precede submission and are rerun on the stored entries. """ import json import sys from collections import Counter from fractions import Fraction from functools import lru_cache from math import comb, factorial import numberdb.sage as numberdb from numberdb._write import to_text from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from sage.rings.complex_arb import ComplexBallField from generate import JacobiSelberg, PARAMETERS, exact_form from check import exactness, measure, agrees_with, names_its_rings CHECK_BITS = 1024 R = RealBallField(CHECK_BITS) C = ComplexBallField(CHECK_BITS) I = C(0, 1) def determinant(a): a = [row[:] for row in a] d = QQ(1) for i in range(len(a)): pivot = next((j for j in range(i, len(a)) if a[j][i] != 0), None) if pivot is None: return QQ(0) if pivot != i: a[i], a[pivot] = a[pivot], a[i] d = -d v = a[i][i] d *= v for j in range(i + 1, len(a)): t = a[j][i]/v for k in range(i + 1, len(a)): a[j][k] -= t*a[i][k] return d def pfaffian(a, field): a = [row[:] for row in a] out = field(1) for k in range(0, len(a), 2): pivot = max(range(k+1, len(a)), key=lambda j: abs(a[k][j])) if pivot != k+1: a[k+1], a[pivot] = a[pivot], a[k+1] for row in a: row[k+1], row[pivot] = row[pivot], row[k+1] out = -out v = a[k][k+1] if field is QQ: assert v != 0 else: assert not v.contains_zero() out *= v for i in range(k+2, len(a)): for j in range(i+1, len(a)): a[i][j] += (a[k+1][i]*a[k][j] - a[k][i]*a[k+1][j])/v a[j][i] = -a[i][j] return out @lru_cache(None) def polynomial_moment(a, b, k): return sum((QQ((-1)**t * comb(b-1,t))/QQ(a+k+t) for t in range(b)), QQ(0)) def even_beta_integral(beta, n, r, s): a, b = (beta//2)*(r+1), (beta//2)*(s+1) if beta == 2: mat = [[polynomial_moment(a,b,i+j) for j in range(n)] for i in range(n)] return QQ(factorial(n))*determinant(mat) assert beta == 4 mat = [[QQ(0) for j in range(2*n)] for i in range(2*n)] for i in range(2*n): for j in range(i+1,2*n): mat[i][j] = (j-i)*polynomial_moment(a,b,i+j-1) mat[j][i] = -mat[i][j] return QQ(factorial(n))*pfaffian(mat,QQ) @lru_cache(None) def fourier(p, q): phase = I**(-p) out = {} for j in range(p+q+1): coeff = sum((-1)**t * comb(p,t) * comb(q,j-t) for t in range(max(0,j-q),min(p,j)+1)) if coeff: out[p+q-2*j] = C(QQ(2*coeff)/QQ(2**(p+q))) * phase return out @lru_cache(None) def exp_integral(k): if k == 0: return C(R.pi()/2) return (I**k - 1)/(I*k) @lru_cache(None) def u_exp_integral(k): if k == 0: return C(R.pi()**2/8) return (R.pi()/2)*I**k/(I*k) + (I**k-1)/(k*k) @lru_cache(None) def triangular_exp(k,l): if l == 0: return u_exp_integral(k) return (exp_integral(k+l)-exp_integral(k))/(I*l) def beta_one_integral(n,r,s): fs = [fourier(2*i+r,s) for i in range(n)] frequencies = sorted(set().union(*(set(f) for f in fs))) # Integral of outer Fourier series against each inner exponential primitive. primitive_pair = [{l:sum((c*triangular_exp(k,l) for k,c in f.items()),C(0)) for l in frequencies} for f in fs] pair = [[sum((c*primitive_pair[i][l] for l,c in fs[j].items()),C(0)) for j in range(n)] for i in range(n)] size = n+n%2 mat = [[C(0) for j in range(size)] for i in range(size)] for i in range(n): for j in range(i+1,n): mat[i][j] = pair[j][i]-pair[i][j] mat[j][i] = -mat[i][j] if n%2: for i,f in enumerate(fs): mat[i][n] = sum((c*exp_integral(k) for k,c in f.items()),C(0)) mat[n][i] = -mat[i][n] value = factorial(n)*pfaffian(mat,C) assert value.real().is_finite() and value.imag().is_finite() and value.imag().contains_zero() real = value.real() assert real > 0 and real.accuracy() >= 400, (n,r,s,str(real)) return real def independent(beta,n,r,s): if beta == 1: return beta_one_integral(n,r,s) return even_beta_integral(beta,n,r,s) def decimal_bounds(text): if '.' not in text and 'e' not in text.lower(): q = QQ(text) return q,q mantissa, _, exponent = text.lower().partition('e') scale = int(exponent or 0)-len(mantissa.partition('.')[2]) center = QQ(Fraction(text)) unit = QQ(10)**scale return center-unit,center+unit def enclosed(text,value): lo,hi = decimal_bounds(text) if lo == hi: return value.__contains__(lo) if hasattr(value, 'rad') else value == lo ball = R(value) return bool(ball >= R(lo) and ball <= R(hi)) def controls(): assert determinant([[QQ(1),QQ(2)],[QQ(3),QQ(4)]]) == -2 assert pfaffian([[QQ(0),QQ(2),QQ(3),QQ(5)],[-QQ(2),QQ(0),QQ(7),QQ(11)], [-QQ(3),-QQ(7),QQ(0),QQ(13)],[-QQ(5),-QQ(11),-QQ(13),QQ(0)]],QQ) == 28 assert even_beta_integral(2,1,0,0)==1 assert even_beta_integral(2,2,0,0)==QQ(1)/6 assert even_beta_integral(4,1,0,0)==QQ(1)/6 # Uniform beta=4 weight can be tested directly with moments a=b=1. uniform4 = [[QQ(j-i)/QQ(i+j) if i!=j else QQ(0) for j in range(4)] for i in range(4)] assert 2*pfaffian(uniform4,QQ)==QQ(1)/15 assert beta_one_integral(1,0,0).overlaps(R.pi()) assert beta_one_integral(2,0,0).__contains__(4) assert beta_one_integral(2,1,1).__contains__(QQ(1)/3) assert not beta_one_integral(2,1,1).__contains__(QQ(1)/6) assert enclosed('3.14159', R.pi()) assert not enclosed('3.14169', R.pi()) return ['determinant=-2','Pfaffian=28','uniform beta=2 n=2 integral=1/6', 'uniform beta=4 n=2 integral=1/15','arcsine n=1 integral=pi', 'arcsine n=2 integral=4','uniform beta=1 n=2 integral=1/3', 'ordered-chamber normalization rejected','shifted decimal rejected'] def run(stored=None): control_results = controls() print('Controls passed.',file=sys.stderr,flush=True) generator=JacobiSelberg() params=list(generator.enumerate()) keys=[tuple(int(p[k]) for k in PARAMETERS) for p in params] assert len(keys)==len(set(keys)) if stored is not None: supplied={tuple(int(rec['params'][k]) for k in PARAMETERS):str(rec['number']) for rec in stored} assert list(supplied)==keys, 'Stored order or selection differs' values={key:generator.value(p,100) for key,p in zip(keys,params)} assert not exactness(values) assert not names_its_rings(__import__('generate').__file__) records=[]; independent_exact={}; exact_values={}; counts=Counter(); accuracies=[]; independent_accuracy=[] for idx,(key,value) in enumerate(values.items()): beta,n,r,s=key q,k=exact_form(*key) assert (q,k)==exact_form(beta,n,s,r) other=independent(*key) if beta in (2,4): independent_exact[key]=other; exact_values[key]=value assert other==value, key else: expected=R(q)*R.pi()**k assert other.overlaps(expected),key assert R(value).overlaps(other),key independent_accuracy.append(int(other.accuracy())) txt = supplied[key] if stored is not None else to_text(value,digits=100) assert enclosed(txt,other), (key,txt) assert enclosed(txt,value), (key,txt) if k == 0: assert txt==str(q), (key,txt) counts['exact']+=1 else: accuracy=int(value.accuracy());accuracies.append(accuracy) assert accuracy>=330 counts['inexact']+=1 counts['beta='+str(beta)]+=1 records.append({'params':dict(zip(PARAMETERS,key)),'number':txt, 'rational_factor':str(q),'pi_power':int(k)}) if (idx+1)%50==0: print('Independent integrals checked:',idx+1,file=sys.stderr,flush=True) assert not agrees_with(exact_values,lambda key:independent_exact[key]) summary={'table':'T563','entries':len(values),'controls':control_results, 'counts':dict(counts),'measure':measure(values), 'longest_serialized':max(len(r['number']) for r in records), 'minimum_generator_accuracy_bits':min(accuracies), 'minimum_independent_beta1_accuracy_bits':min(independent_accuracy), 'check_bits':CHECK_BITS,'generator_bits':numberdb.bits(100,losing=96), 'independent_integrals_checked':len(records),'symmetry_checked':len(records), 'serialized_enclosures_checked':len(records),'stored_readback':stored is not None, 'numberdb_version':numberdb.__version__, 'skill_manifest':'../20261027T223837Z-build-skills/manifest.json', 'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'} return {'summary':summary,'records':records} if __name__=='__main__': print(json.dumps(run(),indent=2))