"""Our analysis for https://factchecked.science/dark-energy-instability/ (Fact-Checked Science, October 31, 2026).

Arithmetic checks on numbers and equations typed from the sources; Python standard library only.

1. Redshift-distance coefficients H0*dL = z + Q z^2 + C z^3 + O(z^4), in exact fractions, for
   (a) flat LambdaCDM, E(z)^2 = Om (1+z)^3 + OL                      (the standard form)
   (b) the form printed in Smoller, Temple & Vogler (2017), eq. (3.39) of the journal version
       ((3.18) of arXiv:1412.4001v1): E(z)^2 = OL (1+z)^2 + Om (1+z)^3, with the flat-space integral (3.38)
   (c) an open dust universe, Om = 0.3, Ok = 0.7, no Lambda (with the sinh of the open geometry)
   (d) Einstein-de Sitter (Om = 1)
   and Q = (1 - q0)/2, the model-independent link between Q and the deceleration parameter
   (Visser 2004, Class. Quantum Grav. 21, 2603, eq. 46).
   Then numerical luminosity distances (Simpson's rule) and distance-modulus differences at z = 0.1, 0.5, 1.
2. The STV ODE of orders n = 1 and n = 2 as printed in Alexander, Temple & Vogler (2026, Proc. R. Soc. A 482:
   20250912, eqs. 2.34-2.35 and 2.48-2.51, the same as in arXiv:2510.14228v1): check that SM = (4/3, 2/3, 40/27,
   2/9) and M = (0, 1, 0, 0) are rest points
   and compute the eigenvalues of the linearisation at SM in the log-time tau = ln t.
3. The meaning of the eigenvalues lambda_An = 2n/3, lambda_Bn = (2n-5)/3 (eq. 2.29): a perturbation of the
   coefficient z_2n of xi^(2n) = (r/t)^(2n), with z = kappa rho r^2, changes the density contrast at a fixed areal
   radius r by t^(lambda - 2n + 2) r^(2n-2) and, following the matter (r grows like t^(2/3) in the background),
   by t^(lambda - (2n-2)/3) at a fixed comoving radius.

Run: python3 cosmography_check.py              (prints the results)
     python3 cosmography_check.py OUTFILE      (writes them to OUTFILE, as cosmography_check.txt was made)
"""
import math
import sys
from fractions import Fraction as F

N = 4  # series terms kept: 1, z, z^2, z^3


def mul(a, b):
    out = [F(0)] * N
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            if i + j < N:
                out[i + j] += x * y
    return out


def pow_series(base, p):
    """(1 + u)^p as a series, where base = [1, u1, u2, ...] (base[0] must be 1); p a Fraction."""
    assert base[0] == 1
    u = [F(0)] + base[1:]
    out = [F(0)] * N
    out[0] = F(1)
    term = [F(1)] + [F(0)] * (N - 1)
    coef = F(1)
    for k in range(1, N):
        term = mul(term, u)
        coef = coef * (p - (k - 1)) / k
        for i in range(N):
            out[i] += coef * term[i]
    return out


def one_plus_z_pow(n):
    return pow_series([F(1), F(1), F(0), F(0)], F(n))


def integrate(s):
    """Integral from 0 to z of a series in y: returns a series in z."""
    out = [F(0)] * N
    for i in range(N - 1):
        out[i + 1] = s[i] / (i + 1)
    return out


def dl_series(terms, ok=F(0)):
    """terms: list of (Omega, power of (1+z)) making E^2; ok: curvature for the sinh/sin correction."""
    e2 = [F(0)] * N
    for om, n in terms:
        s = one_plus_z_pow(n)
        for i in range(N):
            e2[i] += om * s[i]
    inv_e = pow_series(e2, F(-1, 2))
    chi = integrate(inv_e)  # comoving distance times H0
    if ok != 0:
        # sinh(sqrt(ok) chi)/sqrt(ok) = chi + ok chi^3/6 + ...
        chi3 = mul(mul(chi, chi), chi)
        chi = [chi[i] + ok * chi3[i] / 6 for i in range(N)]
    return mul([F(1), F(1), F(0), F(0)], chi)  # (1+z) * chi


def e_func(terms):
    return lambda z: math.sqrt(sum(float(om) * (1 + z) ** n for om, n in terms))


