(* This program calculates the nth fibonacci number
 * using algorithm 1A: naive binary recursion
 *
 * compiled: ocamlopt -ccopt -march=native -o f1a f1a.ml
 * executed: ./f1a n
 *)

(* Naive binary recursion: F(n) = F(n-1) + F(n-2) *)
let rec fib = function
  | 0 -> 0
  | 1 -> 1
  | n -> fib (n-1) + fib (n-2)

(* Function f(n) handles the negative arguments: F(-n) = F(n)*(-1)^(n+1) *)
let f n =
  (if n<0 && -n mod 2=0 then fun n -> -n else fun n -> n)
  (fib (abs n))

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = 
   try int_of_string Sys.argv.(1)
   with _ -> invalid_arg "Usage: f1a <n>"  in
(* calculate f(n) *)
let answer = f n in
(* output *)
  Printf.printf "%dth Fibonacci number is %d\n" n answer

(* This program calculates the nth fibonacci number
 * using algorithm 1B: cached binary recursion
 *
 * compiled: ocamlopt -ccopt -march=native -o f1b nums.cmxa f1b.ml
 * executed: ./f1b n
 *)

open Num

(* same double recursion as in algorith 1A, but reads from cache
 * if there is a value stored already
 *)
let fib n =
  let an = abs n in
  let cache = Array.make (an+1) (Int 0) in
  let rec f = function
          | 0 -> Int 0
          | 1 -> Int 1
          | n -> if cache.(n) =/ Int 0
                 then cache.(n) <- f (n-1) +/ f (n-2);
                 cache.(n) in
  (if n<0 && -n mod 2=0 then minus_num else fun n -> n)
  (f an)

(* entry point *)
let _ =
  (* get n from command line *)
   let n =
     try int_of_string Sys.argv.(1)
     with _ -> invalid_arg "Usage: f1b <n>"  in

  (* calculate f(n) *)
  let answer = fib n in

  (* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer)

(* This program calculates the nth fibonacci number
 * using alrogirhtm 2A: cached linear recursion (as lazy infinite list)
 *
 * compiled: ocamlopt -ccopt -march=native nums.cmxa -o f2a f2a.ml
 * executed: ./f2a n
 * 
 *)
open Num
open Lazy

(* The lazy-evaluated list is head of the list and a promise of the tail.
 * Note that there is no Empty option, meaning this list is infinite *)
type 'a inf_list = Cons of 'a * 'a inf_list lazy_t

(* tail, map2, and nth for this list type *)
let tl l = match l with Cons (_, t) -> force t

let rec map2 f l l' =
    match l, l' with
    Cons (h, t), Cons (h', t')
      -> Cons (f h h', lazy (map2 f (force t) (force t')))

let rec nth lst n = match (lst, n) with
      | (_, n) when n < 0 -> invalid_arg "negative index in nth"
      | (Cons(x, _), 0) -> x
      | (Cons(_, t), n) -> nth (force t) (n - 1)
 
(* Define the infinite lazy-evaluated list of fibonacci numbers *)
let rec fibs = Cons (Int 0, lazy
                (Cons (Int 1, lazy
                  (map2 ( +/ ) fibs (tl fibs)))))

