binomial.iml

(* Binomial Theorem in Imandra. 
   Theorem 44 in the 100 Theorems list.
  
   Grant Passmore, Imandra
 *)

let rec pow (a:int) (n:int) : int =
  if n <= 0 then 1 else a * pow a (n - 1)
[@@measure (Ordinal.of_int n)] [@@by auto]

lemma pow_zero a = pow a 0 = 1 [@@by auto]

lemma pow_succ a n =
  n >= 0 ==> pow a (n + 1) = a * pow a n
[@@by intros @> auto]

lemma pow_succ_l a t =
  t >= 0 ==> a * pow a t = pow a (t + 1)
[@@by intros @> [%use pow_succ a t] @> auto]

(* binomial coefficients via Pascal recursion *)

let rec choose (n:int) (k:int) : int =
  if k < 0 || k > n then 0
  else if k = 0 || k = n then 1
  else choose (n-1) (k-1) + choose (n-1) k
[@@measure (Ordinal.of_int (max 0 n + max 0 k))] [@@by auto]

lemma choose_out_of_range n k =
  (k < 0 || k > n) ==> choose n k = 0
[@@by intros @> auto]

lemma choose_edges n =
  n >= 0 ==> (choose n 0 = 1) && (choose n n = 1)
[@@by auto]

lemma pascal n k =
  (0 < k && k < n) ==> choose n k = choose (n-1) (k-1) + choose (n-1) k
[@@by intros @> simplify () @> auto]

(* binomial sum up to index k (descending) *)

let rec binom_sum (n:int) (k:int) (x:int) (y:int) : int =
  if k < 0 then 0
  else choose n k * pow x (n - k) * pow y k + binom_sum n (k - 1) x y
[@@measure (Ordinal.of_int (max 0 k))] [@@by auto]

(* full sum is binom_sum n n x y *)
let binom_eval (n:int) (x:int) (y:int) : int = binom_sum n n x y

(* shift/combiner lemma: for all k>=0,
     x * binom_sum n k  +  y * binom_sum n (k-1)  =  binom_sum (n+1) k
   This is the key algebra behind the Pascal identity in the binomial step. *)
lemma binom_sum_step n k x y =
  k >= 0 && n >= 0 
   ==> x * binom_sum n k x y + y * binom_sum n (k - 1) x y
             = binom_sum (n + 1) k x y
[@@by
  induction ()
  @>| [
    (* base: k = 0 *)
    intros
    @> simplify ()                          (* expands all 3 binom_sum occurrences at k=0 *)
    (* we need: x*C(n,0) x^n + y*C(n,-1) x^(n-(-1)) y^(-1) (but choose n (-1)=0) *)
    @> [%use choose_out_of_range n (-1)]
    @> [%use choose_edges n]
    @> [%use pow_succ_l x (n - 0)]          (* x * x^(n-0) = x^(n - 0 + 1) *)
    @> auto;

    (* step: assume for k-1, prove for k *)
    intros
    (* separate the head terms and the recursive tails via IH *)
    (* head terms want Pascal and the exponent shifts: *)
    @> [%cases (0 < k && k < n + 1)]
    @>| [
      (* interior case: 0<k<n+1, so 1≤k≤n *)
      intros
      @> [%use pascal (n + 1) k]            (* C(n+1,k)=C(n,k-1)+C(n,k) *)
      @> [%use pow_succ_l x (n - k + 1)]    (* x * x^(n-k+1) = x^(n-k+2) == x^((n+1)-k) *)
      @> [%use pow_succ_l y (k - 1)]        (* y * y^(k-1) = y^k *)
      @> auto;

      (* boundary case(s): k=0 or k=n+1. Both are handled by the unfolded definitions:
         - if k=0, we are the base (proved above);
         - if k=n+1, choose(n,k)=0; the identities reduce using choose_out_of_range. *)
      intros
      @> [%use choose_out_of_range n k]
      @> [%use choose_out_of_range n (k - 1)]
      @> [%use choose_out_of_range (n + 1) k]
      @> auto
    ]
  ]
]

lemma binom_sum_extend_top n x y =
  binom_sum n (n + 1) x y = binom_sum n n x y
[@@by
  simplify ()                             (* unfold binom_sum at k = n+1 *)
  @> [%use choose_out_of_range n (n + 1)] (* choose n (n+1) = 0 *)
  @> auto
]

lemma pred_nonneg_of_pos n =
  (not (n <= 0)) ==> (n - 1) >= 0
[@@by intros @> auto]

lemma mul_add_distr_left x y t =
  ((x + y) * t) = (x * t + y * t)
[@@by auto]

(* Binomial Theorem! *)

theorem binomial_theorem n x y =
  n >= 0 && x >= 0 && y >= 0 
  ==> pow (x + y) n = binom_eval n x y
[@@by
  induction ()
  @>| [ auto;

  (* Inductive step goal:
   ((n >= 0) && x >= 0 && y >= 0) ==> pow (x+y) n = binom_sum n n x y
   Hyp: not (n <= 0)  and  IH: ((n-1) >= 0 && x >= 0 && y >= 0) ==> pow(x+y)(n-1) = binom_sum(n-1)(n-1) x y *)
  intros
  (* 1) ensure (n-1) >= 0 so we can use IH *)
  @> [%use pred_nonneg_of_pos n]

  (* 2) expand power at n via n-1 *)
  @> [%use pow_succ (x + y) (n - 1)]          (* pow(x+y) n = (x+y) * pow(x+y)(n-1) *)

  (* 3) distribute (x+y) over the sum *)
  @> [%use mul_add_distr_left x y (binom_sum (n - 1) (n - 1) x y)]

  (* 4) bump the first partial sum to k = n using extend-top *)
  @> [%use binom_sum_extend_top (n - 1) x y]   (* binom_sum(n-1,n) = binom_sum(n-1,n-1) *)

  (* 5) combine via the combiner at k = n: x*S(n,n) + y*S(n,n-1) = S(n+1,n) *)
  @> [%use binom_sum_step (n - 1) n x y]       (* gives binom_sum n n x y *)

  (* 6) apply IH and boom! *)
  @> auto
 ]
]