Math Ace: the plot thickens

  • Difficulty: 025/100
  • Author: jeremy theler
  • Keywords: PARAMETRIC, DIMENSIONS, TYPE, DEFAULT_ARGUMENT_VALUE, OUTER_STEPS, MIN, MAX, OFFSET, IF, ENDIF, :=, sobol, niederreiter, halton, abs,

But for me, the pursuit of money has always played second fiddle to

the pursuit of beauty. Thank goodness for mathematics and

computer science, whose methods offer many ways by which our left

brains and our right brains can both be happy, simultaneously.

Donald. E. Knuth, Selected Papers on Fun & Games, 2011

Mankind’s primary concern is beauty. Call it love, art or science. Whatever. “Donald E. Knuth”: “http://en.wikipedia.org/wiki/Donald_Knuth” knows it. And it is himself that shows us one of the most aesthetically pleasant examples of mathematical and computer graphics I have ever seen. Back in 1959 he asked its readers to predict what a certain equation meant. Please try to figure out what the result is before turning the page or scrolling down. Believe me, it is worth the shot:

!bt \begin{align*} \left( \Big| \big| 3-|x| \big| - 3 + |x| \Big| + \bigg| \Big| \sqrt{ |9-x^2| } - \big| y-2/3 |x| \big| \Big| - \sqrt{|9-x^2|} + \big| y - 2/3 |x| \big| \bigg| \right) \quad \\ \Bigg[ \big| |x y| + xy \big| + \left( \bigg| \Big| 2 - \big|33 - 3 |x| \big| \Big| - 2 + \big| 33 - 3 |x| \big| \bigg| + \big| 14 - |y| \big| \right) \quad\quad \\ \left( \Big| 16 - |y| - 3 \big|11 - |x| \big| \Big| \Bigg| \bigg| \sqrt{\big| 1 - (11 - |x|)^2 \big|} - \Big| 11 - |y| + 2/3 \big|11 - |x|\big| \Big| \bigg| \right. \quad\quad\quad \\ \left. - \sqrt{\big| 1 - (11 - |x|)^2 \big|} + \Big| 11 - |y| + 2/3 \big|11 - |x|\big| \Big| \Bigg| \right. \quad\quad\quad \\ \left. + \bigg|\Big| 1 - \big|11 - |x|\big| \Big| - 1 + \big|11 - |x|\big|\bigg| \right) \Bigg] &= 0 \end{align*}

!et

The absolute value bars may be dazzling, but the “wasora”: “http://www.talador.com.ar/jeremy/wasora” input below shows the actual implementation of the left hand side of the equation as a function \(f(x,y)\) where the arguments to the abs function are clearly enclosed in parentheses. If this is not love, art and science, I wonder what we are doing in this tiny planet.

mathace.was