def dl_numeric(terms, z, ok=0.0, steps=2000):
    e = e_func(terms)
    h = z / steps
    s = 1 / e(0) + 1 / e(z)
    for i in range(1, steps):
        s += (4 if i % 2 else 2) / e(i * h)
    chi = s * h / 3
    if ok > 0:
        chi = math.sinh(math.sqrt(ok) * chi) / math.sqrt(ok)
    return (1 + z) * chi


def main():
    if len(sys.argv) > 1:
        out = open(sys.argv[1], "w", encoding="utf-8", newline="\n")
    else:
        sys.stdout.reconfigure(encoding="utf-8")
        out = sys.stdout

    def p(*a):
        print(*a, file=out)

    om = F(3, 10)
    cases = {
        "(a) flat LambdaCDM, Om=0.3, OL=0.7: E^2 = Om(1+z)^3 + OL": ([(om, 3), (1 - om, 0)], F(0)),
        "(b) STV 2017 eq. (3.39): E^2 = OL(1+z)^2 + Om(1+z)^3, flat integral": ([(1 - om, 2), (om, 3)], F(0)),
        "(c) open dust, Om=0.3, Ok=0.7, no Lambda (sinh)": ([(om, 3), (1 - om, 2)], 1 - om),
        "(d) Einstein-de Sitter, Om=1": ([(F(1), 3)], F(0)),
    }
    p("1. Redshift-distance coefficients, H0 dL = z + Q z^2 + C z^3 + O(z^4)")
    res = {}
    for name, (terms, ok) in cases.items():
        s = dl_series(terms, ok)
        res[name] = s
        p(f"  {name}")
        p(f"     linear {s[1]}, Q = {s[2]} = {float(s[2]):.4f}, C = {s[3]} = {float(s[3]):.5f}; "
          f"q0 from Q = (1-q0)/2: {1 - 2 * s[2]} = {float(1 - 2 * s[2]):.3f}")
    p("  Paper values (STV 2017 eq. 3.42; ATV 2026 pp. 22-23): Q = 17/40 = 0.425, C = -433/2400 = -0.1804")
    p("  General Om, flat LambdaCDM: Q = 1 - 3 Om/4, C = -(5/4) Om + (9/8) Om^2  (q0 = 3 Om/2 - 1)")
    for omv in (F(1), F(3, 10), F(315, 1000), F(0)):
        s = dl_series([(omv, 3), (1 - omv, 0)])
        s2 = dl_series([(1 - omv, 2), (omv, 3)])
        p(f"     Om = {float(omv):.3f}: flat LCDM Q = {float(s[2]):.4f}, C = {float(s[3]):.5f};"
          f"   STV-2017 form Q = {float(s2[2]):.4f}, C = {float(s2[3]):.5f}")
    p("  So the STV-2017 form gives 0.25 <= Q <= 0.5 (q0 from 1/2 to 0: decelerating to coasting),")
    p("  while flat LambdaCDM gives 0.25 <= Q <= 1 (q0 from 1/2 to -1); they agree only at Om = 1.")
    p("  Cross-check by Visser (2004): a component's density scales as (1+z)^(3(1+w)) (eq. 48), so (1+z)^2 means")
    p("  w = -1/3 and a constant density w = -1; w0 is the density-weighted mean of the w's (eq. 49);")
    p("  q0 = (1 + 3 w0)/2 for a flat universe (eq. 20); Q = (1 - q0)/2 (eq. 46).")
    for label, w_l in (("a cosmological constant, w = -1", F(-1)), ("the 2017 form's (1+z)^2 term, w = -1/3", F(-1, 3))):
        w0 = om * 0 + (1 - om) * w_l
        q0 = (1 + 3 * w0) / 2
        p(f"     Om = 0.3, dark energy as {label}: w0 = {w0} = {float(w0):.4f}, q0 = {q0} = {float(q0):.3f}, "
          f"Q = {(1 - q0) / 2} = {float((1 - q0) / 2):.4f}")
    p("")
    p("  Numerical luminosity distances (H0 dL, Simpson's rule) and distance-modulus differences:")
    for z in (0.1, 0.5, 1.0):
        da = dl_numeric([(0.3, 3), (0.7, 0)], z)
        db = dl_numeric([(0.7, 2), (0.3, 3)], z)
        dc = dl_numeric([(0.3, 3), (0.7, 2)], z, ok=0.7)
        de = dl_numeric([(1.0, 3)], z)
        sa = res["(a) flat LambdaCDM, Om=0.3, OL=0.7: E^2 = Om(1+z)^3 + OL"]
        ser = float(sa[1]) * z + float(sa[2]) * z ** 2 + float(sa[3]) * z ** 3
        p(f"   z = {z}: LCDM {da:.4f} (series {ser:.4f}); STV-2017 form {db:.4f}; open dust {dc:.4f}; EdS {de:.4f}")
        p(f"          mu(LCDM) - mu(STV-2017 form) = {5 * math.log10(da / db):+.3f} mag;"
          f" mu(LCDM) - mu(open dust) = {5 * math.log10(da / dc):+.3f} mag;"
          f" mu(LCDM) - mu(EdS) = {5 * math.log10(da / de):+.3f} mag")
    p("")

    # 2. STV ODE (n = 1, 2) as printed in ATV 2026 eqs. (2.34)-(2.35), (2.48)-(2.51), as in arXiv v1; d/dtau = t d/dt.
    def rhs(u):
        z2, w0, z4, w2 = u
        return [2 * z2 - 3 * z2 * w0,
                -F(1, 6) * z2 + w0 - w0 ** 2,
                F(5, 12) * z2 ** 2 * w0 + 4 * z4 - 5 * w0 * z4 - 5 * z2 * w2,
                -F(1, 24) * z2 ** 2 + F(1, 4) * z2 * w0 ** 2 - F(1, 10) * z4 + 3 * w2 - 4 * w0 * w2]

    sm = [F(4, 3), F(2, 3), F(40, 27), F(2, 9)]
    mm = [F(0), F(1), F(0), F(0)]
    p("2. STV ODE of order 2 (eqs. 2.48-2.51): rest points")
    p(f"   F(SM) = {[str(x) for x in rhs(sm)]}   F(M) = {[str(x) for x in rhs(mm)]}")

    def jac(u):
        z2, w0, z4, w2 = u
        return [[2 - 3 * w0, -3 * z2, 0, 0],
                [-F(1, 6), 1 - 2 * w0, 0, 0],
                [F(5, 6) * z2 * w0 - 5 * w2, F(5, 12) * z2 ** 2 - 5 * z4, 4 - 5 * w0, -5 * z2],
                [-F(1, 12) * z2 + F(1, 4) * w0 ** 2, F(1, 2) * z2 * w0 - 4 * w2, -F(1, 10), 3 - 4 * w0]]

    def eig2(a, b, c, d):
        tr, det = a + d, a * d - b * c
        disc = tr * tr - 4 * det
        r = math.sqrt(float(disc))
        return (float(tr) + r) / 2, (float(tr) - r) / 2, tr, det, disc

    for label, u in (("SM", sm), ("M", mm)):
        j = jac(u)
        # block lower-triangular: eigenvalues are those of the two diagonal 2x2 blocks
        e1 = eig2(j[0][0], j[0][1], j[1][0], j[1][1])
        e2 = eig2(j[2][2], j[2][3], j[3][2], j[3][3])
        p(f"   at {label}: n=1 block eigenvalues {e1[0]:.6f}, {e1[1]:.6f} (trace {e1[2]}, det {e1[3]}, disc {e1[4]});"
          f" n=2 block {e2[0]:.6f}, {e2[1]:.6f} (trace {e2[2]}, det {e2[3]}, disc {e2[4]})")
    p("   Paper: at SM (eqs. 2.29-2.31), lambda_A1 = 2/3, lambda_B1 = -1, lambda_A2 = 4/3, lambda_B2 = -1/3;"
      " at M (eq. 2.28), lambda = -1 (degenerate)")
    p("")
    p("3. Eigenvalues (2.29) and the density contrast of a perturbation of the xi^(2n) coefficient")
    p("   exponent at a fixed areal radius: lambda - 2n + 2 (times r^(2n-2)); at a fixed comoving radius"
      " (r ~ t^(2/3)): lambda - (2n-2)/3")
    for n in range(1, 7):
        la, lb = F(2 * n, 3), F(2 * n - 5, 3)
        p(f"   n = {n}: lambda_A = {la}, lambda_B = {lb}; comoving exponents {la - F(2 * n - 2, 3)} and "
          f"{lb - F(2 * n - 2, 3)}")
    p("   Textbook linear perturbations of dust in Einstein-de Sitter: growing mode t^(2/3), decaying mode t^(-1).")
    if out is not sys.stdout:
        out.close()


if __name__ == "__main__":
    main()