## 2011-12-30

### (One by) Four by Nine

The four nines puzzle asks which positive numbers can be obtained from arithmetic expressions involving four "9" digits and an assortment of operations, minimally "+", "-" (binary and unary), "×" and "/", also frequently "√" and sometimes "!" (factorial). I'll show how to find them by brute force. In this case, I'll forgo using factorials; this means that not every number under 100 can be so obtained. As it is, at 130-odd lines, this is a longish program.

Expressions are labeled trees, where type α is the type of labels and type β is the type of leaves:

```type ('a, 'b) expr =
| Con of 'a * 'b
| Neg of 'a * ('a, 'b) expr
| Sqr of 'a * ('a, 'b) expr
| Add of 'a * ('a, 'b) expr * ('a, 'b) expr
| Sub of 'a * ('a, 'b) expr * ('a, 'b) expr
| Mul of 'a * ('a, 'b) expr * ('a, 'b) expr
| Div of 'a * ('a, 'b) expr * ('a, 'b) expr
```

Later I'll make clear what labels are useful for. For now, the simplest operation on an expression is extracting its `label`:

```let label = function
| Con (l, _    )
| Neg (l, _    )
| Sqr (l, _    )
| Add (l, _, _ )
| Sub (l, _, _ )
| Mul (l, _, _ )
| Div (l, _, _ ) -> l
```

Note that this is not a recursive function; it just extracts the label of the root. Converting expression trees into algebraic syntax is a tedious exercise in unparsing a precedence operator grammar. In this case, since the leaves are of an arbitrary type, `format` needs a conversion function `cnv`:

```let format cnv e =
let buf = Buffer.create 10 in
let rec go level e =
let prf op prec e =
if prec < level then Buffer.add_char buf '(';
go prec e;
if prec < level then Buffer.add_char buf ')'
in
let bin op prec e e' =
if prec < level then Buffer.add_char buf '(';
go prec e;
go prec e';
if prec < level then Buffer.add_char buf ')'
in match e with
| Con (_, x    ) -> Buffer.add_string buf (cnv x)
| Neg (_, e    ) -> prf "-" 10 e
| Sqr (_, e    ) -> prf "\xe2\x88\x9a" 10 e
| Add (_, e, e') -> bin "+"  1 e e'
| Sub (_, e, e') -> bin "-"  2 e e'
| Mul (_, e, e') -> bin "*"  5 e e'
| Div (_, e, e') -> bin "/"  6 e e'
in go 0 e; Buffer.contents buf
```

The inner function `prf` formats prefix operators with precedence `prec`, while `bin` formats binary operators. In both cases, if `op`'s binding power `prec` is less than the current precedence `level`, the whole expression is surrounded by parentheses. Note that I use the UTF-8 representation of the radical sign `U+221A`; who said that OCaml doesn't do Unicode?

If expressions are labeled with their values, they can be evaluated in constant time with `label`; to avoid losing precision, I use rational labels. If any expression is out of range, I use the "fraction" 1/0, `infinity_ratio`, as a sentinel. For this to work with OCaml `ratio`s, I must turn off error checking on null denominators:

```let () = Arith_status.set_error_when_null_denominator false

let infinity_ratio = Ratio.(inverse_ratio (ratio_of_int 0))
```

How to compute the square root of a fraction? If both the numerator and denominator are positive perfect squares, the answer is clear. In any other case, I signal failure via `infinity_ratio`:

```let sqrt_ratio r =
let open Ratio in
let open Big_int in
match sign_ratio r with
| -1 -> infinity_ratio
|  0 -> r
|  1 ->
let r1, r2 = numerator_ratio r, denominator_ratio r in
let q1, q2 = sqrt_big_int   r1, sqrt_big_int     r2 in
let s1, s2 = square_big_int q1, square_big_int   q2 in
if eq_big_int r1 s1 && eq_big_int r2 s2
then div_big_int_ratio q1 (ratio_of_big_int q2)
else infinity_ratio
|  _ -> assert false
```

