// Examples from chapter 4 "Numerical Analysis" of the book
// F# for Scientists
//
// (C) Flying Frog Consultancy Ltd., 2007

#light

#r "FSharp.Powerpack.dll"

module NumericalAnalysis =

  /// 4.1
  module NumberRepresentation =

    /// 4.1.1
    module MachinePrecisionIntegers =
      // These definitions are now only kept in F# for backwards
      // compatibility with OCaml:
      min_int, max_int

      let binary_of_int n =
        [ for i in Sys.word_size - 1 .. -1 .. 0 ->
            if (n >>> i) % 2 = 0 then "0" else "1" ]
        |> String.concat ""

      binary_of_int 11
      binary_of_int(-11)

    /// 4.1.2
    module MachinePrecisionFloatingPointNumbers =
      6.02214e23

      fsi.AddPrinter(sprintf "%0.20g")

      min_float, max_float

      1.0 / 0.0

      log -1.0

      nan <= 3.0
      nan >= 3.0
      nan <> nan

      set [nan; nan; nan]

      log(complex -1.0 0.0)

      epsilon_float

      1.0 + epsilon_float

      1.0 + sqrt epsilon_float

      1.0 / 3.0

      1.0 - 0.9

  /// 4.2
  module Algebra =
    (0.1 + 0.2) + 0.3 = 0.1 + (0.2 + 0.3)

    (0.1 + 0.2) + 0.3, 0.1 + (0.2 + 0.3)

    1.3 + 1e15 - 1e15

    let f_1 x =
      sqrt(1.0 + x) - 1.0

    f_1 1e-15

    let f_2 x =
      x / (1.0 + sqrt(1.0 + x))

    f_2 1e-15

  /// 4.3
  module Interpolation =
    
    [0.1 .. 0.2 .. 0.9]

    [0.1 .. 0.16 .. 0.9]

    let interp x1 xn n =
      [ for i in 1 .. n ->
          x1 + float(i - 1) * (xn - x1) / float(n - 1) ]

    interp 0.1 0.9 6

    [ for i in 0 .. 5 ->
        0.1 + float i * 0.8 / 5.0 ]

  /// 4.4
  module QuadraticSolutions =
    let quadratic a b c =
      let y = sqrt(b * b - 4. * a * c) in
      (-b + y) / (2.0 * a), (-b - y) / (2.0 * a)

    quadratic 1.0 1e9 1.0

    module Robust =
      let quadratic a b c =
        let y = sqrt(b * b - 4. * a * c) in
        let x1 = if b < 0.0 then y - b else -(y + b) / (2.0 * a) in
        x1, c / (x1 * a)

      quadratic 1.0 1e9 1.0

  /// 4.5
  module MeanAndVariance =

    let mean list =
      let sum = List.fold_left ( + ) 0. list in
      sum / float(List.length list)

    mean [1.; 3.; 5.; 7.]

    let variance_aux (m, s, k) x =
      let m' = m + (x - m) / k in
      m', s + (x - m) * (x - m'), k + 1.0

    let variance xs =
      let _, s, n2 = Seq.fold variance_aux (0.0, 0.0, 1.0) xs in
      s / (n2 - 2.0)

    variance [|1.; 3.; 5.; 7.|]

    let standard_deviation x = sqrt(variance x)

  /// 4.6
  module OtherFormsOfArithmetic =
    
    open Math

    /// 4.6.1
    module ArbitraryPrecisionIntegerArithmetic =
    
      123I
      
      BigInt.Factorial 33I
    
    /// 4.6.2
    module ArbitraryPrecisionRationalArithmetic =
      
      123N / 456N
      
      let rec factorial = function
        | 0 -> 1N
        | n -> BigNum.FromInt n * factorial(n - 1)
      
      factorial 33
      
      let rec e = function
        | 0 -> 1N
        | n -> 1N / factorial n + e (n-1)
      
      e 17
      
      float(e 17)
      
      exp 1.0