mpri-funcprog-project/ocaml/NewtonRaphson.ml
2017-09-22 11:32:50 +02:00

152 lines
3 KiB
OCaml

(* A couple abbreviations. *)
type 'a thunk = 'a Lazy.t
let force = Lazy.force
(* The definition of (finite or infinite) lazy lists. *)
type 'a stream =
'a head thunk
and 'a head =
| Nil
| Cons of 'a * 'a stream
(* Calling [tail xs] demands the head of the stream, that is, forces
the computation of the first element of the stream (if there is one). *)
let tail xs =
match force xs with
| Nil -> assert false
| Cons (_, xs) -> xs
(* Newton-Raphson approximation, following Hughes,
"Why functional programming matters", 1990. *)
let next n x =
(x +. n /. x) /. 2.
(* An infinite stream obtained by iterating [f]. *)
(* The following definition, copied almost literally from Hughes'
paper, is correct but somewhat unsatisfactory; can you see why? Can
you fix it? Think about it before reading the solution below. *)
let rec repeat (f : 'a -> 'a) (a : 'a) : 'a stream =
lazy (Cons (a, repeat f (f a)))
(*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*)
(* In the above definition of [repeat], the function call [f a] takes
place when the *first* element of the list is demanded by the consumer.
That's too early -- ideally, this function call should take place only
when the *second* element is demanded, since the result of [f a] is the
second element in the infinite stream [a], [f a], [f (f a)], ... *)
(* This code exhibits the problem: *)
let () =
let x = ref 0 in
let f () = incr x; () in
let xs = repeat f () in
let xs = tail xs in
(* This assertion fails because [x] has been incremented once: *)
assert (!x = 0);
ignore xs
(* This can be fixed in several ways. One solution is to let [repeat] take an
argument of type ['a thunk] instead of ['a]. This approach is in fact the
standard encoding of call-by-need into call-by-value, applied to Hughes'
code. *)
let rec repeat (f : 'a -> 'a) (a : 'a thunk) : 'a stream =
lazy (
Cons (
force a,
repeat f (lazy (f (force a)))
)
)
(* It can also be written as follows. *)
let rec repeat (f : 'a -> 'a) (a : 'a thunk) : 'a stream =
lazy (
let a = force a in
Cons (
a,
repeat f (lazy (f a))
)
)
(* We define a wrapper so [repeat] has the desired type again: *)
let repeat (f : 'a -> 'a) (a : 'a) : 'a stream =
repeat f (lazy a)
(* The problematic code now behaves as desired: *)
let () =
let x = ref 0 in
let f () = incr x; () in
let xs = repeat f () in
(* These assertions succeed: *)
let xs = tail xs in
assert (!x = 0);
let xs = tail xs in
assert (!x = 1);
let xs = tail xs in
assert (!x = 2);
ignore xs
(* Back to Newton-Rapshon. *)
let rec within (eps : float) (xs : float stream) : float =
match force xs with
| Nil -> assert false
| Cons (a, xs) ->
match force xs with
| Nil -> assert false
| Cons (b, _) ->
if abs_float (a /. b -. 1.) <= eps then b
else within eps xs
let sqrt (n : float) : float =
within 0.0001 (repeat (next n) n)
let sqrt2 =
sqrt 2.