- 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.

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.