Low-level but straightforward. Now smart constructors make sure that every expression is correctly labeled with its value:

```let con n    = Con (Ratio.ratio_of_int        n            , n    )
and neg e    = Neg (Ratio.minus_ratio  (label e)           , e    )
and sqr e    = Sqr (      sqrt_ratio   (label e)           , e    )
and sub e e' = Sub (Ratio.sub_ratio    (label e) (label e'), e, e')
and mul e e' = Mul (Ratio.mult_ratio   (label e) (label e'), e, e')
and div e e' = Div (Ratio.div_ratio    (label e) (label e'), e, e')
and abs e    =
let r = label e in
if Ratio.sign_ratio r != -1 then e else
Neg (Ratio.minus_ratio r, e)
```

The stage is set for generating all the expressions. If `nines size` generates all the expressions using `size` nines, the recursion schema it obeys is something like this:

```  if size == 1 then [   9] else
if size == 2 then [  99] @ mix (nines 1) (nines 1) else
if size == 3 then [ 999] @ mix (nines 2) (nines 1) @ mix (nines 1) (nines 2) else
if size == 4 then [9999] @ mix (nines 3) (nines 1) @ mix (nines 2) (nines 2) @ mix (nines 1) (nines 3) else...
```

where `mix` is a hypothetical function that merges two lists of expressions using all the possible binary operators. In each case, the constant 10size - 1 = 99…9 is the base of the recursion, which proceeds by building all binary trees of a given `size` by partitioning `size` in two summands. OK, maybe it is simpler to show the actual code than trying to explain it in English. A small function generates the number having its n digits equal to d:

```let rec rep d n = if n == 0 then 0 else d + 10 * (rep d (pred n))
```

(Yes, the same code can solve the four fours puzzle, or the nine nines puzzle, or…) Another basic function generates a list of integers in a given `range`:

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

(This one is tail recursive just because.) Now, given an expression e, I will count it as valid by adding it to the list es of expressions if it is finite and positive. Furthermore, if it is valid and it has a rational square root, I will count the latter as valid too:

```let with_root e es =
let open Ratio in
let e = abs e in
if null_denominator (label e)
|| sign_ratio (label e) == 0 then es else
let q = sqr e in
if null_denominator (label q) then e :: es else q :: e :: es
```

Finally, the workhorse:

```let rec puzzle digit size =
List.fold_right (fun i ->
let j = size - i in
List.fold_right (fun e0 ->
List.fold_right (fun e1 ->
List.fold_right (fun op ->
with_root (op e0 e1)
) (puzzle digit j)
) (puzzle digit i)
) (range 1 (size - 1))
(with_root (con (rep digit size)) [])
```

It looks more complicated than it is, really; just a list comprehension in a different shape that forces one to read it from the outside in, alas. It all starts with the `con`stant "dd… d", together `with_`(its)`root` if it is valid, as a seed for the list of expressions. Then for each i in the `range` from 1 to `size - 1`, it recursively builds solutions e0 of size i and e1 of size j = `size - i`. For each binary operation op, it builds the expression `op e0 e1` and adds it to the resulting list of solutions, together `with_`(its)`root` if they are valid.

Almost done! The problem is that there could be many possible expressions for a given number; it would be best to find just one exemplar for it. I've written about `group`ing lists in another era:

```let rec group ?(by=(=)) =
let rec filter e l = match l with
| [] -> [], []
| x :: xs ->
if not (by e x) then [], l else
let gs, ys = filter e xs in
x :: gs, ys
in function
| [] -> []
| x :: xs ->
let gs, ys = filter x xs in
(x :: gs) :: group ~by ys
```

A bit of syntax will let me build a filtering pipeline to select the best candidates:

```let (|>) x f = f x
```

Best in this case means that, out of all the expressions evaluating to a given integer, I prefer the one having the shortest representation. Now from all expressions I `filter` those that have integral value, decorate the expression as a string with its value, sort them and group them by value and select the first:

```let fournines =
let cmp (n, s) (n', s') =
let c = Pervasives.compare n n' in
if c != 0 then c else
let c = Pervasives.compare (String.length s) (String.length s') in
if c != 0 then c else
Pervasives.compare s s'
in puzzle 9 4
|> List.filter (fun e ->  Ratio.is_integer_ratio (label e))
|> List.map    (fun e -> (Ratio.int_of_ratio (label e), format string_of_int e))
|> List.sort cmp
|> group ~by:(fun (n, _) (n', _) -> n == n')
|> List.map List.hd
```

Very operational. It only remains to format the list to standard output:

```let () = List.iter (fun (n, s) -> Printf.printf "%4d = %s\n" n s) fournines
```

Behold, in all its glory, the solution to the puzzle:

1 2 1 = 99 / 99 2 = 99 / 9 - 9 3 = (9 + 9 + 9) / 9 4 = 9 / 9 + 9 / √9 5 = 9 - 9 / 9 - √9 6 = 9 * 9 / 9 - √9 7 = 9 - (9 + 9) / 9 8 = 99 / 9 - √9 9 = √(99 - 9 - 9) 10 = (99 - 9) / 9 11 = (9 + 9) / 9 + 9 12 = (9 + 99) / 9 13 = 9 + 9 / 9 + √9 14 = 99 / 9 + √9 15 = 9 + 9 - 9 / √9 17 = 9 + 9 - 9 / 9 18 = 99 - 9 * 9 19 = 9 + 9 + 9 / 9 20 = 9 + 99 / 9 21 = 9 + 9 + 9 / √9 24 = 99 / √9 - 9 26 = 9 * √9 - 9 / 9 27 = 9 * 9 * √9 / 9 28 = 9 * √9 + 9 / 9 30 = (99 - 9) / √9 32 = (99 - √9) / √9 33 = 99 * √9 / 9 34 = (99 + √9) / √9 36 = 9 + 9 + 9 + 9 39 = 9 * √9 + 9 + √9 42 = 9 + 99 / √9 45 = 9 * √9 + 9 + 9 51 = (9 + 9) * √9 - √9 54 = 9 * 9 - 9 * √9 57 = (9 + 9) * √9 + √9 63 = 9 * 9 - 9 - 9 69 = 9 * 9 - 9 - √9 72 = 99 - 9 * √9 75 = 9 * 9 - 9 + √9 78 = 9 * 9 - 9 / √9 80 = 9 * 9 - 9 / 9 81 = 99 - 9 - 9 82 = 9 * 9 + 9 / 9 84 = 9 * 9 + 9 / √9 87 = 99 - 9 - √9 90 = (9 + 9 / 9) * 9 93 = 99 - 9 + √9 96 = 99 - 9 / √9 98 = 99 - 9 / 9 99 = 9 * 99 / 9 100 = 9 / 9 + 99 102 = 9 / √9 + 99 105 = 9 + 99 - √9 108 = 99 + √(9 * 9) 111 = 999 / 9 117 = 9 + 9 + 99 126 = 9 * √9 + 99 135 = (9 + 9 - √9) * 9 144 = (9 + √9) * (9 + √9) 153 = (9 + 9) * 9 - 9 159 = (9 + 9) * 9 - √9 162 = 9 * 9 + 9 * 9 165 = (9 + 9) * 9 + √9 171 = (9 + 9) * 9 + 9 180 = 9 * 9 + 99 189 = (9 + 9 + √9) * 9 198 = 99 + 99 216 = (9 * 9 - 9) * √9 234 = 9 * 9 * √9 - 9 240 = 9 * 9 * √9 - √9 243 = (9 + 9 + 9) * 9 246 = 9 * 9 * √9 + √9 252 = 9 * 9 * √9 + 9 270 = (99 - 9) * √9 288 = 99 * √9 - 9 294 = 99 * √9 - √9 297 = 9 * 99 / √9 300 = 99 * √9 + √9 306 = 9 + 99 * √9 324 = (9 + 99) * √9 333 = 999 / √9 486 = (9 + 9) * 9 * √9 594 = (9 - √9) * 99 648 = (9 * 9 - 9) * 9 702 = (9 * 9 - √9) * 9 720 = 9 * 9 * 9 - 9 726 = 9 * 9 * 9 - √9 729 = 9 * 9 * √(9 * 9) 732 = 9 * 9 * 9 + √9 738 = 9 * 9 * 9 + 9 756 = (9 * 9 + √9) * 9 810 = (99 - 9) * 9 864 = (99 - √9) * 9 882 = 9 * 99 - 9 888 = 9 * 99 - √9 891 = 99 * √(9 * 9) 894 = 9 * 99 + √9 900 = 9 * 99 + 9 918 = (99 + √9) * 9 972 = (9 + 99) * 9 990 = 999 - 9 996 = 999 - √9 1002 = 999 + √9 1008 = 9 + 999 1188 = (9 + √9) * 99 1458 = (9 + 9) * 9 * 9 1782 = (9 + 9) * 99 2187 = 9 * 9 * 9 * √9 2673 = 9 * 99 * √9 2997 = 999 * √9 6561 = 9 * 9 * 9 * 9 8019 = 9 * 9 * 99 8991 = 9 * 999 9801 = 99 * 99 9999 = 9999