(* Replacing this definition with a more direct:
 *
 *  let rec fibs = let rec aux n m =
 *     Cons (n, lazy (aux m (n +/ m))) in lazy (aux (Int 0) (Int 1))
 *
 * did not change the speed or improve memory usage
 * 
 * I've also tried the new camlp5 functional streams, like this:
 *
 * value fibs =
 *    let rec f a b = fstream [: `a; f b (Num.add_num a b) :]
 *       in f (Int 0) (Int 1);
 *
 * value rec nth s =
 *   let next = Option.get (Fstream.next s) in fun
 *     [ 0          -> fst next
 *     | n when n>0 -> nth (snd next) (n-1)
 *     | _ -> invalid_arg "negative argument" ];
 *
 * but this was considerably slower and also had exponential memory use
 *)

(* entry point *)
let _ =
  (* get n from command line *)
   let n =
     try int_of_string Sys.argv.(1)
     with _ -> invalid_arg "Usage: f2a <n>"  in

  (* calculate f(n), taking care of negation *)
  let answer =
      (if n<0 && -n mod 2=0 then minus_num else fun n -> n)
      (nth fibs (abs n)) in

  (* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer);

(* This program calculates the nth fibonacci number
 * using algorithm 2B: linear recursion with accumulator
 *
 * compiled: ocamlopt -ccopt -march=native -o f2b nums.cmxa f2b.ml
 * executed: ./f2b n
 *)
open Num

(* linear recursion: the parameters a and b are accumulators,
 * n is the iteration counter, decremented on every iteration.
 * When it reaches 1, the accumulators contain (n-1)th and nth
 * fibonacci numbers
 *)
let rec fib ?(a=Int 0) ?(b=Int 1) n = match n with
    | 0 -> a
    | 1 -> b
    | n -> fib ~a:b ~b:(a +/ b) (n-1)

(* Function f(n) handles the negative arguments: F(-n) = F(n)*(-1)^(n+1) *)
let f n =
    (if n<0 && -n mod 2=0 then minus_num else fun n -> n)
    (fib (abs n))

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = try int_of_string Sys.argv.(1)
        with _ -> invalid_arg "Usage: f2b <n>"  in
(* calculate f(n) *)
let answer = f n in
(* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer)

(* This program calculates the nth fibonacci number
 * using algorithm 2C: imperative loop with variables
 *
 * compiled: ocamlopt -ccopt -march=native -o f2c nums.cmxa f2c.ml
 * executed: ./f2c n
 *)
open Num

(* Repeat the loop abs(n) times using mutable variables as accumulators *)
let f n =
  let tmp = ref (Int 0) in
  let x1 = ref (Int 0) in
  let x2 = ref (Int 1) in
  for i = 1 to abs n do
    tmp := !x1 +/ !x2;
    x1 := !x2;
    x2 := !tmp;
  done;
  if n < 0 && (-n mod 2) = 0 then minus_num !x1 else !x1

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = try int_of_string Sys.argv.(1)
        with _ -> invalid_arg "Usage: f2c <n>"  in
(* calculate f(n) *)
let answer = f n in
(* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer)

(* This program calculates the nth fibonacci number
 * using algorithm 3A: matrix multiplication
 *
 * compiled: ocamlopt -ccopt -march=native -o f3a nums.cmxa f3a.ml
 * executed: ./f3a n
 *)
open Num
open List

(* Ocaml does not come with matrix multiplication, and I didn't feel like using
 * PsiLAB or another uncommon library just for that, so wrote my own. A matrix
 * in basic ocaml can be a list of lists, an array of arrays, or a bigarray
 * (especially to pass them to C/Fortran libraries) I went with a list of lists,
 * for purity's sake. *)

(* general inner product of two lists, from Harrop *)
let inner base f l1 l2 g =
    fold_left2 (fun accu e1 e2 -> g accu (f e1 e2)) base l1 l2

(* dot product of two lists of Num *)
let dot a b =
    inner (Int 0) ( */ ) a b ( +/ )

(* matrix transposition, from rosettacode.org *)
let rec transpose m =
    if mem [] m then []
    else map hd m :: transpose (map tl m)

(* matrix multiplication for matrices of Num list list *)
let matmult ml mr =
  let m2 = transpose mr in
  map (fun row -> map (fun col -> dot row col) m2) ml

(* general function to raise an object x to an integer power n using the
 * repeated squaring algorithm tail-recursively.
 * additional arguments:
 * fmul: function to multiply two objects of the same type as x
 * one: unit value for this multiplication (1, 1., unit matrix)
 *)
let rec fast_power fmul one x n = match n with
  | n when n < 0 -> invalid_arg "fast_power requires non-negative n"
  | 0 -> one (* "one" is used as product accumulator *)
  | n -> let sq = fmul x x in (* x contains the last square *)
  if n mod 2=1 then fast_power fmul (fmul one x) sq (n/2)
               else fast_power fmul       one    sq (n/2)

(* Now the fibonacci function can be defined. It will raise the (1,1)(1,0)
 * matrix to n-1'st power *)
let rec f = function
  | 0 -> Int 0
  | n -> let m = [[Int 1; Int 1;];
                  [Int 1; Int 0 ]] in
         let one = [[Int 1; Int 0;];
                    [Int 0; Int 1 ]] in
         (if n<0 && -n mod 2=0 then minus_num else fun n -> n)
         (hd (hd (fast_power matmult one m ((abs n)-1))))

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = try int_of_string Sys.argv.(1)
        with _ -> invalid_arg "Usage: f3a <n>"  in
(* calculate f(n) *)
let answer = f n in
(* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer)

(* This program calculates the nth fibonacci number
 * using algorithm 3B: fast recursion
 *
 * compiled: ocamlopt -ccopt -march=native -o f3b nums.cmxa f3b.ml
 * executed: ./f3b n
 *)

open Num
(* calculate a pair of fibonacci numbers according to
 * the recurrent relationship: 
 * F(2n-1) = F(n-1)^2 + F(n)^2
 * F(2n) = (2F(n-1) + F(n))F(n)
 *
 * in this function, k1 = F(n-1), k2 = F(n), k3 = F(n+1)
 *)
let rec f = function
   | n when n<0 -> invalid_arg "negative argument"
   | 0 -> (Int 0, Int 0)
   | 1 -> (Int 1, Int 0)
   | n -> let (k2, k1) = f (n/2) in
          let k3  = k2 +/ k1 in
          let k2s = k2 */ k2 in
          let k1s = k1 */ k1 in
          let k3s = k3 */ k3 in
          if (n mod 2) = 0 then (k1+/k3)*/k2, k1s+/k2s
                           else k3s+/k2s, (k1+/k3)*/k2

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = try int_of_string Sys.argv.(1)
        with _ -> invalid_arg "Usage: f3b <n>"  in
(* calculate f(n) *)
let answer =
    (if n<0 && -n mod 2=0 then minus_num else fun n -> n)
    (fst (f (abs n))) in
(* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n (string_of_num answer)

(* This program calculates the nth fibonacci number
 * using algorithm 3C: Binet's formula with rounding
 *
 * note: for arbitrary precision floating point numbers, I use
 * ML GMP by David Monniaux, an interface to GNU MP library
 *
 * compiled: ocamlopt -I /usr/lib/ocaml/gmp -ccopt -match=native -o f3c
 *           /usr/lib/ocaml/gmp/gmp.cmxa f3c.ml
 * executed: ./f3c n
 *)
open Gmp
let ( +. ), ( *. ), ( -. ) = F.add, F.mul, F.sub

(* find the inverse root of a using the following recurrent equation:
 * y(n+1) = 0.5*y(n)*[3 - a*y(n)^2]
 * It is faster than typical Newton-Raphson equations because it does not use
 * division.
 *)
let invsqrt a =
  let start = F.from_float (1. /. sqrt a ) in
  let f x = F.from_float 0.5 *. x *. (F.from_float 3. -. F.from_float a *. x *. x) in
  let rec loop f x =
      let x' = f x in
      if F.eq x x' !F.default_prec then x
      else loop f x'
  in
  loop f start

(* general function to raise an object to an integer power using the
 * repeated squaring algorithm tail-recursively.
 * additional arguments:
 * fmul: function to multiply two objects of the same type as x
 * one: unit value for this multiplication (1, 1., unit matrix)
 *)
let rec fast_power fmul one x n = match n with
   | n when n < 0 -> invalid_arg "fast_power requires non-negative n"
   | 0 -> one (* "one" is used as product accumulator *)
   | n -> let sq = fmul x x in (* x contains the last square *)
          if n mod 2=1 then fast_power fmul (fmul one x) sq (n/2)
                       else fast_power fmul       one    sq (n/2)

(*
 * This function calculates the nth fibonacci number
 * as floor( phi^n/sqrt(5) + 1/2)
 *)
let fib n =
    F.default_prec := 7*n/10; (* actually n * log_2(phi) = 0.694n *)
    let one = F.from_float 1. in
    let half = F.from_float 0.5 in
    let is5 = invsqrt 5. in
    let s5 = F.div one is5 in
    let phi = (one +. s5) *. half in
    F.floor ( half +. is5 *. fast_power F.mul one phi n )

(* negative argument handling:  F(-n) = F(n)*(-1)^(n+1) *)
let f n =
  (if n<0 && -n mod 2=0 then F.neg else fun n -> n)
  (fib (abs n))

(* entry point *)
let _ =
(* read the integer argument from command line *)
let n = try int_of_string Sys.argv.(1)
        with _ -> invalid_arg "Usage: f3c <n>"  in
(* calculate f(n), note that for zero, F.to_string... returns empty string *)
let answer = F.to_string_exp_base_digits ~base:10 ~digits:n (f n) in
let answer = if n=0 then "0" else fst answer in
(* output *)
  Printf.printf "%dth Fibonacci number is %s\n" n answer