セグメント木を用いた高速な区間積と、セグメント木上の二分探索の検証

Segtree.v

Require Import Program.
From mathcomp Require Import all_ssreflect.

Section Segtree.
  Variable A : Type.
  Variable idm : A.
  Hypothesis mul : Monoid.law idm.
  Local Notation "x * y" := (mul x y).

  Variable table : nat -> nat -> A.
  Hypothesis table_accumulate : forall n i, table n.+1 i = table n i.*2 * table n i.*2.+1.

  Lemma upper_table_product n l : forall r,
    \big[mul/idm]_(l <= i < r) table n.+1 i = \big[mul/idm]_(l.*2 <= i < r.*2) table n i.
  Proof.
    elim => [ | r IHr ].
    - by rewrite double0 !big_geq // leq0n.
    - case (leqP l r) => H.
      + by rewrite doubleS big_nat_recr // IHr table_accumulate Monoid.mulmA !big_nat_recr ?leq_Sdouble ?leq_double.
      + by rewrite !big_geq // leq_double.
  Qed.

  Lemma double_uphalf n : (uphalf n).*2 = odd n + n.
  Proof.
    by rewrite uphalf_half doubleD -(addnn (odd n)) -addnA odd_double_half.
  Qed.

  Lemma double_half n : n./2.*2 = n - odd n.
  Proof.
    by rewrite -{2}(odd_double_half n) addKn.
  Qed.

  Lemma left_segment n l :
    (if odd l then table n l else idm) = \big[mul/idm]_(l <= i < odd l + l) table n i.
  Proof.
    case (odd l).
    - by rewrite big_nat1.
    - by rewrite big_geq.
  Qed.

  Corollary left_segment_acc n l lp :
    (if odd l then lp * table n l else lp) = lp * \big[mul/idm]_(l <= i < odd l + l) table n i.
  Proof.
    move: (left_segment n l) => /(f_equal (mul lp)).
    by rewrite fun_if Monoid.mulm1.
  Qed.

  Lemma right_segment n r :
    (if odd r then table n r.-1 else idm) = \big[mul/idm]_(r - odd r <= i < r) table n i.
  Proof.
    case (@idP (odd r)) => [ /odd_gt0 /prednK {3}<- | ].
    - by rewrite subn1 big_nat1.
    - by rewrite big_geq // subn0.
  Qed.

  Corollary right_segment_acc n r rp :
    (if odd r then table n r.-1 * rp else rp) = \big[mul/idm]_(r - odd r <= i < r) table n i * rp.
  Proof.
    move: (right_segment n r) => /(f_equal (mul^~rp)).
    by rewrite (fun_if (mul^~rp)) Monoid.mul1m.
  Qed.

  Lemma leq_double' l r : (l.*2 <= r) = (l <= r./2).
  Proof.
    rewrite -{1}(odd_double_half r).
    case (odd r).
    - by rewrite leq_Sdouble.
    - by rewrite leq_double.
  Qed.

  (*
    let rec product_rec r n l lp rp =
      if r <= l
      then lp * rp
      else
        product_rec (r / 2) (n + 1) ((l + 1) / 2)
          (if l mod 2 = 0 then lp else lp * table n l)
          (if r mod 2 = 0 then rp else table n (r - 1) * rp)
   *)
  Program Definition product_rec :=
    Fix Wf_nat.lt_wf
      (fun r => forall n l lp rp,
        { p | p = lp * (\big[mul/idm]_(l <= i < r) table n i) * rp })
      (fun r product_rec n l lp rp =>
        ( if r <= l as b return (r <= l) = b -> _
          then fun Hle => exist _ (lp * rp) _
          else fun Hgt =>
            let Hdecr := _ in
            let (p, _) :=
              product_rec r./2 Hdecr n.+1 (uphalf l)
                (if odd l then lp * table n l else lp)
                (if odd r then table n r.-1 * rp else rp) in
            exist _ p _ ) erefl).
  Next Obligation.
    by rewrite big_geq // Monoid.mulm1.
  Qed.
  Next Obligation.
    rewrite -divn2.
    apply /ltP /ltn_Pdiv => //.
    by move: (posnP r) Hgt (leq0n l) => [ ] // -> ->.
  Qed.
  Next Obligation.
    rewrite upper_table_product double_half double_uphalf left_segment_acc right_segment_acc !Monoid.mulmA.
    congr (_ * _).
    rewrite -!Monoid.mulmA.
    congr (_ * _).
    rewrite -!big_cat_nat ?leq_addl ?leq_subr //.
    - apply /(@leq_trans l.+1).
      + by rewrite (leq_add2r l _ 1) leq_b1.
      + by rewrite ltnNge Hgt.
    - case (@idP (odd r)) => Hoddr.
      + case (@idP (odd l)) => Hoddl.
        { rewrite subn1 ltn_predRL.
          case (leqP l.+2 r) => //.
          rewrite ltnS => Hleq.
          have Hsucc : r = l.+1 by apply /anti_leq; rewrite Hleq ltnNge Hgt.
          by move: Hsucc Hoddl Hoddr => -> /= ->. }
        move: Hoddr Hgt => /odd_gt0 /prednK => {1}<- /negbT.
        by rewrite -ltnNge ltnS subn1.
      + apply /(@leq_trans l.+1).
        * by rewrite (leq_add2r l _ 1) leq_b1.
        * by rewrite subn0 ltnNge Hgt.
  Defined.

  (* let product = product_rec r 0 l idm idm *)
  Program Definition product l r : { p | p = \big[mul/idm]_(l <= i < r) table 0 i } :=
    let (p, _) := product_rec r 0 l idm idm in
    exist _ p _.
  Next Obligation.
    by rewrite Monoid.mulm1 Monoid.mul1m.
  Qed.

  Variable P : pred A.

  (*
    let rec upper_bound_rec r n l lp =
      if r <= l
      then (l, lp)
      else
        let lp' = if l mod 2 = 0 then lp else lp * table n l in
        if l mod 2 = 1 && not (P lp')
        then (l, lp)
        else
          let (m, p) = upper_bound_rec (r / 2) (n + 1) ((l + 1) / 2) lp' in
          if r <= m * 2
          then (m * 2, p)
          else
            let p' = p * table n (m * 2) in
            if P p'
            then (m * 2 + 1, p')
            else (m * 2, p)
   *)
  Program Definition upper_bound_rec :=
    Fix Wf_nat.lt_wf
      (fun r => forall n l lp,
        (exists2 m, l <= m <= r &
          forall k, P (lp * \big[mul/idm]_(l <= i < k) table n i) = (k <= m)) ->
        {'(m, p) | [ /\ l <= m <= r,
          p = lp * \big[mul/idm]_(l <= i < m) table n i &
          forall k, P (lp * \big[mul/idm]_(l <= i < k) table n i) = (k <= m) ]})
      (fun r upper_bound_rec n l lp Hex2 =>
        ( if r <= l as b return (r <= l) = b -> _
          then fun Hle => exist _ (l, lp) _
          else fun Hgt =>
            let lp' := if odd l then lp * table n l else lp in
            ( if odd l && ~~ P lp' as b return odd l && ~~ P lp' = b -> _
              then fun Htrue => exist _ (l, lp) _
              else fun Hfalse =>
                let (pair, _) := upper_bound_rec r./2 _ n.+1 (uphalf l) lp' _ in
                ( let (m, p) as pair0 return pair = pair0 -> _ := pair in fun _ =>
                  ( if r <= m.*2 as b return (r <= m.*2) = b -> _
                    then fun Hge => exist _ (m.*2, p) _
                    else fun Hlt =>
                      let p' := p * table n m.*2 in
                      ( if P p' as b return P p' = b -> _
                        then fun HP => exist _ (m.*2.+1, p') _
                        else fun HnP => exist _ (m.*2, p) _ ) erefl ) erefl ) erefl ) erefl ) erefl).
  Next Obligation.
    move: H => /andP [ Hle1 Hle2 ].
    move: (@anti_leq l r) (Hle2) H0 => <-.
    + move => /(conj Hle1) /andP /anti_leq ->.
      by rewrite big_geq // Monoid.mulm1 !leqnn.
    + by rewrite Hle (@leq_trans Hex2).
  Qed.
  Next Obligation.
    move: H => /andP [ Hle1 ? ].
    rewrite leqnn (@leq_trans Hex2) // big_geq // Monoid.mulm1.
    move: Htrue (H0) => /andP [ -> ].
    by rewrite -(@big_nat1 _ idm mul) H0 -ltnNge ltnS => /(conj Hle1) /andP /anti_leq <-.
  Qed.
  Next Obligation.
    apply /ltP.
    by rewrite -divn2 ltn_Pdiv // (@leq_ltn_trans l) // ltnNge Hgt.
  Qed.
  Next Obligation.
    move: H => /andP [ Hle Hlt ].
    exists Hex2./2 => [ | k ].
    - rewrite !uphalf_half half_leq ?andbT ?ltnS //.
      move: (@idP (odd l)) (negbT Hfalse) => [ Hodd /= | ? ? ].
      + rewrite negbK -{1}(@big_nat1 _ idm mul) H0 => /half_leq.
        by rewrite /= uphalf_half Hodd.
      + exact /half_leq.
    - rewrite !uphalf_half upper_table_product doubleD -addnn -addnA odd_double_half.
      move: (@idP (odd l)) Hfalse => [ Hodd /= /negbT | ? ? ].
      + case (leqP k.*2 l) => [ ? | ? ? ].
        { rewrite negbK big_geq.
          - rewrite Monoid.mulm1 => ->.
            by rewrite -(doubleK k) half_leq // (@leq_trans l).
          - exact /leqW. }
        by rewrite -Monoid.mulmA -(@big_nat1 _ idm mul _ (table n)) -big_cat_nat // H0 leq_double'.
      + by rewrite H0 leq_double'.
  Qed.
  Next Obligation.
    move: (H1) H => /andP [ Hge0 Hle0 ] [ /andP [ ] ].
    rewrite -(leq_double (uphalf l)) -(leq_double m) upper_table_product double_half double_uphalf left_segment_acc => Hge' Hle'.
    rewrite -Monoid.mulmA -big_cat_nat ?leq_addl => // -> /(_ m).
    rewrite upper_table_product double_uphalf -Monoid.mulmA -big_cat_nat ?leq_addl // leqnn => ?.
    rewrite (@anti_leq m.*2 Hex2) // -H2.
    by rewrite {2}(@anti_leq m.*2 r) ?Hle0 ?andbT // (leq_trans Hle' (leq_subr _ _)).
  Qed.
  Next Obligation.
    move: (H1) H HP => /andP [ Hge0 Hle0 ] [ /andP [ ] ].
    rewrite -(leq_double (uphalf l)) -(leq_double m) upper_table_product double_half double_uphalf left_segment_acc => Hge' Hle'.
    rewrite -Monoid.mulmA -big_cat_nat ?leq_addl => // -> /(_ Hex2./2).
    rewrite upper_table_product double_uphalf double_half -!Monoid.mulmA => Hub.
    rewrite -(@big_nat1 _ idm mul) -big_cat_nat ?(leq_trans (leq_addl _ _) Hge') // H2 => Hgt'.
    rewrite (@anti_leq m.*2.+1 Hex2) //.
    move: Hub.
    rewrite -big_cat_nat ?leq_addl //.
    - rewrite H2 leq_subr -{1}leq_double double_half leq_subLR => /(@Logic.eq_sym _ _ _) H.
      by move: (leq_add2r m.*2) (leq_b1 (odd Hex2)) andbT => <- /(leq_trans H) -> ->.
    - move: (Hgt').
      by rewrite -(@leq_subRL 1 m.*2 Hex2) ?(leq_ltn_trans (leq0n _) Hgt') => // /(leq_trans Hge') /((@leq_trans _ _ _)^~(leq_sub2l _ (leq_b1 _))).
  Qed.
  Next Obligation.
    move: H1 H (negbT HnP) => /andP [ Hge0 Hle0 ] [ /andP [ ] ].
    rewrite -(leq_double (uphalf l)) -(leq_double m) upper_table_product double_half double_uphalf left_segment_acc => Hge' Hle'.
    rewrite -Monoid.mulmA -big_cat_nat ?leq_addl => // -> /(_ m).
    rewrite leqnn upper_table_product double_uphalf -!Monoid.mulmA => Hub.
    rewrite -(@big_nat1 _ idm mul) -big_cat_nat ?(leq_trans (leq_addl _ _) Hge') //= H2 -leqNgt => Hle''.
    rewrite (@anti_leq m.*2 Hex2) //.
    move: Hub.
    by rewrite -big_cat_nat ?leq_addl // H2 => ->.
  Qed.

  (* let upper_bound l r = upper_bound_rec r 0 l idm *)
  Program Definition upper_bound l r
    (Hex2: exists2 m, l <= m <= r &
      forall k, P (\big[mul/idm]_(l <= i < k) table 0 i) = (k <= m)):
    {'(m, p) | [ /\ l <= m <= r,
      p = \big[mul/idm]_(l <= i < m) table 0 i &
      forall k, P (\big[mul/idm]_(l <= i < k) table 0 i) = (k <= m) ]} :=
    let (pair, _) := upper_bound_rec r 0 l idm _ in
    exist _ pair _.
  Next Obligation.
    exists Hex2 => // ?.
    by rewrite Monoid.mul1m.
  Qed.
  Next Obligation.
    case H => ? -> Hub.
    split => // [ | ? ].
    - exact /Monoid.mul1m.
    - by rewrite -Hub Monoid.mul1m.
  Qed.
End Segtree.