(Code for Lecture 8 *)

(* First some basic fuctions we will use *)
fun cube n = n*n*n 
fun rcube (n:real) = n*n*n 
fun square n = n * n 
fun exp (b,n) = if n = 0 then 1 else b * exp(b, n-1)


(* Here are some well-known sums *)
(* sumInts(a,b) = sum k from k=a to k=b  *)
fun sumInts(a,b) = if (a > b) then 0 else a + sumInts(a+1,b)

(* sumSquare(a,b) = sum k^2 from k=a to k=b  *)
fun sumSquare(a,b) = if (a > b) then 0 else square(a) + sumSquare(a+1,b)

(* sumCubes(a,b) = sum k^3 from k=a to k=b  *)
fun sumCubes(a,b) = if (a > b) then 0 else cube(a) + sumCubes(a+1,b)

(* sumExp(a,b) = sum 2^k from k = a to k = b *)
fun sumExp(a,b) = if (a > b) then 0 else exp(2,a) + sumExp(a+1,b)



(*------------------------------------------------------------------- *)
(* Step 1: 
   We will abstract over the function f (i.e. cube, square, exp etc) 
   to get a general sum function. *)

(* sum: (int -> int) * int * int -> int *)
fun sum(f,a,b) =
    if (a > b) then 0
    else (f a) + sum(f,a+1,b)

(* Now we can get the previous sums as follows *)
fun sumInts'(a,b) = sum(fn x => x, a, b)
fun sumCubes'(a,b) = sum(cube, a, b) 
fun sumSquare'(a,b) = sum(square, a, b) 
fun sumExp'(a,b) = sum(fn x => exp(2,x), a, b) 

(*------------------------------------------------------------------- *)
(* What if we want to abstract over the increment function? *)
(* For example, if we want to sum up all the odd numbers? *)

(* sumOdd: int -> int -> int *)
fun sumOdd (a, b) = 
  let
    fun sum' (a,b) = 
      if a > b then 0
      else a + sum'(a+2, b)
  in  
    if (a mod 2) = 1 then
      (* a was odd *)
      sum'(a,b)
    else 
      (* a was even *)
      sum'(a+1, b)
  end 
    

(* sum': (int -> int) * int * int * (inc -> inc) -> int *)
fun sum' (f, a, b, inc) = 
    if (a > b) then 0
    else (f a) + sum'(f,inc(a),b, inc)

(* Now we can rewrite sumOdd using sum' *)
fun sumOdd' (a,b) = 
  if (a mod 2) = 1 then    
    sum'(fn x => x, a, b, fn x => x + 2)
  else 
    sum'(fn x => x, a+1, b, fn x => x + 2)

(*------------------------------------------------------------------- *)

(* Product in analogy with sum.  *)
fun product(f,lo,hi,inc) =
  if (lo > hi) then 1
  else (f lo) * product(f,inc(lo),hi,inc)

(* Using product to define factorial. *)
fun fact n = product(fn x => x, 1, n, fn x => (x + 1))

(*------------------------------------------------------------------- *)
(* The general series written in the tail-recursive form. *)
(* series: 
     (int * int -> int) *     combiner
     (int -> int) *           f  
     int * int *              lo and hi
     (int -> int) *           inc 
     int                      result 
  -> int 

*)
fun series(comb,f,lo,hi,inc,r) =
  let
    fun series' (lo, r) = 
      if (lo > hi) then r
      else
	series' (inc(lo), comb(r, f(lo)))
  in 
    series'(lo, r)
  end 

fun sumSeries (f,a,b,inc) = series(fn (x,y) => x + y, f, a, b, inc, 0)
fun prodSeries (f,a,b,inc) = series(fn (x,y) => x * y, f, a, b, inc, 1)

(*------------------------------------------------------------------- *)
(* An iterative version of sum modified to deal with real values. *)

fun iter_sum(f, lo:real, hi, inc) =
  let
    fun sum' (lo, r) = 
      if (lo > hi) then r
    else 
      sum'(inc(lo), r + f(lo))
  in 
    sum'(lo, 0.0)
  end


(* The integral of f(x) is the area below the curve f(x) in the interval [a, b]. 
   We will use rectangle method to approximate the integral of f(x) in 
   the interval [a,b], made by summing up a series of small rectangles
   using midpoint approximation.
*)

fun integral(f,lo,hi,dx) =
    dx * iter_sum(f,(lo + (dx / 2.0)), hi, fn x => (x + dx)) 


(* ---------------------------------------------------------- *)

(* halfInterval (f, (a,b) , t) = (a', b') 

Given a function f : real -> real,
      two real numbers a, b s.t. f(a) < 0 and f(b) > 0
  and a real number tolerance t,

 the x for which f(x) = 0 will be between  a' and b',
 and the size of this interval is t.

val halfInterval: (real -> real) * (real * real) * real -> (real * real)
*)
fun abs(x:real) = if (x < 0.0) then ~x else x 

fun halfint(f, a, b, t) =
  let
    val mid = (a + b)/2.0
  in
    if  (abs(f(mid)) < t)
      then mid
    else
      if (f(mid) < 0.0)
        then halfint(f, mid, b, t)
      else halfint(f, a, mid, t)
  end

(*------------------------------------------------------------------- *)
(* Combining higher order functions with recursive datatypes:  *)
(*------------------------------------------------------------------- *)

(* map : ('a -> 'b) -> 'a list -> 'b list 
   
   map f l = l'  where l' is the list we obtain by applying f to
   each element in l

*)

fun map f nil = nil
  | map f (h::t) = (f h)::(map f t)

(* val filter : ('a -> bool) -> 'a list -> 'a list

   filter p l ==> sublist containing the elements of l for which p holds.

   Invariants: none
   Effects: none
*)

fun filter (p : 'a -> bool) (nil : 'a list)  : 'a list = nil
  | filter (p : 'a -> bool) (x::l : 'a list) : 'a list = 
      if p(x) then x::filter p l
      else filter p l;


(* val pos : int list -> int list
       pos(l) ==> sublist of l consisting of strictly positive integers.

   Invariants: none
   Effects: none
*)

val pos : int list -> int list = filter (fn (n:int) => n > 0);

(* Now try: *)

pos [2, ~1, ~4, 3, 0, 8];