The following input defines the left hand side of the mystery equation as a function \(f(x,y)\) and evaluates it by sweeping both \(x\) and \(y\) in the two-dimensional interval \([-13,13]\times[-18,18]\). It does so by sampling \((x,y)\) pairs from a quasi-random number generator. By default, it uses the sobol sequence (see the “GSL’s Quasi-random number generator algorithms”: “http://www.gnu.org/software/gsl/manual/html_node/Quasi_002drandom-number-generator-algorithms.html#Quasi_002drandom-number-generator-algorithms”), but niederreiter, halton and reversehalton may be chosen. Moreover, wasora allows an offset to be given, so finer sweepings of the parametric space can be done without needing to recompute the coarse sequence. The reported running times in the terminal mimic shows the difference between incremental and full computations.

If you are not extremely impatient and curious, please proceed as follows. First, take a look at the proposed input below but do not advance to see its output. Try to run the example by yourself and then fire up gnuplot or your favorite plotting program and draw the result. If you, as me, cannot wait to find out what the equation means, advance two pages. In any case, for sure you will then go back to the previous page and take a new look at the equation from a different point of view. Beautiful, isn’t it?

# Knuth's mystery equation
f(x,y) := {
 ( abs(abs(3-abs(x)) - 3 + abs(x)) + abs(abs(sqrt(abs(9-x^2)) - abs(y-2/3*abs(x))) - sqrt(abs(9-x^2)) + abs(y - 2/3*abs(x)))) *
 ( abs(abs(x*y) + x*y) +
   ( abs(abs(2 - abs(33-3*abs(x))) - 2 + abs(33-3*abs(x))) + abs(14 - abs(y)) ) *
   ( abs(16 - abs(y) - 3*abs(11 - abs(x))) *
      abs(abs(sqrt(abs(1 - (11 - abs(x))^2)) - abs(11 - abs(y) + 2/3*abs(11 - abs(x)))) - sqrt(abs(1 - (11 - abs(x))^2)) + abs(11 - abs(y) + 2/3*abs(11 - abs(x))))
      + abs(abs(1 - abs(11 - abs(x))) - 1 + abs(11 - abs(x))) )
 )
}

DEFAULT_ARGUMENT_VALUE 1 sobol
DEFAULT_ARGUMENT_VALUE 2 16
DEFAULT_ARGUMENT_VALUE 3 0

# the range to sweep
PARAMETRIC x y TYPE $1 MIN -13 -18 MAX 13 18 OUTER_STEPS 2^$2 OFFSET 2^$3-1
                 
IF abs(f(x,y))<1
 PRINT %f x y
ENDIF
$ time wasora mathace.was niederreiter  > mathace-niederreiter.dat

real    0m0.365s
user    0m0.340s
sys 0m0.020s
$ time wasora mathace.was sobol         > mathace-sobol.dat

real    0m0.366s
user    0m0.356s
sys 0m0.008s
$ time wasora mathace.was halton        > mathace-halton.dat

real    0m0.368s
user    0m0.352s
sys 0m0.012s
$ time wasora mathace.was reversehalton > mathace-reversehalton.dat

real    0m0.385s
user    0m0.364s
sys 0m0.020s
$ time wasora mathace.was reversehalton 18 16 > mathace-reversehalton2.dat

real    0m1.158s
user    0m1.096s
sys 0m0.056s
$ time wasora mathace.was reversehalton 19 17 > mathace-reversehalton3.dat

real    0m2.294s
user    0m2.180s
sys 0m0.108s
$ time wasora mathace.was reversehalton 19 > mathace-reversehalton-fine.dat

real    0m2.998s
user    0m2.900s
sys 0m0.092s
$ gnuplot4 mathace.gp
$ 

FIGURE: [025-mathace/mathace-niederreiter] The Math Ace solved with a \(256 \times 256\) niederreiter sequence

FIGURE: [025-mathace/mathace-sobol] The Math Ace solved with a \(256 \times 256\) sobol sequence

FIGURE: [025-mathace/mathace-halton] The Math Ace solved with a \(256 \times 256\) halton sequence

FIGURE: [025-mathace/mathace-reversehalton] The Math Ace solved with a \(256 \times 256\) reversehalton sequence

FIGURE: [025-mathace/mathace-reversehalton2] The Math Ace solved with a \(512 \times 512\) reversehalton sequence, starting from the \(256\times256\) solution

FIGURE: [025-mathace/mathace-reversehalton3] The Math Ace solved with a \(768 \times 768\) reversehalton sequence, starting from both the \(256\times256\) and \(512\times512\) solutions

FIGURE: [025-mathace/mathace-reversehalton-fine] The Math Ace solved with a \(768 \times 768\) reversehalton sequence, starting from scratch

Note that the running time is fairly high due to the fact that wasora’s sometime-to-be-replaced-algebraic parser and evaluator does not work in a tree-based way as proposed by Knuth in is “LR Parser”: “http://en.wikipedia.org/wiki/LR_parser” algorithm of 1965, but in a easier-to-write-but-unefficient multilevel array organization.