Monday, 18 June 2012

Taylor series expansion of cosine that is generic over number type


Taylor series expansion of cos(x) about x=0:
cosine(x) = 1 - x²/2! + x⁴/4! - x⁶/6! + ...
Let's start by writing a float-specific version in F#:
let cosine n x =
  let rec loop i q t c =
    if i=n then c else
      loop (i + 1) (q + 10 + 8*i) (-t * x * x / float q) (c + t)
  loop 0 2 1.0 0.0
For example, we can compute 1M terms of the expansion of x=0.1:
cosine 1000000 0.1
The best way to make this polymorphic in F# is to parameterize the function over the operators it uses and mark it as inline in order to remove the performance overhead of this parameterization:
let inline cosine zero one ofInt ( ~-. ) ( +. ) ( *. ) ( /. ) n x =
  let rec loop i q t c =
    if i=n then c else
      loop (i + 1) (q + 10 + 8*i) (-.t *. x *. x /. ofInt q) (c +. t)
  loop 0 2 one zero
Now we can compute 1M terms using float like this, which is just as fast as before:
cosine 0.0 1.0 float ( ~- ) (+) (*) (/) 1000000 0.1
But we can also do single-precision float:
cosine 0.0f 1.0f float32 ( ~- ) (+) (*) (/) 1000000 0.1f
And arbitrary-precision rational:
cosine 0N 1N BigNum.FromInt (~-) (+) (*) (/) 10 (1N / 10N)
And even symbolic:
type Expr =
  | Int of int
  | Var of string
  | Add of Expr * Expr
  | Mul of Expr * Expr
  | Pow of Expr * Expr

  static member (~-) f = Mul(Int -1, f)
  static member (+) (f, g) = Add(f, g)
  static member (*) (f, g) = Mul(f, g)
  static member (/) (f, g) = Mul(f, Pow(g, Int -1))

cosine (Int 0) (Int 1) Int (~-) (+) (*) (/) 3 (Var "x")
To make it faster, hoist the common subexpression -x*x out of loop.

No comments: