Pi in the Mandelbrot set and computational limits. Part 1: Meet the set
Lately I’ve been learning more about an old mathematical friend: the Mandelbrot set. Amazingly, you can use it to compute the digits of 𝜋! But it doesn’t give up those digits easily. Along the way, it taught me a lesson about computational limits.
But first, have you met the Mandelbrot set?
I can’t remember exactly when I met the Mandelbrot set, but this surprising and intriguing object has been a touchstone throughout my programming life. The textbook for my first programming class was Oh! Pascal!, by Doug Cooper, a delightfully inspiring primer which guides readers from zero knowledge of programming to plotting the Mandelbrot set using the graphics library in Turbo Pascal 6.0. Wonderful.1
If you haven’t yet met the Mandelbrot set, allow me to introduce you. But I’m not just going to insert one of the internet’s myriad high-resolution color renderings of it here. That would be cheating. I hope to show you how to render it yourself, using the Swift programming language and its Standard Library, high school math, and a bunch of asterisks. If you already know how to plot the Mandelbrot set and want to get to the stuff about \(\pi\) and computational limits, skip to part 2 of this series.
Meeting the Mandelbrot set
The Mandelbrot set was discovered in the study of iterated functions—functions whose output is repeatedly fed back into themselves as input. Researchers are interested in the behavior of iterated functions over time. Do the output values grow ever larger—do they “blow up”, to use the jargon? Or do they approach a constant value, or repeat the same sequence endlessly? The jargon for these latter cases is “bounded”.
One example of a function that is iterated in this way is
\[f(x) = x^2 + c\]where \(x\) is initially set to \(0\), and \(c\) is some number that remains the same throughout the iterations. Let’s try some values for \(c\) and see what happens.
For \(c = 0\), we start with \(x = 0\), and compute
\[f(0) = 0^2 + 0 = 0\]Then we take the output, \(0\), plug that back into the function as \(x\), and compute again:
\[f(0) = 0^2 + 0 = 0\]We can see that this is going nowhere. We’ll always get \(0\), so the sequence of values for \(c = 0\) is just
\[0, 0, 0, 0, ...\]We’ve learned \(f\) is bounded for \(c = 0\).
Let’s try \(c = -1\).
\[\begin{align} f(0) &= 0^2 + -1 = -1 \\ f(-1) &= -1^2 + -1 = 1 - 1 = 0 \\ f(0) &= 0^2 + -1 = -1 \\ f(-1) &= -1^2 + -1 = 1 - 1 = 0 \end{align}\]That’s a bit more interesting. The output alternates between \(-1\) and \(0\) forever. The sequence of output values is
\[-1, 0, -1, 0, ...\]\(f\) is bounded for \(c = -1\), too.
How about \(c = 1\)?
\[\begin{align} f(0) &= 0^2 + 1 = 1 \\ f(1) &= 1^2 + 1 = 2 \\ f(2) &= 2^2 + 1 = 5 \\ f(5) &= 5^2 + 1 = 26 \\ f(26) &= 26^2 + 1 = \end{align}\]Okay, now I need a calculator. But didn’t I say we were going to write some Swift code? Let’s start by writing a program to compute this function. It’s a simple function and we should be able to see the result quickly in a Swift playground. Launch Xcode2, choose File -> New -> Playground…, and pick any of the templates. A blank iOS playground is fine.
By the way, as you’re adding code to the playground and running it, you may occasionally encounter an error message in the output area below that reads
error: Couldn't lookup symbols
If that happens, you’ll want to stop the playground. You can do that by clicking the little outline of a square that’s in the divider between the code editor and the output area. It’s labeled “Stop Playground” if you hover your pointer over it. Then you can resume running your playground.
We want to compute \(f(26)\). The function translates almost directly into Swift. We’ve been working with integers so far, so we’ll type x and c as Int.
import UIKit
func f(_ x: Int) -> Int {
x*x + 1
}
print(f(26))
// -> 677
Note that we use _ as x’s argument label so our calls to f() look more mathy: f(26) instead of f(x: 26).
We could already see what was happening with the sequence of output values we computed by hand, and \(f(26)\) emphasizes it: these values are getting big fast.
print(f(677))
// -> 458330
The sequence of output values is
\[1, 2, 5, 26, 677, 458330, ...\]When \(c = 1\), \(f\) isn’t bounded—it blows up.
Furthermore, we can guess \(f\) will blow up even quicker for any values of \(c\) larger than \(1\). Let’s test that.
But our initial Swift implementation of \(f\) is limiting. Even though mathematicians would say \(f\) is a function of \(x\) only, since \(c\) is a constant throughout the iterations, we want to be able to easily compute \(f\) for different values of \(c\). In Swift, you can do more with functions than just define and call them. You can pass functions as arguments to other functions, and return functions from functions. Anything you can do with, say, an Int, or a String, or a type you’ve defined, you can do with a function. We can write a function that takes a value for \(c\), and returns a function \(f\) for that \(c\).
func f(c: Int) -> (Int) -> Int {
{ x in
x*x + c
}
}
This new function that returns a function can still be named f, because Swift can tell that it’s a different type of
function. Our
original f(_:) takes an Int argument and returns an Int, so its type is
(Int) -> Int
Our new function f(c:)
also takes an Int argument, but returns another function of type (Int) -> Int, so the new function’s type is
(Int) -> ((Int) -> Int)
or, without the extra parentheses around the return type,
(Int) -> (Int) -> Int
We can assign the function returned by f(c:) to a variable, and call that variable like a function (because it
is a function!).
let fC2 = f(c: 2)
print(fC2(0))
// -> 2
print(fC2(2))
// -> 6
print(fC2(6))
// -> 38
print(fC2(38))
// -> 1446
print(fC2(1446))
// -> 2090918
print(fC2(2090918))
// -> 4371938082726
We can also pass a value for x directly to the result of f(c:) like so:
print(f(c: 2)(38))
// -> 1446
When \(c = 2\), the sequence of output values of \(f\) is
\[2, 6, 38, 1446, 2090918, 4371938082726, ...\]As we guessed, \(f\) blows up for \(c = 2\) even quicker than for \(c = 1\).
So far we’ve seen that \(f\) is bounded for \(c = 0\) and \(c = -1\), and blows up for \(c \geq 1\). What about when \(c\) is something like \(\frac{1}{2}\) or \(-\frac{1}{3}\) or \(\pi\)?
Now we find that f(c:) is still too limiting, because it can only work with integers. We need a new f(c:) that can work with the numbers between the integers: the real numbers. Swift’s recommended real number type (and the type it automatically infers when you use a real number literal like 3.14) is Double, which is a 64-bit floating-point type. But Swift also has a 32-bit floating-point type. Might we want to use that instead? Or what if we want to go beyond the real numbers?3 How many different versions of f(c:) will we need to write? Swift has a great solution for this problem: using a generic function, we can write one version of f(c:) that works with any kind of number.
One way to think about generic functions is that they allow us to use variables for types as well as values. When declaring a generic function, you specify the type variables you want to use in angle brackets (< >) after the function name, and then you use those type variables where you would normally use type names. For example, here is a generic version of f(c:) that uses a type variable named Number instead of the specific type Int.
func f<Number>(c: Number) -> (Number) -> Number {
{ x in
x*x + c
}
}
This lets us pass a value of any type as c. It could be 1, or 3.14159, or even "hello", but we’re about to deal with that.
This function doesn’t compile. We get an error at the multiplication operator * that says
Protocol 'FloatingPoint' requires that 'Number' conform to 'FloatingPoint'
Swift is inferring that x has the variable type Number, since we’re adding x*x to c, and we’ve said c has the variable type Number. We intend Number to be some kind of number, but we haven’t told Swift that. Right now, Number could be anything: a boolean value, a string, a dictionary. And yet we’re trying to multiply x by itself.
Swift won’t allow just any type of value to be multiplied. We can’t multiply two strings or two dictionaries. Swift is telling us that Number needs to be a type that can be multiplied. To accomplish that, it’s proposing we constrain Number to types that conform to the FloatingPoint protocol4. I don’t know why Swift picks FloatingPoint5, as there are other, more generic protocols in the standard library that provide the ability to multiply.
One example is the protocol Numeric, which specifies the properties and methods that allow a type’s values to be multiplied and added. We can tell Swift that Number should only stand in for types that support multiplication and addition by constraining it to types that conform to Numeric:
func f<Number: Numeric>(c: Number) -> (Number) -> Number {
{ x in
x*x + c
}
}
We can make that definition more concise—and save ourselves typing in the future—by defining a new name for the function
type (Number) -> Number. This is done with Swift’s typealias keyword.
typealias Func<Number> = (Number) -> Number
func f<Number: Numeric>(c: Number) -> Func<Number> {
{ x in
x*x + c
}
}
Now we have an f(c:) to which we can pass any type of value that implements the Numeric protocol; i.e., any
type that supports addition and multiplication. Let’s try it out with some familiar integer values of c to verify it’s
working.
print(f(c: 0)(0))
// -> 0
print(f(c: -1)(-1))
// -> 0
print(f(c: 2)(38))
// -> 1446
Looking good. Now we can try \(c = \frac{1}{2}\).
let fC0_5 = f(c: 0.5)
print(fC0_5(0))
// -> 0.5
print(fC0_5(0.5))
// -> 0.75
print(fC0_5(0.75))
// -> 1.0625
print(fC0_5(1.0625))
// -> 1.62890625
print(fC0_5(1.62890625))
// -> 3.1533355712890625
print(fC0_5(3.1533355712890625))
// -> 10.443525225156918
It happens more slowly than for \(c = 1\), but these values are blowing up. Let’s try \(c = \frac{1}{4}\).
let fC0_25 = f(c: 0.25)
print(fC0_25(0))
// -> 0.25
print(fC0_25(0.25))
// -> 0.3125
print(fC0_25(0.3125))
// -> 0.34765625
print(fC0_25(0.34765625))
// -> 0.3708648681640625
print(fC0_25(0.3708648681640625))
// -> 0.38754075043834746
print(fC0_25(0.38754075043834746))
// -> 0.4001878332503175
Interesting. The sequence of output values is increasing, but very slowly. We need to try more iterations, but this is getting rather tedious. We’ve got a computer—let’s program it compute as many terms of the sequence of output values as we like. The Swift Standard Library defines a Sequence protocol will let us do exactly that, plus much more.
To conform to Sequence, a type need only6 implement the method
mutating func next() -> Element?
next() returns the next value in the sequence. Element is an alias for whatever type of values the sequence contains. This means we can only access the iterated function values in one direction (from the beginning), but that’s fine—that’s the way we compute them. With this one method, our custom sequence type gains access to a huge amount of functionality for free: iteration over the sequence using for-in, finding methods like first(), min(), and max(), and higher-order functions like map() and reduce().
Let’s build it. We’ll define a struct named IteratedFuncValues that is generic with respect to the type of the values, just like f(c:). It will store the function we want to iterate, and how many times we want to iterate.
Why store the desired number of iterations? A function can be iterated infinitely many times—our sequence of iterated function values is infinitely long, at least mathematically. But our code is not running in the limitless space of mathematical thought. It’s running on a computer with a finite amount of memory devoted to each value. I said I had learned something about computational limits, and there’s much more to say about that in part 2, but here’s an appetizer. The largest value an Int can represent is
print(Int.max)
// -> 9223372036854775807
For a Double:
print(Double.greatestFiniteMagnitude)
// -> 1.7976931348623157e+308
Int and Double both use 64 bits to store each value. As we’ll see shortly, it’s astonishing how quickly either of these maximum values can be exceeded by iterating \(f\) with even a relatively small value of \(c\).
Unless we know in advance how many iterations each different \(f\) can bear for each different value type, storing the iteration count in our sequence type isn’t going to prevent all possibility of numeric overflow. But some of the Sequence methods and properties we’ll want to use iterate over the entire sequence, so they’ll definitely crash due to overflow or memory exhaustion if we haven’t set a finite limit to the sequence length. Storing the iteration count in a property is the the most straightforward way I know to do that.7 As with everything on this blog, suggestions for improvements are welcome!
We also need to store the result of the previous iteration. We can initialize this property to the first argument passed to the function, which for our purposes in this program is always \(0\). If we constrain our variable-type number to Numeric again, we can use its handy .zero property.
struct IteratedFuncValues<Number: Numeric> {
let f: Func<Number>
var iterations: Int
var currentValue = Number.zero
}
Now we just need to define next(). We’ll do it in an extension that adopts the Sequence protocol. We could have defined next() in the struct above, but doing it this way is a convention that makes the code more digestible and clarifies what roles different methods play in a custom type.
extension IteratedFuncValues: Sequence, IteratorProtocol {
mutating func next() -> Number? {
if iterations > 0 {
currentValue = f(currentValue)
iterations -= 1
return currentValue
} else {
return nil
}
}
}
Each time next() is called, it passes currentValue to f and updates it with the result.
Note that IteratedFuncValues also declares conformance to IteratorProtocol, which is really the protocol that requires next(). The compiler is doing magic to synthesize the requirements of Sequence from this one method belonging to another protocol. Marvelous.
Note also that next() returns an optional8 Number, indicated by the trailing ?. We use this to limit the number of iterations. After iterations calls to next(), it returns nil.
Let’s test our freshly-baked type with some familiar values of \(c\).
for value in IteratedFuncValues(f: f(c: 0), iterations: 3) {
print(value)
}
// -> 0
// -> 0
// -> 0
IteratedFuncValues(f: f(c: -1), iterations: 3).forEach { value in
print(value)
}
// -> -1
// -> 0
// -> -1
let fC2Values = IteratedFuncValues(f: f(c: 2), iterations: 6)
let valueDescriptions = fC2Values.map(\.description)
print(valueDescriptions.joined(separator: "\n"))
// -> 2
// -> 6
// -> 38
// -> 1446
// -> 2090918
// -> 4371938082726
It’s working! And since IteratedFuncValues is a full-fledged Swift Sequence, we can use it in a variety of different styles: built-in for-in looping, forEach(), and map() are all producing the same results.
The pleasant readability of
fC2Values.map(\.description)
is courtesy of the ability of key path expressions, e.g., \.somePropertyOfAType, to be used as functions. That feature was new in Swift 5.2, which first shipped in Xcode 11.4. If you’re using an earlier version of Swift, you’ll have to say something like
fC2Values.map { $0.description }
Check this out:
print(fC2Values.max()!)
// -> 4371938082726
print(Int.max)
// -> 9223372036854775807
After only six iterations with \(c = 2\), we don’t have enough headroom left in an Int to square it again. But even Double can’t give us many more iterations.
print(IteratedFuncValues(f: f(c: 2.0), iterations: 10).max()!)
// -> 1.781492681120714e+202
print(Double.greatestFiniteMagnitude)
// -> 1.7976931348623157e+308
These are unfathomably huge numbers, vastly greater than the number of atoms in the universe!9
Now let’s see if \(f\) is really bounded when \(c = \frac{1}{4}\).
var fC0_25Values = IteratedFuncValues(f: f(c: 0.25), iterations: 10)
for value in fC0_25Values {
print(value)
}
// -> 0.25
// -> 0.3125
// -> 0.34765625
// -> 0.3708648681640625
// -> 0.38754075043834746
// -> 0.4001878332503175
// -> 0.4101503018815839
// -> 0.4182232701335544
// -> 0.42491070368120404
// -> 0.430549106102856
Hmmm. It’s still not clear. Let’s try a lot more iterations and check the last several of them.
fC0_25Values.iterations = 1000
let lastFive = fC0_25Values.suffix(5)
print(lastFive.map(\.description).joined(separator: "\n"))
// -> 0.4990046581204365
// -> 0.4990056488258937
// -> 0.4990066375601512
// -> 0.4990076243290881
// -> 0.49900860913856027
It looks like the sequence is tending towards \(0.5\), but it’s still slowly increasing, even after a thousand iterations. It seems there are some values of \(c\) for which it can take a lot of iterations to determine the boundedness of \(f\). We’ll explore this further, with surprising results, in part 2.
So far we know (or strongly suspect) \(f\) is bounded for \(c = {-1, 0, \frac{1}{2}}\), and blows up for the cases of \(c > \frac{1}{2}\) that we’ve tested. Let’s test a few \(c < -1\).
IteratedFuncValues(f: f(c: -2), iterations: 4).forEach { print($0) }
// -> -2
// -> 2
// -> 2
// -> 2
That seems weird at first glance, but
\[\begin{align} f(0) &= 0^2 - 2 = -2 \\ f(-2) &= -2^2 - 2 = 4 - 2 = 2 \\ f(2) &= 2^2 - 2 = 4 - 2 = 2 \end{align}\]so it checks out. Let’s try more negative integers.
IteratedFuncValues(f: f(c: -3), iterations: 5).forEach { print($0) }
// -> -3
// -> 6
// -> 33
// -> 1086
// -> 1179393
IteratedFuncValues(f: f(c: -4), iterations: 5).forEach { print($0) }
// -> -4
// -> 12
// -> 140
// -> 19596
// -> 384003212
Each starts out negative, but in the second iteration the squaring of a negative results in a positive of greater magnitude than \(c\), and they blow up quickly from there. It looks like values of \(c \lt -2\) cause \(f\) to blow up increasingly rapidly.
We’ve only checked a few points in the infinite realm of numbers. Let’s get serious. Let’s use the computer to visualize how \(f\) behaves with values of \(c\) spanning some interesting segment of the number line—say, \(-3.0 \leq c \leq 1.0\). We could generate an image that shows this on a number line, but let’s not get bogged down in graphics code. We can crudely but adequately represent \(f\)’s behavior with a string of spaces and asterisks. If \(f\) is bounded for a value of \(c\) in the segment, we can write an asterisk to the string at the character position corresponding with \(c\). If instead \(f\) blows up for that value of \(c\), we can write a space to that position in the string. Thus a range of bounded \(c\) values will look like a string of asterisks such as *****, and a range of unbounded values will be a string of spaces.
But wait: how do we determine whether \(f\) is going to blow up for a given \(c\)? So far, we’ve been iterating anywhere from two to a thousand times and eyeballing the results. Do we teach the computer to eyeball the results? Let’s avoid artificial intelligence problems for now, and see if we can find a simple heuristic to use. Let’s examine the values after many thousands of iterations for several values of \(c\) in the range where we think \(f\) is bounded.
for c in stride(from: -2.0, through: 0.25, by: 0.25) {
let manyValues = IteratedFuncValues(f: f(c: c), iterations: 30000)
let lastTen = manyValues.suffix(10)
print("c = \(c): ..., \(lastTen.map(\.description).joined(separator: ", "))")
}
// -> c = -2.0: ..., 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0
// -> c = -1.75: ..., -0.054946929270315525, -1.746980834963763, 1.3019420377306865, -0.05494693038966769, -1.746980834840753, 1.3019420373008943, -0.054946931508796704, -1.7469808347177676, 1.3019420368711878, -0.05494693262770256
// -> c = -1.5: ..., -1.4897386587095798, 0.719321271253818, -0.9825769087217913, -0.5345426184467287, -1.214264189064115, -0.025562479156467033, -1.4993465596593751, 0.7480401059624042, -0.9404359998717551, -0.6155801301452123
// -> c = -1.25: ..., -1.206116011493828, 0.20471583318178, -1.2080914276446897, 0.20948489754858457, -1.206116077699059, 0.20471599288416242, -1.2080913622574516, 0.20948473956106528, -1.2061161438910326, 0.20471615255457398
// -> c = -1.0: ..., -1.0, 0.0, -1.0, 0.0, -1.0, 0.0, -1.0, 0.0, -1.0, 0.0
// -> c = -0.75: ..., -0.504073637556147, -0.4959097679209142, -0.5040735020806251, -0.4959099045001741, -0.5040733666186282, -0.49591004106576203, -0.5040732311701542, -0.49591017761768025, -0.5040730957352009, -0.495910314155931
// -> c = -0.5: ..., -0.3660254037844387, -0.3660254037844386, -0.3660254037844387, -0.3660254037844386, -0.3660254037844387, -0.3660254037844386, -0.3660254037844387, -0.3660254037844386, -0.3660254037844387, -0.3660254037844386
// -> c = -0.25: ..., -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752, -0.20710678118654752
// -> c = 0.0: ..., 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
// -> c = 0.25: ..., 0.4999666700852357, 0.49996667119611893, 0.4999666723069281, 0.4999666734176632, 0.49996667452832433, 0.4999666756389114, 0.4999666767494244, 0.49996667785986343, 0.49996667897022845, 0.4999666800805195
You may have noticed that it took a while to compute each of those 30,000-element sequences. That’s not because it’s hard for a modern computer to perform that many multiplications and additions; rather, there’s a lot of behind-the-scenes instrumentation happening in a Swift playground to display all the debug output next to the code. This same code runs much faster in an Xcode command-line app project. But in part 2, we’ll make computations that take hours, even in a native command-line app. We should comment out that last block of code so we don’t have to keep waiting for it.
//for c in stride(from: -2.0, through: 0.25, by: 0.25) {
// let manyValues = IteratedFuncValues(f: f(c: c), iterations: 30000)
// let lastTen = manyValues.suffix(10)
// print("c = \(c): ..., \(lastTen.map(\.description).joined(separator: ", "))")
//}
Back to the task at hand. We already knew that for \(c = {-2, -1, 0}\) we have repeating sequences, and we see that confirmed here. And we have further confirmation that for \(c = \frac{1}{4}\) the sequence slowly approaches \(\frac{1}{2}\).
Of the sequences we hadn’t yet seen, \(c = -\frac{1}{4}\) is quite stable as far as the precision of a Double is able to show. \(c = -\frac{1}{2}\) is similar, though if we look closely, it’s alternately increasing and decreasing by \(1\) in its sixteenth digit of precision.
All the other sequences are alternating between at least two completely different numbers, and for \(c = -1.5\), there’s no discernible pattern in the last ten values of the sequence. However, in no case does a sequence contain a value with an absolute value or magnitude greater than \(2\).
This is hardly rigorous, but from our investigation so far, we might reasonably use the following heuristic: if the magnitude of \(f\) is less than or equal to \(2\) after some number of iterations for a value of \(c\), then it’s not going to blow up.
Thirty thousand iterations took too long, though. Let’s check whether a more reasonable number like thirty still supports our heuristic.
for c in stride(from: -1.75, through: -0.25, by: 0.25) {
let biggest = IteratedFuncValues(f: f(c: c), iterations: 30).map(\.magnitude).max()!
print("c = \(c): \(biggest)")
}
// -> c = -1.75: 1.75
// -> c = -1.5: 1.5
// -> c = -1.25: 1.25
// -> c = -1.0: 1.0
// -> c = -0.75: 0.75
// -> c = -0.5: 0.5
// -> c = -0.25: 0.25
None of those sequences contains a value with a magnitude greater than \(2\) in the first thirty iterations.
Using our heuristic and more goodies from the Standard Library, we can write some simple computed properties to test whether \(f\) is bounded for any value of \(c\).
First, let’s give ourselves some flexibility in defining what we mean by “magnitude”. Numeric’s magnitude property is exactly what we want for all the types we’ve worked with so far, but it may not be defined the way we want for other types we might want to use. The concept of the magnitude or absolute value of some number can also be thought of as the length measured from zero to that number. We’ll define a new protocol Lengthy which requires a property length, which has the same type as Numeric’s magnitude.
protocol Lengthy: Numeric {
var length: Magnitude { get }
}
We can provide a default implementation of length that will work with all the Numeric types we’ve been using so far, since they already have a magnitude implementation that meets our needs. Though we still have to tell Swift explictly that Int and Double conform to Lengthy.
extension Lengthy {
var length: Magnitude {
magnitude
}
}
extension Int: Lengthy {}
extension Double: Lengthy {}
Now we have the ability to redefine length for other types, if necessary.10
We’ll also define a convenient property for testing whether a value’s length is \(2\) or less.
extension Lengthy {
var lengthTwoOrLess: Bool {
length <= 2
}
}
Finally, we’ll define a property that lets us concisely ask whether a particular iterated function value sequence is bounded.
extension IteratedFuncValues where Number: Lengthy {
var bounded: Bool {
allSatisfy(\.lengthTwoOrLess)
}
}
print(IteratedFuncValues(f: f(c: 2), iterations: 30).bounded)
// -> false
print(IteratedFuncValues(f: f(c: 0.25), iterations: 30).bounded)
// -> true
Thank you, Swift Standard Library! Not only is allSatisfy(_:) a concise way to implement our heuristic, it’s efficient, too: as soon as its predicate fails, it short-circuits and returns false. We can verify that by checking the number of times lengthTwoOrLess was invoked for the two examples we tried. The playground counts 32 invocations: 30 for \(c = 0.25\), because all of those values were less than \(2\), but only 2 invocations for \(c = 2\), since the second value of that sequence is \(6\).
Now let’s use this to plot the boundedness of \(f\) for \(c\) in \(-3.0 \leq c \leq 1.0\). If we could generate a sequence of values in that range, we could use map() to transform bounded cases into a * and blow-ups into a space. We could then build a string from that Character array and print it to see our visualization.
We’ve already seen that stride() can create a sequence of values similar to what we want. But stride requires us to specify the step between values in the desired range. Since we want to draw a graphic (low-resolution though it may be), it would be nicer to specify how many points we’d like to plot in the desired range, and have the step size computed for us. This is a good opportunity to write another custom type, and this time we’ll adopt Collection, because we’ll find use for its subscript later.
struct PlotPoints<Number> where Number: FloatingPoint {
typealias Index = Int
let start: Number
let end: Number
let pointCount: Index
let step: Number
init(from start: Number, to end: Number, pointCount: Index) {
self.start = start
self.end = end
self.pointCount = pointCount
self.step = (end - start) / Number(pointCount - 1)
}
}
Compared to IteratedFuncValues, PlotPoints has a few more properties, and we’ve added an init() so we can compute
step. To allow that computation, we constrain Number to FloatingPoint instead of Numeric, because we need to be able to perform division on Number values.11
Now we’ll adopt Collection.
extension PlotPoints: Collection {
typealias Element = Number
typealias Indices = Range<Int>
var startIndex: Index { 0 }
var endIndex: Index { pointCount }
var indices: Indices { 0..<pointCount }
subscript(position: Index) -> Element {
start + step * Number(position)
}
func index(after i: Index) -> Index {
i + 1
}
}
Collection conformance requires more work, as promised, but it’s all quite straightforward. We define computed properties for the starting and ending indices, and for a range of all the indices. The subscript does the real work here, computing the value for the point at a requested position in the collection.
Let’s take it for a spin.
let points = PlotPoints(from: -3.0, to: 1.0, pointCount: 100)
print(points.count)
// -> 100
print(points.map(\.description).joined(separator: "\n"))
// -> -3.0
// -> -2.95959595959596
// -> -2.919191919191919
// -> [...]
// -> 0.9191919191919196
// -> 0.9595959595959598
// -> 1.0
Looking good. We’re one convenience function away from our first plot. Let’s factor out our plotting rule into its own function, so we don’t have to duplicate its logic for every plot we want to make.
func plotValue(for predicate: Bool) -> Character {
predicate ? Character("*") : Character(" ")
}
Now we can plot!
let plots = points.map { point in
plotValue(for: IteratedFuncValues(f: f(c: point), iterations: 30).bounded)
}
print(plots.count)
print("first plot!")
print(String(plots))
100
first plot!
********************************************************
That’s somewhat unsatisfying. It’s roughly what we were expecting. It’s 100 characters long. We know \(f\) blows up at the points we’ve tried from \(c = -3.0\) to \(-2.0\), and then is stable at points from \(-2.0\) through probably \(0.25\). We know for sure it blows up again at points from \(0.5\) and beyond. So in the complete range from \(-3.0\) through \(1.0\), we expect to see a string of spaces, followed by a string of asterisks, followed by more spaces. And we do see that. But it would be nice to see some indices on the number line to confirm this. Here’s a function for that:
func renderIndices<Number>(_ indices: [Number], for points: PlotPoints<Number>) -> String where Number: CVarArg, Number: Comparable {
var renderedIndices = ""
var indexIterator = indices.sorted().makeIterator()
var currentIndex = indexIterator.next()
while renderedIndices.count < points.count {
guard let index = currentIndex else { break }
let plotPoint = points[renderedIndices.count]
if plotPoint >= index {
renderedIndices.append(String(format: "|%.2f", plotPoint))
currentIndex = indexIterator.next()
} else {
renderedIndices.append(" ")
}
}
return renderedIndices
}
renderIndices() takes an Array of the indices you’d like rendered, and what points to render them for, and returns a string that will line up nicely under the plot. It considers the provided indices to be suggestions; what it actually renders depends on the the values of the points nearest the suggestions, and how much space is available. It uses PlotPoints’s subscript to find the next point there’s room to write an index for.
This time we need Number to conform to two protocols. It needs Comparable because we’re sorting the provided indices, so the caller doesn’t have to. Very hospitable! And it needs CVarArg because we’re formatting the indices to show only 2 digits after the decimal point. Space is precious in the low-resolution environment of the console!
We’ll choose some interesting (and not necessarily sorted) indices, and print them using renderIndices().
let indices = [-3.0, -2.0, -1.0, 0.0, 0.25, 1.0, -0.75]
print(renderIndices(indices, for: points))
Now our plot makes sense.
********************************************************
|-3.00 |-1.99 |-0.98|-0.74 |0.03 |0.27 |1.00
This matches our expectations, at least to the resolution possible in a 100-point plot. It looks like \(f\) is bounded for all \(c\) in \(-2.0 \leq c \leq 0.25\), and unbounded everywhere else we tested.
I thought we were going to meet the Mandelbrot set?
Wasn’t I talking about meeting the Mandelbrot set a long time ago? Well, it’s so close now. Because the Mandelbrot set is simply the extension of the 1-dimensional plot we just drew to the 2-dimensional complex plane. Somewhat more formally (though I am not a mathematician), the Mandelbrot set is the set of all points \(c\) in the complex plane for which the absolute value of the function
\[f(z) = z^2 + c\]is bounded under iteration. Not coincidentally, this is the same function we’ve been studying. \(z\) is the conventional variable name for complex numbers, like \(x\) is for real numbers.
If you’d like a refresher on complex numbers and the complex plane, I highly recommend 3Blue1Brown’s Complex Number fundamentals lesson on YouTube.12 I’ll just quickly mention that complex numbers are composed of two different kinds of numbers: a real number, and an imaginary number that is some multiple of \(i \equiv \sqrt{-1}\). Even though \(i\) seems to violate the rules of multiplication, it’s extremely useful in math and physics, so it exists. Because a complex number has two distinct parts that can’t be commingled, it can’t be plotted on the number line. It has to be plotted in two dimensions, on a plane. By convention, we plot the real part along the horizontal axis, and the imaginary part along the vertical axis. Happily, you can draw much more interesting pictures on a plane than on a line—for example, the Mandelbrot set.
Let’s plot it! We have almost all the code we need. We can already plot the boundedness of \(f\) for any points \(c\) in the complex plane whose imaginary part is some constant value, because these are just horizontal lines in the plane. For example, we can use PlotPoints to plot the line for all points \(x + 0.5i\), where is \(x\) is in the range specified by the start and end properties.
We just need a new type that collects PlotPoints for as many horizontal lines as we want to plot in the range of imaginary values of \(c\) we want to cover. For example, we could plot 100 horizontal lines, each corresponding to a multiple of \(i\) from \(2.0i\) through \(-2.0i\). Each horizontal line could have 100 points corresponding to values of \(c\) from \(-3.0\) through \(1.0\). Assuming all of the complex points \(c\) for which \(f\) is bounded are found within the rectangle in the complex plane whose upper left corner is at \(-3.0 + 2.0i\) and whose lower right corner is at \(1.0 - 2.0i\), we’d have a 10,000 point plot containing the Mandelbrot set.
First, though, we need to be able to create complex values so we can use them with the generic functions and types we’ve written. A complex number type is provided in the Swift Numerics package, which can be downloaded from GitHub, but it’s not currently possible to use Swift packages in a playground.13 We’ll have to move our code to an Xcode command-line project.
In Xcode, hit cmd-shift-n, and select the macOS Command Line Tool template. I named my project mandelplot, of course. Once the project is created, click on the File menu and choose Swift Packages -> Add Package Dependency…. Type “numerics” in the search box, hit Return, and scroll until you find the “swift-numerics” package whose owner is “apple”. Select it and click Next. The displayed repository should be “git@github.com:apple/swift-numerics.git”. Make sure “Version” is selected. As I write this, the most recent version is 0.0.6, so enter “0.0.6” in the text field to the right of the pop-up menu. “Up to Next Major” is a fine selection for our purposes. Click Next. Check the box next to “ComplexModule” and click Next.
Xcode will download the swift-numerics package. If you have the Project Navigator open (cmd-1), you’ll see a new Swift Package Dependencies section appear, which should contain an entry for the swift-numerics package.
Now we can use swift-numerics’ Complex type. It’s generic over the type of its real and imaginary parts. We’ve been using Double for our work on the real number line. Let’s keep using it in each dimension of the complex plane. To avoid having to type ComplexModule.Complex<Double> a lot14, we’ll alias that to something more succinct. Open main.swift in your newly-created project, and make it look like this:
import ComplexModule
import Foundation
import RealModule
typealias Complex = ComplexModule.Complex<Double>
Now we can compute with complex values. Let’s explore the behavior of \(f\) in the new dimension that’s available to it. We know \(f\)’s behavior on the real number line: it blows up for \(c < -2.0\) and \(c > 0.25\), and it’s bounded between those points. And now we know that’s exactly the behavior of \(f\) for complex values of \(c\) whose imaginary part is \(0\). But what happens to \(f\) in the vast space above and below the real number line?
Addition and multiplication of complex numbers is more…er…complex than for real numbers. Some might also say more tedious. But it’s good to start with some computation by hand so we get a feel for what’s going on, before we let the computer do all the work. The simplest value we can try is \(c = 0 + i\), or just \(i\). Recall that when we iterated for real values of \(c\), we initially set \(z\) to \(0\). For the first iteration of complex \(f\), we set \(z = 0 + 0i\). We add complex numbers by adding the real and imaginary parts separately, and we multiply complex numbers using the FOIL (First, Outer, Inner, Last) algorithm: summing the products of the first, outer, inner, and last components of the multiplicands.
\[\begin{align} m(0 + 0i) &= (0 + 0i)^2 + i = 0 + i \\ m(0 + i) &= (0 + i)^2 + i = (0 + i)(0 + i) + i = 0 + 0i + 0i + i^2 + i = i^2 + i = -1 + i \\ m(-1 + i) &= (-1 + i)^2 + i = (-1 + i)(-1 + i) + i = 1 - i - i + i^2 + i = 1 - 2i -1 + i = 0 - i \\ m(0 - i) &= (0 - i)^2 + i = (0 - i)(0 - i) + i = 0 - 1 + i = -1 + i \\ m(-1 + i) &= 0 - i \end{align}\]So for \(c = i\), \(f\) visits a few different points on the complex plane, then settles into a repeating sequence of \(-1 + i, 0 - i\). Which means \(f\) is bounded for \(c = i\), so \(i\) is the first complex member of the Mandelbrot set we’ve seen!
Whew, that manual computation was a pain. Let’s have IteratedFuncValues do the work from now on. Copy all the code from your playground to main.swift, below the line starting with typealias.
typealias Complex = ComplexModule.Complex<Double>
// Paste all code from playground here.
import UIKit will cause a compile error if it was copied over. Delete that line. Hit cmd-r to build and run, and you’ll see all the output from your playground below, in the debug console.
To be clear, this is not a nice command line program. In a future post I’ll show how to make one of those using Swift Argument Parser. For now, we’re just writing code into an Xcode command line project, instead of a playground, so we have access to the swift-numerics package.
Let’s compute the sequences for some interesting complex \(c\) above and below the real number line. Add the following below the code you pasted in from the playground.
let interestingCValues = [
Complex(0, 2),
Complex(0, 1.1),
Complex(0, 1),
Complex(0, 0.5),
Complex(0, -0.5),
Complex(0, -1),
Complex(0, -1.1),
Complex(0, -2)
]
for c in interestingCValues {
print("c = \(c): \(IteratedFuncValues(f: f(c: c), iterations: 10).map(\.description).joined(separator: ", "))")
}
// -> c = (0.0, 2.0): (0.0, 2.0), (-4.0, 2.0), (12.0, -14.0), (-52.0, -334.0), (-108852.0, 34738.0), (10642029260.0, -7562601550.0), (5.605984456663374e+19, -1.609628539536427e+20), (-2.2766334180066574e+40, -1.8047105147285954e+40), (1.926079678012724e+80, 8.217328535318233e+80), (-6.38147053313017e+161, 3.165445899886102e+161)
// -> c = (0.0, 1.1): (0.0, 1.1), (-1.2100000000000002, 1.1), (0.2541000000000002, -1.5620000000000007), (-2.375277190000002, 0.306191599999999), (5.548188433423746, -0.3545798464992045), (30.656668025233646, -2.8345516061441076), (931.796611601489, -172.69581517990548), (838421.080811364, -321833.65084478585), (599373009932.8123, -539664234764.405), (6.801051875206817e+22, -6.469203534876586e+23)
// -> c = (0.0, 1.0): (0.0, 1.0), (-1.0, 1.0), (0.0, -1.0), (-1.0, 1.0), (0.0, -1.0), (-1.0, 1.0), (0.0, -1.0), (-1.0, 1.0), (0.0, -1.0), (-1.0, 1.0)
// -> c = (0.0, 0.5): (0.0, 0.5), (-0.25, 0.5), (-0.1875, 0.25), (-0.02734375, 0.40625), (-0.1642913818359375, 0.477783203125), (-0.2012851310428232, 0.3430086746811867), (-0.07713924692761764, 0.3619149079359445), (-0.12503193716972322, 0.44416431309988635), (-0.18164895171908027, 0.3889305510229235), (-0.11827063185835415, 0.3587023462303234)
// -> c = (0.0, -0.5): (0.0, -0.5), (-0.25, -0.5), (-0.1875, -0.25), (-0.02734375, -0.40625), (-0.1642913818359375, -0.477783203125), (-0.2012851310428232, -0.3430086746811867), (-0.07713924692761764, -0.3619149079359445), (-0.12503193716972322, -0.44416431309988635), (-0.18164895171908027, -0.3889305510229235), (-0.11827063185835415, -0.3587023462303234)
// -> c = (0.0, -1.0): (0.0, -1.0), (-1.0, -1.0), (0.0, 1.0), (-1.0, -1.0), (0.0, 1.0), (-1.0, -1.0), (0.0, 1.0), (-1.0, -1.0), (0.0, 1.0), (-1.0, -1.0)
// -> c = (0.0, -1.1): (0.0, -1.1), (-1.2100000000000002, -1.1), (0.2541000000000002, 1.5620000000000007), (-2.375277190000002, -0.306191599999999), (5.548188433423746, 0.3545798464992045), (30.656668025233646, 2.8345516061441076), (931.796611601489, 172.69581517990548), (838421.080811364, 321833.65084478585), (599373009932.8123, 539664234764.405), (6.801051875206817e+22, 6.469203534876586e+23)
// -> c = (0.0, -2.0): (0.0, -2.0), (-4.0, -2.0), (12.0, 14.0), (-52.0, 334.0), (-108852.0, -34738.0), (10642029260.0, 7562601550.0), (5.605984456663374e+19, 1.609628539536427e+20), (-2.2766334180066574e+40, 1.8047105147285954e+40), (1.926079678012724e+80, -8.217328535318233e+80), (-6.38147053313017e+161, -3.165445899886102e+161)
\(\pm2i\) are reaching truly dizzying magnitudes, and \(\pm1.1i\) are headed there too. \(-i\) joins \(i\) as a complex value of \(c\) for which iterations of \(f\) settle into a repeating sequence, so it is also a Mandelbrot set member. \(\pm0.5i\) look to be oscillating around some attractor. Perhaps they are also members of the Mandebrot set.
We’ve gotten a sense of the bounds of the Mandelbrot set in the complex plane: we already knew it ranged from \(-2\) through \(0.25\) along the real axis, and now it looks like it doesn’t extend much beyond \(\pm{i}\) along the complex axis. It’s not symmetrical left-to-right, but it might be symmetrical top-to-bottom. Let’s put our hunches to the test and plot it.
Here’s a new collection type whose elements are as many PlotPoints as we specify to sample the rectangle in the complex plane that we want to investigate.
struct ComplexPlotPoints<RealType: Real> {
typealias SomeComplex = ComplexModule.Complex<RealType>
typealias Index = Int
let upperLeft: SomeComplex
let lowerRight: SomeComplex
let horizontalPointCount: Index
let verticalPointCount: Index
let verticalStep: RealType
init(from upperLeft: SomeComplex, to lowerRight: SomeComplex, horizontalPointCount: Index, verticalPointCount: Index) {
self.upperLeft = upperLeft
self.lowerRight = lowerRight
self.horizontalPointCount = horizontalPointCount
self.verticalPointCount = verticalPointCount
verticalStep = (upperLeft.imaginary - lowerRight.imaginary) / RealType((verticalPointCount - 1))
}
}
It looks very similar to the base definition of PlotPoints, with an extra property to store the desired number of points in the imaginary dimension. We’re computing the step size in init() as before, but this time it’s the vertical step.
The number of horizontal and vertical points essentially define a sampling grid. We can sample the rectangle we want to plot with as few or as many points as we like, depending on how much resolution we want or the size of our output device, which in this case is the console. If our rectangle is a square, and we use a 2 x 2 grid of points, our plot of the Mandelbrot set won’t be very interesting. It might look like
**
**
or
*
*
or
*
*
But if we use too many points, it might not fit in the console.
The code to adopt the Collection protocol is also quite familiar. The only significant difference is in subscript(), where we return a collection of plot points, instead of a simple value.
extension ComplexPlotPoints: Collection {
typealias Element = PlotPoints<SomeComplex>
typealias Indices = Range<Int>
var startIndex: Index { 0 }
var endIndex: Index { verticalPointCount }
var indices: Indices { 0..<verticalPointCount }
subscript(position: Index) -> Element {
let imaginaryPart = upperLeft.imaginary - verticalStep * RealType(position)
return PlotPoints(from: SomeComplex(upperLeft.real, imaginaryPart), to: SomeComplex(lowerRight.real, imaginaryPart), pointCount: horizontalPointCount)
}
func index(after i: Index) -> Index {
i + 1
}
}
Unfortunately, this code causes a few compile errors, the key one being
Type 'ComplexPlotPoints<RealType>.SomeComplex' (aka 'Complex<RealType>') does not conform to protocol 'FloatingPoint'
Sorry about that. Up to this point, we haven’t had to edit any previously-written code as we added more plotting capabilities, but we’ll have to now.
The error occurs because we required Number in PlotPoints to conform to FloatingPoint, so we could perform division to compute the step size. We still need to perform division to compute the vertical step size in ComplexPlotPoints, but with complex values, which don’t conform to FloatingPoint. We need a protocol that allows for division of both real and complex numbers. Thankfully, the swift-numerics package provides one: AlgebraicField.
Change the definition of PlotPoints as commented below.
struct PlotPoints<Number> where Number: AlgebraicField { // Change constraint to `AlgebraicField`
typealias Index = Int
let start: Number
let end: Number
let pointCount: Index
let step: Number
init(from start: Number, to end: Number, pointCount: Index) {
self.start = start
self.end = end
self.pointCount = pointCount
self.step = (end - start) / Number(exactly: pointCount - 1)! // Change initializer.
}
}
extension PlotPoints: Collection {
typealias Element = Number
typealias Indices = Range<Int>
var startIndex: Index { 0 }
var endIndex: Index { pointCount }
var indices: Indices { 0..<pointCount }
subscript(position: Index) -> Element {
start + step * Number(exactly: position)! // Change initializer.
}
func index(after i: Index) -> Index {
i + 1
}
}
Complex conforms to Numeric, so it already has a magnitude property that will work for checking the boundedness of a sequence of Complex values. However, there are multiple ways to determine the magnitude or absolute value of a complex number. Complex’s magnitude property returns the absolute value of the real part, or the absolute value of the imaginary part, whichever is larger.
print(Complex(0, -0.25).magnitude)
// -> 0.25
print(Complex(3.14, -3.0).magnitude)
// -> 3.14
This is known in mathematics as the max norm, among other names.15
Another way to find the absolute value of a complex number is the Euclidean norm, which is the length of the vector from the origin of the complex plane at \(0 + 0i\) to the complex value whose absolute value we’re seeking. It’s defined as
\[|a + bi| = \sqrt{a^2 + b^2}\]Complex already has a property that gives the Euclidean norm: the appropriately-named length.
print(Complex(0, -0.25).length)
// -> 0.25
print(Complex(3.14, -3.0).length)
// -> 4.3427640967476
We see that when a complex value has non-zero real and imaginary parts, its length is greater than its magnitude, because both parts contribute to the length.
The max norm and the Euclidean norm are mathematically equivalent for determining boundedness. But it might take an extra iteration for \(f\) to blow up using the max norm instead of the Euclidean norm, in the cases where the Euclidean norm is greater. In part 2, we want to duplicate some results found using the Euclidean norm, so we’ll use it instead of the max norm.
We’ve already defined IteratedFuncValues’ bounded property to use the length property of the Lengthy protocol, so all we have to do is tell Swift that Complex already conforms to Lengthy.
extension ComplexModule.Complex: Lengthy {}
Now we’re using the Euclidean norm to check boundedness of sequences of Complex values. I love it when a plan comes together!
It’s time to meet the Mandelbrot set. Finally! As discussed above, we think the set lives in the rectangle whose upper left corner is at \(-3.0 + 2.0i\) and whose lower right corner is at \(1.0 - 2.0i\). That rectangle is actually a 4 x 4 square. Let’s try sampling it with a grid of 50 x 50 points.
let mandelPoints = ComplexPlotPoints(from: Complex(-3.0, 2.0), to: Complex(1.0, -2.0), horizontalPointCount: 50, verticalPointCount: 50)
let mandelPlot = mandelPoints.map { points in
points.map { point in
plotValue(for: IteratedFuncValues(f: f(c: point), iterations: 30).bounded)
}
}
for row in mandelPlot {
print(String(row))
}
Build and run. If you haven’t seen it before, I’m delighted to introduce you to the Mandelbrot set!
*
***
**
* ******
***********
************
*************
**************
**** **************
********************
*********************
*********************
********************
**** **************
**************
*************
************
***********
* ******
**
***
*
You may find that plot both intriguing and disappointing. How much of each of those you feel will depend on how familiar you are with the Mandelbrot set. For the moment, let’s empty our minds of thoughts of high-resolution color Mandelbrot renderings we may have seen, and just consider the plot we’ve made.
We already expected the plot to be horizontally asymmetric, and it is, but wow. Look how the upper and lower protuberances seem to bend towards the left. The whole thing looks like a winged insect. Who expected a plot of the behavior of a simple iterated function to be so…organic?
But this plot can be improved. For example, it seems we overshot the bounds of the set a bit. There’s a lot of empty space above and below. Let’s zoom in for a closer look. How much can we zoom in? It would be nice to index those axes so we can choose better coordinates. For the real axis, renderIndices() will work unmodified, but we need a real-valued PlotPoints to pass to it. Here’s a method for that.
extension ComplexPlotPoints {
var realPoints: PlotPoints<RealType> {
PlotPoints(from: upperLeft.real, to: lowerRight.real, pointCount: horizontalPointCount)
}
}
For the imaginary axis, a simple solution is to append the corresponding imaginary value to the end of each plotted line. We can do that by using the index of each plot row to access the imaginary value of a point for that row.
for (i, row) in mandelPlot.enumerated() {
print(String(row) + "|\(mandelPoints[i].first!.imaginary)")
}
Now we print the real indices.
let interestingIndices = [-3.0, -2.0, -1.0, -0.75, 0.0, 0.25, 0.5, 1.0]
print(renderIndices(interestingIndices, for: mandelPoints.realPoints))
Build and run.
|2.0
|1.9183673469387754
|1.836734693877551
|1.7551020408163265
|1.6734693877551021
|1.5918367346938775
|1.510204081632653
|1.4285714285714286
|1.3469387755102042
|1.2653061224489797
|1.183673469387755
|1.1020408163265307
|1.0204081632653061
|0.9387755102040818
* |0.8571428571428572
*** |0.7755102040816328
** |0.6938775510204083
* ****** |0.6122448979591837
*********** |0.5306122448979593
************ |0.44897959183673475
************* |0.3673469387755104
************** |0.2857142857142858
**** ************** |0.20408163265306145
******************** |0.12244897959183687
********************* |0.04081632653061229
********************* |-0.04081632653061229
******************** |-0.12244897959183643
**** ************** |-0.204081632653061
************** |-0.2857142857142856
************* |-0.36734693877551017
************ |-0.4489795918367343
*********** |-0.5306122448979589
* ****** |-0.6122448979591835
** |-0.693877551020408
*** |-0.7755102040816326
* |-0.8571428571428568
|-0.9387755102040813
|-1.020408163265306
|-1.1020408163265305
|-1.1836734693877546
|-1.2653061224489792
|-1.3469387755102038
|-1.4285714285714284
|-1.510204081632653
|-1.591836734693877
|-1.6734693877551017
|-1.7551020408163263
|-1.8367346938775508
|-1.9183673469387754
|-1.9999999999999996
|-3.00 |-1.94 |-0.96|-0.47|0.02|0.43|0.84
We can see the set doesn’t extend beyond \(\pm1.0i\). It should fit within a 4 x 2 rectangle with top left coordinate \(-3.0 + 1.0i\).
The real indices are revealing another inadequacy of our current plot. We know the set extends to and includes \(-2.0\) along the real axis. But in this plot the set doesn’t seem to extend past \(-1.5\). Why aren’t we seeing the same result we plotted before for the line of points whose imaginary part is \(0i\)?
The imaginary axis indices tell the story. In the vicinity of \(0i\), we see
******************** |0.12244897959183687
********************* |0.04081632653061229
********************* |-0.04081632653061229
******************** |-0.12244897959183643
We’re not even plotting the line for \(0i\). We start at \(2.0i\), decrement by steps of
\[4.0i \div 49 \approx 0.082i\]and skip right over \(0i\). Intriguingly, the lines we do render are exactly half a step above and below \(0i\), and \(0i\) is exactly halfway between our upper and lower bounds. But we’re rendering an even number of lines, so there can’t be a line between all the others. Floating-point math can be weird, but maybe if we plotted an odd number of lines, the one in the center would be at exactly \(0i\), or close enough.
To zoom in, we’re shrinking the rectangle that is our window onto the complex plane closer to the bounds of the Mandelbrot set. We’re halving the height of the rectangle: switching from a 1:1 to a 2:1 aspect ratio, as discussed above. But we don’t want to sample half as many lines—we want greater resolution, not less. And we want an odd number of lines to try to capture the line with an imaginary coordinate of \(0i\). So we want an odd number of lines greater than 50. Let’s try a grid of 110 x 55.
let moreMandelPoints = ComplexPlotPoints(from: Complex(-3.0, 1.0), to: Complex(1.0, -1.0), horizontalPointCount: 110, verticalPointCount: 55)
let moreMandelPlot = moreMandelPoints.map { points in
points.map { point in
plotValue(for: IteratedFuncValues(f: f(c: point), iterations: 30).bounded)
}
}
for (i, row) in moreMandelPlot.enumerated() {
print(String(row) + "|\(moreMandelPoints[i].first!.imaginary)")
}
print(renderIndices(interestingIndices, for: moreMandelPoints.realPoints))
Build and run.
|1.0
|0.962962962962963
|0.9259259259259259
* |0.8888888888888888
** |0.8518518518518519
***** |0.8148148148148149
****** |0.7777777777777778
****** |0.7407407407407407
***** |0.7037037037037037
**** * |0.6666666666666667
** ************ * |0.6296296296296297
*** *************** * * |0.5925925925925926
********************* *** |0.5555555555555556
*********************** |0.5185185185185186
* ************************ |0.4814814814814815
*************************** |0.4444444444444444
**************************** |0.40740740740740744
***************************** |0.37037037037037046
******************************* |0.33333333333333337
****************************** |0.2962962962962963
** *** * ******************************* |0.2592592592592593
********** ****************************** |0.22222222222222232
********** ******************************* |0.18518518518518523
************ ******************************** |0.14814814814814814
************* ****************************** |0.11111111111111116
** ******************************************** |0.07407407407407418
********************************************** |0.03703703703703709
************************************************************* |0.0
********************************************** |-0.03703703703703698
** ******************************************** |-0.07407407407407396
************* ****************************** |-0.11111111111111116
************ ******************************** |-0.14814814814814814
********** ******************************* |-0.18518518518518512
********** ****************************** |-0.2222222222222221
** *** * ******************************* |-0.2592592592592591
****************************** |-0.2962962962962963
******************************* |-0.33333333333333326
***************************** |-0.37037037037037024
**************************** |-0.40740740740740744
*************************** |-0.4444444444444444
* ************************ |-0.4814814814814814
*********************** |-0.5185185185185184
********************* *** |-0.5555555555555554
*** *************** * * |-0.5925925925925926
** ************ * |-0.6296296296296295
**** * |-0.6666666666666665
***** |-0.7037037037037037
****** |-0.7407407407407407
****** |-0.7777777777777777
***** |-0.8148148148148147
** |-0.8518518518518516
* |-0.8888888888888888
|-0.9259259259259258
|-0.9629629629629628
|-1.0
|-3.00 |-1.97 |-0.98 |-0.72 |0.01 |0.27 |0.52 |1.00
Cool! We can see more detail. The set is looking even more insect-like. And we captured the real number line! It extends to about \(-2.0\) as expected.
At this point, if you hadn’t seen the Mandelbrot set before, you might think you’ve now met it properly. But if you have seen it before, this plot still doesn’t look right. It’s stretched top-to-bottom. The “head” of the insect is supposed to be circular, not an oval. This is because the horizontal spacing of our points is much tighter than their vertical spacing. Our “pixels” aren’t square! That’s what we get for doing our drawing in the console instead of with real graphics APIs.
Our code is flexible enough that we can easily plot what should be a circle of radius \(1\) to demonstrate this aspect ratio problem.
let badCirclePoints = ComplexPlotPoints(from: Complex(-1.1, 1.1), to: Complex(1.1, -1.1), horizontalPointCount: 55, verticalPointCount: 55)
let badCirclePlot = badCirclePoints.map { points in
points.map { point in
plotValue(for: point.length <= 1.0)
}
}
for row in badCirclePlot {
print(String(row))
}
***********
*****************
*********************
*************************
*****************************
*******************************
*********************************
***********************************
*************************************
***************************************
*****************************************
*****************************************
*******************************************
*******************************************
*********************************************
*********************************************
***********************************************
***********************************************
***********************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
*************************************************
***********************************************
***********************************************
***********************************************
*********************************************
*********************************************
*******************************************
*******************************************
*****************************************
*****************************************
***************************************
*************************************
***********************************
*********************************
*******************************
*****************************
*************************
*********************
*****************
***********
This “circle” shows the same vertical stretching as our Mandelbrot plot. We can fix this. It looks like we could fit another line of points between each two we have currently. So the vertical point spacing is roughly twice the horizontal point spacing. We can’t control the vertical spacing of lines in the console, but we can draw twice as many points horizontally to compensate.
let circlePoints = ComplexPlotPoints(from: Complex(-1.1, 1.1), to: Complex(1.1, -1.1), horizontalPointCount: 110, verticalPointCount: 55)
let circlePlot = circlePoints.map { points in
points.map { point in
plotValue(for: point.length <= 1.0)
}
}
for row in circlePlot {
print(String(row))
}
********************
**********************************
********************************************
****************************************************
**********************************************************
**************************************************************
********************************************************************
************************************************************************
****************************************************************************
******************************************************************************
**********************************************************************************
************************************************************************************
**************************************************************************************
****************************************************************************************
******************************************************************************************
********************************************************************************************
**********************************************************************************************
**********************************************************************************************
************************************************************************************************
**************************************************************************************************
**************************************************************************************************
**************************************************************************************************
**************************************************************************************************
****************************************************************************************************
****************************************************************************************************
****************************************************************************************************
**************************************************************************************************
**************************************************************************************************
**************************************************************************************************
**************************************************************************************************
************************************************************************************************
**********************************************************************************************
**********************************************************************************************
********************************************************************************************
******************************************************************************************
****************************************************************************************
**************************************************************************************
************************************************************************************
**********************************************************************************
******************************************************************************
****************************************************************************
************************************************************************
********************************************************************
**************************************************************
**********************************************************
****************************************************
********************************************
**********************************
********************
Perfect!16 Let’s apply the same fix to the Mandelbrot set: we’ll plot 220 x 55 points. We don’t need the indices anymore.
let bestMandelPoints = ComplexPlotPoints(from: Complex(-3.0, 1.0), to: Complex(1.0, -1.0), horizontalPointCount: 220, verticalPointCount: 55)
let bestMandelplot = bestMandelPoints.map { points in
points.map { point in
plotValue(for: IteratedFuncValues(f: f(c: point), iterations: 30).bounded)
}
}
print("Best mandelplot!")
for row in bestMandelplot {
print(String(row))
}
Prepare to meet the Mandelbrot set in about as much glory as the console can muster! Build and run.
*
*
* ***
********** *
* ************
*************
***********
******* *
* * ** ** ********************* ***
***** ******************************* *
******************************************* ******
***********************************************
*************************************************
* ******************************************************
*******************************************************
******************************************************* **
*************************************************************
** * * ************************************************************
***** * **** * ************************************************************
****************** ************************************************************** *
********************** *************************************************************
************************* **************************************************************
*************************** ************************************************************
* * ****************************************************************************************
**********************************************************************************************
****************************************************************************************************************************
**********************************************************************************************
* * ****************************************************************************************
*************************** ************************************************************
************************* **************************************************************
********************** *************************************************************
****************** ************************************************************** *
***** * **** * ************************************************************
** * * ************************************************************
*************************************************************
******************************************************* **
*******************************************************
* ******************************************************
*************************************************
***********************************************
******************************************* ******
***** ******************************* *
* * ** ** ********************* ***
******* *
***********
*************
* ************
********** *
* ***
*
*
Even in our low resolution output device, we can start to see the surprising intricacy of this object. To do better, we’d have to use graphics APIs to draw a bitmap. I plan to do that in a future post.
Remember, this is a plot of the complex values for which a simple iterated function is bounded. Why should it look like this? Why not some basic geometric shape, like a circle? Or even a slightly more complicated shape like a heart? What we see is heart-like, if we squint, but the real Mandelbrot set is infinitely more complicated than a heart. Why should so much complexity arise from repeating a multiplication and an addition?
This mystery, for me, is the enduring appeal of the Mandelbrot set. Which is not to say that it’s quite so mysterious to everyone. Mathematicians have learned a lot about it. And I have learned more about it, too. Remember a long time ago when I said you can use it to compute the digits of \(\pi\)? That’s what I want to talk about in the next part.
Revision History
-
2020-07-06: Corrected errors in, and expanded, the discussion about determining the absolute value of complex numbers.
magnitudeandlengthare equivalent for determining boundedness; we choose to uselengthfor reasons relevant to part 2. Thanks to Steve Canon for the correction. -
2020-07-05: Added
Lengthyprotocol so the boundedness check can use the Euclidean norm for complex values, instead of the ∞-norm thatComplex’smagnitudegives. Changed name of propertytwoOrLessto the more-accuratelengthTwoOrLess. -
2020-06-26: Fixed an infinite loop that occurred in
renderIndices()if the largest index was less than the largest point. Didn’t happen in this post, but it might have done in a future post. -
2020-06-24: Switched from awkward ranged subscript to
suffix()for finding last several values of a sequence. -
2020-06-24: Added
realPointsproperty toComplexPlotPointsfor much more pleasant rendering of real axis indices on complex plots. -
2020-06-23: Added information about required versions of Swift and Xcode due to the use of key-path expressions as functions, and a workaround if that feature isn’t available.
-
2020-06-22: Rewrote all code in a functional and protocol-oriented style. Thanks again to Doug Lovell for encouraging me to do this.
-
2020-06-11: Changed “x” to “c” in
renderBehaviorOfF(),renderXIndices(), top-level code, and prose, and changed a reference to “the x axis” to “the number line”. We’re plotting different values of \(c\) on the real number line and the complex plane, not different values of \(x\) on the \(x\)-axis. Thanks to Doug Lovell for the correction.
-
Specifically Oh! Pascal!: Turbo Pascal 6.0, Third Edition. ISBN 0-393-96077-3 (you know, the one bundled with the 3.5” rather than 5.25” floppy). ↩
-
All the code in this post was developed and tested in Xcode 11.5. I use at least one language feature that requires Swift 5.2, which first shipped in Xcode 11.4. ↩
-
Spoiler alert: we do! ↩
-
Protocols are an important concept in Swift that I won’t explain in detail here. See The Swift Programming Language for an introduction. ↩
-
Just one of the many aspects of Swift generics I want to learn more about. ↩
-
Sequencehas more requirements, but if we implementnext(), the compiler can synthesize the rest of them. ↩ -
Another strategy would be to allow
IteratedFuncValuesto be iterated infinitely with its built-in iterator, and define another, finite iterator that can be used when a specific number of iterations is desired. I think this approach would be more complicated to write and use, but I haven’t tried it yet, as this post needs to be completed at some finite point in the future! ↩ -
See The Swift Programming Language for an introduction to optional values. ↩
-
which is only \(6x10^{79}\), according to Wolfram Alpha. ↩
-
Spoiler alert: it will be! ↩
-
There may be a more generic protocol in the Standard Library that enables division.
FloatingPointis what Xcode suggested, and I didn’t search further. Soon enough, we’ll have reason to seek a different protocol for division, which we’ll find outside of the standard Standard Library. ↩ -
In fact, watch every 3Blue1Brown video! ↩
-
It’s been announced that playgrounds will allow the use of Swift packages in Xcode 12. Cool! ↩
-
See swift-numerics’ README for further discussion on this topic. ↩
-
I’m out of my mathematical depth here. Corrections are welcome. ↩
-
It depends greatly on what font’s being used to display the plot. The corrected circle plot looks perfect to me in the Xcode debug console, but still a bit wonky when I view this post in Safari. ↩