dissection.py
17.2 kB · python · 453 lines
1import sys2import time34GAMMA = 0.577215664901532956# CONSTANTS78def harm(n):9 from math import log10 return log(n) + GAMMA + 1 / (2 * n)1112def phi_q(q):13 from math import pi14 n = -(-(q - 2) // 2)15 return (4 / pi) * q + (2 * q / pi) * harm(n) + (1 - 2 / pi) * (q - 2) + 0.7271617def pb_step3(q, m):18 return m ** 0.5 + phi_q(q) / q1920def psi2(q):21 from math import pi22 p = q // 223 h = harm(p - 1)24 if q % 2 == 0:25 return (q / pi) * (2 * h - 1 + 1 / p) + (1 - 2 / pi) * q / 226 return (q / pi) * (2 * h - 1 + 2 / p) + (1 - 2 / pi) * (q / 2 + 1 / (2 * q))2728def pb_chord(q, e0):29 from math import cos, pi30 sec = 1 / cos(pi * (e0 - (q - 1) / 2) / q)31 return ((4 / pi) * q + psi2(q) + q / 2 - sec / 2) / q3233def gap_hp(q, e0, b, chord):34 import mpmath35 mpmath.mp.dps = 4036 q = mpmath.mpf(q)37 pi = mpmath.pi38 g = mpmath.euler39 if chord:40 p = mpmath.floor(q / 2)41 h = mpmath.log(p - 1) + g + 1 / (2 * (p - 1))42 if int(q) % 2 == 0:43 ps = (q / pi) * (2 * h - 1 + 1 / p) + (1 - 2 / pi) * q / 244 else:45 ps = (q / pi) * (2 * h - 1 + 2 / p) + (1 - 2 / pi) * (q / 2 + 1 / (2 * q))46 sec = 1 / mpmath.cos(pi * (e0 - (q - 1) / 2) / q)47 pb = ((4 / pi) * q + ps + q / 2 - sec / 2) / q48 else:49 n = mpmath.ceil((q - 2) / 2)50 ph = (4 / pi) * q + (2 * q / pi) * (mpmath.log(n) + g + 1 / (2 * n)) + (1 - 2 / pi) * (q - 2) + mpmath.mpf("0.727")51 pb = 1 + ph / q52 return (q - 1) * q ** (-mpmath.mpf(b)) - pb5354def scan(chord, worst, lo, hi, b=0.8):55 import numpy56 from math import pi57 q = numpy.arange(lo, hi + 1, dtype=numpy.float64)58 if chord:59 p = numpy.floor(q / 2)60 h = numpy.log(p - 1) + GAMMA + 1 / (2 * (p - 1))61 ex = numpy.where(q % 2 == 0, 1 / p, 2 / p)62 fl = numpy.where(q % 2 == 0, q / 2, q / 2 + 1 / (2 * q))63 ps = (q / pi) * (2 * h - 1 + ex) + (1 - 2 / pi) * fl64 c = numpy.where(q % 2 == 1, 0.0, 0.5) if worst else (q - 1) / 265 pb = ((4 / pi) * q + ps + q / 2 - 1 / numpy.cos(pi * c / q) / 2) / q66 else:67 n = numpy.ceil((q - 2) / 2)68 pb = 1 + ((4 / pi) * q + (2 * q / pi) * (numpy.log(n) + GAMMA + 1 / (2 * n)) + (1 - 2 / pi) * (q - 2) + 0.727) / q69 gap = (q - 1) * q ** (-b) - pb70 ok = gap > 071 i = numpy.nonzero(ok)[0]72 first = int(q[i[0]])73 return first, int(ok.sum()), bool(ok[i[0]:].all()), float(gap[i[0]:].min()), int(q[i[0] + int(numpy.argmin(gap[i[0]:]))])7475def wall():76 print("THE REGION A WALL: the least base with PB < (base-1) base^(-4/5) at one excluded digit")77 print("form every e_0 e_0 in {0, base-1} scan held up-set least float gap above the wall")78 hi = 20000079 for name, chord, lo in [("step 3, PB_base(1)", False, 17), ("chord, PB'_base(1, e_0)", True, 36)]:80 rows = [scan(chord, w, lo, hi) for w in (True, False)]81 print("%-28s %10d %20d %d..%d %d %d %s %s %.3e at %d, %.3e at %d"82 % (name, rows[0][0], rows[1][0], lo, hi, rows[0][1], rows[1][1], rows[0][2], rows[1][2],83 rows[0][3], rows[0][4], rows[1][3], rows[1][4]))84 print("the worst e_0 is the middle digit at odd base and c = 1/2 at even base, the best e_0 in {0, base-1}")85 print()86 print("THE MARGIN AT EACH WALL, 40 digits, gap = (base-1) base^(-4/5) - PB")87 for name, chord, q, worst in [("step 3", False, 92317, False), ("chord, worst e_0", True, 39363, True),88 ("chord, e_0 = 0", True, 28352, False)]:89 g0 = gap_hp(q - 1, (q - 1) // 2 if worst else 0, 0.8, chord)90 g1 = gap_hp(q, q // 2 if worst else 0, 0.8, chord)91 print("%-18s base %6d gap %+.6e base %6d gap %+.6e" % (name, q - 1, float(g0), q, float(g1)))92 print()93 print("REGION A ASKS alpha_1 < 1/5, at each wall: 1/5 - alpha_1 = alpha - c - 4/5 = log(1 + gap/PB)/log base > 0, truncated down")94 import mpmath95 for name, chord, q, worst in [("step 3", False, 92317, False), ("chord, worst e_0", True, 39363, True),96 ("chord, e_0 = 0", True, 28352, False)]:97 e0 = q // 2 if worst else 098 g = gap_hp(q, e0, 0.8, chord)99 pb = (q - 1) * mpmath.mpf(q) ** (-mpmath.mpf(0.8)) - g100 v = mpmath.log(1 + g / pb) / mpmath.log(q)101 e = int(mpmath.floor(mpmath.log10(v)))102 mant = mpmath.floor(v / mpmath.mpf(10) ** e * 10 ** 4) / 10 ** 4103 print("%-18s base %6d 1/5 - alpha_1 >= %.4fe%03d" % (name, q, float(mant), e))104105 print()106 print("THE UNSHIFTED MASS region A pays: c_k = sum_(a mod base^k) |hat F_k(a/base^k)|, readings of c_k/c_(k-1) against 2(base - 1)")107 import numpy108 for q, e0, k in [(33, 16, 4), (33, 0, 4), (17, 8, 5)]:109 F = [v for v in range(q) if v != e0]110 cs = []111 for kk in (k - 1, k):112 ind = numpy.zeros(q ** kk)113 ind[strings(q, F, kk)] = 1.0114 cs.append(numpy.abs(numpy.fft.fft(ind)).sum())115 print(" base %d, excluded %2d, k = %d: c_k/c_(k-1) = %.4f against 2(base - 1) = %d" % (q, e0, k, cs[1] / cs[0], 2 * (q - 1)))116117# GRID118119def mu_upto(N):120 import numpy121 mu = numpy.ones(N + 1, dtype=numpy.int8)122 mu[0] = 0123 sieve = numpy.ones(N + 1, dtype=bool)124 sieve[:2] = False125 for p in range(2, int(N ** 0.5) + 1):126 if sieve[p]:127 sieve[p * p:: p] = False128 primes = numpy.nonzero(sieve)[0]129 for p in primes:130 mu[p:: p] = -mu[p:: p]131 for p in primes:132 if p * p > N:133 break134 mu[p * p:: p * p] = 0135 return mu136137def strings(q, F, L):138 import numpy139 out = numpy.zeros(1, dtype=numpy.int64)140 d = numpy.array(F, dtype=numpy.int64)141 for i in range(L):142 out = (out[:, None] + d[None, :] * (q ** i)).ravel()143 out.sort()144 return out145146def convergent(num, den, cap):147 import numpy148 n = num.astype(numpy.int64).copy()149 d = numpy.full_like(n, den) if numpy.isscalar(den) else den.astype(numpy.int64).copy()150 cap = numpy.broadcast_to(numpy.asarray(cap, dtype=numpy.int64), n.shape)151 h1 = numpy.ones_like(n)152 h2 = numpy.zeros_like(n)153 k1 = numpy.zeros_like(n)154 k2 = numpy.ones_like(n)155 live = d > 0156 c = numpy.zeros_like(n)157 c[live] = n[live] // d[live]158 h1, h2 = c * h1 + h2, h1159 k1, k2 = c * k1 + k2, k1160 n, d = d, n - c * d161 live = d > 0162 while live.any():163 c = numpy.zeros_like(n)164 c[live] = n[live] // d[live]165 hn = c * h1 + h2166 kn = c * k1 + k2167 go = live & (kn <= cap)168 h2 = numpy.where(go, h1, h2)169 k2 = numpy.where(go, k1, k2)170 h1 = numpy.where(go, hn, h1)171 k1 = numpy.where(go, kn, k1)172 nn = numpy.where(go, d, n)173 dd = numpy.where(go, n - c * d, 0)174 n, d = nn, dd175 live = go & (d > 0)176 return h1, k1177178def smooth(d, q):179 from math import gcd180 if d < 1:181 return False182 while True:183 g = gcd(d, q)184 if g == 1:185 return d == 1186 while d % g == 0:187 d //= g188189GAL_CAP = 10 ** 5190191def gmod(F, th):192 import numpy193 z = numpy.zeros(th.shape, dtype=complex)194 for v in F:195 z += numpy.exp(2j * numpy.pi * v * th)196 return numpy.abs(z)197198def l1_norms(q, F, i1):199 import numpy200 V = q ** i1201 t = numpy.arange(4 * V) / (4 * V)202 gs = [sum(numpy.exp(2j * numpy.pi * v * (q ** j * t % 1.0)) for v in F) for j in range(i1)]203 dg = [sum(2j * numpy.pi * v * numpy.exp(2j * numpy.pi * v * (q ** j * t % 1.0)) for v in F) for j in range(i1)]204 f = numpy.ones(t.size, dtype=complex)205 for g in gs:206 f = f * g207 fp = numpy.zeros(t.size, dtype=complex)208 for j in range(i1):209 term = q ** j * dg[j]210 for jj in range(i1):211 if jj != j:212 term = term * gs[jj]213 fp = fp + term214 return float(numpy.abs(f).mean()), float(numpy.abs(fp).mean())215216def regions_one(q, e0, k, Z):217 import numpy218 from math import log, pi, gcd219 F = [v for v in range(q) if v != e0]220 fill = len(F)221 y = q ** k222 D = numpy.zeros(y)223 D[strings(q, F, k)] = 1.0224 mu = mu_upto(y)[:y].astype(numpy.float64)225 hatF = numpy.conj(numpy.fft.fft(D))226 Sneg = numpy.fft.fft(mu)227 exact = float(mu[D > 0].sum())228 a = numpy.arange(y, dtype=numpy.int64)229 Q = int(y ** 0.6)230 l, d = convergent(a, y, Q)231 h = numpy.abs(a * d - l * y)232 assert (d <= Q).all() and (h * Q <= y).all()233 assert (numpy.gcd(l, d) == 1).all()234 y25 = y ** 0.4235 sm = numpy.array([smooth(int(v), q) for v in range(Q + 1)])236 A = d >= y25237 C = (~A) & (d < Z) & (h < Z)238 B = (~A) & (~C)239 C2 = C & sm[d]240 C1 = C & (~sm[d])241 term = hatF * Sneg / y242 parts = {nm: term[msk].sum().real for nm, msk in [("A", A), ("B", B), ("C1", C1), ("C2", C2)]}243 tot = sum(parts.values())244 print("base %d, excluded %d, level %d, y = %d, Q = y^(3/5) = %d, Z = %d" % (q, e0, k, y, Q, Z))245 print(" exact sum of mu over the strings %d, the four regions sum to %.6f, difference %.2e"246 % (exact, tot, abs(tot - exact)))247 assert abs(tot - exact) < 1e-6 * y248 for nm in ("A", "B", "C1", "C2"):249 msk = {"A": A, "B": B, "C1": C1, "C2": C2}[nm]250 print(" region %-2s points %8d contribution %+14.4f = %+.6f fill^k" % (nm, msk.sum(), parts[nm], parts[nm] / fill ** k))251 assert C1.sum() <= 3 * Z * Z252 nZ = sum(1 for v in range(1, Z) if sm[v])253 assert C2.sum() <= 3 * Z * nZ254 print(" counts: C1 %d <= 3 Z^2 = %d, C2 %d <= 3 Z N_Z = %d" % (C1.sum(), 3 * Z * Z, C2.sum(), 3 * Z * nZ))255 dd = d[C2]256 assert (y % dd == 0).all()257 j = a[C2] - l[C2] * (y // dd)258 j = (j + y // 2) % y - y // 2259 assert (numpy.abs(j) * dd == h[C2]).all() and (numpy.abs(j) * dd < Z).all()260 print(" C2 algebra: every d divides y and a = l y/d + j with |j| d = h < Z at all %d points" % C2.sum())261 ib = numpy.nonzero(B & (h >= 1))[0]262 capq = (2 * y) // h[ib]263 l2, d2 = convergent(a[ib], y, capq)264 lo = d2 * 2 * h[ib] >= y265 hi = d2 * h[ib] <= 2 * y266 diff = (l2 * d[ib] != l[ib] * d2)267 assert lo.all() and hi.all() and diff.all()268 print(" second approximation: y/(2h) <= d' <= 2y/h and l'/d' != l/d at all %d points of B with h >= 1" % ib.size)269 B3 = q * pb_step3(q, 1)270 a1 = log(B3 / fill) / log(q)271 af = numpy.abs(hatF)272 Dc = numpy.floor(numpy.log2(d)).astype(numpy.int64)273 Hc = numpy.where(h > 0, numpy.floor(numpy.log2(numpy.maximum(h, 1))).astype(numpy.int64) + 1, 0)274 key = Dc * 64 + Hc275 order = numpy.argsort(key, kind="stable")276 ks, starts = numpy.unique(key[order], return_index=True)277 sums = numpy.add.reduceat(af[order], starts)278 worst, tested, skipped = 0.0, 0, 0279 gal, top_ok, norms, galn = 0.0, True, {}, 0280 for kk, s, st in zip(ks, sums, starts):281 Dv = 2 ** int(kk // 64)282 hc = int(kk % 64)283 Hv = 0 if hc == 0 else 2 ** (hc - 1)284 V1 = 1285 while V1 < 4 * Dv * Dv:286 V1 *= q287 V2 = 1288 while V2 < 4 * Hv / Dv + 1:289 V2 *= q290 if V1 * V2 > y or 16 * Dv * Hv > y:291 skipped += 1292 continue293 tested += 1294 if Dv == 1 and Hv == 0:295 continue296 bound = (1 + pi * q) * fill ** k * (V1 * V2) ** a1297 worst = max(worst, s / bound)298 mem = order[st: st + (len(order) if st == starts[-1] else starts[numpy.searchsorted(starts, st) + 1] - st)]299 dm = d[mem]300 fk = (l[mem] % dm) * (Q + 1) + dm301 i1 = round(log(V1) / log(q))302 bot = numpy.ones(mem.size)303 for j in range(i1):304 bot = bot * gmod(F, (q ** j * a[mem] % y) / y)305 o2 = numpy.argsort(fk, kind="stable")306 _, s2 = numpy.unique(fk[o2], return_index=True)307 lhs = numpy.maximum.reduceat(bot[o2], s2).sum()308 if V1 <= GAL_CAP:309 if i1 not in norms:310 norms[i1] = l1_norms(q, F, i1)311 rhs = 4 * Dv * Dv * norms[i1][0] + norms[i1][1]312 gal = max(gal, lhs / rhs)313 galn += 1314 grp = numpy.split(a[mem][o2] % V2, s2[1:])315 top_ok = top_ok and all(numpy.unique(g).size == g.size for g in grp)316 assert worst <= 1.0 and gal <= 1.0 and top_ok317 print(" hybrid lemma, alpha_1 <= %.6f from the step 3 constant: %d classes (D, H) meet V_1 V_2 <= y and 16 D H <= y, %d do not"318 % (a1 + 5e-7, tested, skipped))319 print(" past the class of a = 0, where |hat F_k(0)| = fill^k: largest ratio of a class l^1 sum to the lemma %.6f;"320 % (worst + 5e-7))321 print(" the Gallagher step at the %d classes with V_1 <= %d: the fractions' largest |hat F_(i_1)| at their residues, summed,"322 % (galn, GAL_CAP))323 print(" against 4 D^2 ||f||_1 + ||f'||_1 with both norms read on 4 V_1 points, largest ratio %.6f; every fraction's residues"324 " distinct mod V_2" % (gal + 5e-7))325 L = log(y)326 rb = numpy.abs(Sneg[ib]) / (y ** 0.8 + y * L ** 3 * numpy.maximum(d[ib], h[ib]) ** -0.5)327 ra = numpy.abs(Sneg[A]) / y ** 0.8328 print(" readings, not bounds: max |S|/(y^(4/5) + y (log y)^3 max(d,h)^(-1/2)) on B %.3e, max |S|/y^(4/5) on A %.4f"329 % (rb.max(), ra.max()))330 print()331332def lemma_a_prime(q, e0, k, Z, draws, seed):333 import random334 from math import cos, sin, pi, log, floor, gcd335 F = [v for v in range(q) if v != e0]336 fill = len(F)337 y = q ** k338 rng = random.Random(seed)339 cp = 1 - (2 / fill) * (1 - cos(pi / (4 * q)))340 worst, n = -1e9, 0341 while n < draws:342 d = rng.randrange(2, Z)343 if smooth(d, q):344 continue345 l = rng.randrange(d)346 if gcd(l, d) != 1:347 continue348 j = rng.randrange(-Z, Z)349 a = (l * y) // d + j350 h = abs(a * d - l * y)351 if h >= Z:352 continue353 assert (h / (d * y)) < y ** (-2 / 3) / (4 * q * (q - 1))354 lf = 0.0355 for i in range(k):356 ph = (pow(q, i, y) * a % y) / y357 s = abs(sum(complex(cos(2 * pi * v * ph), sin(2 * pi * v * ph)) for v in F))358 lf += log(s / fill)359 md = max(1, floor(log(d / 2) / log(q) + 1e-12) + 1)360 rhs = floor(2 * k / (3 * md)) * log(cp)361 worst = max(worst, lf - rhs)362 n += 1363 assert worst <= 1e-9364 print("perturbed Lemma A' at the grid point, base %d, excluded %d, level %d, d < %d with a prime outside the base, h < %d"365 % (q, e0, k, Z, Z))366 print(" %d seeded draws: log |hat F_k(a/y)|/fill^k minus floor(2k/(3 m_d)) log c' is at most %.4f, never above 0"367 % (draws, worst))368369def regions():370 t = time.time()371 print("FALSIFICATION OF THE DISSECTION at small base and level, every grid point a mod y")372 print("A: d >= y^(2/5); B: d < y^(2/5) and max(d, h) >= Z; C1, C2: d < Z and h < Z, d with or without a prime outside the base")373 print()374 regions_one(10, 5, 6, 16)375 regions_one(10, 0, 6, 16)376 regions_one(5, 2, 9, 12)377 lemma_a_prime(10, 5, 30, 60, 3000, 1009)378 lemma_a_prime(5, 0, 45, 40, 3000, 1013)379 print("runtime %.1f s" % (time.time() - t))380381# BLOCKS382383def blocks_of(x, q, F):384 ds = []385 v = x386 while v:387 ds.append(v % q)388 v //= q389 L = len(ds)390 out = []391 if 0 in F:392 if L >= 2:393 out.append((0, L - 1))394 else:395 for ell in range(1, L):396 out.append((0, ell))397 for j in range(L - 1, -1, -1):398 top = 0399 for i in range(L - 1, j, -1):400 top = top * q + ds[i]401 for f in F:402 if f < ds[j] and not (j == L - 1 and f == 0):403 out.append((top * q + f, j))404 if ds[j] not in F:405 break406 tail = all(v in F for v in ds)407 return out, tail, L408409def blocks():410 import numpy411 import random412 t = time.time()413 print("OFF THE POWERS OF THE BASE: S_F below x as blocks P base^k + D_k, at most fill + 1 per scale")414 N = 2 * 10 ** 6415 mu = mu_upto(N)416 rng = random.Random(7)417 for q, e0 in [(10, 0), (10, 9), (5, 0), (7, 3), (4, 1)]:418 F = [v for v in range(q) if v != e0]419 fill = len(F)420 memb = numpy.zeros(N + 1, dtype=bool)421 L = 1422 while q ** L <= N:423 L += 1424 for ell in range(1, L + 1):425 s = strings(q, F, ell)426 s = s[s >= q ** (ell - 1)] if 0 in F else s427 memb[s[s <= N]] = True428 memb[0] = False429 cum = numpy.cumsum(numpy.where(memb, mu, 0).astype(numpy.int64))430 cnt = numpy.cumsum(memb.astype(numpy.int64))431 worst = 0432 for x in [rng.randrange(1, N) for _ in range(400)] + [q ** e - 1 for e in range(1, L) if q ** e <= N]:433 bl, tail, Lx = blocks_of(x, q, F)434 per = {}435 tot = 0436 for P, kk in bl:437 per[kk] = per.get(kk, 0) + 1438 s = strings(q, F, kk) + P * q ** kk439 tot += int(mu[s[s > 0]].sum())440 if tail:441 tot += int(mu[x])442 assert tot == cum[x], (q, e0, x)443 assert max(per.values(), default=0) <= fill + 1444 assert cnt[x] >= fill ** (Lx - 1) - 1445 worst = max(worst, max(per.values(), default=0))446 print(" base %2d excluded %d: 400 random x below %d and every base^e - 1, the block sums meet M_F(x) exactly;"447 % (q, e0, N))448 print(" at most %d blocks at one scale against fill + 1 = %d, and A_F(x) >= fill^(L-1) - 1 at every x" % (worst, fill + 1))449 print("runtime %.1f s" % (time.time() - t))450451if __name__ == "__main__":452 verb = sys.argv[1] if len(sys.argv) > 1 else "wall"453 {"wall": wall, "regions": regions, "blocks": blocks}[verb]()