Bisection_Method_Code.sml

structure Product_Type : sig
  val apsnd : ('a -> 'b) -> 'c * 'a -> 'c * 'b
  val fst : 'a * 'b -> 'a
  val snd : 'a * 'b -> 'b
end = struct

fun apsnd f (x, y) = (x, f y);

fun fst (x1, x2) = x1;

fun snd (x1, x2) = x2;

end; (*struct Product_Type*)

structure Arith : sig
  datatype int = Int_of_integer of IntInf.int
  type nat
  datatype num = One | Bit0 of num | Bit1 of num
  val one_int : int
  val one_nat : nat
  val less_int : int -> int -> bool
  val plus_int : int -> int -> int
  val zero_int : int
  val plus_nat : nat -> nat -> nat
  val zero_nat : nat
  val equal_int : int -> int -> bool
  val minus_int : int -> int -> int
  val times_int : int -> int -> int
  val uminus_int : int -> int
  val divide_int : int -> int -> int
  val modulo_integer : IntInf.int -> IntInf.int -> IntInf.int
end = struct

datatype int = Int_of_integer of IntInf.int;

datatype nat = Nat of IntInf.int;

datatype num = One | Bit0 of num | Bit1 of num;

val one_int : int = Int_of_integer (1 : IntInf.int);

val one_nat : nat = Nat (1 : IntInf.int);

fun integer_of_int (Int_of_integer k) = k;

fun less_int k l = IntInf.< (integer_of_int k, integer_of_int l);

fun plus_int k l =
  Int_of_integer (IntInf.+ (integer_of_int k, integer_of_int l));

val zero_int : int = Int_of_integer (0 : IntInf.int);

fun integer_of_nat (Nat x) = x;

fun plus_nat m n = Nat (IntInf.+ (integer_of_nat m, integer_of_nat n));

val zero_nat : nat = Nat (0 : IntInf.int);

fun divmod_integer k l =
  (if ((k : IntInf.int) = (0 : IntInf.int))
    then ((0 : IntInf.int), (0 : IntInf.int))
    else (if IntInf.< ((0 : IntInf.int), l)
           then (if IntInf.< ((0 : IntInf.int), k)
                  then IntInf.divMod (IntInf.abs k, IntInf.abs l)
                  else let
                         val (r, s) =
                           IntInf.divMod (IntInf.abs k, IntInf.abs l);
                       in
                         (if ((s : IntInf.int) = (0 : IntInf.int))
                           then (IntInf.~ r, (0 : IntInf.int))
                           else (IntInf.- (IntInf.~ r, (1 : IntInf.int)),
                                  IntInf.- (l, s)))
                       end)
           else (if ((l : IntInf.int) = (0 : IntInf.int))
                  then ((0 : IntInf.int), k)
                  else Product_Type.apsnd IntInf.~
                         (if IntInf.< (k, (0 : IntInf.int))
                           then IntInf.divMod (IntInf.abs k, IntInf.abs l)
                           else let
                                  val (r, s) =
                                    IntInf.divMod (IntInf.abs k, IntInf.abs l);
                                in
                                  (if ((s : IntInf.int) = (0 : IntInf.int))
                                    then (IntInf.~ r, (0 : IntInf.int))
                                    else (IntInf.- (IntInf.~
              r, (1 : IntInf.int)),
   IntInf.- (IntInf.~ l, s)))
                                end))));

fun equal_int k l = (((integer_of_int k) : IntInf.int) = (integer_of_int l));

fun minus_int k l =
  Int_of_integer (IntInf.- (integer_of_int k, integer_of_int l));

fun times_int k l =
  Int_of_integer (IntInf.* (integer_of_int k, integer_of_int l));

fun uminus_int k = Int_of_integer (IntInf.~ (integer_of_int k));

fun divide_integer k l = Product_Type.fst (divmod_integer k l);

fun divide_int k l =
  Int_of_integer (divide_integer (integer_of_int k) (integer_of_int l));

fun modulo_integer k l = Product_Type.snd (divmod_integer k l);