(The table is built out of the actual output to Mac OS Terminal. The Unicode characters are printed perfectly.) Using factorials to fill in the gaps is left as an exercise to the reader (it is not simple).

## 2011-12-29

### Vose's Alias Method

This is a note regarding the implementation of Vose's Method to construct the alias tables for non-uniform discrete probability sampling presented by Keith Schwartz (as of 2011-12-29). Mr. Schwartz otherwise very informative write-up contains an error in the presentation of the Method. Step 5 of the algorithm fails if the `Small` list is nonempty but the `Large` list is. This can happen either:

1. if the initial probabilities sum to less than 1, or
2. if updating in step 5.5 the `Large` probability results in cancellation so that pg results in an approximation less than its true value

The first case can be mitigated by scaling by n / ∑ pj in step 3. The second case can be mitigated by replacing 5.5 by "Set pg := (pg + pl) - 1." In any case, there is no proof that the algorithm as shown is correct, which is not the case for Vose's presentation. In the interest of correctness and completeness, here is an OCaml implementation of Vose's Method:

```type alias = { n : int; prob : float array; alias : int array; }

let alias pa =
let rec split pa n j (small, large as part) =
if j == n then part else
if pa.(j) > 1.
then split pa n (succ j) (     small, j :: large)
else split pa n (succ j) (j :: small,      large)
in
let rec init r pa part = match part with
| j :: small, k :: large ->
r.prob.(j)  <- pa.(j);
r.alias.(j) <- k;
pa.(k)      <- (pa.(k) +. pa.(j)) -. 1.;
if pa.(k) > 1.
then init r pa (small, k :: large)
else init r pa (k :: small, large)
| j :: small, []         ->
r.prob.(j) <- 1.;
init r pa (small, [])
| []        , k :: large ->
r.prob.(k) <- 1.;
init r pa ([], large)
| []        , []         -> r
in
let n  = Array.length pa in
if n == 0 then invalid_arg "alias" else
let sc = float n /. Array.fold_left (fun s p ->
if p < 0. then invalid_arg "alias" else
s +. p) 0. pa in
let sa = Array.map (( *. ) sc) pa in
let r  = { n; prob = Array.make n 0.; alias = Array.make n (-1); } in
init r sa (split sa n 0 ([], []))

let choose r =
let p, e = modf (Random.float (float r.n)) in
let j = truncate e in
if p <= r.prob.(j) then j else r.alias.(j)
```

