(* Demo: A formally verified implementation of binary search *)

Require Import NArith.
Open Scope positive_scope.

(* binary search: Find an index i in sorted array A[b..(b+len-1)] satisfying A[i]=x. *)
Fixpoint binsearch (A: positive -> positive) b len x :=
  match len with
  | 1 => b (* array length is 1, so return the only element *)
  | n~1 => (* array length is odd; split at (n-1)/2 *)
      match A(b+n) ?= x with
      | Eq => b+n
      | Lt => binsearch A (b+n+1) n x
      | Gt => binsearch A b n x
      end
  | n~0 => (* array length is even; split at n/2 *)
      match A(b+n) ?= x with
      | Eq => b+n
      | Lt => binsearch A (b+n) n x
      | Gt => binsearch A b n x
      end
  end.

Compute binsearch (fun n => n*2) 1 100 2.
Compute binsearch (fun n => n*2) 1 100 14.
Compute binsearch (fun n => n*2) 1 100 84.
Compute binsearch (fun n => n*2) 50 100 200.

Definition sorted A : Prop :=
  forall i j, (i <= j) -> (A(i) <= A(j)).

Definition member (x:positive) A b len : Prop :=
  exists i, (b <= i < b+len) /\ A(i)=x.

Theorem binsearch_correct:
  forall A x b len,
    sorted A /\ member x A b len -> A(binsearch A b len x) = x.
Proof.
  intros.
  assert (SRT: forall i j, A i < A j -> i < j).
    intros. apply Pos.lt_nle. intro. revert H0. apply Pos.le_nlt, H, H1.
  apply proj2 in H. revert b H. induction len; intros; destruct H as [j [[LO HI] EQ]]; simpl;
  only 1-2: (destruct (_ ?= _) eqn:CMP; [apply Pos.compare_eq|apply IHlen; exists j; repeat split ..]); try assumption.
    rewrite Pos.add_1_r. apply Pos.le_succ_l, SRT. apply -> Pos.compare_lt_iff. rewrite EQ. assumption.
    rewrite <- !Pos.add_assoc, (Pos.add_comm 1), (Pos.add_assoc len), Pos.add_1_r, Pos.add_diag, <- Pos.xI_succ_xO. assumption.
    apply SRT. apply -> Pos.compare_gt_iff. rewrite EQ. assumption.
    apply Pos.lt_le_incl, SRT. apply -> Pos.compare_lt_iff. rewrite EQ. assumption.
    rewrite <- Pos.add_assoc, Pos.add_diag. assumption.
    apply SRT. apply Pos.compare_gt_iff. rewrite EQ. assumption.
    replace b with j. assumption. apply Pos.le_antisym.
      apply Pos.lt_succ_r. rewrite <- (Pos.add_1_r b). assumption.
      assumption.
Qed.