Showing posts with label Maths. Show all posts
Showing posts with label Maths. Show all posts

2012-08-06

Merge Right

The so-called master-transaction update is one of the, if not the defining algorithms of the discipline formerly known as "data processing". Given two sorted files of records with increasing keys, the process applies each record in the transaction file to each record of the the master file and outputs the result, if any, to the updated master file in one pass over each input. The same algorithm can compute the union, intersection or difference of sorted sequences. For instance, the union of two sets represented as sorted lists of unique elements is:

let union       =
  let rec go l r = match l, r with
  | [], xs | xs, []  -> xs
  | x :: xs, y :: ys ->
    match compare x y with
    | -1 -> x :: go       xs (y :: ys)
    |  0 -> x :: go       xs       ys
    |  1 -> y :: go (x :: xs)      ys
    |  _ -> assert false
  in go

Intersection is:

let inter       =
  let rec go l r = match l, r with
  | [], _  | _, []   -> []
  | x :: xs, y :: ys ->
    match compare x y with
    | -1 ->      go       xs (y :: ys)
    |  0 -> x :: go       xs       ys
    |  1 ->      go (x :: xs)      ys
    |  _ -> assert false
  in go

while difference is:

let diff       =
  let rec go l r = match l, r with
  | [], _  | _, []   -> l
  | x :: xs, y :: ys ->
    match compare x y with
    | -1 -> x :: go       xs (y :: ys)
    |  0 ->      go       xs       ys
    |  1 ->      go (x :: xs)      ys
    |  _ -> assert false
  in go

And so on. The three functions use the same underlying merge schemata; what varies is the operation to perform in each of the five possible cases:

  1. The left sequence is nil
  2. The right sequence is nil
  3. The next element in the left sequence is less than the next element in the right sequence
  4. The next element in the left sequence is equal to the next element in the right sequence
  5. The next element in the left sequence is greater than the next element in the right sequence

The question is, then, how many set operations can the merge algorithm implement? These five cases partition both input sequences in disjoint sets such that the output sequence is the natural merge of zero or more of them. For example, processing the sets { 1, 3, 4, 6, 7, 8 } ⋈ { 2, 3, 5, 6, 8, 9 } results in the following steps:

LNRNLTEQGTArguments
1 { 1, 3, 4, 6, 7, 8 } ⋈ { 2, 3, 5, 6, 8, 9, 10 }
2{ 3, 4, 6, 7, 8 } ⋈ { 2, 3, 5, 6, 8, 9, 10 }
3 { 3, 4, 6, 7, 8 } ⋈ { 3, 5, 6, 8, 9, 10 }
4 { 4, 6, 7, 8 } ⋈ { 5, 6, 8, 9, 10 }
5{ 6, 7, 8 } ⋈ { 5, 6, 8, 9, 10 }
6 { 6, 7, 8 } ⋈ { 6, 8, 9, 10 }
7 { 7, 8 } ⋈ { 8, 9, 10 }
8 { 8 } ⋈ { 9, 10 }
9,10∅ ⋈ { 9, 10 }

Abstracting away the operations to perform in each of these five cases we have the following schema:

let merge ln rn lt eq gt : 'a list -> 'a list -> 'a list =
  let rec go l r = match l, r with
  | [], ys -> ln ys
  | xs, [] -> rn xs
  | x :: xs, y :: ys ->
    match compare x y with
    | -1 -> lt x (go       xs (y :: ys))
    |  0 -> eq x (go       xs       ys )
    |  1 -> gt y (go (x :: xs)      ys )
    |  _ -> assert false
  in go

Both ln and rn must decide what to do with the remaining list and so have type α list → α list, while lt, eq and gt must decide what to do with the element in consideration and so have type α → α list → α list; thus the type of merge is (α list → α list) → (α list → α list) → (α → α list → α list) → (α → α list → α list) → (α → α list → α list) → α list → α list → α list. The operations on the remainder either pass it unchanged or return nil, while the operations on the next element either add it to the output sequence or not:

let self   xs =      xs
and null   _  =      []
and cons x xs = x :: xs
and tail _ xs =      xs

(some of these have well-known names in functional programming, but here I choose to use these neat, four-letter ones.) With the proviso that the output sequence must be increasing these four functions exhaust the possibilities by parametricity; otherwise, duplications and rearrangements would satisfy the parametric signature. Now union, intersection and difference are simply:

let union l r = merge self self  cons cons cons l r
and inter l r = merge null null  tail cons tail l r
and diff  l r = merge null self  cons tail tail l r

It is obvious that the question I posed above is answered as 25 = 32 possible set operations obtainable by varying each of the five operations. The next question is, then, what are these 32 operations? Let me characterize each of the five sets ln, rn, lt, eq and gt. The easiest one is eq, as it obviously is the intersection of both sets:

eq(A, B) = A ∩ B

By substitution in merge it is possible to show that ln(A, B) = rn(B, A) and vice versa; hence just one set expression suffices. The merge ends with rn for every element in A that is greater than every element in B, as the latter were included in the comparison sets; and conversely for ln. Hence:

ln(A, B) = B / A ≡ ⟨S y : B : ⟨∀ x : A : y > x⟩⟩
rn(A, B) = A / B ≡ ⟨S x : A : ⟨∀ y : B : x > y⟩⟩

(read A / B as "A over B"). All three sets are pairwise disjoint, since A / B ⊆ A and A / B ∩ B = ∅, and conversely, by construction.

The two remaining sets are also symmetric in the sense that lt(A, B) = gt(B, A) but are more difficult to characterize. My first attempt was to think of each element in A as being processed in turn and put into lt(A, B) just when strictly less than all the elements in B against which it could be matched, namely lt(A, B) = ⟨S x : A : x < ⟨min y : B : x ≤ y⟩⟩. The condition can be simplified with a bit of equational reasoning:

   x∈A ∧ x < ⟨min y : B : x ≤ y⟩
≡ { GLB }
   x∈A ∧ ⟨∀ y : y∈B ∧ x ≤ y : x < y⟩⟩
≡ { Trading }
   x∈A ∧ ⟨∀ y : y∈B : x > y ∨ x < y⟩⟩
≡ { Trichotomy }
   x∈A ∧ ⟨∀ y : y∈B : x ≠ y⟩⟩
≡
   x∈A ∧ x∉B

In other words, A − B. The problem is that, since the quantification over an empty set is trivially true, this set is too big as it includes the respective remainder; that is to say A / B ⊆ A − B as I showed above. To preserve disjointness I define:

lt(A, B) = A − B − A / B
gt(A, B) = B − A − B / A

In a Venn diagram, these five sets are:

Venn diagram of the five component sets

So by including or excluding one of the five components depending on the function passed to each of the five operations, the 32 set operations achievable by merge are:

Venn diagrams for all 32 set operations

Or in tabular form:

NlnrnlteqgtFunction
0selfselfconsconsconsA∪B
1selfselfconsconstailA ∪B/A
2selfselfconstailconsA∆B
3selfselfconstailtailA−B∪B/A
4selfselftailconsconsB ∪A/B
5selfselftailconstailA∩B∪B/A∪A/B
6selfselftailtailconsB−A∪A/B
7selfselftailtailtail B/A∪A/B
8selfnullconsconsconsA∪B−A/B
9selfnullconsconstailA ∪B/A−A/B
10selfnullconstailconsA∆B−A/B
11selfnullconstailtailA−B∪B/A−A/B
12selfnulltailconsconsB
13selfnulltailconstailA∩B∪B/A
14selfnulltailtailconsB−A
15selfnulltailtailtail B/A
16nullselfconsconsconsA∪B−B/A
17nullselfconsconstailA
18nullselfconstailconsA∆B−B/A
19nullselfconstailtailA−B
20nullselftailconsconsB −B/A∪A/B
21nullselftailconstailA∩B∪A/B
22nullselftailtailconsB−A−B/A∪A/B
23nullselftailtailtail A/B
24nullnullconsconsconsA∪B−B/A−A/B
25nullnullconsconstailA −A/B
26nullnullconstailconsA∆B−B/A−A/B
27nullnullconstailtailA−B−A/B
28nullnulltailconsconsB −B/A
29nullnulltailconstailA∩B
30nullnulltailtailconsB−A−B/A
31nullnulltailtailtail∅

Arguably, besides the traditional five set operations A ∪ B, A ∩ B, A − B, B − A and A ∆ B, only the remainders A / B, B / A and perhaps A / B ∪ B / A = A ⊔ B, the join of A and B (not to be confused with the relational operation), are independently useful. These three are obscure, and as far as I know have no name, although I'd love to hear if there is literature about them. This might mean that this exhaustive catalog of set merges is rather pointless, but at least now I know for sure.

2012-07-17

An Odd Lemma

While proving that every monad is an applicative functor, I extracted the following derivation as a lemma:

  fmap f ∘ (λh. fmap h x)
≡ { defn. ∘, β-reduction }
  λg. fmap f (fmap g x)
≡ { defn. ∘ }
  λg. (fmap f ∘ fmap g) x
≡ { Functor }
  λg. fmap (f ∘ g) x
≡ { abstract f ∘ g }
  λg. (λh. fmap h x) (f ∘ g)
≡ { defn. ∘, η-contraction }
  (λh. fmap h x) ∘ (f ∘)

for all f, x. This is the Yoneda Lemma in a special form. The term λh. fmap h x is the natural transformation between an arbitrary functor F and the functor Hom(X, —), where X is fixed by the type of x. Dan Piponi calls it the hardest trivial thing in mathematics. I didn't recognize it immediately (my category-fu is nonexistent), but the striking symmetry made me curious and I started investigating.

2010-05-15

Observable sharing and executable proofs

Every rational number has an infinite repeating decimal expansion. In particular, every integer admits two decimal expansions: one with and one without a "tail of nines". This means that 1 = 0.9999…, a provable fact that many people still choose to disbelieve. What do I mean here by "provable"? I mean "provable by finitistic reasoning", that is, by a purely algebraic proof free of notions of continuity or convergence. The proof in this particular (and particularly popular) case goes like this:

    x     = 0.9999…
≡ { product by base is shift of decimal point }
  10x     = 9.9999…
≡ { subtraction of finitely-expanded integer }
  10x - 9 = 0.9999…
≡ { identity of decimal expansion }
  10x - 9 = x
≡ { algebra }
   9x     = 9
≡ { cancellation of multiplication }
    x     = 1

The last two steps are elementary consequences of the fact that integers form a ring. The first three steps require justification. By definition, a decimal expansion of a positive real number x is a (possibly semi-infinite) sequence of digits 0 ≤ di < 10 such that for all i ≥ N, di = 0, and

x = ⟨∑ i :: 10i⋅di⟩

and the decimal point separates the digits corresponding to 100 and 10-1, namely d0 from d-1.

Lemma 0 Any decimal number x can be expressed as the sum of an integer and a (possibly infinite) decimal part

Proof: By definition, x has a finite number of decimal digits to the left of the decimal point, namely for 0 ≤ i < N, hence:

  x
= { definition }
  ⟨∑ i :: 10i⋅di ⟩
= { associativity of addition, a finite number N of times }
  ⟨∑ i : 0 ≤ i < N : 10i⋅di⟩ + ⟨∑ i : 0 > i : 10i⋅di⟩
= { naming the first term }
  n + ⟨∑ i : 0 > i : 10i⋅di⟩

This justifies representing a positive real number as a pair of an integer and a possibly infinite purely fractional decimal expansion:

type dec = Dec of int * int list

A sufficient condition for two numbers to be equal is that their decimal representations are the same:

let eq (Dec (e, ds)) (Dec (e', ds')) = e = e' && ds == ds'

Here, since the decimal expansions could be infinite, I use physical equality on lists. This is a stronger condition than stated above. Now, the proof for the first step:

Lemma 1 Multiplying a positive real number x by 10 shifts the decimal point in its decimal expansion one place to the right

Proof:

  10⋅x
= { lemma 0 }
  10⋅(n + ⟨∑ i : 0 > i : 10i⋅di⟩)
= { product finitely distributes over sum }
  10⋅n + 10⋅⟨∑ i : 0 > i : 10i⋅di⟩
= { product formally distributes over sum }
  10⋅n + ⟨∑ i : 0 > i : 10i+1⋅di⟩
= { substitute i := i-1 }
  10⋅n + ⟨∑ i : 1 > i : 10i⋅di-1⟩
= { lemma 0 }
  (10⋅n + d0) + ⟨∑ i : 0 > i : 10i⋅di-1⟩

Note that the third step is not strictly speaking finitist but formal: it requires taking the distributivity of product over a summation as axiomatic if induction is not accepted. This justifies the following:

let mul10 = function
| Dec (e, []     ) -> Dec (10 * e    , [])
| Dec (e, d :: ds) -> Dec (10 * e + d, ds)

The proof for the second step:

Lemma 2 Subtracting an integer k from a positive real number x only affects its integer part

Proof:

  x - k
= { lemma 0 }
  (n + ⟨∑ i : 0 > i : 10i⋅di⟩) - k
= { commutativity and associativity of sum }
  (n - k) + ⟨∑ i : 0 > i : 10i⋅di⟩

This justifies the following:

let subi (Dec (e, ds)) n =
  if e < n then failwith "subi" else
  Dec (e - n, ds)

Now the proof for the third step is purely executable and relies on physical equality of infinite lists of digits represented as circular lists:

let nines = let rec ds = 9 :: ds in Dec (0, ds)

Of course, for more general repeating expansions, equality on heads won't suffice, but this case is nicely handled by a one-liner:

let lemma = assert (eq nines (subi (mul10 nines) 9))

The purely verbal mathematical statement of this truth is that, since the first two steps operate on a finite prefix of the infinite expansion, the inifinite suffix is unchanged by the operations and hence identical to itself. In other words, the infinite summation is a purely formal object that remains unchanged from step to step.

Will this proof finally lay the "controversy" to rest, persuading even the staunchest finitist?

Edit: Thanks to id that alerted me in my confusion between structural and physical equality.

2009-06-05

Euler's Power Riffs

Ed Sandifer is an exegete of Euler's methods. In his latest "How Euler Did It", he shows the method through which he arrived at closed polynomials for the sums of n-th powers of the first x integers:

In general, Euler has us multiply the terms of Σxn by (n + 1)/(n + 2) x, (n + 1)/(n + 1) x, (n + 1)/n x, … (n + 1)/2 x and (n + 1)/1 x, respectively, and add the appropriate linear term and to get the formula for Σxn+1.

where Σxn is Euler's notation for the sum in question. The algorithm is very easily translated to a short OCaml function. First, I'll need a way to generate an integer range:

let range i j =
  let rec go l j =
    if j < i then l else
    go (j :: l) (pred j)
  in go [] j

Simple, eager, tail-call; almost imperative, definitely dirty. With that, Sandifer's rendition of Euler's method is simply:

open Num

let rec eusum n =
  if n = 0 then [Int 0; Int 1] else
  let l = eusum (pred n)
  (* n/1 n/2 n/3 ... n/(n+1) *)
  and k = List.map (fun d -> Int n // Int d) (range 1 (n + 1)) in
  (* the lowest coefficient is that of x^1 and always 0 *)
  let _ :: xs = List.map2 ( */ ) k l in
  (* normalize the sum of coefficients to 1 *)
  let x = List.fold_left ( -/ ) (Int 1) xs in
  Int 0 :: x :: xs

I want to clarify a couple of subtleties in the implementation. First, I use lists of Nums representing dense polynomials, least-order coefficient first. The base case is that for the sum of x ones (with exponent n = 0), namely x itself. As the article points out, the constant term is always zero, because a sum of zero powers is zero.

Then, the general case is built out of Σxn (called l) and the terms with which to multiply it (called k, in reverse order to that reported by Sandifer, which follows Euler in putting the highest power first). Both polynomials are multiplied term by term with map2; however, the degree is not yet raised at this stage. Also, the lowest term is always zero by the previous observation and I disregard it with an irrefutable pattern match (I could've used List.tl, but I prefer not to, in general). A simple fold_left computes the coefficient α for the linear term, and the polynomial is completed with it and the zero linear term, thus raising the degree by one (one term discarded, two added).

That's it, seven lines of code. Regrettably, it is not of much use as it is:

# eusum 4 ;;
- : Num.num list =
[Int 0; Ratio <abstr>; Int 0; Ratio <abstr>; Ratio <abstr>; Ratio <abstr>]

Bummer, abstract data types. It turns out, pretty-printing polynomials is surprisingly tedious to get right:

let pp_print_poly pp l =
  let pow = ref 0
  and any = ref false in
  Format.pp_open_box pp 1;
  List.iter (fun k ->
    if k <>/ Int 0 then begin
      if !any then Format.pp_print_space pp ();
      Format.pp_open_hbox pp ();
      (* Print coefficient *)
      if k </ Int 0 then Format.pp_print_char pp '-' else
      if !any then Format.pp_print_char pp '+';
      if !any then Format.pp_print_space pp ();
      let n = Num.abs_num k in
      if n <>/ Int 1 || !pow = 0 then Format.pp_print_string pp (Num.string_of_num n);
      (* Print formal power *)
      if n <>/ Int 1 && !pow > 0 then Format.pp_print_char pp '*';
      if !pow > 0 then Format.pp_print_char pp 'x';
      if !pow > 1 then Format.fprintf pp "^%d" !pow;
      Format.pp_close_box pp ();
      any := true
    end;
    incr pow
  ) l;
  Format.pp_close_box pp ()

I follow some very general rules for writing pretty printers. First, always accept the pretty printing channel pp as a parameter for maximum generality (you can then use it with #install_printer in the top-level). Second, always box everything. In this case I use a vertical-or-horizontal box around the entire polynomial, with one space for continuation lines, and a horizontal-only box around monomials.

As it is evident here, the usual rules for polynomials make the function logic-heavy: the first term uses the sign of the coefficient, but subsequent monomials are added or subtracted with the operation sign offset with spaces, and the coefficient is taken in absolute value. If the coefficient or the variable are 1 they are not shown, unless it is the constant term. If both the coefficient and the power are shown, they are separated by the multiplication sign.

With this it is easy to recreate the table in page 2 of the paper:

# let () =
    Format.open_box 1;
    List.iter (fun i ->
      pp_print_poly Format.std_formatter (eusum i);
      Format.print_newline ()) (range 0 8);
    Format.close_box ()
  ;;
x
1/2*x + 1/2*x^2
1/6*x + 1/2*x^2 + 1/3*x^3
1/4*x^2 + 1/2*x^3 + 1/4*x^4
-1/30*x + 1/3*x^3 + 1/2*x^4 + 1/5*x^5
-1/12*x^2 + 5/12*x^4 + 1/2*x^5 + 1/6*x^6
1/42*x - 1/6*x^3 + 1/2*x^5 + 1/2*x^6 + 1/7*x^7
1/12*x^2 - 7/24*x^4 + 7/12*x^6 + 1/2*x^7 + 1/8*x^8
-1/30*x + 2/9*x^3 - 7/15*x^5 + 2/3*x^7 + 1/2*x^8 + 1/9*x^9

As noted by Euler, the coefficient of the linear term is a Bernoulli number. There are more efficient ways to compute them, but in a pinch this serves:

let bernoulli n = let _ :: b :: _ = eusum n in b

2009-01-09

Why so Series

Notational Notions has an entry raising to the challenge of showing with power series that (sin x)² + (cos x)² = 1. I jestingly but unfairly called the author a cheater for using what to me was an obvious appeal to Euler's formula (I made a mistake in naming it as De Moivre's). I reasoned thus:

  (sin x)² + (cos x)²
= { Euler: e^ix - e^(-ix) = 2i sin x }
  [(e^ix - e^(-ix))/(2i)]² + [(e^ix + e^(-ix))/2]²
= { Binomial, twice }
 -(e^2ix - 2 e^ix e^(-ix) + e^(-2ix))/4 + (e^2ix + 2 e^ix e^(-ix) + e^(-2ix))/4
= { Algebra }
  (-e^2ix + 2 e^ix e^(-ix) - e^(-2ix) + e^2ix + 2 e^ix e^(-ix) + e^(-2ix))/4
= { Algebra throughout }
  (2 e^ix e^(-ix) + 2 e^ix e^(-ix))/4
= { Algebra }
  e^ix e^(-ix)
= { Negative power is inverse }
  1

The author defended the choice of purely imaginary exponentials as the most natural representation for the real series (1, 0, -1, 0, …). Meeting the challenge of showing a purely real derivation incurred in an ugly and possibly unavoidable case analysis. I conceded, but tried my hand at it. It is well-known that the series development for the sine and cosine are:

sin x = x - x^3/3! + x^5/5! - x^7/7! + …

cos x = 1 - x^2/2! + x^4/4! - x^6/6! + …

Squaring the first series:

 (sin x)²
= { Definition }
  (x - x^3/3! + x^5/5! - x^7/7! + …)²
= { Distribution }
      x^2/1!1! - x^4/3!1! + x^6/5!1! - x^8/7!1! + …
               - x^4/1!3! + x^6/3!3! - x^8/5!3! + x^10/3!7! - …
                          + x^6/1!5! - x^8/3!5! + x^10/5!5! - x^12/7!5! + …
                                     - x^8/1!7! + x^10/7!3! - x^12/7!5! + x^14/7!7! - …

(the use of pre-formatted text is crucial here to see what's going on). Squaring the second:

 (cos x)²
= { Definition }
 (1 - x^2/2! + x^4/4! - x^6/6! + …)²
= { Distribution }
  1 - x^2/2!0! + x^4/4!0! - x^6/6!0! + …
    - x^2/0!2! + x^4/2!2! - x^6/4!2! + x^8/6!2! - …
               + x^4/0!4! - x^6/2!4! + x^8/4!4! - x^10/6!4! + …
                          - x^6/0!6! + x^8/2!6! - x^10/4!6! + x^12/6!6! - …

Note that the terms are set so that columns line up by exponent. Note also that the factorial dividends add up to the exponent. Now:

 (sin x)² + (cos x)²
=
      x^2/1!1! - x^4/3!1! + x^6/5!1! - x^8/7!1! + …
               - x^4/1!3! + x^6/3!3! - x^8/5!3! + x^10/3!7! - …
                          + x^6/1!5! - x^8/3!5! + x^10/5!5! - x^12/7!5! + …
                                     - x^8/1!7! + x^10/7!3! - x^12/7!5! + x^14/7!7! - …
+
  1 - x^2/2!0! + x^4/4!0! - x^6/6!0! + …
    - x^2/0!2! + x^4/2!2! - x^6/4!2! + x^8/6!2! - …
               + x^4/0!4! - x^6/2!4! + x^8/4!4! - x^10/6!4! + …
                          - x^6/0!6! + x^8/2!6! - x^10/4!6! + x^12/6!6! - …
= { Associativity of addition }
  1 - x^2/2!0! + x^4/4!0! - x^6/6!0! + …
    + x^2/1!1! - x^4/3!1! + x^6/5!1! - x^8/7!1! + …
    - x^2/0!2! + x^4/2!2! - x^6/4!2! + x^8/6!2! - …
               - x^4/1!3! + x^6/3!3! - x^8/5!3! + x^10/3!7! - …
               + x^4/0!4! - x^6/2!4! + x^8/4!4! - x^10/6!4! + …
                          + x^6/1!5! - x^8/3!5! + x^10/5!5! - x^12/7!5! + …
                          - x^6/0!6! + x^8/2!6! - x^10/4!6! + x^12/6!6! - …
                                     - x^8/1!7! + x^10/7!3! - x^12/7!5! + x^14/7!7! - …
= { Binomial theorem }
  1 - (x - x)^2/2! + (x - x)^4/4! - (x - x)^6/6! + (x - x)^8/8! - …
=
  1

Of course there are two fishy manouvers going on here. The first is the use of ellipses to denote the series. This is a notational convenience that allows me to show the full expansions for easy pattern matching of terms. These ellipses actually encode a rigorous inductive law defining each term, and could be replaced by summations. The second is the use of associativity to freely reorder the terms of an infinite series without regard to convergence, of which I can't say anything as I'm no mathematician but which I believe it can be made airtight. In any case, you can take this as a purely symbolic manipulation.

2009-01-05

Ring Laws

This is a really basic derivation of elementary facts about rings, just for the record. I follow Milne's definition of a ring (R, +, 0, -, ⋅ 1) as a commutative group (R, +, 0, -) equipped with an additional operation ⋅ and a distinguished 1 ∈ R such that (R, ⋅, 1) is a monoid, and ⋅ distributes over + both on the left and on the right: for all a, b, c ∈ R

(a + b)⋅c = a⋅c + b⋅c
a⋅(b + c) = a⋅b + a⋅c

Lemma for any a, b ∈ R, a⋅0 = 0 = 0⋅b

   a⋅0 = 0
⇐ { Cancellation on the left }
   a⋅b + a⋅0 = a⋅b + 0
≡ { Distributivity on the right }
   a⋅(b + 0) = a⋅b + 0
≡ { Group neutral, twice }
   a⋅b = a⋅b
≡
   true

The other direction is analogous:

   0⋅b = 0
⇐ { Cancellation on the right }
   0⋅b + a⋅b = 0 + a⋅b
≡ { Distributivity on the right }
   (0 + a)⋅b = 0 + a⋅b
≡ { Group neutral, twice }
   a⋅b = a⋅b
≡
   true

Theorem for any a, b ∈ R, (-a)⋅b = -(a⋅b) = a⋅(-b)

   (-a)⋅b = -(a⋅b)
⇐ { Cancellation on the right }
   (-a)⋅b + a⋅b = -(a⋅b) + a⋅b
≡ { Distributivity on the right, group inverse }
   (-a + a)⋅b = 0
≡ { Group inverse }
   0⋅b = 0
≡ { Lemma }
   true

The other direction is analogous:

   a⋅(-b) = -(a⋅b)
⇐ { Cancellation on the left }
   a⋅b + a⋅(-b) = a⋅b + -(a⋅b)
≡ { Distributivity on the left, group inverse }
   a⋅(b + -b) = 0
≡ { Group inverse }
   a⋅0 = 0
≡ { Lemma }
   true

Also, the following is immediate:

Theorem If 1 = 0, then R = {0}

Proof: note that, for arbitrary a, b ∈ R:

   a
= { Definition }
   a⋅1
= { Hypothesis }
   a⋅0
= { Lemma }
   0
= { Lemma }
   0⋅b
= { Hypothesis }
   1⋅b
= { Definition }
   b

which shows that there is at most one element in R.

2008-08-31

Just draw something on the screen!

Nostalgia for times past sometimes leads to frustration and anger. Yes, times are harder for would-be pixel artists, but nothing a little work can't solve. Let's do fractals in OCaml!

Suppose you've already installed the wonderful findlib (vielen Dank, Herr Stolpmann!). Suppose, furthermore, you've succeded in compiling OcamlSDL. This might look like a heavy prerequisite to fulfill before showing pretty graphics, but it's a job you have to do only once. Then, something like plotting Newton fractals requires less than ninety lines of code, and the result is this:

I find natural to abstract iteration patterns into higher-order functions. In this case, centering and scaling the pixel plane over the complex plane (the so-called screen-space transformation) is tedious but repetitive code that could be reused. Furthermore, the details of striding, color conversion and such are the same from program to program, where the only thing that varies is the actual calculation performed on each point:

let iter_surface ?(mag=1.0) f g surface =
  let info   = Sdlvideo.surface_info surface
  and pixels = Sdlvideo.pixel_data_32 surface in
  let size   = min info.Sdlvideo.w info.Sdlvideo.h in
  let scale  = mag /. float size
  and cx, cy = info.Sdlvideo.w / 2, info.Sdlvideo.h / 2 in
  for i = 0 to info.Sdlvideo.h - 1 do
    let b = i * info.Sdlvideo.pitch / 4
    and y = float (cy - i) *. scale in
    for j = 0 to info.Sdlvideo.w - 1 do
      let x = float (j - cx) *. scale in
      let c = f x y in
      pixels.{b + j} <- Sdlvideo.map_RGB surface c
    done;
    g i
  done

The magnification mag specifies how many real units the shortest dimension spans. The outer loop iterates with i for every scanline, whose base address into the pixels array is b. The inner loop iterates with j the index to successive pixels, and from both i and j the coordinates x, y are computed. With those the function f is called, and the result c (an RGB triple of type Sdlvideo.color) is mapped and assigned to the pixel. At the end of every scanline, the function g is called so that the main program can display progress.

The simplest SDL event loop waits for the user to quit the application:

let rec event_loop () =
  match Sdlevent.wait_event () with
  | Sdlevent.QUIT -> ()
  | _ -> event_loop ()

Anything more complicated is easy to accomodate. Now, there are a myriad schemes for color rendering of mathematical data, but using the color wheel is very common, especially since the complex argument maps naturally to a hue value. Converting HSV colors to RGB is somewhat involved, but the following code is entirely reusable:

let quomod d x =
  let f, e = modf (float d *. (x -. floor x)) in truncate e, f

let byte_of_frac x =
  let mask x = lnot (x asr 30) land 0xff in
  let n = truncate (ldexp x 8) in
  (n land (mask n))  lor  mask (n - 256)

let hsv_to_rgb (h, s, v) =
  let epsilon = 1e-3 in
  let w = byte_of_frac v in 
  if abs_float s < epsilon then (w, w, w) else
  let (e, f) = quomod 6 h in
  let x = byte_of_frac (v *. (1. -. s))
  and y = byte_of_frac (v *. (1. -. s *. f))
  and z = byte_of_frac (v *. (1. -. s *. (1. -. f))) in
  match e with
  | 0 -> (w, z, x)
  | 1 -> (y, w, x)
  | 2 -> (x, w, z)
  | 3 -> (x, y, w)
  | 4 -> (z, x, w)
  | 5 -> (w, x, y)
  | _ -> assert false

quomod scales the fractional part of x by the integer d, and splits the result into an entier part e and a fractional part f. This clasifies the hue into a sextant on which the fractional part linearly interpolates. To map a real number in the interval [0… 1] into [0… 256), byte_of_frac scales the fraction by 28 and clips the result to a byte by constructing the appropriate bit-mask.

The code so far doesn't incorporate any specific functionality, and could be shipped off to a standard module to use again and again. The particular application that I had in mind, Newton fractals, are the plot on the complex plane of the basins of attraction to the complex roots of a polynomial, as found by applying the Newton method.

Polynomials and their derivatives are especially easy to evaluate at a point in one pass using Horner's method, when represented as the list of their coefficients with the independent term at the head, a so-called dense representation of polynomials:

open Complex

let eval p x =
  List.fold_left (fun (y, dy) p ->
    (add p (mul y x), add y (mul dy x)))
    (zero, zero) p

With that, the Newton iteration is very simple:

let nr_root ?(max_iter=32) ?(epsilon=1e-6) p x =
  let rec go (x, i as res) =
    if i = max_iter then res else
    let y, dy = eval p x in
    if norm y < epsilon then res else
    let dx = div y dy in
    go (sub x dx, succ i)
  in go (x, 0)

The root finder returns the number of iterations required, up to a maximum and/or specified precision, together with the closest estimate to the root. This allows me to "color" the fractal according to the time needed to converge to a root. With the chosen representation, the polynomial zn - 1 (which represents, by definition, the n-th roots of unity) is straightforward to construct:

let pow n x =
  let rec go i l =
    if i = 0 then l else go (pred i) (zero :: l)
  in go n [x]

let nth_roots n = (neg one) :: pow (n - 1) one

These twenty lines are all the mathematical machinery required. It only remains to put all the pieces together:

let () =
  let max_iter = 128 in
  Sdl.init [`EVERYTHING];
  let screen = Sdlvideo.set_video_mode ~w:640 ~h:480 ~bpp:32 [`SWSURFACE] in
  Sdlwm.set_caption ~title:"nrfractal" ~icon:"";
  let p = nth_roots 3 in
  iter_surface ~mag:2.0 (fun x y ->
    let (z, n) = nr_root ~max_iter p { re = x; im = y } in
    let hue = Complex.arg z *. 0.159154943091895346 +. 0.5
    and vlu = float (max_iter - n) /. float max_iter in
    hsv_to_rgb (hue, 1., vlu) )
    (fun i -> if i land 7 = 0 then Sdlvideo.flip screen)
    screen;
  Sdlvideo.flip screen;
  event_loop ();
  Sdl.quit ()

The iteration limit max_iter controls the fractal detail. I first initialize the SDL library and create a big enough window, with a descriptive title. Then, I build a polynomial p representing the third roots of unity, and invoke iter_surface with a function that finds p's root starting from the complex value corresponding to the coordinates of the point to be plotted. I convert the complex angle of the returned approximation z to a hue value (by dividing by 2π and normalizing to the range [0… 1]); analogously, the number of iterations determines the vividness, or value of the color. This is converted to an RGB triple and returned. Every eight scanlines the screen is refreshed with the contents of the backing buffer. Once every point is plotted, the event loop is entered until the user decides to quit.

So, not counting the generic, reusable code, less than 40 lines of code suffice to just draw something on the screen. Have fun!

2008-07-09

The Zero-One Principle

Continuing my treatment of comparator networks, I've experimenting with the discovery of minimal sorting networks using evolutionary programming, about which I plan to write in the future. One of the things needed for that is a means of quickly evaluating if a given network actually sorts all sequences of a given length or not. There is a very powerful theoretical result, called the zero-one principle that aids in identifying sorting networks from general comparator networks. It says:

(The zero-one principle): A comparator network N of width n sorts all sequences of length m if it sorts all 2n binary sequences of length n.

The proof is not difficult, and is a good showcase for the equational reasoning I'm so fond of. First, some definitions.

Given an ordered set (A, ≤), the minimum of elements a, a′ ∈ A is the unique element a↓a′ such that:

a↓a′ ≤ a ∧ a↓a′ ≤ a′

Dually, the maximum of a, a′ ∈ A is the unique element a↑a′ such that:

a ≤ a↑a′ ∧ a′ ≤ a↑a′

A monotonic function from ordered set (A, ≤) to ordered set (B, ⊑) is a function that preserves orderings, that is:

⟨∀a, a′ ∈ A : a ≤ a′ : f.a ⊑ f.a′⟩

With that, I can prove the first

Lemma: monotonic functions preserve minima (dually, maxima)

Let f be monotonic on A. Then ∀a, a′ ∈ A, by definition of minimum:

   a↓a′ ≤ a ∧ a↓a′ ≤ a′
⇒ { monotonicity of f }
   f.(a↓a′) ≤ f.a ∧ f.(a↓a′) ≤ f.a′
≡ {definition of ↓ }
   f.(a↓a′) = f.a ↓ f.a′

And dually for maxima.

Recall that an comparator box is a function (x ↕ y) = (x↓y, x↑y). Applied to a sequences, a comparator [i ↕ j] with 0 ≤ i ≤ j < n sorts elements i and j of [l0, l1,… ln-1] into [l0,… li-1, li↓lj, li+1,… lj-1, li↑lj, lj+1,… ln-1]. A comparator network N is a sequence of such boxes, and as a sequence, it has a monoidal structure under composition. For a function f : A → B, the pointwise extension map f : An → Bn maps a sequence [l0, l1,… ln-1]∈An to the sequence [f.l0, f.l1,… f.ln-1]∈Bn. Hence, the following

Lemma: comparator networks commute with the pointwise extension of monotonic functions

Let f be monotonic on A. By induction on the structure of N, we have for the base case a single comparator box [i ↕ j]:

   ([i ↕ j] ∘ map f).l
= { defn. pointwise extension }
   [i ↕ j].[f.l0, f.l1,… f.ln-1]
= { defn. comparator box }
   [f.l0,… f.li-1, f.li ↓ f.lj, f.li+1,… f.lj-1, f.li ↑ f.lj, f.lj+1,… f.ln-1]
= { lemma }
   [f.l0,… f.li-1, f.(li↓lj), f.li+1,… f.lj-1, f.(li↑lj), f.lj+1,… f.ln-1]
= { defn. pointwise extension }
   map f.[l0,… li-1, li↓lj, li+1,… lj-1, li↑lj, lj+1,… ln-1]
= { defn. comparator box }
  (map f ∘ [i ↕ j]).l

For the inductive case, a comparator network N can be decomposed as the composition of comparator networks NL and NR, such that:

   N ∘ map f
= { compositionality of network }
   (NL ∘ NR) ∘ map f
= { associtativity of composition, base case }
   NL ∘ (map f ∘ NR)
= { associtativity of composition, base case }
   map f ∘ (NL ∘ NR)
= { compositionality of network }
   map f ∘ N

Now I can prove the

Theorem: The zero-one principle as stated above is valid.

Suppose there is a sequence l∈An that is not sorted by N, that is, for m = N.l there is an 0 ≤ k < n-1 such that mk ≰ mk+1. For all c∈A, define f : A → {0, 1} as f.c = 1 if mk ≤ c, and f.c = 0 otherwise. Now f is monotonic:

   f.c ≤ f.c′
⇐ { range of f }
   f.c = 0 ∨ f.c′ = 1
≡ { defn. f }
   mk ≰ c ∧ mk ≤ c′
≡
   mk ≤ c ⇒ mk ≤ c′
⇐ { transitivity ≤ }
    c ≤ c′

Furthermore, f.mk = 1 and f.mk+1 = 0 which means that (map f ∘ N).l is unsorted. But since f is monotonic, by the second Lemma, (N ∘ map f).l is unsorted too. In other words, if there's a sequence over An that N doesn't sort, then there is a sequence over 2n that N doesn't sort either; by contrapositive the theorem follows.

In the next post I'll show how this translates in an efficient test for comparator networks.

2008-06-09

Random Experiments

I often find it is much more rewarding to tackle an algorithmic or mathematical problem from an experiential perspective rather than from a formal one. Trying to get some test data for an earlier program, I thought about a couple of ways of generating random increasing sequences.

Intent on generating n random points 0 = x0 ≤ x1 ≤ … ≤ xn-1 = 1 on the unit interval, I thought the most efficient way to do it was to at each step randomly split the remaining interval. This suggested itself:

rnd0[n_Integer?Positive] := 
 Append[NestList[# + RandomReal[1 - #]&, 0., n - 2], 1.]

(NestList generates a list of n+1 elements.) The result was less than satisfying:

ListLinePlot[rnd0[10], PlotRange -> All]

Not good; in fact, more points show that it indeed wasn't what I wanted:

ListLinePlot[rnd0[100], PlotRange -> All]

It is apparent that the not after long a pretty large chunk of the interval is carved off, and the remaining points compete for less and less room, giving rise to a knee. Well, back to the drawing board. It then occurred to me that instead of dividing the remainder I could step away from the last point a random distance. The problem with this is that I couldn't know in advance how far I'd end; hence, I needed a way to normalize the resulting random walk:

scale[l_List] := With[{h = l[[-1]] - l[[1]]}, (l - l[[1]])/h]

The sequence could be generated simply with:

rnd1[n_Integer?Positive] := 
 scale[Accumulate[Array[RandomReal[1]&, n]]]

(Accumulate gives the list of partial running sums.) To see how I fared, I plotted some points:

ListLinePlot[rnd1[10], PlotRange -> All]

Not bad; however, 100 points showed that the result was a bit too smooth:

ListLinePlot[rnd1[100], PlotRange -> All]

What to do now? Sometimes the most obvious solution glares so intently to you that you can't help but notice it the last, if at all. Since the difference between a random sequence and a monotone random sequence is order, why not try to sort a plain old bunch of random points?

rnd2[n_Integer?Positive] := scale[Sort[Array[RandomReal[1] &, n]]]

(I used scale to ensure that the first and last points where 0 and 1 respectively. I could have added them to a sequence n-1 points long just as well.) Few points showed an encouraging picture:

ListLinePlot[rnd2[10], PlotRange -> All]

Indeed, a few big jumps nicely interspersed with smaller increments. More points still maintain the picture of satisfying roughness:

ListLinePlot[rnd2[100], PlotRange -> All]

Now I wasn't satisfied with this. Is there a way to generate a just-rough-enough sequence like this without having to do an O(n log n) sorting pass? Time to investigate a bit what exactly is the origin of the roughness. One thing I found that distracted me from seeing this quality was the monotonicity. Maybe by looking at the detrended sequences could shed some light into this. The detrending is simply achieved by subtracting the line between the first and last points:

norm[l_List] := 
 With[{h = l[[-1]] - l[[1]], n = Length[l] - 1}, 
  l - l[[1]] - Table[i*h/n, {i, 0, n}]]

(This is more general than I needed, as by construction the trend is the line y = i/n.) I didn't feel that investing analysis on rnd0 was necessary. This is what rnd1 looks like detrended:

ListLinePlot[norm[rnd1[100]], PlotRange -> All]

And this is rnd2:

ListLinePlot[norm[rnd2[100]], PlotRange -> All]

If there was a difference I couldn't see it. But maybe there is some other way to look at the data. Maybe this roughness is a quality of the differences, the jumps between values in the sequence. Now the difference themselves would also be random, albeit with a different underlying distribution. And that's exactly the right idea; since I'm after all accumulating differences, what should their underlying distribution be to achieve the right degree of roughness? I needed first a way to compute histograms of the differences:

hist[k_Integer?Positive, l_List] := With[{d = Differences[l]},
  With[{hi = Max[d], lo = Min[d]}, 
   With[{st = k/(hi - lo)}, 
    Module[{h = Array[0&, k + 1]}, 
     Do[h[[Floor[st (x - lo)] + 1]]++, {x, d}]; h]]]]

(The fact that With binds in parallel, like ML's let is not very practical, as Mathematica's syntax is a bit heavy.) The function hist counts frequencies in k bins. Now, after seeing this:

ListLinePlot[hist[50, rnd1[10000]], PlotRange -> All]

it struck me as obvious: by construction, rnd1 differences are uniformly distributed. The mean 200 is simply 10,000 points divided by 50 bins. The last bin is empty as, since it's evident from the code above, I used a last value as a sentinel for the maximum possible value. Here is rnd2:

ListLinePlot[hist[50, rnd2[10000]], PlotRange -> All]

Now I was getting somewhere. The histogram looked like an exponential distribution. I could have read the article and tried to determine λ to within some confidence interval, in effect testing the hypothesis that the differences were indeed exponentially distributed; or I could try fitting an exponential to the data to see if it would lead me to somewhere. This is what I finally tried:

j = hist[100, rnd2[10000]];

FindFit[j, a k Exp[-k x], {a, k}, x]
{a -> 10420.3, k -> 0.0990637}

A reasonable fit: the a⋅e-λ is approximately 10,000, the number of points. To really see how good the fit was I needed, however, a picture:

Plot[Evaluate[ a k Exp[-k x] /. %], {x, 0, 50}, 
 Epilog -> MapIndexed[Point[{#2[[1]], #1}]&, j]]

Very tight fit, indeed. With this, I could directly generate a monotone sequence with the right distribution of jumps, in one pass without having to sort:

rnd3[n_Integer?Positive] := 
 scale[Accumulate[Array[-Log[1 - RandomReal[1]]/n&, n]]]

(I used the standard algorithm to generate an exponentially distributed variate.) Graphically, there was no appreciable qualitative difference between rnd2 and rnd3:

ListLinePlot[rnd3[100], PlotRange -> All]

And the histogram, as I expected it would be by construction, was likewise similar:

ListLinePlot[hist[50, rnd3[10000]], PlotRange -> All]

I was initially quite satisfied after this hands-on approach to discovering the right algorithm. This satisfaction would prove to be ephemeral, as I realized none of this explained the why. Maybe it's time to consult The Art of Computer Programming, but that would be for another post.

2008-06-05

A Circular Triviality

(Minor corrections) This is my own take at the corollary to the Intermediate Value theorem on the unit circle. It is not substantially different to the proof linked to, so its value is zero everywhere (maybe except for a countable subset). In particular, this is a special case of the much more general Borsuk-Ulam theorem.

Consider f a continuous function with domain on the unit circle. Then there exists at least a pair of points x0, x1 ∈ [0…2π) such that f(x0) = f(x0 + π) and f(x1) = f(x1 + π).

This is a purely existential statement, so it would be somewhat pointless (!) to strive to identify (or construct) a value x0 witnessing to the theorem. The only obvious possibility is to work with f by itself; that is, by focusing on the range and not on the domain of f. I must head towards something to which I can apply the IVT, and this lets me fix a number of things.

First, f is real-valued, as required by the IVT. Second, f is periodic, since if there were an x for which f(x) ≠ f(x + 2π), the difference |f(x + 2π) - f(x + ε)| would be large for |x + 2π - x - ε| = |ε| arbitrarily small on the unit circle. Third, since I see that the hypothesis mentions f(x) and f(x + π) being equal, I'm forced to consider the difference g(x) = f(x) - f(x + π). Fourth, g is continuous as f is.

Now, I have a choice in stating the proof; or rather, I have two different ways to show the existence of an interval satisfying the IVT.

The first one is to consider an arbitrary x ∈ [0…2π), and to see that g(x + π) = f(x + π) - f(x + 2π) = f(x + π) - f(x) = -g(x). Since -|g(x)| ≤ 0 ≤ |g(x)| and g is continuous, there must be some x0 ∈ [-|g(x)| … |g(x)|] such that g(x0) = 0.

The other possibility is to note that there must be a pair of values x+ and x- such that g(x+) ≥ 0 and g(x-) ≤ 0, or else f would be monotonic, hence discontinuous. Since g(x-) ≤ 0 ≤ g(x+) and g is continuous, there exists an x0 ∈ [x- … x+] such that g(x0) = 0.

Both things amount to the same; namely, there is a zero of g in the unit circle. As the domain is periodic, there must be another zero in the interval "across π" with respect to the original one (that is, the one having the endpoints listed in the opposite order, modulo 2π), call it x1. But this is precisely the posited pair of points that make f(x0) = f(x0 + π) and f(x1) = f(x1 + π).

Note that both choices are not equally potent: the first interval [-|g(x)| … |g(x)|] has the endearing property of being a function of an arbitrary x on the unit circle. This not only means that I can trade two existentials for one universal, making the theorem more constructive; it also permits a direct application of a root-finding method to hone in the zero.

The usual formulation of this is that there exists a pair of antipodal points where the function takes the same value. We see there is more than one, as the inventors of the Gömböc show. Also, as opposed to the proof by contradiction stated there, this is a direct proof.

2008-05-04

Romantic, or Enlightened?

It is not much of a confession to say that I'm of a rather dogmatic disposition regarding the need for strict, clear and unambiguous formalisms —let this serve— in the endeavors of formal thought. Mathematics are a testament to this, but, somewhat inexplicably, no other human science has gone as far. Computer Science, its little sibling, should have been expected to follow its steps; instead, it diverges more and more from its path.

What's more, the tendency to cast "programming" (in the most extended of its senses: every thing that pertains to computing, to use a single umbrella term) in formal languages is often derided as "academicism". This is a feud that is being waged for two thousand and five hundred years now, and a resolution is nowhere to be expected. On the one hand, Platonists continue to peer in the gloom for the shapes that make the shadows (pity they never prefigured the Necker Cube); on the other, Cynics bark their laughter at them while at the same time correctly pointing out that so much is lost in the translation from idea to actuality that they might as well build castles in the sand.

And yet I personally loathe Platonism as a general philosophy, and would rather side with Naturalism. For me, the way out of this inconsistency is not to regard algebra as operating on thought entities but as a mere shorthand for the operations to be actually carried out in the real world. I'm not so much a Platonist as I am a poet: words not only have meaning, that is, representational power in the real world, but grammar is inescapable. And that is a good thing, I think.

This autobiographical irrelevance is prompted by a couple of blog posts that resonated deeply with my convictions and I find them oddly fitting in a web of thought that has lately come to preoccupy me somewhat. I've written on Reddit something that summarizes quite pithily what I think in general of the art of programming:

[…] it violates the fundamental heuristic of using a stricter object (in the mathematical sense) even when a laxer one would do. An array is isomorphic to a finite map from an initial natural interval [0..N), so I can't really see a way to justify on stylistic grounds using a partial function (a cascade of ifs) instead of a finite map.

(emphasis added.) Apart from the specifics I was objecting to, I was merely echoing Dijkstra's emphasis on strengthening the antecedents as an effective way to carve out, and thus reduce the complexity associated with an unbounded state- (problem-) space. This was not lost on someone that replied with what I surmise was intended to be faint praise:

You and Dijkstra would have gotten along great.

I felt quite honored, but that's not to the point. My conviction on the generative power of restrictions comes from not only the practice of programming, but also from many years of Literary Workshop attendance. I was more apt to thrive with baroquely artificial writing assignments than with "open-ended" ones. This is the Oulipo lesson, too (and it is not fortuitous that it numbered a mathematician) and von Trier's Dogme 95's, and I'm sure there must be many other examples of "creating by restricting". The other side of this coin is what I think is the Romantic imperative of giving free rein to the creative process. The opening paragraph of the Wikipedia article is, for me, revelatory not so much in what it says but in how much, how viscerally, violently I disagree with what it expounds:

Romanticism […] was partly a revolt against aristocratic social and political norms of the Age of Enlightenment and a reaction against the scientific rationalization of nature, and was embodied most strongly in the visual arts, music, and literature.

I cannot think of a better expression of the Romantic mind than this:

Want to know who are the damned creative magicals? Philosophers. I hate their guts.

Analyzing this would throw me on a long limb, and I am probably misrepresenting the quote anyway; suffice to say that I place myself on the antipode of that thought. And fortunately there are still some thinkers that would accompany me in that lonely spot. This is a wonderful sample of that kind of rational, analytic, articulate thought that I feel is much needed in the discipline of programming:

But the power of languages is to abstract these things away, so that we can concentrate on the algebra while the machine does the arithmetic. Yet, I think we owe the acceptance of null largely to the fact that people are used to sloppy thinking. They’re used to treating nothing as if it were a thing, and using concepts of which they only know the approximate meaning.

(again, emphasis mine.) And this conclusion is reached by unwaveringly adhering to such fundamental logical principles as the Law of Identity and the Law of the Excluded Middle. Computers are Aristotelian, whether we like it or not.

Following a link from that article, I found another gem of the Euclidean approach to computing, that I leave for you to ponder. Personally, I find exciting to think about the result of embarking on a path that begins by recognizing that a value is unique, eternal, immutable, and is not fixed in time or space and continuing from there.

By the way, it is Romantics I blame for artificially distinguishing between "art" and "discipline": as of the XVIII Century, "art", as in Greek techne sufficed for —or rather there was no distinction between— both.

2008-04-11

Fixnum Figurations

The polytopic numbers P.k.n are a kind of figurate numbers giving the number of lattice points on a k-dimensional simplex (a right triangle, tetrahedron, …) x0 + x1 + … + xk-1 ≤ n of integer side n. For k = 2 we have the Pythagorean triangular numbers 1, 3, 6, 10…, the last one called the tetractys. In general, the way to calculate P.k.n is to stack n simplices of dimension k - 1 with decreasing sides:

P.0.n = 1
P.k.n = ⟨∑ i : 0 ≤ i < n : P.(k-1).i⟩

This can be put in closed form by taking antidifferences:

P.k.n = n↑k / k!

where n↑k = n⋅(n + 1)⋅…⋅(n + k - 1) is the rising factorial power. For positive, integer k, n↑k = (n + k - 1)↓k, which is the falling factorial power, defined analogously. Since by definition the binomial coefficient C.n.k = n↓k / k!, we have that:

P.k.n = C.(n + k - 1).k

and the problem of calculating polytopic number can be reduced to the calculation of binomial coefficients. Mark-Jason Dominus has posted a while ago about the best method to compute binomial coefficients of machine-size precision without integer overflow. I'll weave his insight in the formal derivation of a program satisfying these requirements.

What is the largest binomial coefficient we can compute? Let B be the bit width of a positive machine integer (in OCaml, B = 30). For a given n, the largest binomial coefficient is C.n.⌊k/2⌋, so that:


   C.n.k < 2B
⇐
   C.(2⋅k).k < 2B
≡
   (2⋅k)↓k / k! < 2B
≡
   ⟨∏ i : 0 ≤ i < k : 2⋅k - i⟩/⟨∏ i : 1 ≤ i ≤ k : i⟩ < 2B
≡ { Changing the index of summation }
   ⟨∏ i : 1 ≤ i ≤ k : k + i⟩/⟨∏ i : 1 ≤ i ≤ k : i⟩ < 2B
≡ { Algebra }
   ⟨∏ i : 1 ≤ i ≤ k : (k + i)/i⟩ < 2B
⇐ { Taking logarithms }
   ⟨∑ i : 1 ≤ i ≤ k : lg.((k + i)/i)⟩ < B

Define S.k = ⟨∑ i : 1 ≤ i ≤ k : lg.((k + i)/i)⟩. Then S.0 = 0, and S.(k + 1) - S.k = 1 + lg.(2⋅k + 1) - lg.(k + 1) (the calculations are tedious but straightforward.) Hence, the following program is immediate:

let maxcomb b =
  let rec iter s k =
    let d = log (float (4*k + 2) /. float (k + 1)) in
    if s < d then k else iter (s -. d) (k + 1)
  in iter (log 2. *. float b) 0

With this, we can determine that for B = 30, k must be at most 16. This is a bit too strong, since C.33.15 < 230 < C.33.16. I'm now set to derive the program to compute binomial coefficients. But first, I'll need this later (foreshadowing!):

let rem m n =
  let r = m mod n in
  if r < 0 then r + abs n else r

let rec gcd m n =
  if n = 0 then m else gcd n (rem m n)

(cf Daan Leijen, Division and Modulus for Computer Scientists.) The program must fulfill the following specification:

let binomial n k =
  (* 0 ≤ k ≤ n ≤ 32 *)
  let r = ref ... in
  (* Binomial *)
  !r

However, with the stated conditions, C.n.k = C.n.(n - k), and I can take advantage of that to minimize work by forcing k ≤ n - k:

let rec binomial n k =
  if k > n - k then binomial n (n - k) else
  (* 0 ≤ k ≤ n - k < n ≤ 32 *)
  let r = ref ... in
  (* Binomial *)
  !r

From the definition of binomial coefficient:

P ≡ r = ⟨∏ i : 0 ≤ i < k : n - i⟩/⟨∏ i : 0 ≤ i < k : i + 1⟩

Changing the constant k by a variable h:

P.h ≡ r = ⟨∏ i : 0 ≤ i < h : n - i⟩/⟨∏  i : 0 ≤ i < h : i + 1⟩

with P.k ≡ P. The initialization:

  let r, h = ref 1, ref 0 in ...

establishes P.0. We have:

let rec binomial n k =
  if k > n - k then binomial n (n - k) else
  (* 0 ≤ k ≤ n - k < n ≤ 32 *)
  let r, h = ref 1, ref 0 in
  (* P.h *)
  while !h != k do
    (* increment h under invariance of P.h *)
  done;
  (* P.k *)
  !r

Now:

   P.(h+1)
≡
   r = ⟨∏ i : 0 ≤ i < h+1 : n - i⟩/⟨∏ i : 0 ≤ i < h+1 : i + 1⟩
≡ { Splitting the range }
  r = (⟨∏ i : 0 ≤ i < h : n - i⟩⋅(n - h)) / (⟨∏ i : 0 ≤ i < h : i + 1>⋅(h+1))

which is restored with:

  r := !r * (n - !h) / (!h + 1);
  h := !h + 1

The assignment to r might unfortunately cause overflow. But here's the trick: if (n-h) / (h+1) = a/b in lowest terms, then r⋅(n-h) / (h+1) = (r/g)⋅a / (b/g), where g = (r, b). This is the most simplifying we can do while computing the binomial coefficient stepwise. Hence:

let rec binomial n k =
  if k > n - k then binomial n (n - k) else
  (* 0 ≤ k ≤ n - k < n ≤ 32 *)
  let r, h = ref 1, ref 0 in
  (* P.h *)
  while !h != k do
    let f = gcd (n - !h) (!h + 1) in
    let a, b = (n - !h) / f, (!h + 1) / f in
    let g = gcd !r b in
    r := (!r / g) * a / (b / g);
    incr h
    (* P.(h+1) *)
  done;
  (* P.k *)
  !r

is the desired program. With this, computing polytopic numbers is just a matter of defining:

let polytopic k n =
  if n = 0 then 0 else
  binomial (n + k - 1) k