end; (*struct Arith*)

structure GCD : sig
  val gcd_int : Arith.int -> Arith.int -> Arith.int
end = struct

fun gcd_integer k l =
  IntInf.abs
    (if ((l : IntInf.int) = (0 : IntInf.int)) then k
      else gcd_integer l (Arith.modulo_integer (IntInf.abs k) (IntInf.abs l)));

fun gcd_int (Arith.Int_of_integer x) (Arith.Int_of_integer y) =
  Arith.Int_of_integer (gcd_integer x y);

end; (*struct GCD*)

structure Rat : sig
  type rat
  val of_int : Arith.int -> rat
  val less_rat : rat -> rat -> bool
  val plus_rat : rat -> rat -> rat
  val zero_rat : rat
  val minus_rat : rat -> rat -> rat
  val times_rat : rat -> rat -> rat
  val divide_rat : rat -> rat -> rat
end = struct

datatype rat = Frct of (Arith.int * Arith.int);

fun of_int a = Frct (a, Arith.one_int);

fun normalize p =
  (if Arith.less_int Arith.zero_int (Product_Type.snd p)
    then let
           val a = GCD.gcd_int (Product_Type.fst p) (Product_Type.snd p);
         in
           (Arith.divide_int (Product_Type.fst p) a,
             Arith.divide_int (Product_Type.snd p) a)
         end
    else (if Arith.equal_int (Product_Type.snd p) Arith.zero_int
           then (Arith.zero_int, Arith.one_int)
           else let
                  val a =
                    Arith.uminus_int
                      (GCD.gcd_int (Product_Type.fst p) (Product_Type.snd p));
                in
                  (Arith.divide_int (Product_Type.fst p) a,
                    Arith.divide_int (Product_Type.snd p) a)
                end));

fun quotient_of (Frct x) = x;

fun less_rat p q = let
                     val (a, c) = quotient_of p;
                     val (b, d) = quotient_of q;
                   in
                     Arith.less_int (Arith.times_int a d) (Arith.times_int c b)
                   end;

fun plus_rat p q =
  Frct let
         val (a, c) = quotient_of p;
         val (b, d) = quotient_of q;
       in
         normalize
           (Arith.plus_int (Arith.times_int a d) (Arith.times_int b c),
             Arith.times_int c d)
       end;

val zero_rat : rat = Frct (Arith.zero_int, Arith.one_int);

fun minus_rat p q =
  Frct let
         val (a, c) = quotient_of p;
         val (b, d) = quotient_of q;
       in
         normalize
           (Arith.minus_int (Arith.times_int a d) (Arith.times_int b c),
             Arith.times_int c d)
       end;

fun times_rat p q = Frct let
                           val (a, c) = quotient_of p;
                           val (b, d) = quotient_of q;
                         in
                           normalize (Arith.times_int a b, Arith.times_int c d)
                         end;

fun divide_rat p q =
  Frct let
         val (a, c) = quotient_of p;
         val (b, d) = quotient_of q;
       in
         normalize (Arith.times_int a d, Arith.times_int c b)
       end;

end; (*struct Rat*)

structure Set : sig
  type 'a set
end = struct

datatype 'a set = Set of 'a list | Coset of 'a list;

end; (*struct Set*)

structure List : sig
  val map : ('a -> 'b) -> 'a list -> 'b list
end = struct

fun map f [] = []
  | map f (x21 :: x22) = f x21 :: map f x22;

end; (*struct List*)

structure Real : sig
  datatype real = Ratreal of Rat.rat
  val less_real : real -> real -> bool
  val plus_real : real -> real -> real
  val zero_real : real
  val minus_real : real -> real -> real
  val times_real : real -> real -> real
  val divide_real : real -> real -> real
end = struct

datatype real = Ratreal of Rat.rat;

fun less_real (Ratreal x) (Ratreal y) = Rat.less_rat x y;

fun plus_real (Ratreal x) (Ratreal y) = Ratreal (Rat.plus_rat x y);

val zero_real : real = Ratreal Rat.zero_rat;

fun minus_real (Ratreal x) (Ratreal y) = Ratreal (Rat.minus_rat x y);

fun times_real (Ratreal x) (Ratreal y) = Ratreal (Rat.times_rat x y);

fun divide_real (Ratreal x) (Ratreal y) = Ratreal (Rat.divide_rat x y);

end; (*struct Real*)

structure Option : sig
  val map_option : ('a -> 'b) -> 'a option -> 'b option
end = struct

fun map_option f NONE = NONE
  | map_option f (SOME x2) = SOME (f x2);

end; (*struct Option*)

structure ITree_Countable_Nondeterminism : sig
  type pndet
end = struct

datatype pndet = Pndetev_C of Arith.nat;

end; (*struct ITree_Countable_Nondeterminism*)

structure Partial_Inj : sig
  type ('a, 'b) pinj
end = struct

datatype ('a, 'b) pinj = Pinj_of_alist of ('a * 'b) list;

end; (*struct Partial_Inj*)

structure Partial_Fun : sig
  type ('a, 'b) pfun
  val map_pfun : ('a -> 'b) -> ('c, 'a) pfun -> ('c, 'b) pfun
end = struct

datatype ('a, 'b) pfun = Pfun_of_alist of ('a * 'b) list |
  Pfun_of_map of ('a -> 'b option) | Pfun_of_pinj of ('a, 'b) Partial_Inj.pinj |
  Pfun_entries of 'a Set.set * ('a -> 'b);

fun map_pfun f (Pfun_of_map g) = Pfun_of_map (fn x => Option.map_option f (g x))
  | map_pfun f (Pfun_of_alist m) =
    Pfun_of_alist (List.map (fn (k, v) => (k, f v)) m);

end; (*struct Partial_Fun*)

structure Interaction_Trees : sig
  datatype ('a, 'b) itree = Ret of 'b | Sil of ('a, 'b) itree |
    Vis of ('a, ('a, 'b) itree) Partial_Fun.pfun
  val bind_itree : ('a, 'b) itree -> ('b -> ('a, 'c) itree) -> ('a, 'c) itree
  val kcomp_itree :
    ('a -> ('b, 'c) itree) -> ('c -> ('b, 'd) itree) -> 'a -> ('b, 'd) itree
end = struct

datatype ('a, 'b) itree = Ret of 'b | Sil of ('a, 'b) itree |
  Vis of ('a, ('a, 'b) itree) Partial_Fun.pfun;

fun bind_itree (Vis t) k = Vis (Partial_Fun.map_pfun (fn x => bind_itree x k) t)
  | bind_itree (Sil t) k = Sil (bind_itree t k)
  | bind_itree (Ret v) k = k v;

fun kcomp_itree p q = (fn s => bind_itree (p s) q);

end; (*struct Interaction_Trees*)

structure ITree_Iteration : sig
  val iterate :
    ('a -> bool) ->
      ('a -> ('b, 'a) Interaction_Trees.itree) ->
        'a -> ('b, 'a) Interaction_Trees.itree
end = struct

fun iterate b p s =
  (if b s
    then Interaction_Trees.bind_itree (p s)
           (Interaction_Trees.Sil o iterate b p)
    else Interaction_Trees.Ret s);

end; (*struct ITree_Iteration*)

structure ITree_Circus : sig
  val assigns : ('a -> 'b) -> 'a -> ('c, 'b) Interaction_Trees.itree
  val cond_itree :
    ('a -> ('b, 'c) Interaction_Trees.itree) ->
      ('a -> bool) ->
        ('a -> ('b, 'c) Interaction_Trees.itree) ->
          'a -> ('b, 'c) Interaction_Trees.itree
end = struct

fun assigns sigma = (fn s => Interaction_Trees.Ret (sigma s));

fun cond_itree p b q = (fn s => (if b s then p s else q s));

end; (*struct ITree_Circus*)

structure Expressions : sig
  val sexp : ('a -> 'b) -> 'a -> 'b
end = struct

fun sexp x = x;

end; (*struct Expressions*)

structure Lens_Laws : sig
  datatype ('a, 'b, 'c) lens_ext =
    Lens_ext of ('b -> 'a) * ('b -> 'a -> 'b) * 'c
  val lens_get : ('a, 'b, 'c) lens_ext -> 'b -> 'a
  val lens_put : ('a, 'b, 'c) lens_ext -> 'b -> 'a -> 'b
end = struct

datatype ('a, 'b, 'c) lens_ext = Lens_ext of ('b -> 'a) * ('b -> 'a -> 'b) * 'c;

fun lens_get (Lens_ext (lens_get, lens_put, more)) = lens_get;

fun lens_put (Lens_ext (lens_get, lens_put, more)) = lens_put;

end; (*struct Lens_Laws*)

structure Bisection : sig
  type 'a state_ext
  val bisection :
    (Real.real -> Real.real) * (Real.real * (Real.real * Real.real)) ->
      unit state_ext ->
        (ITree_Countable_Nondeterminism.pndet, unit state_ext)
          Interaction_Trees.itree
end = struct

datatype 'a state_ext =
  State_ext of
    Arith.nat * Real.real * Real.real * Real.real * Real.real * Real.real *
      Real.real * Real.real * 'a;

fun fa_v_update fa_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_va fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v,
        more);

fun fa_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = fa_v;

val fa : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (fa_v, (fn sigma => fn u => fa_v_update (fn _ => u) sigma), ());

fun fb_v_update fb_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_va fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v,
        more);

fun fb_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = fb_v;

val fb : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (fb_v, (fn sigma => fn u => fb_v_update (fn _ => u) sigma), ());

fun upper_v_update upper_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_v, lower_v, upper_va upper_v, xmid_v, ymid_v, root_v,
        more);

fun upper_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = upper_v;

val upper : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (upper_v, (fn sigma => fn u => upper_v_update (fn _ => u) sigma), ());

fun lower_v_update lower_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_v, lower_va lower_v, upper_v, xmid_v, ymid_v, root_v,
        more);

fun lower_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = lower_v;

val lower : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (lower_v, (fn sigma => fn u => lower_v_update (fn _ => u) sigma), ());

fun ymid_v_update ymid_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_va ymid_v, root_v,
        more);

fun ymid_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = ymid_v;

val ymid : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (ymid_v, (fn sigma => fn u => ymid_v_update (fn _ => u) sigma), ());

fun xmid_v_update xmid_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_va xmid_v, ymid_v, root_v,
        more);

fun xmid_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = xmid_v;

val xmid : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (xmid_v, (fn sigma => fn u => xmid_v_update (fn _ => u) sigma), ());

fun root_v_update root_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_va root_v,
        more);

fun root_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = root_v;

val root : (Real.real, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (root_v, (fn sigma => fn u => root_v_update (fn _ => u) sigma), ());

fun iter_v_update iter_va
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = State_ext
      (iter_va iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v,
        more);

fun iter_v
  (State_ext
    (iter_v, fa_v, fb_v, lower_v, upper_v, xmid_v, ymid_v, root_v, more))
  = iter_v;

val iter : (Arith.nat, 'a state_ext, unit) Lens_Laws.lens_ext =
  Lens_Laws.Lens_ext
    (iter_v, (fn sigma => fn u => iter_v_update (fn _ => u) sigma), ());

fun bisection x =
  (fn (f, (a, (b, tol))) =>
    Interaction_Trees.kcomp_itree
      (ITree_Circus.assigns
        (fn s =>
          Lens_Laws.lens_put iter s
            (Expressions.sexp (fn _ => Arith.zero_nat) s)))
      (Interaction_Trees.kcomp_itree
        (ITree_Circus.assigns
          (fn s => Lens_Laws.lens_put fa s (Expressions.sexp (fn _ => f a) s)))
        (Interaction_Trees.kcomp_itree
          (ITree_Circus.assigns
            (fn s =>
              Lens_Laws.lens_put fb s (Expressions.sexp (fn _ => f b) s)))
          (Interaction_Trees.kcomp_itree
            (ITree_Circus.assigns
              (fn s =>
                Lens_Laws.lens_put lower s (Expressions.sexp (fn _ => a) s)))
            (Interaction_Trees.kcomp_itree
              (ITree_Circus.assigns
                (fn s =>
                  Lens_Laws.lens_put upper s (Expressions.sexp (fn _ => b) s)))
              (Interaction_Trees.kcomp_itree
                (ITree_Circus.assigns
                  (fn s =>
                    Lens_Laws.lens_put xmid s
                      (Expressions.sexp (Lens_Laws.lens_get lower) s)))
                (Interaction_Trees.kcomp_itree
                  (ITree_Circus.assigns
                    (fn s =>
                      Lens_Laws.lens_put ymid s
                        (Expressions.sexp
                          (fn sa => f (Lens_Laws.lens_get xmid sa)) s)))
                  (Interaction_Trees.kcomp_itree
                    (ITree_Iteration.iterate
                      (Expressions.sexp
                        (fn s =>
                          Real.less_real tol
                            (Real.minus_real (Lens_Laws.lens_get upper s)
                              (Lens_Laws.lens_get lower s))))
                      (Interaction_Trees.kcomp_itree
                        (ITree_Circus.assigns
                          (fn s =>
                            Lens_Laws.lens_put iter s
                              (Expressions.sexp
                                (fn sa =>
                                  Arith.plus_nat (Lens_Laws.lens_get iter sa)
                                    Arith.one_nat)
                                s)))
                        (Interaction_Trees.kcomp_itree
                          (ITree_Circus.assigns
                            (fn s =>
                              Lens_Laws.lens_put xmid s
                                (Expressions.sexp
                                  (fn sa =>
                                    Real.divide_real
                                      (Real.plus_real
(Lens_Laws.lens_get lower sa) (Lens_Laws.lens_get upper sa))
                                      (Real.Ratreal
(Rat.of_int (Arith.Int_of_integer (2 : IntInf.int)))))
                                  s)))
                          (Interaction_Trees.kcomp_itree
                            (ITree_Circus.assigns
                              (fn s =>
                                Lens_Laws.lens_put ymid s
                                  (Expressions.sexp
                                    (fn sa => f (Lens_Laws.lens_get xmid sa))
                                    s)))
                            (ITree_Circus.cond_itree
                              (Interaction_Trees.kcomp_itree
                                (ITree_Circus.assigns
                                  (fn s =>
                                    Lens_Laws.lens_put lower s
                                      (Expressions.sexp
(Lens_Laws.lens_get xmid) s)))
                                (ITree_Circus.assigns
                                  (fn s =>
                                    Lens_Laws.lens_put fa s
                                      (Expressions.sexp
(Lens_Laws.lens_get ymid) s))))
                              (Expressions.sexp
                                (fn s =>
                                  Real.less_real Real.zero_real
                                    (Real.times_real (Lens_Laws.lens_get fa s)
                                      (Lens_Laws.lens_get ymid s))))
                              (Interaction_Trees.kcomp_itree
                                (ITree_Circus.assigns
                                  (fn s =>
                                    Lens_Laws.lens_put upper s
                                      (Expressions.sexp
(Lens_Laws.lens_get xmid) s)))
                                (ITree_Circus.assigns
                                  (fn s =>
                                    Lens_Laws.lens_put fb s
                                      (Expressions.sexp
(Lens_Laws.lens_get ymid) s)))))))))
                    (ITree_Circus.assigns
                      (fn s =>
                        Lens_Laws.lens_put root s
                          (Expressions.sexp
                            (fn sa =>
                              Real.divide_real
                                (Real.plus_real (Lens_Laws.lens_get lower sa)
                                  (Lens_Laws.lens_get upper sa))
                                (Real.Ratreal
                                  (Rat.of_int
                                    (Arith.Int_of_integer (2 : IntInf.int)))))
                            s)))))))))))
    x;

end; (*struct Bisection*)