"""Independent checks for centroid distances -- numberdb.org/T469. Run beside generate.py and check.py with: sage -python checks.py sage -python checks.py --stored For a draft, the second command needs NUMBERDB_API_KEY. Set NUMBERDB_KEY_FROM_STDIN=1 to read it on stdin instead; that also selects stored verification for a runner that cannot pass command-line arguments. The generator evaluates closed forms in Arb. This file instead integrates radial densities for polygons and balls and polar coordinates for rectangles. For the cube, splitting into six face pyramids and integrating along rays reduces the volume integral to one smooth angular integral. All quadratures are checked first at powers 0 and 2, whose answers are independently known. They are repeated at 130 and 170 decimal digits; no quadrature error estimate is used as a proof. Certified bounds come from the generator's Arb balls. A seeded Monte Carlo check draws points directly in Cartesian coordinates, or by area-preserving triangular sampling, and tests centroids and second moments before the mean. Those low-precision samples provide a separate check of the measure and the size convention. They do not generate stored digits. """ import os import sys import json from decimal import Decimal from pathlib import Path import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField import mpmath as mp import numpy as np from generate import CentroidDistance, compute, configure_stdin from check import exactness, measure, names_its_rings # Independent enumeration and geometry, deliberately not imported from generate. REGIONS = ('unit-interval', 'unit-square', 'rectangle-1x2', 'equilateral-triangle', 'unit-disk', 'regular-pentagon', 'regular-hexagon', 'regular-heptagon', 'regular-octagon', 'unit-cube', 'unit-ball-3', 'unit-circle', 'unit-sphere-2') N_SIDES = {'equilateral-triangle': 3, 'regular-pentagon': 5, 'regular-hexagon': 6, 'regular-heptagon': 7, 'regular-octagon': 8} RATIONALS = {'unit-interval': QQ(1)/4, 'unit-disk': QQ(2)/3, 'unit-ball-3': QQ(3)/4, 'unit-circle': QQ(1), 'unit-sphere-2': QQ(1)} def second_moment(region, scale=1): if region in N_SIDES: n = N_SIDES[region] radius = 1/mp.sqrt(3) if n == 3 else mp.mpf(1) h = radius*mp.cos(mp.pi/n) answer = (radius**2+2*h**2)/6 else: answer = {'unit-interval': mp.mpf(1)/12, 'unit-square': mp.mpf(1)/6, 'rectangle-1x2': mp.mpf(5)/12, 'unit-disk': mp.mpf(1)/2, 'unit-cube': mp.mpf(1)/4, 'unit-ball-3': mp.mpf(3)/5, 'unit-circle': mp.mpf(1), 'unit-sphere-2': mp.mpf(1)}[region] return scale**2*answer def moment(region, p, scale=1): """Integrate the uniform geometric measure, without the closed forms.""" scale = mp.mpf(scale) if region in N_SIDES: n = N_SIDES[region] theta = mp.pi/n radius = scale/mp.sqrt(3) if n == 3 else scale h = radius*mp.cos(theta) area = n*radius**2*mp.sin(2*theta)/2 # Angular measure at radius r is 2*pi below h, then # 2*pi-2*n*acos(h/r). Substitute r=h*sec(t) in the outer annulus. inner = 2*mp.pi*h**(p+2)/(p+2) outer = mp.quad(lambda t: 2*(mp.pi-n*t)*(h/mp.cos(t))**(p+2)*mp.tan(t), [0, theta]) return (inner+outer)/area if region in ('unit-square', 'rectangle-1x2'): a = scale/2 b = scale/2 if region == 'unit-square' else scale theta = mp.atan(b/a) return (mp.quad(lambda t: (a/mp.cos(t))**(p+2), [0, theta])+ mp.quad(lambda t: (b/mp.sin(t))**(p+2), [theta, mp.pi/2])) / ((p+2)*a*b) if region == 'unit-cube': angular = mp.quad(lambda t: (1+mp.sec(t)**2)**(mp.mpf(p+2)/2)-1, [0, mp.pi/4]) return 6*(scale/2)**p*angular/((p+3)*(p+2)) if region == 'unit-interval': return mp.quad(lambda r: 2*r**p/scale, [0, scale/2]) if region in ('unit-disk', 'unit-ball-3'): d = 2 if region == 'unit-disk' else 3 return mp.quad(lambda r: d*r**(p+d-1)/scale**d, [0, scale]) if region == 'unit-circle': return mp.quad(lambda t: (scale*mp.sqrt(mp.cos(t)**2+mp.sin(t)**2))**p, [0, 2*mp.pi])/(2*mp.pi) if region == 'unit-sphere-2': return mp.quad(lambda u: (scale*mp.sqrt(u*u+(1-u*u)))**p, [-1, 1])/2 raise ValueError(region) def rational_to_mp(q): return mp.mpf(int(q.numerator()))/int(q.denominator()) def ball_mid(value): return rational_to_mp(QQ(value.mid())) def contained_written(text, estimate): text = str(text) if '/' in text or ('.' not in text and 'e' not in text.lower()): return abs(estimate-rational_to_mp(QQ(text))) < mp.mpf('1e-120') decimal = Decimal(text) center = mp.mpf(text) ulp = mp.power(10, decimal.as_tuple().exponent) return center-ulp <= estimate <= center+ulp def decimal_qq(text): sign, digits, exponent = Decimal(text).as_tuple() integer = int(''.join(map(str, digits))) * (-1 if sign else 1) return QQ(integer)*QQ(10)**exponent def strict_stored_ball(text, value): if '/' in text or ('.' not in text and 'e' not in text.lower()): return value == QQ(text) center = decimal_qq(text) ulp = QQ(10)**Decimal(text).as_tuple().exponent field = value.parent() return field(center-ulp) < value and value < field(center+ulp) def simulate(means): rng = np.random.default_rng(46920261024) size = 300000 rows = [] for region in REGIONS: if region in N_SIDES: n = N_SIDES[region] radius = 1/np.sqrt(3) if n == 3 else 1 sector = rng.integers(n, size=size) angle1 = 2*np.pi*sector/n angle2 = 2*np.pi*(sector+1)/n a = radius*np.stack([np.cos(angle1), np.sin(angle1)], axis=1) b = radius*np.stack([np.cos(angle2), np.sin(angle2)], axis=1) u, v = rng.random((2, size)) rho = np.sqrt(u) points = rho[:, None]*((1-v[:, None])*a+v[:, None]*b) elif region in ('unit-interval', 'unit-square', 'rectangle-1x2', 'unit-cube'): dim = {'unit-interval': 1, 'unit-square': 2, 'rectangle-1x2': 2, 'unit-cube': 3}[region] points = rng.random((size, dim))-0.5 if region == 'rectangle-1x2': points[:, 1] *= 2 else: dim = 2 if region in ('unit-disk', 'unit-circle') else 3 points = rng.normal(size=(size, dim)) points /= np.linalg.norm(points, axis=1)[:, None] if region in ('unit-disk', 'unit-ball-3'): points *= rng.random(size=size)[:, None]**(1.0/dim) distances = np.linalg.norm(points, axis=1) centroid_z = np.abs(points.mean(axis=0))/(points.std(axis=0, ddof=1)/np.sqrt(size)) assert np.max(centroid_z) < 6, (region, 'centroid control', centroid_z) if region in ('unit-circle', 'unit-sphere-2'): assert np.max(np.abs(distances-1)) < 1e-14 mean_z = second_z = 0.0 else: second_z = abs(np.mean(distances**2)-float(second_moment(region)))/(np.std(distances**2, ddof=1)/np.sqrt(size)) assert second_z < 6, (region, 'second-moment control', second_z) mean_z = abs(distances.mean()-float(means[region]))/(distances.std(ddof=1)/np.sqrt(size)) assert mean_z < 6, (region, 'mean', mean_z) rows.append({'region': region, 'mean_z': float(mean_z), 'centroid_max_z': float(np.max(centroid_z)), 'second_moment_z': float(second_z)}) print('MONTE_CARLO', json.dumps({'samples_per_region': size, 'seed':46920261024, 'rows': rows})) def main(): configure_stdin() generator = CentroidDistance() assert tuple(p['region'] for p in generator.enumerate()) == REGIONS values = {r: generator.value({'region': r}, 100) for r in REGIONS} assert not exactness(values), exactness(values) assert not names_its_rings(str(Path(__file__).with_name('generate.py'))) print('MEASUREMENT', json.dumps(measure(values))) numerical = {} for digits in (130, 170): mp.mp.dps = digits estimates = {} for region in REGIONS: assert abs(moment(region, 0)-1) < mp.mpf('1e-120'), (region, 'mass control') assert abs(moment(region, 2)-second_moment(region)) < mp.mpf('1e-120'), (region, 'second-moment control') value = moment(region, 1) assert 0 < value <= mp.sqrt(second_moment(region))+mp.mpf('1e-120') assert abs(moment(region, 1, mp.mpf(3)/2)-mp.mpf(3)*value/2) < mp.mpf('1e-120'), (region, 'scaling') if region in numerical: assert abs(numerical[region]-value) < mp.mpf('1e-120'), (region, 'precision agreement') estimates[region] = value numerical = estimates print('QUADRATURE', digits, 'digits: all 13 means, normalization, second moments and scaling passed', flush=True) report = [] for region in REGIONS: value = values[region]['number'] if region in RATIONALS: assert value == RATIONALS[region] center = rational_to_mp(value) radius = mp.mpf(0) else: assert value.is_finite() center = ball_mid(value) radius = rational_to_mp(QQ(value.rad())) assert 0 < radius/abs(center) < mp.mpf('1e-120') assert abs(numerical[region]-center) < mp.mpf('1e-120'), (region, 'closed form versus quadrature') report.append({'region':region, 'mean':mp.nstr(numerical[region], 125), 'arb_radius':mp.nstr(radius, 10), 'quadrature_difference':mp.nstr(abs(numerical[region]-center),10)}) # The simplified expressions printed in the triangle and hexagon comments. assert abs(numerical['equilateral-triangle']-(mp.sqrt(3)/9+mp.log(2+mp.sqrt(3))/18)) < mp.mpf('1e-120') assert abs(numerical['regular-hexagon']-(mp.mpf(1)/3+mp.log(3)/4)) < mp.mpf('1e-120') print('EVIDENCE', json.dumps(report), flush=True) simulate(numerical) if '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1': document = numberdb.table('T469') numbers = document['Numbers'] assert tuple(numbers) == REGIONS for region in REGIONS: text = str(numbers[region]['number']) assert contained_written(text, numerical[region]), (region, 'quadrature outside stored interval') assert strict_stored_ball(text, values[region]['number']), (region, 'Arb enclosure outside stored interval') verification = generator.verify(sample=None) assert verification.ok, str(verification) print('VERIFY', verification) print('STORED: all 13 read-back entries contain their independent quadrature values and full certified balls.') print('PASS') if __name__ == '__main__': main()