verify.py
12.8 kB · python · 357 lines
1import cmath2import math3from decimal import Decimal, getcontext4from fractions import Fraction56# CONSTANTS78GAMMA = 0.57721566490153286060651209008240243104219GUARD = 1e-1210RELGUARD = 1e-911DPREC = 6012DGAMMA = Decimal("0.57721566490153286060651209008240243104215933593992")13DPI = Decimal("3.14159265358979323846264338327950288419716939937511")1415getcontext().prec = DPREC1617def same(got, want, what):18 assert got == want, "%s: got %r want %r" % (what, got, want)1920def atmost(got, want, what):21 assert got <= want, "%s: got %.12g want at most %.12g" % (what, got, want)2223def atleast(got, want, what):24 assert got >= want, "%s: got %.12g want at least %.12g" % (what, got, want)2526def near(got, want, tol, what):27 assert abs(got - want) <= tol, "%s: got %.12g want %.12g" % (what, got, want)2829def below(got, want, what):30 assert got < want, "%s: got %.12g want below %.12g" % (what, got, want)3132def sci(mantissa, exponent):33 whole, frac = mantissa.split(".")34 return int(whole + frac) * 10 ** (exponent - len(frac))3536# THE FORMULAS3738def harmonic(n):39 return math.log(n) + GAMMA + 1.0 / (2 * n)4041def phi(q):42 n = -((-(q - 2)) // 2)43 return (4 / math.pi) * q + (2 * q / math.pi) * harmonic(n) + (1 - 2 / math.pi) * (q - 2) + 0.7274445def proved(q, m=1):46 return math.sqrt(m) + phi(q) / q4748def alpha(q, m=1):49 return math.log(q - m) / math.log(q)5051def cost(q, m=1):52 return math.log(proved(q, m)) / math.log(q)5354def saving(q, m=1):55 a = alpha(q, m)56 return (a - 0.75 - cost(q, m)) / a5758def gap(q, m=1, b=0.75):59 return (q - m) * q ** (-b) - proved(q, m)6061def savingfree(q, m=1):62 return math.log1p(gap(q, m) / proved(q, m)) / (alpha(q, m) * math.log(q))6364def down(x, d):65 return math.floor(x * 10 ** d + GUARD) / 10 ** d6667def up(x, d):68 return math.ceil(x * 10 ** d - GUARD) / 10 ** d6970def row(q, m=1):71 return ("%.6f" % down(alpha(q, m), 6), "%.5f" % up(cost(q, m), 5), "%.5f" % down(saving(q, m), 5))7273# THE FORMULAS IN DECIMAL7475def dphi(q):76 q = Decimal(q)77 n = -((-(int(q) - 2)) // 2)78 hn = Decimal(n).ln() + DGAMMA + 1 / (2 * Decimal(n))79 return (4 / DPI) * q + (2 * q / DPI) * hn + (1 - 2 / DPI) * (q - 2) + Decimal("0.727")8081def dproved(q, m=1):82 return Decimal(m).sqrt() + dphi(q) / Decimal(q)8384def dgap(q, b, m=1):85 q = Decimal(q)86 b = Decimal(b.numerator) / Decimal(b.denominator)87 return (q - m) * (-b * q.ln()).exp() - dproved(int(q), m)8889def dmono(q, b):90 q = Decimal(q)91 b = Decimal(b.numerator) / Decimal(b.denominator)92 return (1 - b) * (q - 2) * (-b * (q + 1).ln()).exp()9394def bisect(test, lo, hi):95 while lo < hi:96 mid = (lo + hi) // 297 if test(mid):98 hi = mid99 else:100 lo = mid + 1101 return lo102103def floorbase(b):104 return bisect(lambda q: dmono(q, b) >= Decimal("1.291"), 40, 10 ** 90)105106def wallbase(b, start):107 return bisect(lambda q: dgap(q, b) > 0, start, 10 ** 90)108109# THE LADDER EXPONENTS110111def bh(a):112 if a < Fraction(11, 20):113 return a + Fraction(1, 4)114 if a < Fraction(3, 5):115 return Fraction(4, 5)116 return (a + 1) / 2117118def zhang(a):119 return (8 * a - 7 * a * a) / (4 - 2 * a)120121def bexp(a):122 v = bh(a)123 if a <= Fraction(4, 7):124 v = min(v, zhang(a))125 return v126127RUNGS = [128 (Fraction(1, 2), Fraction(3, 4), "both", 3690, 723),129 (Fraction(13, 25), Fraction(1417, 1850), "Zhang", 8578, 1486),130 (Fraction(11, 20), Fraction(913, 1160), "Zhang", 33547, 4754),131 (Fraction(4, 7), Fraction(4, 5), "both", 92317, 11221),132 (Fraction(3, 5), Fraction(4, 5), "BH", 92317, 11221),133 (Fraction(2, 3), Fraction(5, 6), "BH", 3107080, 216023),134 (Fraction(3, 4), Fraction(7, 8), "BH", 6939524168, 129458304),135 (Fraction(4, 5), Fraction(9, 10), "BH", None, 128606353005),136 (Fraction(9, 10), Fraction(19, 20), "BH", None, None),137 (Fraction(19, 20), Fraction(39, 40), "BH", None, None),138]139140CAPPED = {141 Fraction(4, 5): (sci("3.09358", 13), None),142 Fraction(9, 10): (sci("3.23663", 34), sci("1.73431", 28)),143 Fraction(19, 20): (sci("9.24614", 83), sci("3.30712", 68)),144}145146# THE TABLE147148TABLE = [149 (1000, "0.999855", "0.28087", "-0.03102", False),150 (2000, "0.999934", "0.26335", "-0.01342", False),151 (3000, "0.999958", "0.25430", "-0.00434", False),152 (3689, "0.999966", "0.24997", "-0.00001", False),153 (3690, "0.999967", "0.24997", "0.00000", True),154 (5000, "0.999976", "0.24393", "0.00605", True),155 (10 ** 4, "0.999989", "0.23141", "0.01858", True),156 (10 ** 5, "0.999999", "0.19906", "0.05094", True),157 (10 ** 6, "0.999999", "0.17589", "0.07411", True),158 (10 ** 9, "0.999999", "0.13305", "0.11695", True),159]160161MULTI = [162 (10 ** 4, 6, "0.999934", "0.24865", "0.00129"),163 (10 ** 5, 78, "0.999932", "0.24972", "0.00022"),164 (10 ** 6, 451, "0.999967", "0.24994", "0.00003"),165]166167BUDGET = [168 (Fraction(1, 2), 1971),169 (Fraction(13, 25), 1002),170 (Fraction(11, 20), 365),171 (Fraction(4, 7), 176),172 (Fraction(3, 5), 176),173 (Fraction(2, 3), 8),174]175176def check_table():177 for q, a, c, d, closes in TABLE:178 same(row(q), (a, c, d), "table row q = %d" % q)179 same(gap(q) > 0, closes, "table gap sign q = %d" % q)180 for q in range(3, 20000):181 same(0.75 + cost(q) < alpha(q), gap(q) > 0, "test equivalence q = %d" % q)182183def check_wall():184 for q in range(3, 3690):185 below(gap(q), 0.0, "gap below the wall q = %d" % q)186 atleast(gap(3690), 3.752213034e-4, "gap at the wall")187 atmost(gap(3689), -1.533059397e-4, "gap under the wall")188 near(gap(3690), 3.752213034e-4, RELGUARD * 3.752213034e-4, "gap 3690 digits")189 near(gap(3689), -1.533059397e-4, RELGUARD * 1.533059397e-4, "gap 3689 digits")190 atleast(savingfree(3690), 5.863425182e-6, "delta at the wall")191 atmost(savingfree(3689), -2.395807653e-6, "delta under the wall")192 near(savingfree(3690), 5.863425182e-6, RELGUARD * 5.863425182e-6, "delta 3690 digits")193 near(savingfree(3689), -2.395807653e-6, RELGUARD * 2.395807653e-6, "delta 3689 digits")194 near(savingfree(3690), saving(3690), 1e-12, "the two forms of delta agree")195 steps = [gap(q + 1) - gap(q) for q in range(3690, 100000)]196 same(len(steps), 96310, "step count of the sweep")197 same(sum(1 for s in steps if s <= 0), 0, "steps that fail to rise")198 atleast(min(steps), 3.172e-5, "smallest step of the sweep")199 same(3690 + steps.index(min(steps)), 99998, "where the smallest step sits")200201def widest(q, b=0.75):202 m = 1203 while gap(q, m, b) > 0:204 m += 1205 return m - 1206207def check_multi():208 for q, m, a, c, d in MULTI:209 same(widest(q), m, "widest digit set at q = %d" % q)210 same(row(q, m), (a, c, d), "multi row q = %d" % q)211 for a, m in BUDGET:212 same(widest(10 ** 7, float(bexp(a))), m, "ladder budget at a = %s" % a)213214def check_exponents():215 for (a, b, source, _, _) in RUNGS:216 same(bexp(a), b, "b(a) at a = %s" % a)217 if source == "both":218 same(bh(a), zhang(a), "the tables agree at a = %s" % a)219 elif source == "Zhang":220 assert zhang(a) < bh(a), "Zhang is not smaller at a = %s" % a221 for den in range(1, 201):222 for num in range(1, den):223 a = Fraction(num, den)224 if not (Fraction(1, 2) < a <= Fraction(4, 7)):225 continue226 z = zhang(a)227 same(z - (a + Fraction(1, 4)), -5 * (a - Fraction(1, 2)) * (a - Fraction(2, 5)) / (4 - 2 * a), "first factorisation")228 same(z - Fraction(4, 5), -7 * (a - Fraction(4, 7)) * (a - Fraction(4, 5)) / (4 - 2 * a), "second factorisation")229 same(z - Fraction(3, 4), -7 * (a - Fraction(1, 2)) * (a - Fraction(6, 7)) / (4 - 2 * a), "third factorisation")230 assert z > Fraction(3, 4), "Zhang below three quarters at a = %s" % a231 if a < Fraction(11, 20):232 assert z < a + Fraction(1, 4), "Zhang not smaller at a = %s" % a233 elif a < Fraction(4, 7):234 assert z < Fraction(4, 5), "Zhang not smaller at a = %s" % a235 for den in range(1, 21):236 for num in range(1, den):237 a = Fraction(num, den)238 if Fraction(1, 2) <= a < 1:239 assert bexp(a) >= Fraction(3, 4), "b(a) below three quarters at a = %s" % a240241def check_ladder():242 for (a, b, _, wall, floor) in RUNGS:243 got = floorbase(b)244 if floor is not None:245 same(got, floor, "Q(b) at a = %s" % a)246 assert dmono(floor, b) >= Decimal("1.291"), "mono at Q, a = %s" % a247 assert dmono(floor - 1, b) < Decimal("1.291"), "mono at Q - 1, a = %s" % a248 else:249 atmost(got, CAPPED[a][1], "Q(b) upper bound at a = %s" % a)250 for q in range(3, 3690):251 below(gap(q, 1, float(b)), 0.0, "no rung closes at q = %d, a = %s" % (q, a))252 if wall is not None:253 assert dgap(wall, b) > 0, "gap at the wall, a = %s" % a254 assert dgap(wall - 1, b) < 0, "gap under the wall, a = %s" % a255 same(wallbase(b, got), wall, "q_0(a) by bisection at a = %s" % a)256 if wall <= 10 ** 5:257 scan = next(q for q in range(3, wall + 1) if gap(q, 1, float(b)) > 0)258 same(scan, wall, "q_0(a) by exhaustive scan at a = %s" % a)259 else:260 assert dgap(CAPPED[a][0], b) > 0, "printed upper bound is not a wall at a = %s" % a261 logs = [math.ceil(100 * math.log10(w if w is not None else CAPPED[a][0])) / 100 for (a, _, _, w, _) in RUNGS]262 same(["%.2f" % v for v in logs], ["3.57", "3.94", "4.53", "4.97", "4.97", "6.50", "9.85", "13.50", "34.52", "83.97"], "the log10 trend")263264# THE GENERAL FLOOR265266def lowproved(q):267 return 1 + 4 / math.pi + (2 / math.pi) * (math.log((q - 2) / 2) + GAMMA) + (1 - 2 / math.pi) * (q - 2) / q268269def major(q, b):270 return q ** (1 - b) - lowproved(q)271272def check_floor():273 near(1.291 - (2 / math.pi) * math.log(1.291 * 4) - 2.544 * 0.25, -0.39014598, 1e-8, "the bracket at u = 1/4")274 atmost(4 * (1.291 - (2 / math.pi) * math.log(1.291 * 4) - 2.544 * 0.25), -1.56, "four times the bracket")275 atmost(2 * 723 ** -0.75, 0.015, "the minimality slack")276 kappa = 1 + 4 / math.pi + (2 / math.pi) * (GAMMA - math.log(1448 / 721)) + (1 - 2 / math.pi) * 721 / 723277 atleast(kappa, 0.015 + 2.544, "the assembled constant")278 same(floorbase(Fraction(1417, 1850)), 1486, "Q at the Zhang rung")279 grid = [Fraction(750 + 5 * i, 1000) for i in range(46)]280 for b in grid:281 q = floorbase(b)282 atleast(q, 723, "Q(b) below the base floor at b = %s" % b)283 atmost(float(dgap(q, b)), -1.56, "gap at Q(b) for b = %s" % b)284 if q >= 3690:285 atmost(float(major(q, float(b))) - float(dgap(q, b)), 0.004, "majorant slack at b = %s" % b)286 atmost(major(3690, float(b)), -0.95, "majorant at 3690 for b = %s" % b)287 assert b >= Fraction(1417, 1850), "a low b with a high floor at b = %s" % b288 else:289 atmost(q, 3689, "empty range expected at b = %s" % b)290291# THE ELEMENTARY INEQUALITIES292293def check_elementary():294 for i in range(20001):295 v = i / 20000296 atmost(math.sin(math.pi * v), 4 * v * (1 - v) + 1e-15, "sine against the parabola at v = %.6f" % v)297 for i in range(1, 20001):298 x = math.pi / 2 * i / 20000299 atmost(1 / math.sin(x), 1 / x + 1 - 2 / math.pi + 1e-12, "cosecant bound at x = %.6f" % x)300301def gridsum(q, t, kernel):302 total = 0.0303 for r in range(q):304 u = (t + r) / q305 d = abs(u - round(u))306 if d < 1e-15:307 total += q308 else:309 total += abs(math.sin(math.pi * q * d)) / math.sin(math.pi * d) if kernel else 0.0310 return total311312def check_kernel():313 for q in [3, 5, 7, 10, 37, 100, 257, 1000]:314 worst = max(gridsum(q, i / 2000.0, True) for i in range(2001))315 atmost(worst, phi(q), "the kernel bound at q = %d" % q)316 atmost(worst / phi(q), 0.92, "slack of the kernel bound at q = %d" % q)317318def symbolsum(q, digits, t, power):319 total = 0.0320 for r in range(q):321 u = (t + r) / q322 value = abs(sum(cmath.exp(2j * math.pi * d * u) for d in digits))323 total += value ** power324 return total325326def check_floor_identity():327 for q in range(3, 10):328 for mask in range(1, 1 << q):329 digits = [d for d in range(q) if mask >> d & 1]330 k = len(digits)331 for i in range(20):332 t = i / 20.0333 near(symbolsum(q, digits, t, 2), q * k, 1e-8 * q * k, "the l2 identity at q = %d" % q)334 atleast(symbolsum(q, digits, t, 1), q - 1e-9, "the l1 floor at q = %d" % q)335336def check_shape():337 q = 10 ** 12338 near(cost(q), (math.log(math.log(q)) + math.log(2 / math.pi)) / math.log(q), 0.01, "the asymptotic shape")339 atleast(saving(10 ** 9), 0.11695, "the saving at a large base")340341# DOOR342343def main():344 check_table()345 check_wall()346 check_multi()347 check_exponents()348 check_ladder()349 check_floor()350 check_elementary()351 check_kernel()352 check_floor_identity()353 check_shape()354 print("sparse-mertens-under-grh: every printed number recomputed and asserted")355356if __name__ == "__main__":357 main()