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 =
      Buffer.add_string buf op;
      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;
      Buffer.add_char   buf ' ';
      Buffer.add_string buf op;
      Buffer.add_char   buf ' ';
      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 ratios, 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 add e e' = Add (Ratio.add_ratio    (label e) (label e'), 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)
        ) [add; sub; mul; div]
      ) (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 constant "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 grouping 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:


11 =99 / 99
22 =99 / 9 - 9
33 =(9 + 9 + 9) / 9
44 =9 / 9 + 9 / √9
55 =9 - 9 / 9 - √9
66 =9 * 9 / 9 - √9
77 =9 - (9 + 9) / 9
88 =99 / 9 - √9
99 =√(99 - 9 - 9)
1010 =(99 - 9) / 9
1111 =(9 + 9) / 9 + 9
1212 =(9 + 99) / 9
1313 =9 + 9 / 9 + √9
1414 =99 / 9 + √9
1515 =9 + 9 - 9 / √9
1617 =9 + 9 - 9 / 9
1718 =99 - 9 * 9
1819 =9 + 9 + 9 / 9
1920 =9 + 99 / 9
2021 =9 + 9 + 9 / √9
2124 =99 / √9 - 9
2226 =9 * √9 - 9 / 9
2327 =9 * 9 * √9 / 9
2428 =9 * √9 + 9 / 9
2530 =(99 - 9) / √9
2632 =(99 - √9) / √9
2733 =99 * √9 / 9
2834 =(99 + √9) / √9
2936 =9 + 9 + 9 + 9
3039 =9 * √9 + 9 + √9
3142 =9 + 99 / √9
3245 =9 * √9 + 9 + 9
3351 =(9 + 9) * √9 - √9
3454 =9 * 9 - 9 * √9
3557 =(9 + 9) * √9 + √9
3663 =9 * 9 - 9 - 9
3769 =9 * 9 - 9 - √9
3872 =99 - 9 * √9
3975 =9 * 9 - 9 + √9
4078 =9 * 9 - 9 / √9
4180 =9 * 9 - 9 / 9
4281 =99 - 9 - 9
4382 =9 * 9 + 9 / 9
4484 =9 * 9 + 9 / √9
4587 =99 - 9 - √9
4690 =(9 + 9 / 9) * 9
4793 =99 - 9 + √9
4896 =99 - 9 / √9
4998 =99 - 9 / 9
5099 =9 * 99 / 9
51100 =9 / 9 + 99
52102 =9 / √9 + 99
53105 =9 + 99 - √9
54108 =99 + √(9 * 9)
55111 =999 / 9
56117 =9 + 9 + 99
57126 =9 * √9 + 99
58135 =(9 + 9 - √9) * 9
59144 =(9 + √9) * (9 + √9)
60153 =(9 + 9) * 9 - 9
61159 =(9 + 9) * 9 - √9
62162 =9 * 9 + 9 * 9
63165 =(9 + 9) * 9 + √9
64171 =(9 + 9) * 9 + 9
65180 =9 * 9 + 99
66189 =(9 + 9 + √9) * 9
67198 =99 + 99
68216 =(9 * 9 - 9) * √9
69234 =9 * 9 * √9 - 9
70240 =9 * 9 * √9 - √9
71243 =(9 + 9 + 9) * 9
72246 =9 * 9 * √9 + √9
73252 =9 * 9 * √9 + 9
74270 =(99 - 9) * √9
75288 =99 * √9 - 9
76294 =99 * √9 - √9
77297 =9 * 99 / √9
78300 =99 * √9 + √9
79306 =9 + 99 * √9
80324 =(9 + 99) * √9
81333 =999 / √9
82486 =(9 + 9) * 9 * √9
83594 =(9 - √9) * 99
84648 =(9 * 9 - 9) * 9
85702 =(9 * 9 - √9) * 9
86720 =9 * 9 * 9 - 9
87726 =9 * 9 * 9 - √9
88729 =9 * 9 * √(9 * 9)
89732 =9 * 9 * 9 + √9
90738 =9 * 9 * 9 + 9
91756 =(9 * 9 + √9) * 9
92810 =(99 - 9) * 9
93864 =(99 - √9) * 9
94882 =9 * 99 - 9
95888 =9 * 99 - √9
96891 =99 * √(9 * 9)
97894 =9 * 99 + √9
98900 =9 * 99 + 9
99918 =(99 + √9) * 9
100972 =(9 + 99) * 9
101990 =999 - 9
102996 =999 - √9
1031002 =999 + √9
1041008 =9 + 999
1051188 =(9 + √9) * 99
1061458 =(9 + 9) * 9 * 9
1071782 =(9 + 9) * 99
1082187 =9 * 9 * 9 * √9
1092673 =9 * 99 * √9
1102997 =999 * √9
1116561 =9 * 9 * 9 * 9
1128019 =9 * 9 * 99
1138991 =9 * 999
1149801 =99 * 99
1159999 =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
  Buffer.add_string buf "l = Map[N[\
(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.add_string buf "};";
  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.