packages feed

egison-5.1.0: sample/math/geometry/thurston.egi

---
--- Calculation of the WCS Invariant on the Thurston Example (Section 4)
---

declare symbol θ₁, θ₂, θ₃, θ₄, κ, p

def x~i := [| θ₁, θ₂, θ₃, θ₄ |]~i

def β := `(1 + θ₂ - θ₂^2)

def g_i_j :=
  [|[| 1, 0, 0, 0 |],
    [| 0, 1, 0, 0 |],
    [| 0, 0, κ / (sqrt β), (-1 * θ₂ * κ) / (sqrt β)  |],
    [| 0, 0, (-1 * θ₂ * κ) / (sqrt β), (`(1 + θ₂) * κ) / (sqrt β) |]|]

def g~i~j :=
  [|[| 1, 0, 0, 0 |],
    [| 0, 1, 0, 0 |],
    [| 0, 0, `(1 + θ₂) / (κ * (sqrt β)), θ₂ / ((sqrt β) * κ) |],
    [| 0, 0, θ₂ / ((sqrt β) * κ), 1 / ((sqrt β) * κ) |]|]

-- Sanity check on the metric definition.
assertEqual "g_3_3"
  g_3_3
  (κ / sqrt `(- θ₂^2 + θ₂ + 1))
assertEqual "g~3~3"
  g~3~3
  (`(1 + θ₂) / (κ * sqrt `(- θ₂^2 + θ₂ + 1)))