Since this algorithm is recursive, it might not be immediately obvious how to translate it to more mainstream languages. A perhaps more faithful rendition in rather low-level Java is the following:

```import java.util.Random;

public class Vose {
private final Random random;
private final int limit;
private final double[] prob;
private final int[] alias;

public Vose(final double[] pa) {
this(pa, new Random());
}

public int getLimit() {
return limit;
}

public Vose(final double[] pa, final Random random) {
final int limit = pa.length;
if (limit == 0)
throw new IllegalArgumentException("Vose");
double sum = 0;
for (int j = 0; j != limit; j++) {
if (pa[j] < 0)
throw new IllegalArgumentException("Vose");
sum += pa[j];
}
final double scale = (double) limit / sum;
final double[] sa = new double[limit];
for (int j = 0; j != limit; j++)
sa[j] = pa[j] * scale;

this.random = random;
this.limit  = limit;
this.prob   = new double[limit];
this.alias  = new int[limit];
init(sa);
}

private void init(final double[] sa) {
final int[] small = new int[limit];
final int[] large = new int[limit];
int ns = 0;
int nl = 0;
for (int j = 0; j != limit; j++)
if (sa[j] > 1)
large[nl++] = j;
else
small[ns++] = j;
while (ns != 0 && nl != 0) {
final int j = small[--ns];
final int k = large[--nl];
prob[j]  = sa[j];
alias[j] = k;
sa[k]    = (sa[k] + sa[j]) - 1; // sic
if (sa[k] > 1)
large[nl++] = k;
else
small[ns++] = k;
}
while (ns != 0)
prob[small[--ns]] = 1;
while (nl != 0)
prob[large[--nl]] = 1;
}

public int choose() {
final double u = limit * random.nextDouble();
final int    j = (int) Math.floor(u);
final double p = u - (double) j;
return p <= prob[j] ? j : alias[j];
}
}
```

In all generality, it is sufficient that the "probabilities" be non-negative, since dividing the values by the mean normalizes them so that the conditions of the Method are satisfied.

## References

1. Michael D. Vose. A Linear Algorithm For Generating Random Numbers With a Given Distribution, IEEE TRANSACTIONS ON SOFTWARE ENGINEERING VOL. 17, NO. 9 SEPTEMBER 1991. Online.

## 2011-12-21

### A Better (Gauss) Error Function

