#!/usr/bin/env python3 """ test_346 -- HOMOGENEOUS THREE-METRICS AS HYPERSURFACES OF R^4: THE THEOREM, THE LEFT-INVARIANT CLASSIFICATION, AND THE THURSTON GEOMETRIES AS CHECKS For MB ยง4 (the companion of MA). Conventions as in MA: K_ab = -N . d_a d_b x, so the unit sphere with outward normal has K = g; M = -G = adj K; Cod(K)_abc = nabla_a K_bc - nabla_b K_ac; [nabla_a, nabla_b] V^c = R^c_dab V^d, R_ab = R^c_acb. Left-invariant metrics are handled in an orthonormal frame with [e_i, e_j] = c_ij^k e_k, the structure constants in the Ellis-MacCallum form c_ij^k = eps_ijl n^lk + a_i delta_j^k - a_j delta_i^k (n symmetric; Jacobi iff n a = 0), and the Levi-Civita connection from Koszul's formula, = (c_ijk - c_jki + c_kij)/2. This engine is written here and does not use Milnor's formula; P1 compares the two. PRE-REGISTERED: P1 (the engine; the Milnor-frame table) In a Milnor frame (n = diag(l1, l2, l3), a = 0) the engine's Ricci tensor is diagonal with entries 2 mu_j mu_k, mu_i = (l1 + l2 + l3)/2 - l_i, and M_i = mu_i l_i - mu_j mu_k (symbolic). For (1,1,1), (1,1,-1), (1,0,0), (1,-1,0): M = Id/4; diag(1/4, 1/4, -7/4) with det -7/64; eigenvalues (-3/4, 1/4, 1/4) with det -3/64; diag(-1, -1, 1) with det +1. Bianchi V (n = 0, a = (a,0,0)) has M = -a^2 Id, det -a^6. Jacobi holds for a random n with n a = 0 and fails for one with n a != 0. P2 (Sol^3 in coordinates, with positive controls) For ds^2 = e^{2z}dx^2 + e^{-2z}dy^2 + dz^2, a coordinate pipeline (Christoffels, Riemann, Ricci, M, K~ = sqrt(det M) M^{-1}, Cod) gives: the four Christoffel symbols G^x_xz = 1, G^y_yz = -1, G^z_xx = -e^{2z}, G^z_yy = e^{-2z} and no others; Ric = diag(0, 0, -2); M = diag(-1,-1,1); K~_ab = diag(-e^{2z}, -e^{-2z}, 1); the nonzero Codazzi components with a < b are exactly C_xzx = 2e^{2z} and C_yzy = -2e^{-2z}; the sum of squares over all ordered index triples is 16 cosh 4z; the invariant norm g^aa' g^bb' g^cc' C_abc C_a'b'c' is 16. The Lie engine on (1,-1,0) gives constant orthonormal components of modulus 2 and the same invariant norm 16. Positive controls for the coordinate pipeline: (a) on the unit S^3 in hyperspherical coordinates it gives Ric = 2 gamma, M = Id and a residual that vanishes after expand_trig and simplification, and is below 1e-30 at five random points; (b) on the realizable non-round metric (r^2 + r'^2) dchi^2 + r^2 sin^2 chi (dtheta^2 + sin^2 theta dphi^2), r = 1 + cos(2 chi)/10 (a convex radial graph), the Ricci tensor is diagonal to 1e-30, K~ has eigenvalues k_mer, k_sph, k_sph (test_345's formulas) and the residual is below 1e-30, at three points. P3 (Lemma 4.1, the Minkowski-Hsiung formulas) With s = x . N, H1 = sigma1/3, H2 = sigma2/3, on the radial graphs r = 1 + eps cos 2chi for eps = 1/10 (convex) and 3/10 (indefinite): int (1 - H1 s) dV and int (H1 - H2 s) dV vanish to 1e-25 (40 digits), while int |1 - H1 s| dV > 1e-3 (the integrands are not zero). Symbolic: 18 (H1^2 - H2) = sum_{i 0, H = -(s1+s2+s3)(-s1+s2+s3)(s1-s2+s3) (s1+s2-s3) with s_i = sqrt(l_i); on H = 0 with l_i > 0 (l3 = (s1 + s2)^2), det M < 0. M scalar with l_i > 0 forces l1 = l2 = l3. P6 (E(2) and E(1,1)) For l = (l1, l2, 0): R = -(l1 - l2)^2/2 and det M = -(l1 - l2)^4 (3l1 + l2)(l1 + 3l2)/64. So det M < 0 for l1, l2 > 0, l1 != l2 (non-flat E(2)); for l1 > 0 > l2 (E(1,1), Sol), det M > 0 iff -3 l1 < l2 < -l1/3, and rank M = 2 on the two boundary lines (checked at l = (1, -3, 0) and (3, -1, 0)). P7 (Corollary 4.6, the degenerate left-invariant metrics) Unimodular: any two of M_1, M_2, M_3 vanish together only where M = 0 (all common real zeros of each pair, over the three affine charts, have M = 0). Class B in the Schuecking frame (a = (a,0,0), n = diag(0, n2, n3), a != 0): det M = -PQST/64 with test_215's P, Q, S, T; rank M <= 1 iff (n2, n3) = (-a, a) or (a, -a), and there M has eigenvalues (0, 0, -4a^2); its unit eigenvector v for -4a^2 is parallel (nabla_{e_i} v = 0 for all i), so the metric is locally H^2(-4a^2) x R. P8 (class B: test_215's algebra from the engine) Every nonzero component of Cod(adj M) in the Schuecking frame is divisible by (n2 - n3); C_231 = -(n2 - n3)^2 (n2 + n3) P Q / 8; on n3 = -n2, C_122 = 4 a n2^2 (n2^2 - a^2) and det M = -(a^2 - n2^2)^3; at n2 = n3, det M = -a^6. P9 (Theorem 4.5 and Corollary 4.6, sampled exactly) Random left-invariant metrics in exact rational arithmetic, in randomly rotated orthonormal frames (rational Cayley rotations): 200 unimodular, 200 generic class B (n = P n0 P with P the projection off a), 100 class B near the Schuecking line n3 = -n2 with a^2 < n2^2. Jacobi holds for all; at least 20 unimodular and 20 class B samples have det M > 0; every sample with det M > 0 and M not scalar has Cod(adj M) != 0 exactly. Near-round control: for l = (1, 1 + d, 1 + 2d), |Cod(adj M)|^2 / d^2 at d = 10^-3 and 10^-6 agree to 1%, and Cod(adj M) = 0 at d = 0. P10 (Proposition 4.7, the product geometries) H^2 x R (metric (dx^2 + dy^2)/y^2 + dt^2): M has rank 1 with image span(d_t), and d_t is parallel. For symmetric K of rank 2 (generic, symbolic), K adj K = 0 and im(adj K) = ker K; rank adj K is 3, 1, 0 for rank K = 3, 2, 1 (exact random examples). The round cylinder S^2(r) x R: K = diag(1/r, 1/r, 0) in the product frame has adj K = M and Cod(K) = 0. FIRST RUN, AND THE CHANGES: The first run crashed before any check (sympy cannot multiply by a bool; int() added). The second passed every check up to P2's control (a) and was stopped after eleven minutes inside control (b): the symbolic chain for the radial graph was too slow to evaluate, and a third attempt with the profile carried as jet symbols was stopped after ten. Control (b) is re-posed on the Flamm metric of MA Example 3.3 instead, (1 - 2m/rho)^{-1} drho^2 + rho^2 (dtheta^2 + sin^2 theta dphi^2): realizable, not homogeneous, det M > 0 and M indefinite, as for Sol^3. Registered: M = diag(2m, -m, -m)/rho^3, K~ = -diag(k_mer, k_sph, k_sph) with MA's k_mer = -sqrt(2m)/(2 rho^{3/2}), k_sph = sqrt(2m)/rho^{3/2} (the sign is sign det K, MA Theorem A(3)), and Cod(K~) = 0 exactly. Nothing else changed. """ import sys import random from fractions import Fraction as Fr from itertools import product import sympy as sp import mpmath as mp FAILURES = [] def check(name, ok, detail=""): print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else "")) if not ok: FAILURES.append(name) R3 = range(3) # ---------------------------------------------------------------------------------------- the Lie-algebra engine def eps3(i, j, k): return ((i - j) * (j - k) * (k - i)) // 2 def structure(n, a): """c[i][j][k] = <[e_i, e_j], e_k> for c_ij^k = eps_ijl n^lk + a_i delta_jk - a_j delta_ik.""" return [[[sum(eps3(i, j, l) * n[l][k] for l in R3) + a[i] * int(j == k) - a[j] * int(i == k) for k in R3] for j in R3] for i in R3] def rotate(c, Q): return [[[sum(Q[i][p] * Q[j][q] * Q[k][s] * c[p][q][s] for p in R3 for q in R3 for s in R3) for k in R3] for j in R3] for i in R3] def jacobi_ok(c): return all(sum(c[i][j][l] * c[l][k][m] + c[j][k][l] * c[l][i][m] + c[k][i][l] * c[l][j][m] for l in R3) == 0 for i, j, k, m in product(R3, repeat=4)) def connection(c, half): return [[[half * (c[i][j][k] - c[j][k][i] + c[k][i][j]) for k in R3] for j in R3] for i in R3] def ricci(c, G): def Rvec(i, j, l): # R(e_i, e_j) e_l, components return [sum(G[j][l][m] * G[i][m][q] - G[i][l][m] * G[j][m][q] - c[i][j][m] * G[m][l][q] for m in R3) for q in R3] return [[sum(Rvec(cc, b, a)[cc] for cc in R3) for b in R3] for a in R3] def Mmat(Ric, half): R = sum(Ric[i][i] for i in R3) return [[half * R * int(i == j) - Ric[i][j] for j in R3] for i in R3] def adj3(A): def cof(i, j): r = [x for x in R3 if x != i] s = [y for y in R3 if y != j] return (-1) ** (i + j) * (A[r[0]][s[0]] * A[r[1]][s[1]] - A[r[0]][s[1]] * A[r[1]][s[0]]) return [[cof(j, i) for j in R3] for i in R3] def det3(A): return sum(A[0][j] * adj3(A)[j][0] for j in R3) def cod(G, K): def dK(i, j, k): return -sum(G[i][j][m] * K[m][k] + G[i][k][m] * K[j][m] for m in R3) return {(i, j, k): dK(i, j, k) - dK(j, i, k) for i in R3 for j in R3 for k in R3 if i < j} def lie(n, a, half=sp.Rational(1, 2), Q=None): c = structure(n, a) if Q is not None: c = rotate(c, Q) G = connection(c, half) Ric = ricci(c, G) M = Mmat(Ric, half) return c, G, Ric, M def milnor(l1, l2, l3, half=sp.Rational(1, 2)): z = 0 return lie([[l1, z, z], [z, l2, z], [z, z, l3]], [z, z, z], half) def symM(M): return sp.Matrix(3, 3, lambda i, j: sp.expand(M[i][j])) # ------------------------------------------------------------------------------------------------------------- P1 print("P1 the engine against Milnor's formula; the Milnor-frame table") l1, l2, l3 = sp.symbols("l1 l2 l3", real=True) c, G, Ric, M = milnor(l1, l2, l3) lam = [l1, l2, l3] mu = [(l1 + l2 + l3) / 2 - l for l in lam] ok = all(sp.expand(Ric[i][j]) == 0 for i in R3 for j in R3 if i != j) ok &= all(sp.expand(Ric[i][i] - 2 * mu[(i + 1) % 3] * mu[(i + 2) % 3]) == 0 for i in R3) ok &= all(sp.expand(M[i][i] - (mu[i] * lam[i] - mu[(i + 1) % 3] * mu[(i + 2) % 3])) == 0 for i in R3) check("Milnor frame: Ric = diag(2 mu_j mu_k), M_i = mu_i l_i - mu_j mu_k (symbolic)", ok) q = sp.Rational table = {(1, 1, 1): ([q(1, 4)] * 3, q(1, 64)), (1, 1, -1): ([q(1, 4), q(1, 4), q(-7, 4)], q(-7, 64)), (1, 0, 0): ([q(-3, 4), q(1, 4), q(1, 4)], q(-3, 64)), (1, -1, 0): ([-1, -1, 1], 1)} det_t = [] ok = True for lm, (eig, dt) in table.items(): Mm = symM(milnor(*[sp.Integer(v) for v in lm])[3]) ok &= Mm.is_diagonal() and sorted(Mm.diagonal()) == sorted(eig) and Mm.det() == dt det_t.append(f"{lm}: {list(Mm.diagonal())}, det {Mm.det()}") check("S^3 Id/4; SL2 det -7/64; Nil (-3/4, 1/4, 1/4) det -3/64; Sol diag(-1,-1,1) det +1", ok, "; ".join(det_t)) a_ = sp.symbols("a", positive=True) MV = symM(lie([[0] * 3 for _ in R3], [a_, 0, 0])[3]) check("Bianchi V (H^3): M = -a^2 Id, det M = -a^6", MV == -a_**2 * sp.eye(3) and sp.expand(MV.det() + a_**6) == 0, f"M = {list(MV.diagonal())}") random.seed(346) av = [Fr(1), Fr(2), Fr(-1)] P = [[Fr(int(i == j)) - av[i] * av[j] / sum(x * x for x in av) for j in R3] for i in R3] n0 = [[Fr(random.randint(-5, 5)) for _ in R3] for _ in R3] n0 = [[n0[i][j] + n0[j][i] for j in R3] for i in R3] nP = [[sum(P[i][k] * n0[k][l] * P[l][j] for k in R3 for l in R3) for j in R3] for i in R3] okJ = jacobi_ok(structure(nP, av)) and not jacobi_ok(structure(n0, av)) check("Jacobi holds iff n a = 0 (one random example of each)", okJ) # ------------------------------------------------------------------------------------------------------------- P2 print("P2 Sol^3 in coordinates; the Lie engine; positive controls") def coord_pipeline(g, X, D=None, post=sp.expand): """D(expr, a) is the derivative along coordinate a (default: the partial derivative).""" D = D or (lambda e, a: sp.diff(e, X[a])) gi = g.inv() Gam = [[[post(sum(gi[cc, d] * (D(g[d, b], a) + D(g[d, a], b) - D(g[a, b], d)) for d in R3) / 2) for b in R3] for a in R3] for cc in R3] # Gam[c][a][b] def Riem(cc, d, a, b): return (D(Gam[cc][b][d], a) - D(Gam[cc][a][d], b) + sum(Gam[cc][a][e] * Gam[e][b][d] - Gam[cc][b][e] * Gam[e][a][d] for e in R3)) Ric = sp.Matrix(3, 3, lambda d, b: sum(Riem(cc, d, cc, b) for cc in R3)) return Gam, Ric, gi def codazzi_coord(Gam, Kl, X, D=None): D = D or (lambda e, a: sp.diff(e, X[a])) def dK(a, b, cc): return D(Kl[b, cc], a) - sum(Gam[d][a][b] * Kl[d, cc] + Gam[d][a][cc] * Kl[b, d] for d in R3) return {(a, b, cc): dK(a, b, cc) - dK(b, a, cc) for a in R3 for b in R3 for cc in R3 if a != b} x, y, z = sp.symbols("x y z", real=True) X = [x, y, z] g = sp.diag(sp.exp(2 * z), sp.exp(-2 * z), 1) Gam, Ric, gi = coord_pipeline(g, X) want_G = {(0, 0, 2): 1, (0, 2, 0): 1, (1, 1, 2): -1, (1, 2, 1): -1, (2, 0, 0): -sp.exp(2 * z), (2, 1, 1): sp.exp(-2 * z)} okG = all(sp.simplify(Gam[cc][a][b] - want_G.get((cc, a, b), 0)) == 0 for cc in R3 for a in R3 for b in R3) Ric = Ric.applyfunc(sp.simplify) Rs = sp.simplify(sum(gi[i, j] * Ric[i, j] for i in R3 for j in R3)) Mmix = (sp.Rational(1, 2) * Rs * sp.eye(3) - gi * Ric).applyfunc(sp.simplify) Kmix = (sp.sqrt(Mmix.det()) * Mmix.inv()).applyfunc(sp.simplify) Kl = (g * Kmix).applyfunc(sp.simplify) C = {k: sp.simplify(v) for k, v in codazzi_coord(Gam, Kl, X).items()} nz = {k: v for k, v in C.items() if v != 0 and k[0] < k[1]} want_C = {(0, 2, 0): 2 * sp.exp(2 * z), (1, 2, 1): -2 * sp.exp(-2 * z)} sumsq = sp.simplify((sum(v**2 for v in C.values()) - 16 * sp.cosh(4 * z)).rewrite(sp.exp)) inv = sp.simplify(sum(gi[a, a] * gi[b, b] * gi[cc, cc] * v**2 for (a, b, cc), v in C.items())) check("Christoffels: exactly G^x_xz = 1, G^y_yz = -1, G^z_xx = -e^2z, G^z_yy = e^-2z", okG) check("Ric = diag(0,0,-2), M = diag(-1,-1,1), K~_ab = diag(-e^2z, -e^-2z, 1)", Ric == sp.diag(0, 0, -2) and Mmix == sp.diag(-1, -1, 1) and sp.simplify(Kl - sp.diag(-sp.exp(2 * z), -sp.exp(-2 * z), 1)) == sp.zeros(3, 3)) check("nonzero Codazzi (a < b): exactly C_xzx = 2e^2z, C_yzy = -2e^-2z; sum of squares 16 cosh 4z; invariant norm 16", set(nz) == set(want_C) and all(sp.simplify(nz[k] - want_C[k]) == 0 for k in want_C) and sumsq == 0 and inv == 16, f"components {nz}; invariant norm {inv}") c, G, Ric_, M_ = milnor(sp.Integer(1), sp.Integer(-1), sp.Integer(0)) Ks = symM(M_) Kt = (sp.sqrt(Ks.det()) * Ks.inv()) CL = cod(G, [[Kt[i, j] for j in R3] for i in R3]) mods = sorted({abs(v) for v in CL.values() if v != 0}) check("Lie engine on (1,-1,0): orthonormal components of modulus 2, invariant norm 16", mods == [2] and 2 * sum(v**2 for v in CL.values()) == 16, f"moduli {mods}") # positive control (a): the unit S^3 in hyperspherical coordinates ch, th, ph = sp.symbols("chi theta phi", positive=True) Y = [ch, th, ph] gS = sp.diag(1, sp.sin(ch)**2, sp.sin(ch)**2 * sp.sin(th)**2) GamS, RicS, giS = coord_pipeline(gS, Y) RicS = RicS.applyfunc(sp.simplify) MS = (sp.Rational(1, 2) * sp.simplify(sum(giS[i, j] * RicS[i, j] for i in R3 for j in R3)) * sp.eye(3) - giS * RicS).applyfunc(sp.simplify) KS = (gS * sp.sqrt(MS.det()) * MS.inv()).applyfunc(sp.simplify) CS = codazzi_coord(GamS, KS, Y) okS = RicS == (2 * gS).applyfunc(sp.simplify) and MS == sp.eye(3) okS &= all(sp.simplify(sp.expand_trig(v)) == 0 for v in CS.values()) random.seed(3461) pts = [(sp.Rational(random.randint(200, 2900), 1000), sp.Rational(random.randint(200, 2900), 1000), sp.Rational(random.randint(0, 6000), 1000)) for _ in range(5)] worstS = max(abs(v.evalf(40, subs={ch: p[0], th: p[1], ph: p[2]})) for v in CS.values() for p in pts) okS &= worstS < 1e-30 check("control (a), unit S^3: Ric = 2 gamma, M = Id, residual zero (expand_trig) and < 1e-30 at 5 points", okS, f"max |C| = {sp.N(worstS, 3)}") # positive control (b): the Flamm metric of MA Example 3.3, realizable, not homogeneous, det M > 0, M indefinite rho, mm = sp.symbols("rho m", positive=True) Yf = [rho, th, ph] gF = sp.diag(1 / (1 - 2 * mm / rho), rho**2, rho**2 * sp.sin(th)**2) GamF, RicF, giF = coord_pipeline(gF, Yf) RicF = RicF.applyfunc(sp.simplify) MF = (sp.Rational(1, 2) * sp.simplify(sum(giF[i, j] * RicF[i, j] for i in R3 for j in R3)) * sp.eye(3) - giF * RicF).applyfunc(sp.simplify) KFmix = (sp.sqrt(MF.det()) * MF.inv()).applyfunc(sp.simplify) CF = codazzi_coord(GamF, (gF * KFmix).applyfunc(sp.simplify), Yf) kmerF = -sp.sqrt(2 * mm) / (2 * rho**sp.Rational(3, 2)) ksphF = sp.sqrt(2 * mm) / rho**sp.Rational(3, 2) okF = MF == sp.diag(2 * mm / rho**3, -mm / rho**3, -mm / rho**3) okF &= (KFmix + sp.diag(kmerF, ksphF, ksphF)).applyfunc(sp.simplify) == sp.zeros(3, 3) okF &= all(sp.simplify(v) == 0 for v in CF.values()) check("control (b), Flamm metric: M = (2m, -m, -m)/rho^3, K~ = -diag(k_mer, k_sph, k_sph) (MA Ex. 3.3), Cod K~ = 0", okF, f"M = {list(MF.diagonal())}") # ------------------------------------------------------------------------------------------------------------- P3 print("P3 the Minkowski-Hsiung formulas on radial graphs") mp.mp.dps = 40 def mink(epsv): def parts(t): r = 1 + epsv * mp.cos(2 * t) r1 = -2 * epsv * mp.sin(2 * t) r2 = -4 * epsv * mp.cos(2 * t) w = mp.sqrt(r**2 + r1**2) km = (r**2 + 2 * r1**2 - r * r2) / w**3 rcot = -4 * epsv * mp.cos(t)**2 # r' cot chi, smooth at the poles ks = (r - rcot) / (r * w) H1 = (km + 2 * ks) / 3 H2 = (2 * km * ks + ks**2) / 3 s = r**2 / w dV = 4 * mp.pi * w * r**2 * mp.sin(t)**2 return H1, H2, s, dV I1 = mp.quad(lambda t: (1 - parts(t)[0] * parts(t)[2]) * parts(t)[3], [0, mp.pi / 2, mp.pi]) I2 = mp.quad(lambda t: (parts(t)[0] - parts(t)[1] * parts(t)[2]) * parts(t)[3], [0, mp.pi / 2, mp.pi]) A1 = mp.quad(lambda t: abs(1 - parts(t)[0] * parts(t)[2]) * parts(t)[3], [0, mp.pi / 4, mp.pi / 2, 3 * mp.pi / 4, mp.pi]) return I1, I2, A1 ok3 = True det3_ = [] for ev in (mp.mpf(1) / 10, mp.mpf(3) / 10): I1, I2, A1 = mink(ev) ok3 &= abs(I1) < mp.mpf("1e-25") and abs(I2) < mp.mpf("1e-25") and A1 > mp.mpf("1e-3") det3_.append(f"eps={mp.nstr(ev, 2)}: {mp.nstr(I1, 3)}, {mp.nstr(I2, 3)}; int|1-H1 s| = {mp.nstr(A1, 5)}") check("int (1 - H1 s) dV = int (H1 - H2 s) dV = 0 on a convex and an indefinite graph; integrands not zero", ok3, "; ".join(det3_)) k1, k2, k3 = sp.symbols("k1:4", real=True) H1s = (k1 + k2 + k3) / 3 H2s = (k1 * k2 + k1 * k3 + k2 * k3) / 3 check("18 (H1^2 - H2) = sum_{i 0 for v in (l1, l2))] ok5 &= len(pos) >= 1 and all(sp.simplify(s_.get(l1, l1) - l3) == 0 and sp.simplify(s_.get(l2, l2) - l3) == 0 for s_ in pos) check("Cod linear in k; minors +-H/4; kernel k1=k2=k3 where H != 0; R = -H/2; Heron; det M < 0 on H = 0; " "M scalar => l1 = l2 = l3", ok5, f"det M on H = 0: {dH}") # ------------------------------------------------------------------------------------------------------------- P6 print("P6 E(2) and E(1,1)") M6 = symM(milnor(l1, l2, sp.Integer(0))[3]) ok6 = sp.expand(2 * sum(M6[i, i] for i in R3) + (l1 - l2)**2 / 2) == 0 ok6 &= sp.expand(M6.det() + (l1 - l2)**4 * (3 * l1 + l2) * (l1 + 3 * l2) / 64) == 0 ok6 &= symM(milnor(1, -3, 0)[3]).rank() == 2 and symM(milnor(3, -1, 0)[3]).rank() == 2 ok6 &= symM(milnor(1, -1, 0)[3]).det() > 0 and symM(milnor(1, -4, 0)[3]).det() < 0 ok6 &= symM(milnor(2, 1, 0)[3]).det() < 0 check("R = -(l1-l2)^2/2; det M = -(l1-l2)^4 (3l1+l2)(l1+3l2)/64; rank 2 on the boundary lines", ok6, f"det M = {sp.factor(M6.det())}") # ------------------------------------------------------------------------------------------------------------- P7 print("P7 the degenerate left-invariant metrics") L = [sp.Symbol("L1"), sp.Symbol("L2"), sp.Symbol("L3")] MU = symM(milnor(*L)[3]) ok7 = True for i, j in [(0, 1), (0, 2), (1, 2)]: for fix in L: others = [v for v in L if v != fix] for s_ in sp.solve([MU[i, i].subs(fix, 1), MU[j, j].subs(fix, 1)], others, dict=True): if all(sp.im(s_[v]) == 0 for v in others): ok7 &= sp.simplify(MU.subs(fix, 1).subs(s_)) == sp.zeros(3, 3) check("unimodular: two M-eigenvalues vanish together only where M = 0 (flat)", ok7) a, n2, n3 = sp.symbols("a n2 n3", real=True) cB, GB, RicB, MBl = lie([[0, 0, 0], [0, n2, 0], [0, 0, n3]], [a, 0, 0]) MB = symM(MBl) Pp, Qp = 2 * a - n2 + n3, 2 * a + n2 - n3 Sp_ = 4 * a**2 - 3 * n2**2 + 2 * n2 * n3 + n3**2 Tp = 4 * a**2 + n2**2 + 2 * n2 * n3 - 3 * n3**2 ok7b = sp.expand(MB.det() + Pp * Qp * Sp_ * Tp / 64) == 0 mins = [sp.expand(MB.extract(list(rw), list(cl)).det()) for rw in [(0, 1), (0, 2), (1, 2)] for cl in [(0, 1), (0, 2), (1, 2)]] solB = sp.solve([m for m in mins if m != 0], [n2, n3], dict=True) ok7b &= sorted([(s_[n2], s_[n3]) for s_ in solB], key=str) == sorted([(-a, a), (a, -a)], key=str) okp = True for s_ in solB: Ms = MB.subs(s_) ok_e = sorted(Ms.eigenvals().items(), key=str) == sorted({0: 2, -4 * a**2: 1}.items(), key=str) v = (Ms - (-4 * a**2) * sp.eye(3)).nullspace()[0] v = v / sp.sqrt(v.dot(v)) Gs = [[[sp.expand(GB[i][j][k].subs(s_)) for k in R3] for j in R3] for i in R3] nab = [[sp.simplify(sum(v[j] * Gs[i][j][k] for j in R3)) for k in R3] for i in R3] okp &= ok_e and all(e == 0 for row in nab for e in row) check("class B: det M = -PQST/64; rank M <= 1 iff (n2, n3) = (-a, a), (a, -a); there M ~ (0,0,-4a^2) and the " "-4a^2 eigenvector is parallel (H^2 x R)", ok7b and okp, f"solutions {solB}") # ------------------------------------------------------------------------------------------------------------- P8 print("P8 class B: test_215's algebra from the engine") CB = {k: sp.factor(v) for k, v in cod(GB, adj3([[MB[i, j] for j in R3] for i in R3])).items()} nzB = {k: v for k, v in CB.items() if v != 0} ok8 = all(sp.rem(sp.expand(v), n2 - n3, n2) == 0 for v in nzB.values()) ok8 &= sp.expand(CB[(1, 2, 0)] + (n2 - n3)**2 * (n2 + n3) * Pp * Qp / 8) == 0 ok8 &= sp.expand(CB[(0, 1, 1)].subs(n3, -n2) - 4 * a * n2**2 * (n2**2 - a**2)) == 0 ok8 &= sp.expand(MB.det().subs(n3, -n2) + (a**2 - n2**2)**3) == 0 ok8 &= sp.expand(MB.det().subs(n3, n2) + a**6) == 0 check("all components divisible by (n2-n3); C_231 = -(n2-n3)^2(n2+n3)PQ/8; on n3=-n2 C_122 = 4a n2^2(n2^2-a^2), " "det M = -(a^2-n2^2)^3; det M = -a^6 at n2 = n3", ok8, f"{len(nzB)} nonzero components") # ------------------------------------------------------------------------------------------------------------- P9 print("P9 random left-invariant metrics, exact") random.seed(34609) HALF = Fr(1, 2) def rnd(lo=-9, hi=9, den=4): return Fr(random.randint(lo, hi), random.randint(1, den)) def cayley(): s = [rnd(-3, 3, 3) for _ in range(3)] S = [[0, -s[2], s[1]], [s[2], 0, -s[0]], [-s[1], s[0], 0]] IpS = [[Fr(int(i == j)) + S[i][j] for j in R3] for i in R3] ImS = [[Fr(int(i == j)) - S[i][j] for j in R3] for i in R3] dI = det3(IpS) inv = [[adj3(IpS)[i][j] / dI for j in R3] for i in R3] return [[sum(ImS[i][k] * inv[k][j] for k in R3) for j in R3] for i in R3] def sample(kind): if kind == "uni": n = [[rnd() if i == j else Fr(0) for j in R3] for i in R3] a = [Fr(0)] * 3 elif kind == "B": a = [rnd() for _ in R3] while all(v == 0 for v in a): a = [rnd() for _ in R3] Pm = [[Fr(int(i == j)) - a[i] * a[j] / sum(v * v for v in a) for j in R3] for i in R3] m0 = [[rnd() for _ in R3] for _ in R3] m0 = [[m0[i][j] + m0[j][i] for j in R3] for i in R3] n = [[sum(Pm[i][k] * m0[k][l] * Pm[l][j] for k in R3 for l in R3) for j in R3] for i in R3] else: # near the Schuecking line n3 = -n2, a^2 < n2^2 nn2 = rnd(1, 9, 2) * random.choice([1, -1]) aa = abs(nn2) * Fr(random.randint(1, 9), 10) * random.choice([1, -1]) t = Fr(random.randint(-5, 5), 50) n = [[Fr(0)] * 3, [Fr(0), nn2, Fr(0)], [Fr(0), Fr(0), -nn2 + t]] a = [aa, Fr(0), Fr(0)] Q = cayley() c = rotate(structure(n, a), Q) G = connection(c, HALF) M = Mmat(ricci(c, G), HALF) return c, G, M counts = {"uni": [0, 0], "B": [0, 0], "line": [0, 0]} ok9 = True for kind, N in (("uni", 200), ("B", 200), ("line", 100)): for _ in range(N): c, G, M = sample(kind) ok9 &= jacobi_ok(c) if det3(M) > 0: counts[kind][0] += 1 scalar = all(M[i][j] == (M[0][0] if i == j else 0) for i in R3 for j in R3) Cv = cod(G, adj3(M)) if not scalar: ok9 &= any(v != 0 for v in Cv.values()) else: counts[kind][1] += 1 ok9 &= counts["uni"][0] >= 20 and counts["B"][0] + counts["line"][0] >= 20 check("Jacobi for all; every det M > 0, non-scalar sample has Cod(adj M) != 0 exactly; >= 20 of each kind with " "det M > 0", ok9, f"det M > 0 (scalar): {counts}") def near_round(d): c, G, Ric, M = milnor(Fr(1), 1 + d, 1 + 2 * d, HALF) return sum(v * v for v in cod(G, adj3(M)).values()) r3, r6 = near_round(Fr(1, 10**3)) / Fr(1, 10**6), near_round(Fr(1, 10**6)) / Fr(1, 10**12) check("near-round: |Cod|^2/d^2 agrees at d = 1e-3, 1e-6 to 1%; zero at d = 0", abs(r3 / r6 - 1) < Fr(1, 100) and near_round(Fr(0)) == 0, f"{float(r3):.6f} vs {float(r6):.6f}") # ------------------------------------------------------------------------------------------------------------ P10 print("P10 the product geometries") u, v_, t = sp.symbols("u v t", real=True) gH = sp.diag(1 / v_**2, 1 / v_**2, 1) GamH, RicH, giH = coord_pipeline(gH, [u, v_, t]) RicH = RicH.applyfunc(sp.simplify) MH = (sp.Rational(1, 2) * sp.simplify(sum(giH[i, j] * RicH[i, j] for i in R3 for j in R3)) * sp.eye(3) - giH * RicH).applyfunc(sp.simplify) ok10 = MH == sp.diag(0, 0, -1) and MH.rank() == 1 ok10 &= all(sp.simplify(GamH[cc][a_i][2]) == 0 for cc in R3 for a_i in R3) # nabla_a d_t = Gam^c_{a t} d_c = 0 p_, q_, w_ = sp.symbols("p q w", real=True) R2 = sp.Matrix([[sp.cos(w_), -sp.sin(w_), 0], [sp.sin(w_), sp.cos(w_), 0], [0, 0, 1]]) Kr2 = R2 * sp.diag(p_, q_, 0) * R2.T Ad = Kr2.adjugate() ok10 &= (Kr2 * Ad).applyfunc(sp.simplify) == sp.zeros(3, 3) ok10 &= (Ad * sp.Matrix([0, 0, 1]) - p_ * q_ * sp.Matrix([0, 0, 1])).applyfunc(sp.simplify) == sp.zeros(3, 1) ok10 &= (Kr2 * sp.Matrix([0, 0, 1])).applyfunc(sp.simplify) == sp.zeros(3, 1) random.seed(34610) rk = [] for rank_ in (3, 2, 1): B_ = sp.Matrix(3, rank_, lambda i, j: random.randint(-5, 5)) Kq = B_ * sp.diag(*[random.choice([1, -1]) * random.randint(1, 4) for _ in range(rank_)]) * B_.T rk.append((Kq.rank(), Kq.adjugate().rank())) ok10 &= rk == [(3, 3), (2, 1), (1, 0)] rc = sp.symbols("rho", positive=True) gC = sp.diag(rc**2, rc**2 * sp.sin(th)**2, 1) GamC, RicC, giC = coord_pipeline(gC, [th, ph, t]) RicC = RicC.applyfunc(sp.simplify) MC = (sp.Rational(1, 2) * sp.simplify(sum(giC[i, j] * RicC[i, j] for i in R3 for j in R3)) * sp.eye(3) - giC * RicC).applyfunc(sp.simplify) KCmix = sp.diag(1 / rc, 1 / rc, 0) ok10 &= (KCmix.adjugate() - MC).applyfunc(sp.simplify) == sp.zeros(3, 3) CC = codazzi_coord(GamC, (gC * KCmix).applyfunc(sp.simplify), [th, ph, t]) ok10 &= all(sp.simplify(e) == 0 for e in CC.values()) check("H^2 x R: M rank 1, image d_t, d_t parallel; rank-2 K: K adj K = 0, im adj K = ker K; rank adj = 3, 1, 0; " "cylinder: adj K = M, Cod K = 0", ok10, f"(rank K, rank adj K) = {rk}") print() if FAILURES: print(f"test_346: {len(FAILURES)} FAILED") sys.exit(1) print("test_346: ALL PASS")