assertEqual "β"
  β
  `(- θ₂^2 + θ₂ + 1)

def Γ~c_a_b := withSymbols [e]
  (1 / 2) * g~c~e . (∂/∂ g_b_e x~a + ∂/∂ g_a_e x~b - ∂/∂ g_a_b x~e)

def R_i_j_k~l := withSymbols [a]
  ∂/∂ Γ~l_j_k x~i - ∂/∂ Γ~l_i_k x~j + Γ~l_i_a . Γ~a_j_k - Γ~l_j_a . Γ~a_i_k

def R_i_j_k_l := withSymbols [a] R_i_j_k~a . g_a_l

def J_a_b :=
  [|[| 0, 1, 0, 0 |],
    [| -1, 0, 0, 0 |],
    [| 0, 0, 0, κ |],
    [| 0, 0, -1 * κ, 0 |]|]

def J_a~c := withSymbols [b] J_a_b . g~b~c

-- The original EMR Thurston example used an inline ∇J definition (rather
-- than going through an abstracted variadic `∇_c T~(a_1)...~(a_r)_(b_1)..._(b_k)`).
-- The abstracted form caused type-inference errors (the variadic-index
-- pattern doesn't unify with the concrete monomorphic Tensor types),
-- preventing the file from compiling. Inline form matches the original
-- paper's expected S value (same simplified result).
--
-- Sign convention follows the EMR original (`+` Γ terms). The standard
-- textbook ∇ has `-` for lower indices and `+` for upper, but for the
-- WCS invariant the absolute sign cancels out in the product.
def ∇J_m_a_b :=
  withSymbols [n]
    ∂/∂ J_a_b x~m + Γ~n_m_a . J_n_b + Γ~n_m_b . J_a_n

def ∇J~m_a_b :=
  withSymbols [t]
    ∇J_t_a_b . g~t~m

def ∇J_m~a_b :=
  withSymbols [t]
    ∇J_m_t_b . g~t~a

def ∇J_m_a~b :=
  withSymbols [t]
    ∇J_m_a_t . g~t~b

def δ :=
  generateTensor
    (\match as list integer with
       | [$n, #n] -> 1
       | [_, _] -> 0)
    [5, 5]

def R'{_i_j}_k~l : Tensor MathValue :=
  generateTensor
    (\match as list integer with
       | [#1, #1, _, _] -> 0
       | [_, _, #1, #1] -> 0
       | [#1, $b, #1, $d] -> -1 * p^2 * δ~(b - 1)_(d - 1)
       | [$a, #1, #1, $d] ->      p^2 * δ~(a - 1)_(d - 1)
       | [#1, $b, $c, #1] ->      p^2 * g_(b - 1)_(c - 1)
       | [$a, #1, $c, #1] -> -1 * p^2 * g_(a - 1)_(c - 1)
       | [#1, $b, $c, $d] -> -1 * p * ∇J_(b - 1)_(c - 1)~(d - 1)
       | [$a, #1, $c, $d] ->      p * ∇J_(a - 1)_(c - 1)~(d - 1)
       | [$a, $b, #1, $d] -> -1 * p * ∇J~(d - 1)_(a - 1)_(b - 1)
       | [$a, $b, $c, #1] ->      p * ∇J_(c - 1)_(a - 1)_(b - 1)
       | [$a, $b, $c, $d] -> R_(a - 1)_(b - 1)_(c - 1)~(d - 1)
                             + -1 * p^2 * J_(b - 1)_(c - 1) * J_(a - 1)~(d - 1)
                             +      p^2 * J_(a - 1)_(c - 1) * J_(b - 1)~(d - 1)
                             +  2 * p^2 * J_(a - 1)_(b - 1) * J_(c - 1)~(d - 1))
    [5, 5, 5, 5]

def S :=
  withSymbols [i, j, k]
    let (es, os) := evenAndOddPermutations 5 in
      sum (map (\σ -> R'_(σ 1)_j_1~i . R'_(σ 2)_(σ 3)_k~j . R'_(σ 4)_(σ 5)_i~k) es) -
      sum (map (\σ -> R'_(σ 1)_j_1~i . R'_(σ 2)_(σ 3)_k~j . R'_(σ 4)_(σ 5)_i~k) os)

-- WCS (Wodzicki-Chern-Simons) invariant on the Thurston example (Section 4
-- of "Diffeomorphism Groups of Circle Bundles over Integral Symplectic
-- Manifolds"). The Wolfram-simplified form of S, with β = 1+θ₂-θ₂², is:
--   S = p² κ (-25 - 640 p² β² + 3072 p⁴ β⁴) / (16 β⁴)
-- See <https://github.com/egisatoshi/EMR-Paper-Computation>.
--
-- The identity needs the defining relations of the quoted atoms
-- (β = 1+θ₂-θ₂² and `(1+θ₂) = 1+θ₂), which a structural comparison
-- cannot see.  Clear the Laurent denominators and reduce the difference
-- modulo the Groebner basis of the quote relations, applied once at the
-- point of comparison (declaring the relations as always-on rules would
-- blow up the intermediate computation).
def sExpected := p^2 * κ * (- 25 - 640 * p^2 * β^2 + 3072 * p^4 * β^4) / (16 * β^4)
def quoteGb := groebnerBasis ['(1 + θ₂ - θ₂^2 - β), '(1 + θ₂ - `(1 + θ₂))]

assertEqual "WCS invariant S"
  (polyNF quoteGb (16 * β^8 * (S - sExpected)))
  0

-- The same identity through plain rational arithmetic: expanding the
-- quotes makes both sides rational functions of θ₂, which polynomial
-- GCD reduces.
assertEqual "WCS invariant S (expanded)"
  ((expandAll S) = (expandAll sExpected))
  True

-- For reference, the raw (pre-reduction) result of the original EMR
-- computation was:
-- (1536 p^6 κ Sqrt[(1 + θ₂ - θ₂^2)]^16 - 1536 p^6 θ₂^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^14 - 576 p^4 (1 + θ₂) κ Sqrt[(1 + θ₂ - θ₂^2)]^12 + 1536 p^6 (1 + θ₂) κ Sqrt[(1 + θ₂ - θ₂^2)]^14 + 8 p^2 (1 - 2 θ₂)^4 (1 + θ₂)^3 θ₂^2 κ - 88 p^2 (1 - 2 θ₂)^2 (1 + θ₂)^2 θ₂^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ + 48 p^2 (1 - 2 θ₂)^3 (1 + θ₂)^2 θ₂^3 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ - 12 p^2 (1 - 2 θ₂)^4 (1 + θ₂)^2 θ₂^4 κ - 24 p^2 (1 - 2 θ₂)^3 (1 + θ₂)^2 θ₂^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^2 + 288 p^4 (1 - 2 θ₂)^2 (1 + θ₂)^2 θ₂^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^6 - 160 p^2 (1 - 2 θ₂) (1 + θ₂) θ₂^3 Sqrt[(1 + θ₂ - θ₂^2)]^6 κ + 128 p^2 (1 - 2 θ₂)^2 (1 + θ₂) θ₂^4 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ - 48 p^2 (1 - 2 θ₂)^3 (1 + θ₂) θ₂^5 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ - 80 p^2 (1 - 2 θ₂)^2 (1 + θ₂) θ₂^3 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ + 768 p^4 (1 - 2 θ₂) (1 + θ₂) θ₂^3 Sqrt[(1 + θ₂ - θ₂^2)]^8 κ + 8 p^2 (1 - 2 θ₂)^4 (1 + θ₂) θ₂^6 κ + 24 p^2 (1 - 2 θ₂)^3 (1 + θ₂) θ₂^4 κ Sqrt[(1 + θ₂ - θ₂^2)]^2 - 288 p^4 (1 - 2 θ₂)^2 (1 + θ₂) θ₂^4 κ Sqrt[(1 + θ₂ - θ₂^2)]^6 + 112 p^2 (1 - 2 θ₂) (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^2 κ + 20 p^2 (1 - 2 θ₂)^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^2 κ - 64 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^4 κ + 96 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^5 (1 - 2 θ₂) κ - 56 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^6 (1 - 2 θ₂)^2 κ - 80 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^4 (1 - 2 θ₂) κ + 384 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^10 θ₂^4 κ + 16 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^2 θ₂^7 (1 - 2 θ₂)^3 κ + 40 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^5 (1 - 2 θ₂)^2 κ - 384 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^5 (1 - 2 θ₂) κ + 32 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^3 κ + 24 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^3 (1 - 2 θ₂) κ - 448 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^10 θ₂^3 κ - 2 p^2 (1 - 2 θ₂)^4 θ₂^8 κ - 8 p^2 (1 - 2 θ₂)^3 θ₂^6 κ Sqrt[(1 + θ₂ - θ₂^2)]^2 + 96 p^4 (1 - 2 θ₂)^2 θ₂^6 κ Sqrt[(1 + θ₂ - θ₂^2)]^6 - 10 p^2 (1 - 2 θ₂)^2 θ₂^4 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ - 16 p^2 (1 + θ₂)^3 (1 - 2 θ₂)^3 Sqrt[(1 + θ₂ - θ₂^2)]^2 θ₂ κ + 64 p^2 (1 + θ₂)^2 (1 - 2 θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂ κ + 40 p^2 (1 + θ₂)^2 (1 - 2 θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂ κ - 384 p^4 (1 + θ₂)^2 (1 - 2 θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^8 κ θ₂ + 96 p^2 (1 + θ₂) θ₂^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 κ - 32 p^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂ κ - 24 p^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂ (1 - 2 θ₂) κ + 448 p^4 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^10 κ θ₂ - 32 p^2 (1 - 2 θ₂) (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 κ - 10 p^2 (1 - 2 θ₂)^2 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ + 8 p^2 (1 - 2 θ₂)^3 (1 + θ₂)^3 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ + 16 p^2 (1 - 2 θ₂)^2 (1 + θ₂)^3 κ Sqrt[(1 + θ₂ - θ₂^2)]^4 - 2 p^2 (1 - 2 θ₂)^4 (1 + θ₂)^4 κ - 96 p^4 (1 - 2 θ₂)^2 (1 + θ₂)^3 κ Sqrt[(1 + θ₂ - θ₂^2)]^6 + 4 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 (1 + θ₂) (1 - 2 θ₂) κ + 8 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 (1 + θ₂) κ - 112 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^10 κ (1 + θ₂) - 8 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^2 κ - 4 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^2 (1 - 2 θ₂) κ + 112 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^10 θ₂^2 κ - 48 p^2 (1 - 2 θ₂)^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂ κ + 40 p^2 (1 - 2 θ₂)^3 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^2 κ + 12 p^2 (1 - 2 θ₂)^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 κ - 20 p^2 (1 - 2 θ₂)^3 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ + 112 p^2 (1 - 2 θ₂)^2 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^2 κ + 32 p^2 (1 - 2 θ₂) (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^8 κ - 32 p^2 (1 - 2 θ₂)^2 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 κ + 40 p^2 (1 - 2 θ₂)^3 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂ κ - 21 p^2 (1 - 2 θ₂)^4 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^2 θ₂^2 κ + 7 p^2 (1 - 2 θ₂)^4 (1 + θ₂)^3 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ - 80 p^2 (1 - 2 θ₂)^3 (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^3 κ + 21 p^2 (1 - 2 θ₂)^4 (1 + θ₂) θ₂^4 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ - 32 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^2 (1 - 2 θ₂) κ + 48 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^3 (1 - 2 θ₂)^2 κ - 16 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^10 θ₂^2 κ + 64 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 θ₂^3 (1 - 2 θ₂) κ - 80 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 θ₂^4 (1 - 2 θ₂)^2 κ + 40 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^4 θ₂^5 (1 - 2 θ₂)^3 κ - 20 p^2 (1 - 2 θ₂)^3 θ₂^4 Sqrt[(1 + θ₂ - θ₂^2)]^4 κ - 12 p^2 (1 - 2 θ₂)^2 θ₂^2 Sqrt[(1 + θ₂ - θ₂^2)]^6 κ - 7 p^2 (1 - 2 θ₂)^4 θ₂^6 Sqrt[(1 + θ₂ - θ₂^2)]^2 κ - 64 p^4 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^10 κ - 32 p^2 (1 + θ₂)^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 κ + 16 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^10 (1 + θ₂) κ - 64 p^2 Sqrt[(1 + θ₂ - θ₂^2)]^8 (1 - 2 θ₂) θ₂ (1 + θ₂) κ + 576 p^4 θ₂ Sqrt[(1 + θ₂ - θ₂^2)]^12 κ - 144 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^12 κ - 224 p^4 (1 - 2 θ₂)^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^10 (1 + θ₂) - 384 p^4 κ (1 - 2 θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^12 θ₂ + 224 p^4 (1 - 2 θ₂)^2 θ₂^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^10 + 192 p^4 (1 - 2 θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^12 κ + 128 p^4 Sqrt[(1 + θ₂ - θ₂^2)]^14 κ - 320 p^4 θ₂^2 Sqrt[(1 + θ₂ - θ₂^2)]^10 κ (1 + θ₂) + 128 p^4 (1 - 2 θ₂) (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^10 κ - 128 p^4 θ₂^2 (1 - 2 θ₂) κ Sqrt[(1 + θ₂ - θ₂^2)]^10 + 192 p^4 (1 + θ₂)^2 (1 - 2 θ₂) κ Sqrt[(1 + θ₂ - θ₂^2)]^8 - 384 p^4 θ₂^2 (1 - 2 θ₂) (1 + θ₂) κ Sqrt[(1 + θ₂ - θ₂^2)]^8 + 192 p^4 θ₂^4 Sqrt[(1 + θ₂ - θ₂^2)]^8 κ (1 - 2 θ₂) - 256 p^4 θ₂ Sqrt[(1 + θ₂ - θ₂^2)]^10 κ (1 - 2 θ₂) (1 + θ₂) + 256 p^4 θ₂^3 Sqrt[(1 + θ₂ - θ₂^2)]^10 κ (1 - 2 θ₂) + 128 p^4 θ₂^2 (1 - 2 θ₂)^2 κ (1 + θ₂) Sqrt[(1 + θ₂ - θ₂^2)]^8 - 64 p^4 θ₂^4 (1 - 2 θ₂)^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^8 - 64 p^4 (1 - 2 θ₂)^2 (1 + θ₂)^2 κ Sqrt[(1 + θ₂ - θ₂^2)]^8)/(16 Sqrt[(1 + θ₂ - θ₂^2)]^16)
-- The above result is simplified using the Wolfam language as follows:
-- (p^2 (-25 - 640 p^2 (1 + θ₂- θ₂^2)^2 + 3072 p^4 (1 + θ₂ - θ₂^2)^4) κ) / (16 (1 + θ₂ - θ₂^2)^4)