Nice!

The back half of my oven is ~25F too hot. Gotta rotate early and often.

No, sorry. I've never played the game. My friend and a few other people told me about the community, I just relayed what I heard.

Okay, technical details.

**Calculating roots**

There are multiple numerical methods to find complex roots of polynomials. I'm aware about Durand-Kerner method (https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method) and Dandelin-Graeffe method (https://en.wikipedia.org/wiki/Graeffe's_method), where "aware" means "I know that they exist". I did not implement root finding myself; instead I've used GNU scientific library (http://www.gnu.org/software/gsl/) and built-in function of it called gsl_poly_complex_solve (http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations). The documentation claims that they use "balanced-QR reduction of the companion matrix" and I have no idea what is it and was not aware about it existence before. Since I wrote my code in OCaml, I've used OCaml bindings for GSL (from https://github.com/mmottl/gsl-ocaml).

**Rendering**

After the roots are calculated, I transform them to pixel coordinates, making sure that the whole thing scaled to the given image size (I must admit that I crop out some purely real roots, nothing interesting there). Then, according to pixel coordinates, I write 'FFFFFF' (i.e. "white") bytes to the corresponding positions in an array. This array backs up Cairo (http://cairographics.org/) surface of the 20000x20000 size, and Cairo PNG module is used to write the result to a file. Again, I needed OCaml bindings (https://github.com/Chris00/ocaml-cairo).

This approach limits the resolution of the image; for example, Cairo fails to create surface of 40000x40000 size. See below for future plans.

**Performance**

Most of the time is taken by calculating the roots and scaling them to fill the image of the given size; contribution of drawing and writing png file is noticeable but insignificant. On Macbook Pro with 3GHz Intel Core i7 the whole process takes about 40 minutes; drawing the picture and writing it to file takes about 2-3 minutes.

**Criticism**

- I preferred convenience in writing to optimization, so there's a lot of work for garbage collector when it kicks in -- the pauses are quite noticeable.

- The whole thing is written in a pretty dumb way, it runs in one thread only, so no multi-CPU multi-core benefits.

- Everything is kept in memory until image rendered. As as result I cannot do degree 25 or more, as I'm running out of memory (I have 16 Gb of RAM). Memory usage definitely can be optimized. See "convenience over optimization" point above.

**Future plans**

I have started to rewrite everything in a parallel fashion. I also plan to do rendering in chunks (like "top left quarter image, top right quarter image, bottom left quarter image, bottom right quarter image") -- this would allow me to have much higher resolution. When this is done, I'll rent an hour from some Amazon high-memory high-CPU server and will do much better resolution and higher degree.

**Minor disappointment**

I have expected that drawing the picture of degree 23 on top of the picture of degree 24 would be interesting. Not so. Here is the degree 23 (in green) on top of degree 24 (in cyan):

Green is basically distributed everywhere -- more in the middle of the ring, but also in all interesting fractal-like parts.

Here is degree 22 (in orange) on top of degree 23 (in green) on top of degree 24 (in cyan):

Same story. Adding more degrees creates the mess of colored pixels.

But what's clear (from bigger pictures which I do not share) is that my resolution is too rough to see how exactly degree 23 roots overlap degree 24 roots in interesting areas. Stay tuned :)

- I cannot do degree 25 or more, as I'm running out of memory

The image appears to have horizontal and vertical symmetry, but I haven't been able to have a good look at it without my browser overheating.

If it is symmetric, could you calculate one quadrant and then reflect it to build the whole image?

Well, the picture *is* symmetric. Vertical symmetry comes from the fact that if z is a root of a polynomial with real coefficients then the complex conjugate z' is also a root of the same polynomial (if z = a + bi then the complex conjugate z' is defined as z = a - bi). The horizontal symmetry comes from the fact that if z is a root of a polynomial, then -z is a root of polynomial which can be constructed by negating the sign of each coefficient next to each odd degree. As I'm looking at all polynomials with coefficients +1 or -1, for every root z of some polynomial I have a root -z for some other polynomial.

This indeed gives some possibilities:

- I can draw only 1/4 of the picture, and construct the result out of it

- I can store only 1/4 of all the roots (only top right quadrant, for example)

- I can do calculations for only 1/2 of all the polynomials (but I need some logic to keep track if I need to find the roots of the polynomial or not -- would be simple, though)

- I cannot calculate only half of the roots, because, as I've mentioned above, vertically symmetric roots are coming in pairs in the same polynomial, and the method used by GSL library for finding the roots does not produce roots one by one, instead calculating the whole set of roots at the same time.

All in all, these optimizations would save about 50% of the calculation time and about 75% of memory, also I'd be able to double the resolution. I should be able to do the power 25, and, maybe 26. But not 27. And I dream about doing 32, and increasing the resolution of the picture ten-fold. So I'll rather at first focus on doing the whole thing in more parallel manner, with results not being stored in memory -- but yes, symmetry-based optimization is definitely one of the things to do after that. Thanks for suggestion!

Eh, not really a math puzzle, more of an exercise in vigilance. I always feel indignant after seeing this sort of thing: it seems deliberately designed to insult the player's intelligence.

The thing is, life is full of these kinds of problems. Advertisements, speeches, statistics, news stories, and so on. And unlike this little bugger, _those_ problems aren't harmless. They _are_ designed to make a real fool out of you, usually to your cost. So there is a lesson to be learned here. It's just not about math.

So a couple of years ago I was staying in LA with a friend. I was wondering if the tap water in LA is safe to drink, so I did a search from our hotel room, and found this: "Yes, it's fine. Unless you are staying at the Cecil hotel of course." What do you know... we were staying at the Hotel Cecil. I was pretty confused until I looked up the story.

I drank the tap water. I was fine.

Fantastic! I just found a couple of papers in seconds, each of which took me several minutes to find before -- and I have access to a university database.

Here's the study: Sex beyond the genitalia: The human brain mosaic

The findings are interesting, but the conclusion seems unwarranted (even sensational.) Sure, there may be a large overlap between the male distribution and the female distribution for each sexually dimorphic structural feature, but that doesn't imply that you can't identify the sex of a given brain with very high accuracy by considering all such features, thus categorizing brains into two very distinct groups.

- but the conclusion seems unwarranted (even sensational.)

Of course. It's a modern study about gender.

Thumbs up for Cryptonomicon and Anathem.

Another sf recommendation: Wool by Hugh Howey. Not quite free, but well worth it.

Yeah, it's tough, I don't really know how to go about it. But here's at least a demonstration for three denominations, {1, 2, 5}: (n is the total amount to add up to)

`def onetwofive(n):`

numfives = math.floor(n/5.)+1

`return int( numfives * (n % 5 + n + 3)/4 - (numfives % 2) * (n % 2 - 0.5) / 2. )`

Disgusting, I know! But it works: I double checked it against your code, it gives the same result up to n=5000.

An estimate for four denominations {1, 2, 5, 10} is n^3/600; for five denominations {1, 2, 5, 10, 25}, n^4/60000. I got these by summing crudely and dropping low order terms. They seem to be asymptotic: the estimate for five denominations is accurate to within 5% of the true value for n > 2000. That might seem silly but count_change does get quite slow for n >> 10^6 (that's $10,000.00), while the approximation is good to within 0.01%.

edit: You could write a program that returns an estimator function for a given set of denominations. See Faulhaber's formula.

I'm impressed you managed a solution for 3 coins.

> count_change does get quite slow for n >> 10^6 (that's $10,000.00)

From a pure math angle, both memo & dynamic programming versions are linear in time and space. As you said though, those numbers are growing exponentially meaning the space will grow exponentially too.

You're right. I hadn't quite thought it through, but it would work well: eg. for 8 coin values, if you need to calculate the number of ways to make 2.00, you need to make at most 200*8=1600 calls to the recursive function. So you get O(n) efficiency.

But it still seems like you should be able to do better: O(1) would be possible with a closed-form formula.

edit: Oh, I noticed you actually implemented it with memoization. Great!

- O(1) would be possible with a closed-form formula.

That's way beyond my maths ability. Do you think you could find such a formula?

I know there's Binet's formula for Fibonacci, but it's impractical for larger numbers as the floating point errors start to affec tthe results (on most machines anyway.)

Yeah, it's tough, I don't really know how to go about it. But here's at least a demonstration for three denominations, {1, 2, 5}: (n is the total amount to add up to)

`def onetwofive(n):`

numfives = math.floor(n/5.)+1

`return int( numfives * (n % 5 + n + 3)/4 - (numfives % 2) * (n % 2 - 0.5) / 2. )`

Disgusting, I know! But it works: I double checked it against your code, it gives the same result up to n=5000.

An estimate for four denominations {1, 2, 5, 10} is n^3/600; for five denominations {1, 2, 5, 10, 25}, n^4/60000. I got these by summing crudely and dropping low order terms. They seem to be asymptotic: the estimate for five denominations is accurate to within 5% of the true value for n > 2000. That might seem silly but count_change does get quite slow for n >> 10^6 (that's $10,000.00), while the approximation is good to within 0.01%.

edit: You could write a program that returns an estimator function for a given set of denominations. See Faulhaber's formula.

I'm impressed you managed a solution for 3 coins.

> count_change does get quite slow for n >> 10^6 (that's $10,000.00)

From a pure math angle, both memo & dynamic programming versions are linear in time and space. As you said though, those numbers are growing exponentially meaning the space will grow exponentially too.

I remember reading SICP. It's a lot of fun.

The recursive algorithm reminds me of something I was just looking at the other day, the Chu-Vandermonde identity.

One obvious way to speed it up is memoization. But I wonder if that would be as effective in this case as for Fibonacci.

Another way would be to just handle the penny case by just outputting 1 straight away if amount >= 0, and handle the penny-and-nickel case by outputting 1+floor(amount/5). I'd imagine that just doing this would speed up the program quite a bit.

But a straight-up, non-recursive calculation in one fell swoop? I'm not sure how to do that.

- But I wonder if that would be as effective in this case as for Fibonacci.

Sure, why not? It's still **tree recursion** after all. You just have to use an extra procedure that takes the change amount and how many coin values you have.

e.g. For the UK coin values are 1, 2, 5, 10, 20, 50, 100 and 200 so there are 8 coin values.

So instead of

` count_change(amount)`

you'd use

`count_change(amount)`

`cc(amount, num_coins)`

and you memoize cc with indexes of amount and num_coins

You're right. I hadn't quite thought it through, but it would work well: eg. for 8 coin values, if you need to calculate the number of ways to make 2.00, you need to make at most 200*8=1600 calls to the recursive function. So you get O(n) efficiency.

But it still seems like you should be able to do better: O(1) would be possible with a closed-form formula.

edit: Oh, I noticed you actually implemented it with memoization. Great!

- O(1) would be possible with a closed-form formula.

That's way beyond my maths ability. Do you think you could find such a formula?

I know there's Binet's formula for Fibonacci, but it's impractical for larger numbers as the floating point errors start to affec tthe results (on most machines anyway.)

Yeah, it's tough, I don't really know how to go about it. But here's at least a demonstration for three denominations, {1, 2, 5}: (n is the total amount to add up to)

`def onetwofive(n):`

numfives = math.floor(n/5.)+1

`return int( numfives * (n % 5 + n + 3)/4 - (numfives % 2) * (n % 2 - 0.5) / 2. )`

Disgusting, I know! But it works: I double checked it against your code, it gives the same result up to n=5000.

An estimate for four denominations {1, 2, 5, 10} is n^3/600; for five denominations {1, 2, 5, 10, 25}, n^4/60000. I got these by summing crudely and dropping low order terms. They seem to be asymptotic: the estimate for five denominations is accurate to within 5% of the true value for n > 2000. That might seem silly but count_change does get quite slow for n >> 10^6 (that's $10,000.00), while the approximation is good to within 0.01%.

edit: You could write a program that returns an estimator function for a given set of denominations. See Faulhaber's formula.

I'm impressed you managed a solution for 3 coins.

> count_change does get quite slow for n >> 10^6 (that's $10,000.00)

From a pure math angle, both memo & dynamic programming versions are linear in time and space. As you said though, those numbers are growing exponentially meaning the space will grow exponentially too.

Sorry, but I still only have a very hazy idea of what you're trying to say! It seems like a very slippery argument and I'm not sure that it holds water. Could you make it more rigorous?

The first claim is that each black square must be bombed before its adjacent white squares. I take this to mean that if S is a guaranteed-tank-destroying sequence of firing coordinates, and b is some particular black square with w some particular neighbor of b, then the first occurence of b in S must fall before the last occurence of w. (That is, b may be targeted after w as many times as one pleases, as long as w is targeted at least once after the first targeting of b.)

The next statement is: "bombing all blacks and then all whites is the same as bombing them in some other pattern that meets the same requirements and outputs the same number(for instance, bombing all whites in a row then bombing all blacks in a row.)" And here I'm completely lost: What does it mean, precisely, for two firing sequences to be "the same"? Does it mean that either they both are guaranteed to kill the tank or they both are not, and that they are also of the same length?~

well, keep in mind the guaranteed tank-destroying sequence is actually white-black-white, I was simplifying just so I didn't have to type that much. And yes, you got that part exactly right.

As for your question, I suppose for two sequences to be the same means that they both a)destroy the tanks and b)hit the same squares the same amount of times. Thus they'd have to follow the same W-B-W rules and could be morphed into each other so that a more computationally intensive process becomes static or vice-versa. Forgive me if this seems unimportant, I tend to think in terms of algorithms.

- if the tank starts on a black square, you must bomb all black squares and then all white squares. If the tank starts in a white square, you must bomb all white squares and then all black squares and then all white squares. Therefore the bombing run must be W-B-W.

In fact if the tank starts on a black square, you must merely bomb that ONE square before all of its neighbours. You say that it's necessary to bomb all the black squares first, use that claim as proof for a lower bound, then admit in the next paragraph that in fact bombing all the black squares first is unnecessary. So I still don't see how you've shown that it's impossible to destroy the tank in fewer shots. Perhaps by some clever interleaving of black and white squares you might save a missile or two?

What I'm trying to say is that all black squares have to be bombed before their adjacent white squares, so bombing all blacks and then all whites is the same as bombing them in some other pattern that meets the same requirements and outputs the same number(for instance, bombing all whites in a row then bombing all blacks in a row). The checkerboard is the smallest that meets the requirements because it pairs every black with one adjacent white (plus the extra one), so if you added one more black it would require that two pairs of blacks exist.

EDIT:

in fact, I suppose you could use that pairing scheme to always bomb any configuration of squares in at most nSquares*3/2 + remainders

Sorry, but I still only have a very hazy idea of what you're trying to say! It seems like a very slippery argument and I'm not sure that it holds water. Could you make it more rigorous?

The first claim is that each black square must be bombed before its adjacent white squares. I take this to mean that if S is a guaranteed-tank-destroying sequence of firing coordinates, and b is some particular black square with w some particular neighbor of b, then the first occurence of b in S must fall before the last occurence of w. (That is, b may be targeted after w as many times as one pleases, as long as w is targeted at least once after the first targeting of b.)

The next statement is: "bombing all blacks and then all whites is the same as bombing them in some other pattern that meets the same requirements and outputs the same number(for instance, bombing all whites in a row then bombing all blacks in a row.)" And here I'm completely lost: What does it mean, precisely, for two firing sequences to be "the same"? Does it mean that either they both are guaranteed to kill the tank or they both are not, and that they are also of the same length?~

well, keep in mind the guaranteed tank-destroying sequence is actually white-black-white, I was simplifying just so I didn't have to type that much. And yes, you got that part exactly right.

As for your question, I suppose for two sequences to be the same means that they both a)destroy the tanks and b)hit the same squares the same amount of times. Thus they'd have to follow the same W-B-W rules and could be morphed into each other so that a more computationally intensive process becomes static or vice-versa. Forgive me if this seems unimportant, I tend to think in terms of algorithms.