(* Ciaran Dunne
   Newton's method in SML, for polynomials.

   In Newton's method, whilst being more efficient, is more difficult as it
   contains the problem of differentating the function f

   To avoid some hefty lexical analysis of functions (to split them up           into terms, then differentiate based on specific rules), we will use polynomials, and our notation of them will be in the form of a list of coefficients. For example:
 [1.0, 2.0, 3.0] <=> 1 + 2x + 3x^2

 First, we create a function /poly/ to create an actual function for our list of
 coefficients, that we can input a value in and get the corresponding value for
 f(x). *)

fun poly [] x = 0.0
  | poly (hd::tl) x = hd + x*(poly tl x);

(* For efficiency and to avoid using powers of x, I have used
   Horner's method to write the polynomial. So:
   1 + 2x + 3x^2 = 1 + x(2 + x(3))

   If polynomials exists solely as a list of coefficients,
   then to differentiate, we must multiply each element in the list by it's
   corresponding x's power, then shift the list along by one.

   Why do we need the shift?
   [1.0, 2.0, 3.0] => 1 + 2x + 3x^2
   bad_diff [1.0, 2.0, 3.0] => [0.0, 2.0, 6.0] => 2x + 6x^2
   Which is quite obviously not the differentiation of our original function. *)

fun diff L = let fun mult n [] = []
                   | mult n (hd::tl) = (n*hd)::mult (n+1.0) tl
               in tl (mult_elements 0.0 L) end;

(* Now finally to compose our Newton function.
   We make three local declarations f is the function for L, f' the
   differentiated function of f. x is defined recursively by the definition of
   the newton method. *)

fun newton _ x0 0 = x0
  | newton L x0 i = let val x  = newton L x0 (i-1)
            val f  = poly L
                val f' = poly (diff L)
            in x - (f x)/(f' x) end;


newton [~2.0, 0.0, 1.0] 1.6 4;

(* val it = 1.414213562 : real

   Lovely! It converges to a nice approximation of root 2 in just 4 iterations!
   wowee! *)