Oct 16 2011

Pi is wrong! Long live Tau!

Category: Daily Events,Mathematics

At one point or another, we’ve all had a feeling that something is not quite right in the world. It’s a huge relief, therefore, to discover someone else who shares your suspicion. (I’m also surprised that it’s taken me this long to stumble on this!)

It has always baffled me why we define π to be the ratio of the circumference of a circle to its diameter, when it should clearly be the ratio of the circumference to its radius. This would make π become the constant 6.2831853..., or 2 times the current definition of π.

Why should we do this? And what effect would this have?

Well, for starters, this would remove an unnecessary factor of 2 from a vast number of equations in modern physics and engineering.

Most importantly, however, this would greatly improve the intuitive significance of π for students of math and physics. π is supposed to be the “circle constant,” a constant that embodies a very deep relationship between angles, radii, arc lengths, and periodic functions.

The definition of a circle is the set of points in a plane that are a certain distance (the radius) from the center. The circumference of the circle is the arc length that these points trace out. The circle constant, therefore, should be the ratio of the circumference to the radius.

To avoid confusion, we’ll use the symbol tau (τ) to be our new circle constant (as advocated by Michael Hartl, from the Greek τόρνος, meaning “turn”), and make it equal to 6.283..., or 2π.

In high school trigonometry class, students are required to make the painful transition from degrees to radians. And what’s the definition of a radian? It’s the ratio of the length of an arc (a partial circumference) to its radius! Our intuition should tell us that the ratio of a full circumference to the radius should be the circle constant.

Instead, students are taught that a full rotation is 2π radians, and that the sine and cosine functions have a period of 2π. This is intuitively clunky and fails to illustrate the true beauty of the circle constant that π is supposed to be. This is surely part of the reason that so many students fail to grasp these relationships and end up hating mathematics. A full rotation should be τ radians! The period of the sine and cosine functions should be τ!

But... wouldn’t we have to rewrite all of our textbooks and scientific papers that make use of π?

Yes, we would. And, in doing so, we would make them much easier to understand! You can read the Tau Manifesto website to see examples of the beautiful simplifications that τ would bring to mathematics, so I won’t repeat them here. You can also read the original opinion piece by Bob Palais that explores this subject.

It’s not particularly surprising that the ancient Greeks used the diameter of a circle (instead of the radius) in their definition of π, since the diameter is easier to measure, and also because they couldn’t have foreseen the ubiquity of this constant in virtually all sciences.

However, it’s a little unfortunate that someone like Euler, Leibniz, or Bernoulli didn’t pave the way for redefining π to be 6.283..., thus missing the opportunity to simplify mathematics for generations to come.

Aside from all the aesthetic improvements this would bring, considering how vitally important it is for more of our high school students (and beyond) to understand and appreciate mathematics, we need all the “optimizations” we can get to make mathematics more palatable for them. This surely has to be an optimization to consider seriously!

From now on, I'm a firm believer in tauism! Are you?

Aug 14 2009

The Math Book: Get it Now!

Category: Daily Events,Mathematics,Programming

Cliff Pickover, the prolific author of more than forty popular science and mathematics books, has outdone himself with his latest compilation: The Math Book. This is a collection of 250 "milestones" of mathematics throughout history, complete with breathtaking glossy color illustrations for each entry (a first for his books), as well as insightful descriptions that explain the history and the significance of each of these marvels of mathematics.

This book is especially significant in one other way: it contains my artwork! The book's entry on Knight's Tours (p. 186) familiarizes the reader with the history of this problem, dating all the way back to Euler in 1759. And, alongside the article, Pickover displays a 30x30 knight's tour that was solved by my neural network knight's tour implementation. For the picture in the book, I used a modified version of the program that generated a sufficiently hi-res image. That particular knight's tour took about 3 days for my computer to generate.

I'm deeply grateful to have one of my creations published in a book by someone as influential as Cliff Pickover. Of course, it's all of the 250 entries in the book that make it an incredibly fascinating stroll through the history of mathematics. As mentioned elsewhere, this book definitely has bestseller potential, and could easily be one of Pickover's best works. Buy the book now!

Apr 29 2008

Binomial Coefficients and Stirling Numbers in C#

Category: Mathematics,Programming

As long as I'm on a roll with my posts on number theory in C#, I thought I'd briefly discuss how to generate binomial coefficients, and then move on to Stirling numbers of the second kind, both of which are extremely useful in combinatorics and finite calculus.

Some prerequisites for calculating binomial coefficients include a standard factorial function, nothing special:

1
2
3
4
5
6
7

long factorial(long n)
{
    if (n == 0) return 1;
    long t = n;
    while(n-- > 2) t *= n;
    return t;
}

Note that I use longs whenever possible because, despite the performance hit from 64-bit operations, it's worth it to be able to work with numbers that are double the magnitude of ints.

We will also need a function for calculating a falling power of a number, which is defined as: \(x^{\underline{n}} = x(x-1)(x-2)\ldots(x-(n-1))\). You'll see why in a moment.

1
2
3
4
5
6

long fallingPower(long n, long p)
{
    long t = 1;
    for (long i = 0; i < p; i++) t *= n--;
    return t;
}

Notice that the above function handles the case \(n^{\underline{0}}\), returning 1. However, it does not handle the case of p > n, which would give an incorrect result, but this is not necessary for our purposes.

Binomial Coefficients

Recall that the definition for the binomial coefficient is \({n \choose k} = \frac{n!}{k!(n-k)!}\). However, using this exact formula to compute binomial coefficients is a bit naive. If we use falling powers (sometimes called falling factorials), the above formula easily reduces to: \(\frac{n^{\underline{k}}}{k!}\).

We can improve the algorithm a bit more by adding the condition:

\[{n \choose k} =
\begin{cases}
\frac{n^{\underline{k}}}{k!} \quad \mbox{if } k \leq \lfloor n/2 \rfloor,\\
\frac{n^{\underline{n-k}}}{(n-k)!} \quad \mbox{if } k > \lfloor n/2 \rfloor.
\end{cases}
\]

Not only is this algorithm faster, but it can also handle larger coefficients than the original formula, since neither the falling power nor the factorial ever gets larger than n/2.
The code for this is straightforward:

1
2
3
4
5
6

long binomialCoeff(long n, long k)
{
    if ((k < 0) || (k > n)) return 0;
    k = (k > n / 2) ? n - k : k;
    return fallingPower(n, k) / factorial(k);
}

However, this is still not as optimal as it can be. The most optimal approach would be to accumulate the falling power while dividing by each factor of the factorial in place. This would minimize the chance of overflow errors, and allow for even larger coefficients to be calculated. The disadvantage of this algorithm is the necessary use of floating-point math:

1
2
3
4
5
6
7

long binomialCoeff(long n, long k) {
    if ((k < 0) || (k > n)) return 0;
    k = (k > n / 2) ? n - k : k;
    double a = 1;
    for (long i = 1; i <= k; i++) a = (a * (n-k+i)) / i;
    return a + 0.5;
}

Stirling Numbers

Stirling numbers (of the second kind) are useful for, among other things, enumerating the coefficients of the falling-power expansion of a regular power. For example, how would we express x3 in terms of falling powers of x? That is, how do we arrive at an equation of the form \(x^3 = ax^{\underline 3} + bx^{\underline 2} + cx^{\underline 1}\)? Well, we could just solve the equation directly, but this would get unwieldy for higher powers. A neat way of doing this involves the use of Stirling numbers of the second kind, \(\left\{\begin{matrix} n \\ k \end{matrix}\right\}\). A useful theorem for computing these numbers is

\[\left\{\begin{matrix} n \\ k \end{matrix}\right\}=\frac{1}{k!} \sum_{i=0}^{k}(-1)^i{k \choose i}(k-i)^n\]

1
2
3
4
5
6
7
8
9
10
11

long stirling(long n, long k)
{
    long sum = 0, neg = 1;
    for (long i = 0; i <= k; i++)
    {
        sum += neg * binomialCoeff(k, i) * pow(k - i, n);
        neg = -neg;
    }
    sum /= factorial(k);
    return sum;
}

With the Stirling numbers in hand, we can now obtain the coefficients for falling power expansions:

\[x^m = \sum_{k=0}^{m}\left\{\begin{matrix} m \\ k \end{matrix}\right\}x^{\underline k}\]

Apr 28 2008

Project Euler -- Problem 197

Category: Mathematics,Programming

I never thought that solving random math problems can be addictive, but Project Euler does exactly that. Not only does it make you flex your math muscles, it also challenges you to take your programming language of choice to its limits. So it's a total win-win: you brush up on your math skills, and broaden your programming repertoire at the same time.

Once I discovered Project Euler, I couldn't pull myself away from the computer until I solved as many problems as I could. As of this writing I've solved 187 of their 211 problems.

This will be the first of a series of posts about particular Euler problems, and how easily they can be solved with C# and a bit of common-sense math. The problems I discuss will not be in any particular order, mind you.

Problem 197 has to do with finding the nth term of a particular recursively defined sequence: \(u_{n+1} = f(u_n)\), with \(u_0 = -1, f(x) = \lfloor 2^{30.403243784-x^2}\rfloor \cdot 10^{-9}\).

Of course, as with most Euler problems, the value for n is set ridiculously high, presumably to eliminate the possibility of brute-forcing the problem (within the lifetime of the universe).

Fortunately, it takes little more than a superficial examination to see that this problem is actually quite simple in disguise. If we look at the first few terms of the sequence, we can already guess that this sequence appears to be "converging" to a function that oscillates between two values, approximately 0.681 and 1.029 (I won't give precise numbers, since that would give away the solution).

This means that all we need to do is go far enough into the sequence that the deviation of the oscillations is less than the desired precision asked by the problem (10-9). And it so happens that we don't need to go out far at all. The sequence actually settles on its two oscillatory values as early as the 1000th term (probably even earlier)! Therefore, the sum of the 1000th and 1001st term will be equivalent to the sum of the 1012th and (1012+1)st term, which is what the problem asks for!

The code to do this is elementary. I accomplished it with a mere 5 lines of C# code. Can you do better?

Apr 23 2008

Elementary Number Theory in C#

Category: Mathematics,Programming

I thought I'd post a few code snippets in C# that have to do with basic number theory, since I use them in my programs from time to time. These snippets are by no means optimized, and the use of C# pretty much precludes their use in high-performance applications. Still, for relatively small arguments, these routines run surprisingly fast.

Prime numbers

To generate a list of prime numbers, we use the familiar Sieve of Eratosthenes. We can write one routine to generate the sieve:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

bool[] GetPrimeSieve(long upTo)
{
    long sieveSize = upTo + 1;
    bool[] sieve = new bool[sieveSize];
    Array.Clear(sieve, 0, (int)sieveSize);
    sieve[0] = true;
    sieve[1] = true;
    long p, max = (long)Math.Sqrt(sieveSize) + 1;
    for (long i = 2; i <= max; i++)
    {
        if (sieve[i]) continue;
        p = i + i;
        while (p < sieveSize) { sieve[p] = true; p += i; }
    }
    return sieve;
}

The above function returns an array of booleans (the size of the given parameter), each of which is false if the array index is a prime number, or true if it's not a prime number. To get an actual list of prime numbers, we can use a function such as this:

1
2
3
4
5
6
7
8
9
10
11
12

long[] GetPrimesUpTo(long upTo)
{
    if (upTo < 2) return null;
    bool[] sieve = GetPrimeSieve(upTo);
    long[] primes = new long[upTo + 1];
 
    long index = 0;
    for (long i = 2; i <= upTo; i++) if (!sieve[i]) primes[index++] = i;
 
    Array.Resize(ref primes, (int)index);
    return primes;
}

The above function returns an actual list of primes less than or equal to the specified high limit. This means that we can easily compute the prime-counting function \(\pi(x)\) for any integer x by counting the number of elements in the array returned by the function. We can write a similar function to generate a list of composites:

1
2
3
4
5
6
7
8
9
10
11
12
13

long[] GetCompositesUpTo(long upTo)
{
    if (upTo < 2) return null;
    bool[] sieve = GetPrimeSieve(upTo);
    long[] composites = new long[upTo + 1];
 
    long index = 0;
    for (long i = 2; i <= upTo; i++)
        if (sieve[i]) composites[index++] = i;
 
    Array.Resize(ref composites, (int)index);
    return composites;
}

As for determining if a single certain number is prime (without having to generate a giant sieve), we can simply use a function that attempts to factor the number. If the number happens to be divisible by an integer greater than 1 and less than or equal to its square root, then the number is not prime:

1
2
3
4
5
6
7
8

bool IsPrime(long n)
{
    long max = (long)Math.Sqrt(n);
    for (long i = 2; i <= max; i++)
        if (n % i == 0)
            return false;
    return true;
}

Greatest Common Divisor and Totient

To obtain the greatest common divisor (GCD) of two numbers, we use the usual Euclidean algorithm:

1
2
3
4
5
6
7
8
9
10
11

long gcd(long a, long b)
{
    long temp;
    while (b != 0)
    {
        temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

There is a slightly different binary algorithm for computing the GCD which is theoretically more efficient, but in practice (at least in C#) it's actually slightly less efficient than the algorithm above:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

ulong gcdb(ulong a, ulong b)
{
    int shift;
    if (a == 0 || b == 0) return a | b;
 
    for (shift = 0; ((a | b) & 1) == 0; ++shift) { a >>= 1; b >>= 1; }
    while ((a & 1) == 0) a >>= 1;
 
    do {
        while ((b & 1) == 0) b >>= 1;
 
        if (a < b) {
            b -= a;
        } else {
            ulong temp = a - b;
            a = b;
            b = temp;
        }
        b >>= 1;
    } while (b != 0);
    return a << shift;
}

With the GCD readily available, determining whether two numbers are coprime is just a matter of telling whether or not their GCD is equal to 1. Also, calculating the LCM (least common multiple) of two numbers becomes trivial:

1
2
3
4
5
6
7
8
9
10

long lcm(long a, long b)
{
    long ret = 0, temp = gcd(a, b);
    if (temp != 0)
    {
        if (b > a) ret = (b / temp) * a;
        else ret = (a / temp) * b;
    }
    return ret;
}

Also using the GCD algorithm, it becomes easy to calculate Euler's totient function for a certain number, since the totient function is simply the number of integers less than or equal to n that are coprime to n:

1
2
3
4
5
6
7

long eulerTotient(long n)
{
    long sum = 0;
    for (long i = 1; i <= n; i++)
        if (gcd(i, n) == 1) sum++;
    return sum;
}

The above is a really naive algorithm. A much more efficient algorithm (one that is usually given in textbooks) is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

int eulerTotient(int n)
{
    primes = GetPrimesUpTo(n+1);    //this can be precalculated beforehand
    int numPrimes = primes.Length;
 
    int totient = n;
    int currentNum = n, temp, p, prevP = 0;
    for (int i = 0; i < numPrimes; i++)
    {
        p = (int)primes[i];
        if (p > currentNum) break;
        temp = currentNum / p;
        if (temp * p == currentNum)
        {
            currentNum = temp;
            i--;
            if (prevP != p) { prevP = p; totient -= (totient / p); }
        }
    }
    return totient;
}

Jan 24 2007

Hyperbolic Tessellations

Category: Mathematics,Windows Programming

A tessellation refers to a uniform tiling of a plane with polygons, such that an equal number of identical polygons meet at each vertex. For example, the tiles in a bathroom, the squares of linoleum on an office floor, or the honeycomb pattern in a bees' nest are all tessellations of the Euclidean plane.

Hyperbolic TessellationHowever, tessellations are also possible on non-Euclidean spaces, such as the elliptic plane (like the stitching pattern on a soccer ball), and the hyperbolic plane (like... nothing you'd find around the house). In fact, the Euclidean plane has only three regular tessellations (with squares, hexagons, and triangles), while the hyperbolic plane can be tessellated in infinitely many ways.

Since we do not exist in hyperbolic space, we cannot truly "see" hyperbolic tessellations. We can only "represent" them in Euclidean form. A common way of doing this is on the Poincaré disk, which is a finite circle that represents the boundary of the (infinite) hyperbolic plane that is contained inside. The image on the right is a hyperbolic tessellation drawn on the Poincaré disk.

Since tessellations of the hyperbolic plane are especially interesting and mesmerizing to look at, I wrote a small program that generates them, with a great deal of configurable options.

Download the program! | Browse source code repository

Tessellation Application

Using the Program

The program allows you to create an unlimited number of tessellations by selecting "File -> New" from the menu.

When viewing a tessellation, click-and-drag inside of it to shift its position within hyperbolic space. You can also click-and drag with the right mouse button to manipulate the truncation of the tessellation.

Tessellation ControlsThe "Tessellation Controls" window allows you to change the settings for the tessellation that is currently active.

  • p and q -- The numbers specified by p and q refer to the Schläfli symbol {p,q} of the tessellation. The Schläfli symbol is a simple way of classifying tessellations where p is the number of sides in each polygon, and q is the number of polygons that meet at each vertex.
  • Max. Vertices -- This specifies the number of vertices that will be drawn (how far the tessellation will extend towards the disk boundary). More vertices will take exponentially longer to draw. Also, with more vertices, clicking-and-dragging the tessellation will become slower. With a good screen resolution, 10000 vertices fills up the Poincaré disk almost completely.
  • Model -- Select "Poincaré" to draw the tessellation on a Poincaré-style disk, and "Klein" to draw on a Klein-style disk. The Klein disk is similar to the Poincaré disk, except the Klein disk transforms hyperbolic space so that a line between two points appears as a straight line, instead of a circle arc, which is what the Poincaré disk gives.
  • Quality -- Select "Low" to display simple straight lines between vertices (and let the tessellation be drawn much faster). Select "High" to draw actual curved lines between vertices. This will slow down drawing considerably.
  • Truncation -- This is a list of predefined levels of truncation for the tessellation. Select from this list to apply a certain truncation. You can also do free-form truncation by clicking-and-dragging on the tessellation with the right mouse button.
  • Colors -- This allows you to select different colors for each of the components of the tessellation, and to enable or disable drawing of each component.
  • Driver -- This selects what functions the program will use to draw the tessellation. Select "OpenGL" to use OpenGL technology, or "Windows GDI" to use plain Windows functions. In most cases, selecting OpenGL will enable the images to be drawn considerably faster, especially with Antialiasing enabled. However, OpenGL is not supported on some (very) old graphics cards. Also, if you create multiple tessellations, only one can be drawn with OpenGL at any given time.
  • Antialias lines -- Check this box to draw "smooth" lines.
  • Line Thickness -- This specifies the thickness (in pixels) of the lines that make up the tessellation.
  • Advanced -- These options are mostly experimental and will not be discussed here.

Gallery

Here's a brief collection of images created using this program. Click on an image to view a larger version.

{7,3} tessellations, with various truncation:
none, (0,1,0), (0,.5,.5), runcinated, omnitruncated, and snub.
No Truncation(0,1,0) Truncation(0,.5,.5) Truncationruncinatedomnitruncatedsnub

Links

  • This program borrows a substantial amount of code from Don Hatch's page, which has an exhaustive gallery of {p,q} permutations and truncations, as well as a tessellation Java applet.
  • David E. Joyce's tessellation page at Clark University.

Jan 05 2007

Ulam's Prime Number Spiral

Category: Mathematics,Windows Programming

There is an infinite number of prime numbers, and yet the prime numbers themselves do not display any apparent pattern, nor does any formula exist that generates prime numbers. In fact, Legendre proved that there cannot be an algebraic function which always gives primes.

Prime SpiralHowever, prime numbers do exhibit a curious phenomenon when arranged in a spiral along with other consecutive integers, as in the figure to the right (in the figure, prime numbers are highlighted in white, twin primes are green, and Mersenne primes are red).

The Phenomenon

It was first noticed by the physicist Stanisław Ulam in 1963, when he got bored in a meeting and started doodling spirals of numbers. He noticed that, if he makes a spiral of consecutive integers, and circles only the prime numbers, strange diagonal "lines" of prime numbers emerge.

This is quite surprising, since we would intuitively expect a random distribution of prime numbers. However, these diagonal segments occur on an impressively large scale, and arbitrarily far from the center of the spiral. The following image is a spiral containing about 4000 primes, and next to it is the same image with some of the diagonal paths highlighted.
Prime Spiral

Application

Prime Spiral Application
To explore this phenomenon on a large scale, I wrote a small program that generates arbitrarily large spirals, with configurable coloring and other options.

Download the program! | Browse the source code repository

The program generates a spiral based on the total number of integers that you specify. It also allows you to specify the colors to use for the background, prime numbers, twin primes, and Mersenne primes.

In addition, the program allows you to use a custom polynomial (up to degree-2) for generating the spiral. By default, the polynomial is set to

f(n) = n

so the spiral will have the normal sequence 0, 1, 2, 3, etc. However, as an example, suppose you want the spiral to consist only of odd integers (1, 3, 5, 7, etc). You can simply set the polynomial to be

f(n) = 2n + 1

Odd Prime Spiral
by writing "2" in the text box next to "n", and "1" in the last of the three text boxes. The image to the right shows a spiral constructed from odd integers only. Notice the prominent "patterns," this time extending vertically and horizontally.

Conclusions

Ultimately, all that these patterns show is that certain polynomials are more "likely" to generate prime numbers than others. In fact, we can speculate that these kinds of patterns would emerge if we arrange the integers in any ordered design, not just a spiral. Even if we arrange the integers in a simple table, we would still see occasional "streaks" of prime numbers similar to the ones seen in the spiral.

The existence of "prime-generating" polynomials was known since the time of Euler, who discovered a polynomial that gives 40 consecutive prime numbers, namely

f(n) = n2 - n + 41

The reason why some polynomials generate more primes than others is still not known.

Extreme Spirals

Using my program, you can generate extremely large spirals, limited only by the amount of memory on your computer. Here are some fairly large ones that I generated:

Links