Tuesday, 8 October 2013

Fun with infinite power series

This blog post was inspired by Doug McIlroy's excellent Power Serious article about processing infinite power series using infinite lazy lists in Haskell. Although Doug's Haskell code is significantly more elegant than our F# in parts, we add the ability to visualize power series as typeset mathematics including some example approximations to well-known constants using our F# for Visualization library.

Lazy lists are not as elegant in F# as they are in Haskell (where lazy evaluation is the norm) and, in fact, do not even appear in the F# standard library. However, rather than using lazy lists of coefficients to represent infinite series we can use functions that map from the integer index to the rational factor. The first factor in an infinite series is the constant and appear at index zero and so on.

Arbitrary-precision rational numbers are provided by the F# PowerPack:

#r "FSharp.PowerPack"

We shall make extensive use of this type so let us define an alias:

type Q = BigRational

We can begin by defining some fundamental infinite series. The first represents the number one by returning one for the factor at index zero and zero for all other factors:

let one = function 0 -> 1N | _ -> 0N

The function f(x)=x is represented by a series with a one at index one and zeroes at every other index:

let x = function 1 -> 1N | _ -> 0N

A constant is represented as follows:

let constant k = function 0 -> k | _ -> 0N

The following function skips the first factor and moves the rest along so rest(x+2x²+3x³)=1+2x+3x² which is equivalent to dividing by x when the first factor was zero:

let rest f n = f(n+1)

The xTimes function moves the factors down one place and adds a zero at the beginning, which is equivalent to multiplying by x:

let xTimes f n = if n=0 then 0N else f(n-1)

Addition and subtraction of sequences are simply element-wise:

let (+.) f g n : Q = f n + g n

let (-.) f g n : Q = f n - g n

Prepending an element onto the beginning of a lazy list is now replaced by a function that returns a given constant at index zero followed by a sequence:

let cons(x, xs) n = if n=0 then x else xs(n-1)

The derivative of a series is remarkably simple:

let differentiate f n = Q.FromInt(n+1) * f(n+1)

Similarly for integration:

let integrate f n =
  if n=0 then 0N else f(n-1) / Q.FromInt n

As we shall see, the ability to integrate and differentiate power series allows us to derive a wide variety of interesting power series.

The following function scales each factor by another factor:

let scale x xs n = x * xs n

The product of two power series can be expressed as follows:

let rec ( *. ) f g n = f 0 * g n + xTimes(rest f *. g) n

The quotient of one power series divided by another also has a remarkably simple representation:

let rec (/.) f g n =
  let rec qs n =
    cons(f 0 / g 0, scale (1N / g 0) (rest f -. qs *. rest g)) n
  qs n

Composition allows us the calculate f(g(x)) where both f and g are power series:

let rec ( >>. ) f g n = cons(f 0, rest g *. (rest f >>. g)) n

Reversion gives the inverse function of the given power series provided it has a zero constant at the beginning and also has a remarkably simple representation:

let revert f =
  let rec rs n = cons(0N, one /. (rest f >>. rs)) n

Now we're ready to compute the power series of some familiar functions. The power series of ln(1+x) is given by the integral of 1/(1+x):

let logOnePlusX = integrate(one /. (x +. one))

The power series for exponential may now be obtained by reverting the power series for logarithm because they are inverse functions of each other:

let exp = revert logOnePlusX +. one

A similar identity holds for arctan:

let atan = integrate(one /. (one +. x *. x))

Reverting this function recovers the power series for tan:

let tan = revert atan

Amazingly, the power series of sine and cosine can be defined in terms of each other using integration:

let rec sin n = integrate cos n
and cos n = (one -. integrate sin) n

Finally, arc sine can be obtained by reverting the series for sine:

let asin = revert sin

Note, however, that cosine cannot be reverted because it has a non-zero constant.

In order to visualize our power series we begin by loading the F# for Visualization library and opening the relevant namespaces:

#r "FSharpForVisualization"

open FlyingFrog.FSharpForVisualization
open FlyingFrog.Graphics.Typeset
open FlyingFrog.Graphics.Typeset.Math

The following function finds the first six non-zero factors out of the first twelve:

let initialTerms series =
  seq { for i in 0..12 -> i, series i }
  |> Seq.filter (fun (_, q) -> q <> 0N)
  |> Seq.truncate 6
  |> List.ofSeq

Although F# allows arbitrary-precision rational literals in expressions it does not allow them in patterns. Therefore, we define the following active patterns to recognise important rational numbers:

let (|One|_|) q = if q = 1N then Some() else None
let (|MinusOne|_|) q = if q = -1N then Some() else None

We shall classify the sign of the factor in order to determine whether it should be preceded by a minus sign:

type Sign = Negative | NonNegative

The sign of a rational is then given by:

let signOf q =
  if q<0N then Negative else NonNegative

The following function typesets a single term in a power series where the factor is q and the exponent is i:

let typesetTerm = function
  | 0, q -> math q
  | 1, (MinusOne | One) -> math "x"
  | 1, q -> Row[math(abs q); math "x"]
  | n, (MinusOne | One) -> Super(math "x", math n)
  | n, q -> Row[math(abs q); Super(math "x", math n)]

Finally, a power series may be typeset by extracting the first few non-zero terms, finding their signs and typesetting each term before combining them with plus or minus signs between them:

let rec typesetSeries series =
  let rec loop = function
    | [] -> math "+ …"
    | (Negative, x)::xs -> Row[math "-"; x; loop xs]
    | (NonNegative, x)::xs -> Row[math "+"; x; loop xs]
  let terms =
    [ for i, q in initialTerms series ->
        signOf q, typesetTerm(i, q) ]
  match terms with
  | [] -> math "…"
  | (Negative, x)::xs -> Row[math "-"; x; loop xs]
 | (NonNegative, x)::xs -> Row[x; loop xs]

As each visualization will be spawned in its own window it is useful to augment each power series with the mathematical notation for the expression that it expresses. The following draw function does this:

let draw series lbl =
  Row[lbl; math "="; typesetSeries series]
  |> view

For example, our power series for logarithm may be visualized as follows:

Row[Text "ln"; math "(1+x)"]
|> draw logOnePlusX

Here are the other power series we calculated:

Super(math "e", math "x")
|> draw exp
Row[Super(Text "tan", math "-1"); math "(x)"]
|> draw atan
Row[Text "tan"; math "(x)"]
|> draw tan
Row[Text "sin"; math "(x)"]
|> draw sin
Row[Text "cos"; math "(x)"]
|> draw cos
Row[Super(Text "sin", math "-1"); math "(x)"]
|> draw asin

We can also evaluate some approximations by summing the first few terms of some of our power series:

let eval n x (a : int -> BigRational) =
  seq { for i in 0..n-1 -> a i * pown x i }
  |> Seq.sum

Again, it is useful to encapsulate the conventional mathematical notation for a value along with its known decimal representation, a calculated rational approximation and its decimal representation:

let displayApprox lbl (exact: float) (q: Q) =
  Row[lbl; math "="; math exact; Symbols.EqualTilde;
      math q; math "="; math(float q)]
  |> view

We can easily calculate a good approximation to ln(1.1):

eval 10 (1N / 10N) logOnePlusX
|> displayApprox (Row[Text "ln"; math "(1.1)"]) (log 1.1)

Machin's formula for π gives us quite good rational approximations:

16N * eval 8 (1N / 5N) atan - 4N * eval 2 (1N / 239N) atan
|> displayApprox Greek.pi System.Math.PI

Evaluating e¹ gives an approximation to e:

displayApprox (math "e") System.Math.E (eval 10 1N exp)

To create visualizations like these within the comfort of Visual Studio, order a copy of our F# for Visualization library.

No comments: