"""x211scan.py -- fast float scan of X_211(D) and X_1111(D). Means <.> are D-independent: tabulate once per admissible tuple, then sweep D through the coefficients c_q(D) = mu(q) * m_q(floor(D/q)). Support prefilters (necessary conditions, proved in the paper): X_211 term (q; r 1 else 60 # ---------- mu, squarefree, prime sets ---------- mu = np.ones(DMAX + 1, dtype=np.int64) primes = [] is_comp = np.zeros(DMAX + 1, dtype=bool) spf = np.zeros(DMAX + 1, dtype=np.int64) for i in range(2, DMAX + 1): if not is_comp[i]: primes.append(i); spf[i] = i for j in range(2 * i, DMAX + 1, i): is_comp[j] = True if spf[j] == 0: spf[j] = i def factor(n): s = [] while n > 1: p = spf[n]; s.append(p); e = 0 while n % p == 0: n //= p; e += 1 if e > 1: return None # not squarefree return s PS = {} # prime set per squarefree q>1 MU = {1: 1} for q in range(2, DMAX + 1): f = factor(q) if f is None: MU[q] = 0; continue PS[q] = frozenset(f); MU[q] = (-1) ** len(f) SQF = sorted(PS.keys()) # squarefree 2..DMAX # ---------- sigma_q step values on (j,j+1), j mod q ---------- def sigvals(q): out = np.zeros(q) divs = [t for t in range(1, q + 1) if q % t == 0 and MU.get(t, 0) != 0] for j in range(q): x = j + 0.5; v = 0.0 for t in divs: y = t * x / q v += MU[t] / t * ((y - np.floor(y)) - 0.5) out[j] = v return out S = {q: sigvals(q) for q in SQF} # ---------- coefficients c_q(D) = mu(q)*m_q(floor(D/q)) ---------- # m_q(Y) = sum_{m<=Y,(m,q)=1} mu(m)/m ; tabulate partial sums per q Carr = np.zeros((DMAX + 1, DMAX + 1)) # Carr[q, D] for q in SQF: Ymax = DMAX // q msum = 0.0; mq_of_Y = np.zeros(Ymax + 1) for m in range(1, Ymax + 1): if MU.get(m, 0) != 0 and gcd(m, q) == 1: msum += MU[m] / m mq_of_Y[m] = msum for D in range(q, DMAX + 1): Carr[q, D] = MU[q] * mq_of_Y[D // q] # ---------- generic mean of products of step functions ---------- def mean_prod(spec): """spec = list of (q, power); exact over L = lcm.""" L = lcm(*[q for q, e in spec]) n = np.arange(L) p = np.ones(L) for q, e in spec: p *= S[q][n % q] ** e return p.mean(), L # ---------- enumerate X211 terms ---------- t0 = time.time() terms211 = [] # (q, r, s, mean) Lmax211 = 0 for ir in range(len(SQF)): r = SQF[ir] for js in range(ir + 1, len(SQF)): s = SQF[js] g = gcd(r, s); r1 = r // g; s1 = s // g core = r1 * s1 if core > DMAX: continue for q in range(core, DMAX + 1, core): if MU.get(q, 0) == 0 or q == r or q == s: continue m, L = mean_prod([(q, 2), (r, 1), (s, 1)]) Lmax211 = max(Lmax211, L) if m != 0.0: terms211.append((q, r, s, m)) print(f"X211: {len(terms211)} nonzero terms, Lmax={Lmax211}, {time.time()-t0:.1f}s", flush=True) # ---------- enumerate X1111 terms ---------- t0 = time.time() terms1111 = [] Lmax1111 = 0 nsq = len(SQF) for iq in range(nsq): q = SQF[iq]; Pq = PS[q] for ir in range(iq + 1, nsq): r = SQF[ir]; Pr = PS[r] for js in range(ir + 1, nsq): s = SQF[js]; Psx = PS[s] # primes appearing exactly once in {q,r,s} once = (Pq ^ Pr ^ Psx) - (Pq & Pr & Psx) U = 1 for p in once: U *= p if U > DMAX: continue allowed = Pq | Pr | Psx t = U if U > s else ((s // U) + 1) * U while t <= DMAX: if MU.get(t, 0) != 0 and t != q and t != r and t != s and PS[t] <= allowed: m, L = mean_prod([(q, 1), (r, 1), (s, 1), (t, 1)]) Lmax1111 = max(Lmax1111, L) if m != 0.0: terms1111.append((q, r, s, t, m)) t += U print(f"X1111: {len(terms1111)} nonzero terms, Lmax={Lmax1111}, {time.time()-t0:.1f}s", flush=True) # ---------- enumerate X31 terms (r | q by lem:supp(i)) ---------- t0 = time.time() terms31 = [] for q in SQF: for r in SQF: if r != q and q % r == 0: m, L = mean_prod([(q, 3), (r, 1)]) if m != 0.0: terms31.append((q, r, m)) print(f"X31: {len(terms31)} nonzero terms, {time.time()-t0:.1f}s", flush=True) # ---------- enumerate X22 terms (odd common prime by lem:supp(iv)) ---------- t0 = time.time() terms22 = [] V = {q: float((S[q] ** 2).mean()) for q in SQF} for iq in range(len(SQF)): q = SQF[iq] for jr in range(iq + 1, len(SQF)): r = SQF[jr] g = gcd(q, r) if g == 1 or (g == 2): continue m, L = mean_prod([(q, 2), (r, 2)]) cov = m - V[q] * V[r] if cov != 0.0: terms22.append((q, r, cov)) print(f"X22: {len(terms22)} nonzero terms, {time.time()-t0:.1f}s", flush=True) # ---------- sweep D ---------- def sweep(): X211 = np.zeros(DMAX + 1); X1111 = np.zeros(DMAX + 1) X31 = np.zeros(DMAX + 1); X22 = np.zeros(DMAX + 1) A211 = np.zeros(DMAX + 1); A1111 = np.zeros(DMAX + 1) # absolute-value sums if terms31: A3 = np.array([(q, r) for q, r, _ in terms31], dtype=np.int64) M3 = np.array([m for *_, m in terms31]) mx3 = A3.max(axis=1) for D in range(2, DMAX + 1): sel = mx3 <= D if sel.any(): a = A3[sel]; c = Carr[:, D] X31[D] = 4.0 * np.sum(c[a[:, 0]] ** 3 * c[a[:, 1]] * M3[sel]) if terms22: A2 = np.array([(q, r) for q, r, _ in terms22], dtype=np.int64) M2 = np.array([m for *_, m in terms22]) mx2 = A2.max(axis=1) for D in range(2, DMAX + 1): sel = mx2 <= D if sel.any(): a = A2[sel]; c = Carr[:, D] X22[D] = 6.0 * np.sum(c[a[:, 0]] ** 2 * c[a[:, 1]] ** 2 * M2[sel]) if terms211: A = np.array([(q, r, s) for q, r, s, _ in terms211], dtype=np.int64) M = np.array([m for *_, m in terms211]) mx = A.max(axis=1) for D in range(2, DMAX + 1): sel = mx <= D if not sel.any(): continue a = A[sel]; mm = M[sel] c = Carr[:, D] contrib = c[a[:, 0]] ** 2 * c[a[:, 1]] * c[a[:, 2]] * mm X211[D] = 12.0 * contrib.sum() A211[D] = 12.0 * np.abs(contrib).sum() if terms1111: B = np.array([(q, r, s, t) for q, r, s, t, _ in terms1111], dtype=np.int64) N = np.array([m for *_, m in terms1111]) mx = B.max(axis=1) for D in range(2, DMAX + 1): sel = mx <= D if not sel.any(): continue b = B[sel]; nn = N[sel] c = Carr[:, D] contrib = c[b[:, 0]] * c[b[:, 1]] * c[b[:, 2]] * c[b[:, 3]] * nn X1111[D] = 24.0 * contrib.sum() A1111[D] = 24.0 * np.abs(contrib).sum() return X31, X22, X211, X1111, A211, A1111 X31, X22, X211, X1111, A211, A1111 = sweep() out = {"DMAX": DMAX, "X31": X31.tolist(), "X22": X22.tolist(), "X211": X211.tolist(), "X1111": X1111.tolist(), "A211": A211.tolist(), "A1111": A1111.tolist(), "n31": len(terms31), "n22": len(terms22), "n211": len(terms211), "n1111": len(terms1111)} with open(f"x_sectors_{DMAX}.json", "w") as f: json.dump(out, f) np.savez(f"x_terms_{DMAX}.npz", t211=np.array([(q, r, s) for q, r, s, _ in terms211], dtype=np.int64) if terms211 else np.zeros((0,3),dtype=np.int64), m211=np.array([m for *_, m in terms211]), t1111=np.array([(q, r, s, t) for q, r, s, t, _ in terms1111], dtype=np.int64) if terms1111 else np.zeros((0,4),dtype=np.int64), m1111=np.array([m for *_, m in terms1111]), Carr=Carr) print("\n D X31 X22 X211 X1111 X211+X1111 |211| |1111|") for D in [6, 10, 14, 15, 21, 22] + list(range(30, DMAX + 1, max(1, DMAX // 15))): if D <= DMAX: print(f"{D:4d} {X31[D]:+.6f} {X22[D]:+.6f} {X211[D]:+.6f} {X1111[D]:+.6f} {X211[D]+X1111[D]:+.6f} {A211[D]:10.4f} {A1111[D]:10.4f}")