If you do statistics you know of `erf` and `erfc`; if you work in OCaml you surely miss them. It is not very difficult to port the canonical implementation given by Numerical Recipes (which I won't show and not just for licensing reasons); if you Google for a coefficient you'll see that this approximation is, indeed, ubiquitous. There exists a better approximation in the literature, one that is more than 40 years old. Since I find it inconceivable that it is not more widely used (again, Google for a coefficient), I present a straightforward implementation of it.

This approximation boasts a relative error of 10-19.5 in the range |x| ≤ 0.46875 = 15/32, of 10-18.59 in the range 0.46875 ≤ |x| ≤ 4, and of 10-18.24 in the range 4 ≤ |x|. The only difficulty is to accurately transcribe the coefficient tables into code; I've double-checked them and ran some tests to assess the proper functioning of the code:

```let r0 = [|
3.20937_75891_38469_47256_2e+03, 2.84423_68334_39170_62227_3e+03;
3.77485_23768_53020_20813_7e+02, 1.28261_65260_77372_27564_5e+03;
1.13864_15415_10501_55649_5e+02, 2.44024_63793_44441_73305_6e+02;
3.16112_37438_70565_59694_7e+00, 2.36012_90952_34412_09349_9e+01;
1.85777_70618_46031_52673_0e-01, 1.00000_00000_00000_00000_0e+00;
|]

let r1 = [|
1.23033_93547_97997_25272e+03, 1.23033_93548_03749_42043e+03;
2.05107_83778_26071_46532e+03, 3.43936_76741_43721_63696e+03;
1.71204_76126_34070_58314e+03, 4.36261_90901_43247_15820e+03;
8.81952_22124_17690_90411e+02, 3.29079_92357_33459_62678e+03;
2.98635_13819_74001_31132e+02, 1.62138_95745_66690_18874e+03;
6.61191_90637_14162_94775e+01, 5.37181_10186_20098_57509e+02;
8.88314_97943_88375_94118e+00, 1.17693_95089_13124_99305e+02;
5.64188_49698_86700_89180e-01, 1.57449_26110_70983_47253e+01;
2.15311_53547_44038_46343e-08, 1.00000_00000_00000_00000e+00;
|]

let r2 = [|
-6.58749_16152_98378_03157e-04, 2.33520_49762_68691_85443e-03;
-1.60837_85148_74227_66278e-02, 6.05183_41312_44131_91178e-02;
-1.25781_72611_12292_46204e-01, 5.27905_10295_14284_12248e-01;
-3.60344_89994_98044_39429e-01, 1.87295_28499_23460_47209e+00;
-3.05326_63496_12323_44035e-01, 2.56852_01922_89822_42072e+00;
-1.63153_87137_30209_78498e-02, 1.00000_00000_00000_00000e+00;
|]
```

The tables are laid out as pairs of coefficients (pi, qi) for the numerator and denominator polynomials, respectively, in ascending order of monomial power (that is, the independent coefficients are in the first position of the arrays). The rational functions can be evaluated with a suitable variation of Horner's schema:

```let horner2 r x =
let n = Array.length r in
let s = ref 0.
and t = ref 0. in
for i = n - 1 downto 0 do
let p, q = Array.unsafe_get r i in
s := !s *. x +. p;
t := !t *. x +. q
done;
!s /. !t
```

For clarity I haven't inlined Horner's recursion in the following code; this is not very difficult to do, so I've left it as an exercise for the implementor:

```let iqpi = 5.64189_58354_77562_86948_1e-01

let erfc x =
let z  = abs_float x in
let z2 = z *. z in
let y =
if z < 0.46875 then   1. -.         z   *. horner2 r0 z2 else
if z < 4.      then         exp (-. z2) *. horner2 r1 z  else
let z'  = 1. /. z  in
let z'2 = z' *. z' in       exp (-. z2) *. z' *. (iqpi +. z'2 *. horner2 r2 z'2)
in if x < 0. then 2. -. y else y
```

As explained in the paper, `erfc x` must be everywhere equal to `1. -. erf x`:

```let erf x = 1. -. erfc x
```

If you have Mathematica™ lying around, you might want to output a list of points to test the accuracy of the approximation:

```let mma f x0 x1 n =
let buf = Buffer.create 256 in
(1 - 2 BitAnd[1, BitShiftRight[#, 63]])\
2^(BitAnd[16^^7ff, BitShiftRight[#, 52]] - 1023)\
(1 +  2^-52 BitAnd[16^^fffffffffffff, #]),\
\$MachinePrecision]&, #, {2}]&@{";
for i = 0 to n do
let x = (x0 *. float (n - i) +. x1 *. float i) /. float n in
let y = f x in
if i != 0 then Buffer.add_char buf ',';
Printf.bprintf buf "{16^^%Lx,16^^%Lx}"
(Int64.bits_of_float x)
(Int64.bits_of_float y)
done;
Buffer.contents buf
```

In order to exactly marshal the IEEE 754 double values I convert them to rational numbers equivalent to machine floats on Mathematica's end. My informal tests show that around the cut-off points Mathematica evaluates `erf` to the same exact values this code does (that is, the absolute error is zero). When evaluated to sufficient precision, the maximum absolute error among 5,000 equispaced points in the interval [3.99, 4.01] is 5.5437×10-17.

To mathematical libraries implementors everywhere, I urge you: please steal this code!

## References

1. W. J. Cody. Rational Chebyshev Approximations for the Error Function, Mathematics of Computation Vol. 23, No. 107 (Jul., 1969), pp. 631-63. Online.