I've wrote a small program to calculate all (complex) roots of the polynomials of the 24th degree, with all coefficients being either +1 or -1, and to draw the result on the complex plane.
The (very impressive in my opinion) result can be seen here: 
(Warning! 20000x20000 image, over 40 Mb, I've heard that at least Firefox is incapable of showing it)
Technical details, if you care: the program is in OCaml (my general purpose language of choice), roots calculation is done with help of GSL, drawing is done with Cairo.
So I finally finished my little program, it's worth about 25% of one homework assignment (the prof keeps us busy). Written in Python. I was tasked with calculating the complex roots of the function z^n - 1 = 0, with z of the form z = a*x + i*b*y. Here, a and b are coefficients, both of which I varied from -1 to +1 divided into 1001 grid points. In other words, the coordinates of a point on the complex plane from x=-1 to +1 and y=-i to +i were inputted for z. We were only required to map to degree n=3 and n=4, but I can go as high as n=9 before my computer starts protesting. The current hang-up seems to be Python's unwillingness to use scientific notation for the imaginary component. Anything bigger than about 10^16*i gives me errors. The Newton Method was my prescribed workhorse for root solving. So I set up some reasonable constraints for "converging" on a solution, and plotted two things for each value of n: number of iterations it took to converge, and the final solution (converted to a color). As mentioned, both are a function of starting position on the complex plane. Each pair of 1001x1001 plots took about 20 minutes to generate. You must be a wizard to do n=24 for 20,000x20,000 on any timescale less than several days!
Absolutely fantastic. Never seen anything like this! You say degrees 22 and 23 are not so interesting compared to 24. Has the abundance of factors something to do with it? Would 60 degrees be oh so delicious? I know you don't have the power to do 60 degrees (yet), but I just want to plant to seed. Going to spend some time thinking through this. Thank you!
Nice! I have to write a "basin of attraction" program in Python within the next 8 days, using Newton's Method to solve for roots, and I'll post the output here with all relevant details.
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 :)
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?I cannot do degree 25 or more, as I'm running out of memory
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!
That's really something. I didn't know what to expect, but it wasn't that. Can you talk us through this, in Fisher-Price math? Let me see how far I can get on my own. A complex number, like 2 + 3i, has a real (2) and imaginary (3i) part. The real part is simple enough, and we can pretend the i in the imaginary part is just another variable like x, except that when two i's get multiplied together they change into a -1. You wrote a program to solve an equation, actually a lot of equations. Let's start with one equation. A simple polynomial might have the form ax² + bx + c = something, you didn't say what but I'll guess it's zero. You restricted the coefficients, a, b, and c to the values 1 and -1. Each of the three coefficients can take two values, so we should get 2 × 2 × 2 = 8 equations. 1x² + 1x - 1 = 0 1x² - 1x + 1 = 0 1x² - 1x - 1 = 0 -1x² + 1x + 1 = 0 -1x² + 1x - 1 = 0 -1x² - 1x + 1 = 0 -1x² - 1x - 1 = 0 To find solutions, we ask what values x can take that give true equations. Considering real numbers first, we might be lazy and plot all eight equations on a graph and see where they cross the x-axis. We find four values of x that work, each time yielding zero in two different polynomials. The highest is near x = 1.618, where the fourth and fifth polynomials both equal zero. In other words, so or = (1 ± √(1 - 4(1)(-1)))/2 = (1 ± √5)/2 which, in agreement with the graph, is about -0.618 and 1.618, the latter of which is phi, the golden ratio. I thought that number looked familiar. So we could plot these values on a number line, but it would be two boring dots. So let's consider complex numbers. Now x can have two parts, one a real number like before and the other imaginary. For the first equation, 1x² + 1x + 1 = 0, we need the imaginary part or it's hopeless, because the sum is always positive if x is real. If x is complex, the equation takes the form 1(a + b)² + (a + b) + 1 = 0, and when we solve it we will change any b² to a -1 and hopefully the b coefficient will be zero at the end because we're trying to get to zero and there's no i in zero. So a² + 2ab - 1 + a + b + 1 = 0 a² + 2ab + a + b = 0 a² + a + 2ab + b = 0 a² + a + b(2a + 1) = 0 Maybe now we just assume b is zero so that part drops out, and a(a + 1) = 0 and a = 0 or -1. But those don't satisfy the first equation, so maybe it was a bad assumption and this equation has no complex roots. How about the second one, 1x² + 1x - 1 = 0. It has two real roots, near x = -1.618 and x = 0.618. If x is complex, it takes the form a² + ab + b² + a + b - 1 = 0 a² + ab - 1 + a + b - 1 = 0 a² + ab + a + b - 2 = 0 a² + a + ab + b - 2 = 0 a² + a + b(a + 1) - 2 = 0 b(a + 1) = 2 - a² - a b = (2 - a² - a) / (a + 1) and cheating again with the graphing tool I find that when a = -2 or 1 the b equals zero. So that gives me x = (-2 + 0i) and x = (1 + 0i), but these are both real. Is my Fisher-Price arithmetic off, or did I get off track by assuming you would set the polynomials to zero, or perhaps do things only get interesting in the higher degree polynomials? 1x² + 1x + 1 = 0 1x² - 1x - 1 = -1x² + 1x + 1 2x² - 2x - 2 = 0
Now we turn to the quadratic formula to solve for x: x² - x - 1 = 0 x = (-b ± √(b² - 4ac)) / 2a a² + 2ab + b² + a + b + 1 = 0 a² + a = 0 1(a + b)² + 1(a + b) - 1 = 0
Side comment first: yes, root of the polynomial is something that makes it to be equal to zero, so there is 0 in the right side of every equation. Square equations You're completely right up to the point where you're trying to find complex roots of a square equation. It is done using exactly the same way as for real roots: So if your equation is 1x² + 1x + 1 = 0, then a=1, b=1, c=1 and roots are Since √(-3) = i√3 (because i√3 x i√3 = (i x i) x (√3 x √3) = -1 x 3 = -3), the roots are -1/2 ± (√3/2)i It is possible to find complex roots your way You can find the roots the way you calculate things, but be more careful with imaginary part; let's represent x not as 'a + b' as you do but as 'a + bi', where both a and b are real numbers, and i is, well, square root of -1. Then we will get: a² + 2abi + (bi)² + a + bi + 1 = 0 => a² + 2abi - b² + a + bi + 1 = 0 => (a² - b² + a + 1) + (2ab + b)i = 0 Now look at the number at the left side: it has real part (a² - b² + a + 1) and imaginary part (2ab + b)i. For this number to be zero, both parts should be zero (you cannot make zero out of real number by adding imaginary part to it). So we got that I'm going to look at the imaginary part first. 2ab + b = 0 => (2a + 1)b = 0 Now we have two numbers (2a + 1 and b) such that their product is zero, so one of those numbers is zero. It cannot be b (try to substitute b = 0 -- the equation turns out to a² + a + 1 = 0, which has no real solutions, and a was a real number). Hence Now, when we know a, let's look at the real part. 1/4 - b² - 1/2 + 1 = 0 => b² = 3/4 => b = ±√3/2 We substituted x = a + bi, so we can get our two roots now: -1/2 ± (√3/2)i General situation There is a theorem that the number of (complex) roots a polynomial has is equal to its degree. Well, the whole thing is a bit less straightforward as some roots can repeat themselves, like in x² - 2x + 1 = 0 -- there is only one solution, x = 1, but since the same equation can be represented as (x - 1)(x - 1) = 0, the solution x = 1 is, well, counted twice -- once for each "x - 1". I'm not going deeper here, but what's important to know: if you take random polynomial of degree n, almost always it will have exactly n different complex roots. It is not possible to calculate exact roots for general polynomial of degree over 4 -- not that we don't know how, it's proven theorem that the formula cannot exist. About my picture So, I have equations which look like ..x^24 + ..x^23 + ...and so on... + ..x^2 + ..x + .. = 0, where every '..' is either 1 or -1. I lose nothing assuming that x^24 always has 1 (think about it), so there are 24 coefficients left. It gives me 2^24 = 16777216 equations. I find approximate solutions for each of those, and each has exactly 24 solutions. In this way I'm getting 2^24 x 24 = 402653184 complex numbers. For each such number, if it is a + bi, I put a point with coordinates (a, b). And that's how I've got my picture! x = (-b ± √(b² - 4ac)) / 2a (-b ± √(b² - 4ac)) / 2a = (-1 ± √(1 - 4)) / 2 = (-1 ± √(-3)) / 2 1(a + bi)² + 1(a + bi) + 1 = 0 => a² - b² + a + 1 = 0 and (2ab + b)i = 0 (2ab + b)i = 0 => 2a + 1 = 0 => a = -1/2 a² - b² + a + 1 = 0 =>