The Math Less Traveled

Combinatorial proofs

Posted on February 17, 2012 by Brent

Continuing from a previous post, we found that if we begin with $n$ th powers of $(n+1)$ consecutive integers and then repeatedly take successive differences, it seems like we always end up with the factorial of $n$ , that is, $n! = 1 \cdot 2 \cdot 3 \cdot \dots \cdot n$ . We then derived an algebraic expression for the result of the iterative difference procedure. So the goal now is to prove that

$\displaystyle \sum_{i=0}^n (-1)^{n-i} \binom n i (k+i)^n = n!$

that is,

$\displaystyle (-1)^n \binom n 0 k^n + (-1)^{n-1} \binom n 1 (k+1)^n + \dots + \binom n n (k+n)^n = n!$

Now, it’s possible (probable, in fact) that this can be proved by purely algebraic means. If you come up with such a proof I would love to see it! But I must confess that I spent several hours banging my head against it (algebraically speaking) without making any progress. Eventually I turned to the idea of a combinatorial proof.

What do I mean by that? Combinatorics is the subfield of mathematics concerned with counting. The essence of a combinatorial proof is to show that two different expressions are just two different ways of counting the same set of objects—and must therefore be equal. I’ve described some combinatorial proofs before, in counting the number of ways to distribute cookies.

As another simple example, consider the binomial coefficient identity

$\displaystyle \binom n k = \binom {n-1}{k} + \binom{n-1}{k-1}$

It’s certainly possible to prove this algebraically, by expanding out the binomial coefficients (using $\binom n k = \frac{n!}{k!(n-k)!}$ ), but we can give a more elegant proof, based on the fact that $\binom n k$ is the number of ways to choose a subset of $k$ things out of a set of $n$ things. For example, here are the $\binom 5 3 = 10$ ways to choose three things out a set of five:

Consider the first element of the size- $n$ set. Every subset of size $k$ either includes this first element, or it doesn’t. The number of size- $k$ subsets which do not include the first element is $\binom{n-1}{k}$ , since that’s the number of ways to choose $k$ things from the remaining $n-1$ elements. The number of size- $k$ subsets which do include the first element is $\binom{n-1}{k-1}$ , because they correspond to choosing $k-1$ of the remaining $n-1$ things. Therefore $\binom n k = \binom{n-1}{k} + \binom{n-1}{k-1}$ .

Here’s an illustration of how this works in the particular case when $n = 5$ and $k = 3$ :

Notice how the ten subsets from above have been split into two groups: the first group, on the left, are those that don’t include the first element; you can see that each of them corresponds to one of the $\binom 4 3 = 4$ size- $3$ subsets of the remaining four elements. The second group, on the right, are those that do include the first element; each corresponds to one of the $\binom 4 2 = 6$ size- $2$ subsets of the remaining four elements.

So that’s the idea of a combinatorial proof. And we want to do something similar for the identity we are trying to prove—although, of course, it’s going to be a bit more difficult!

You might have fun trying to think about what a combinatorial proof of our target equation might look like; although if you don’t have much experience with combinatorics you may have trouble coming up with what sorts of things the two sides of the equation might be counting! That’s what I’ll talk about in my next post.

Posted in combinatorics, pictures, proof | 7 Comments

Differences of powers of consecutive integers, part II

Posted on February 16, 2012 by Brent

If you spent some time playing around with the procedure from Differences of powers of consecutive integers (namely, raise $n+1$ consecutive integers to the $n$ th power, and repeatedly take pairwise differences until reaching a single number) you probably noticed the curious fact that it always seems to result in a factorial—in the factorial of $n$ , to be precise.

For example:

$3^2, 4^2, 5^2 = 9,16,25 \to 7,9 \to 2$

$196^2, 197^2, 198^2 = 38416,38809,39204 \to 393,395 \to 2$

$7^3, 8^3, 9^3, 10^3 = 343,512,729,1000 \to 169,217,271 \to 48,54 \to 6$

$0^4,1^4,2^4,3^4,4^4 = 0,1,16,81,256 \to 1,15,65,175 \to 14,50,110 \to 36,60 \to 24$

Several commenters figured this out as well. Does this always happen? If so, can we prove it?

Let’s start by thinking about what happens when we do the successive-differencing procedure. If we start with the list $a,b$ , then we get $-a + b$ . (I want to keep the letters in order, which is why I wrote $-a + b$ instead of $b - a$ . Instead of subtracting the first value from the second, we can think of it as adding the negation of the first value to the second.) If we start with $a,b,c$ , we get

$a,b,c \to (-a+b), (-b+c) \to (a - 2b + c)$ .

(The negation of $(-a + b)$ is $(a - b)$ ; adding this to $(-b + c)$ yields $(a - 2b + c)$ .) From $a,b,c,d$ we get

$a,b,c,d \to \\ (-a+b), (-b+c), (-c+d) \to \\ (a - 2b + c), (b - 2c + d) \to \\ (-a + 3b - 3c + d).$

Do you see any patterns yet? Let’s do one more. From the above calculation we can already see that doing four iterations on $a,b,c,d,e$ will result in $(-a + 3b - 3c + d), (-b + 3c - 3d + e)$ (do you see why?). Doing one final iteration gives us

$a - 4b + 6c - 4d + e$ .

Hmm. Let’s make a table.

$n$	Result
1	$-a + b$
2	$a - 2b + c$
3	$-a + 3b - 3c + d$
4	$a - 4b + 6c - 4d + e$
5	$-a + 5b - 10c + 10d - 5e + f$

I included one more row (which you can verify if you like). Now do you see a pattern? The coefficients seem to be taken from Pascal’s triangle, but with alternating signs!

In fact, it’s actually not too hard to see why this happens. At each step we take two offset copies of the previous row (by "offset" I mean that the letters are shifted by one, like $(a - 2b + c)$ and $(b - 2c + d)$ ) and add the negation of the first to the second. Since the signs are alternating, we really end up adding absolute values of the coefficients. Probably the best way to really see this is through an example:

$\begin{array}{cccccccccc} -( & -a & + & 3b & - & 3c & + & d & ) & \\ & & & -b & + & 3c & - & 3d & + & e \\ \hline & a & - & 4b & + & 6c & - & 4d & + & e \end{array}$

Notice how we flip all the signs in the first row, so that they match the signs in the second row. But this is exactly how Pascal’s triangle is generated—each row is the sum of the previous row with itself, offset by one.

Now, in the real problem, we don’t start with $a,b,c,d,\dots$ , but with the $n$ th powers of $n+1$ consecutive integers. Let’s call the first integer $k$ , so the sequence of consecutive integers is $k, k+1, k+2, \dots, k+n$ . Given this, we can now write down an expression for what we end up with after doing the iterated difference procedure:

$\displaystyle \sum_{i=0}^n (-1)^{n-i} \binom{n}{i} (k+i)^n$

Let’s break this down a bit. We know that we get a sum of $n+1$ terms; that’s the $\sum_{i=0}^n$ part (you can read more about sigma notation here). We’ll use $i$ to index the terms. We also know that the terms alternate sign, so we need to include $(-1)$ raised to some power involving $i$ ; the last term is always positive, so we use $(-1)^{n-i}$ , which is equal to $1$ when $i = n$ . The binomial coefficient $\binom n i$ gives us the $i$ th entry on the $n$ th row of Pascal’s triangle. Finally, of course, $(k+i)^n$ is the $n$ th power of one of the integers we started with.

The claim, therefore, is that

$\displaystyle \sum_{i=0}^n (-1)^{n-i} \binom n i (k+i)^n = n!$

And I will prove it to you, with pretty pictures, as promised!

Posted in arithmetic, iteration, pascal's triangle | Tagged binomial coefficients, consecutive, difference, integers, powers | 1 Comment

17×17 4-coloring with no monochromatic rectangles

Posted on February 9, 2012 by Brent

Quick, what’s special about the following picture?

As just announced by Bill Gasarch, this is a $17 \times 17$ grid which has been four-colored (that is, each point in the grid has been assigned one of four colors) in such a way that there are no monochromatic rectangles, that is, no four grid points forming the corners of an axis-aligned rectangle are all of the same color. For example, if we change the top-left grid point to red, we can see several monochromatic rectangles pop up:

Or here’s another version where I randomly picked a grid point in the middle, changed its color, and sure enough, more monochromatic rectangles result:

As you can try verifying for yourself (and as I also verified using a computer program), there are no such monochromatic rectangles in the four-coloring at the top of this post! (If you want to play with the four-coloring yourself, here it is in a simple data format.)

For several years no one knew if this was possible, and Bill had offered a prize of $289 (that’s $\$17^2$ , of course) to anyone who could find such a four-coloring. The prize will be collected by Bernd Steinbach and Christian Posthoff—you can find more details in Bill’s post. No one yet knows exactly how they found their four-coloring, but apparently they will be presenting a paper about it in May. I’ll try to write more about it then (if I understand it at all)!

If you’re interested in reading more about the history and math behind this problem (and to get some intuition for why it is difficult), take a look at these posts by Brian Hayes on his blog, bit-player: The 17×17 challenge and 17 x 17: A nonprogress report. ~~I also remember seeing~~ Here’s a fun interactive applet where you can play around with the problem, created by Martin Schweitzer.

Posted in open problems, pattern, people, pictures | Tagged 17x17, four-coloring, graph, grid, monochromatic, rectangles | 4 Comments

Book review: Nine Algorithms that Changed the Future

Posted on February 4, 2012 by Brent

Nine Algorithms that Changed the Future: the Ingenious Ideas that Drive Today’s Computers, by John MacCormick. Princeton University Press, 2012.

I’m often wary of books written for general audiences on technical topics. It’s quite difficult to write in a way that is both accessible to a wide audience and technically accurate. Many such books end up sacrificing accuracy in the name of accessibility, trying to convey just the “intuition” or “general sense” of some topic, but often end up giving people the wrong idea instead.

I was quite happy to find, therefore, that John MacCormick nails it: “9 Algorithms that Changed the Future” is technically right on the money, but manages to explain things in ways that are both understandable and fun. Want to understand how Google ranks search results? Or how Amazon manages to never lose or mess up your order information, even though they get hundreds of thousands of orders each day and (as we all know) networks and hard drives are unreliable? Ever wonder how you can order something over the internet without your credit card number being stolen? Or how “zip” is able to make your files smaller, seemingly by magic? Even if you have never wondered about these things, perhaps I have made you wonder about them just now. And that’s exactly the point of this book: there are quite a few ingenious algorithmic ideas that most of us rely on every day that we rarely—or never—even stop to wonder about.

For example, I actually learned something new: I knew about public-key cryptography but had never really known much about Diffie-Hellman key exchange, which is what allows your web browser to talk to, say, Amazon’s servers securely even though they have never communicated before. It’s like having a secret conversation in code with a pen-pal whom you’ve never met, even though lots of people are reading your mail. How can you ever agree on a secret code in the first place without the people reading your mail finding out (and hence being able to read all your subsequent coded messages)? Sounds impossible, doesn’t it? But it turns out that it is possible, with some clever ideas, which MacCormick skillfully explains using a fun metaphor about mixing colors of paint.

Each chapter starts out very simply, gradually building up more complex examples until you reach a full understanding of the algorithm being explained. Along the way MacCormick introduces the “tricks”—the clever, central ideas—that make each algorithm work. The writing is excellent: clear, precise, and fun. I highly recommend this book to anyone curious about the ingenious mathematical and algorithmic ideas underlying some of today’s most ubiquitous technology.

Posted in books, computation, review | Tagged algorithms, history, John MacCormick | Leave a comment

Computing with decadic numbers

Posted on January 30, 2012 by Brent

[This is the ninth, and, I think, final in a series of posts on the decadic numbers (previous posts: A curiosity, An invitation to a funny number system, What does "close to" mean?, The decadic metric, Infinite decadic numbers, More fun with infinite decadic numbers, A self-square number, u-tube).]

In a previous post, we found a decadic number

$u = \dots 57423423230896109004106619977392256259918212890625$

with the curious property that it is its own square, even though it is obviously not zero or one. We then derived a more efficient algorithm for generating the digits of $u$ . Here’s some Haskell code (explained in the previous post) which implements the more efficient algorithm, which I include here just so that this post will be a valid literate Haskell file in its entirety.

> {-# LANGUAGE TypeSynonymInstances
>            , FlexibleInstances
>   #-}
>
> module Decadic2 where
>
> import Control.Monad.State
>
> -- State for incrementally constructing u_n.
> -- Invariant: curT = 10^n; un^2 = pn*curT + un
> data UState = UState { pn     :: Integer
>                      , un     :: Integer
>                      , curT   :: Integer
>                      }
>   deriving Show
>
> -- u_1 = 5;  5^2 = 25 = 2*10 + 5
> initUState = UState 2 5 10
>
> uStep :: State UState Int
> uStep = do
>   u <- gets un
>   p <- gets pn
>   t <- gets curT
>
>   let d   = p `mod` 10      -- next digit
>       u'  = d * t + u       -- prepend the next digit to u
>       p'  = (p + 2*d*u + d*d*t) `div` 10   -- see above proof
>
>   put (UState p' u' (10*t)) -- record the new values
>
>   return $ fromInteger d    -- return the new digit
>
> type Decadic = [Int]
>
> u :: Decadic
> u = 5 : evalState (sequence $ repeat uStep) initUState

To round things out, I’d like to show off some of the cool things we can do with this. First, as we know, it’s possible to do arithmetic with decadic numbers. So let’s implement it!

Addition of decadic numbers is done just like addition of the usual decimal numbers: we add corresponding places (i.e., line up the numbers one under the other and then add in columns).

> plus :: Decadic -> Decadic -> Decadic

First, we have special cases for zero, represented by the empty list of digits: in those cases we just return the other number unchanged.

> plus [] n2 = n2
> plus n1 [] = n1

Next, to add a decadic number whose first digit is x to a decadic number whose first digit is y, we just add x and y and then continue adding recursively.

> plus (x:xs) (y:ys) = (x+y) : plus xs ys

Of course, we’re not done: this doesn’t do any carrying. Instead of modifying our plus function to do carrying, we just write a function normalize which makes sure every place in a decadic number is between $0$ and $9$ ; it will come in handy for more than just addition.

> normalize :: Decadic -> Decadic

The normalize function simply calls a recursive helper function normalize' which keeps track of the current "carry". The starting carry, of course, is zero.

> normalize = normalize' 0

To normalize zero (the empty list) when the current carry is zero, just return the empty list.

>   where normalize' 0 [] = []

With a nonzero carry and the empty list, we simply extend the list with a special zero digit and continue normalizing.

>         normalize' carry [] = normalize' carry [0]

In the general case, we add the current carry to the next digit x, and compute the quotient and remainder when dividing this sum by ten. The remainder is the next digit d, and the quotient is the new carry which gets passed along recursively.

>         normalize' carry (x:xs) = d : normalize' carry' xs
>           where (carry', d) = (carry + x) `divMod` 10

And now for multiplication, which is based on the observation that zero times anything is zero, and in the general case

$(a + 10b)(c + 10d) = ac + 10(a \cdot d + b \cdot (c + 10d))$ .

> mul :: Decadic -> Decadic -> Decadic
> mul [] _ = []
> mul _ [] = []
> mul (x:xs) (y:ys) = x*y : (map (x*) ys + (xs * (y:ys)))

Finally, we declare Decadic to be an instance of the Num class, which allows us to use decadic numbers in the same ways that we can use other numeric types:

> instance Num Decadic where

To add or multiply decadic numbers, use the plus and mul functions and then normalize.

>   n1 + n2 = normalize (plus n1 n2)
>   n1 * n2 = normalize (mul n1 n2)

To negate a decadic number, subtract the last digit from 10 and the rest of the digits from 9.

>   negate [] = []
>   negate (x:xs) = normalize $ (10-x) : negate' xs
>     where negate' []     = repeat 9
>           negate' (x:xs) = (9-x) : negate' xs

Finally, to convert an integer into a decadic number, put the integer into a list of one element and normalize.

>   fromInteger = normalize . (:[]) . fromInteger

So, let’s try it! We’ll want a way to display decadic numbers:

> showDecadic :: Decadic -> IO ()
> showDecadic d = putStrLn . dots $ digits
>   where d'   = take 31 d
>         dots | length d' <= 30 = id
>              | otherwise       = ("..." ++)
>         digits =  concat . reverse . map show . take 30 $ d'

Normal decimal integers can also be used as decadic numbers:

*Decadic2> showDecadic 7
7

Here’s $u$ :

*Decadic2> showDecadic u
...106619977392256259918212890625

And here’s $u^2$ ; it had better be the same as $u$ !

*Decadic2> showDecadic (u^2)
...106619977392256259918212890625

Well, looks like it’s the same for the first 30 digits at least! We can also compute $1-u$ . Remember, if $u^2 = u$ then $(1-u)^2 = 1 - 2u + u^2 = 1 - 2u + u = 1 - u$ , so $1-u$ should be another self-square number. Remember how we thought there might be a self-square number ending in $6$ ? Well, this is it!

*Decadic2> showDecadic (1-u)
...893380022607743740081787109376
*Decadic2> showDecadic ((1-u)^2)
...893380022607743740081787109376

Finally, we can check that $u (1 - u) = u - u^2 = u - u = 0$ :

*Decadic2> showDecadic (u * (1-u))
...000000000000000000000000000000

If you recall, this is in some sense the fundamental reason why the decadic numbers act so funny, because it has zero divisors: pairs of numbers (like $u$ and $1-u$ ), neither of which is zero, whose product is nonetheless zero.

Now, if you remember, from even further back, what got us into this whole decadic mess in the first place:

In that first post, I said

I managed to extend this pattern for a few more digits before I got bored. Does it continue forever or does it eventually stop? Is there any deeper mathematical explanation lurking behind this supposed “curiosity”? What’s so special about $f(x) = 2x^2 - 1$ ? Do patterns like this exist for other functions?

Well, by this point I hope it’s clear that there is indeed a deeper mathematical explanation lurking! The equation

$x = 2x^2 - 1$

admits the solutions $x = 1$ and $x = -1/2$ , but does it admit any other decadic solutions? Notice that given $(x - a)(x - b) = 0$ , which has $x = a$ and $x = b$ as solutions, then $u a + (1-u)b$ (and $(1-u)a + ub$ ) are also solutions:

$(ua + (1-u)b - a)(ua + (1-u)b - b) = ((u-1)a + (1-u)b)(ua + ub) = 0$ .

So in this case we get

$\displaystyle u - (1-u)/2 = \frac{3u - 1}{2}$

as a solution (the other solution is not a decadic integer).

To implement it, we need a way to halve decadic numbers (I’ll let you work out what’s going on here):

> halve :: Decadic -> Decadic
> halve [] = []
> halve t@(s:_)
>   | odd s     = error "foo"
>   | otherwise = halve' t
>   where
>     halve' [] = []
>     halve' [x] = [x `div` 2]
>     halve' (x:x':xs) = (x `div` 2 + adj) : halve' (x':xs)
>       where adj | odd x'    = 5
>                 | otherwise = 0

And now we can define

> q = halve (3*u - 1)

*Decadic2> showDecadic q
...159929966088384389877319335937
*Decadic2> showDecadic (2*q^2 - 1)
...159929966088384389877319335937

Woohoo! This clearly shows that the pattern does, in fact, continue forever. It also shows us that $f(x) = 2x^2 - 1$ is not particularly special: any quadratic function that factors as $(x - a)(x - b)$ , at the very least, will lead to a pattern like this, and probably lots of other equations do too.

If you’re interested in reading more, here’s where I got some of my information. For example, you can read about how there is another number $v = \dots 04103263499879186432$ , defined by starting with $2$ and iteratively raising to the fifth power (just as we defined $u$ by starting with $5$ and successively squaring), such that $v^5 = v$ . It even seems that the author of that page, Gérard Michon, has recently added a discussion of this very problem, prompted by my blog posts! Isn’t the internet great?

Posted in arithmetic, programming | Tagged computing, decadic, numbers | 2 Comments

Differences of powers of consecutive integers

Posted on January 29, 2012 by Brent

Patrick Vennebush of Math Jokes 4 Mathy Folks recently wrote about the following procedure that yields surprising results. Choose some positive integer $n$ . Now, starting with $n+1$ consecutive integers, raise each integer to the $n$ th power. Then take pairwise differences by subtracting the first number from the second, the second from the third, and so on, resulting in a list of only $n$ numbers. Do the same thing again, resulting in $n-1$ numbers, and repeat until you are left with a single number.

For example, suppose we choose $n=4$ , and start with the five consecutive integers $23, 24, 25, 26, 27$ . We raise them all to the fourth power, giving us

$279841,331776,390625,456976,531441$

Now we take pairwise differences: $331776 - 279841 = 51935$ , then $390625 - 331776 = 58849$ , and so on, and we get the new list

$51935,58849,66351,74465$

Repeating the difference procedure gives

$6914,7502,8114$

$588,612$

$24$

OK, so we get $24$ . So what?

Well, if you try enough examples, you may notice a surprising pattern. I’ll let you play with it for a while. Over the course of a few future posts I’ll explain the pattern and prove that it always holds—but the proof will be a really cool combinatorial one, with pretty pictures!

Posted in arithmetic, pattern | Tagged consecutive, difference, integers, powers, surprising | 12 Comments

u-tube

Posted on January 10, 2012 by Brent

[This is the eighth in a series of posts on the decadic numbers (previous posts: A curiosity, An invitation to a funny number system, What does "close to" mean?, The decadic metric, Infinite decadic numbers, More fun with infinite decadic numbers, A self-square number).]

In my previous post, we found a decadic number

$u = \dots 57423423230896109004106619977392256259918212890625$

with the curious property that it is its own square.

We also discovered an algorithm for computing $u$ : starting with $u_1 = 5$ , we square each $u_n$ and take the remainder $\bmod{10^{n+1}}$ (that is, we keep only the last $n+1$ digits) to give us $u_{n+1}$ . Then $u$ is defined as the limit of the $u_n$ as $n \to \infty$ .

(As an aside, notice why we are allowed to define $u$ as a limit: because of the funny decadic metric, the distance between successive $u_n$ is actually getting smaller by a factor of $1/10$ each time; they are converging on something, and that something is $u$ .)

So,

$\begin{array}{rcl} u_2 & = & u_1^2 \bmod{100} = 5^2 \bmod{100} = 25 \\ u_3 & = & u_2^2 \bmod{1000} = 25^2 \bmod{1000} = 625 \\ u_4 & = & u_3^2 \bmod{10000} = 625^2 \bmod{10000} = 390625 \bmod{10000} = 625 \\ u_5 & = & \dots \end{array}$

We can pretty easily turn this into a computer program to compute $u$ to as many digits as we like. As usual, I’m using my favorite programming language, Haskell. It doesn’t matter if you don’t know Haskell, I’ll try to explain things well enough that you can still follow along anyway.

First we just need to import something we’ll need later.

> module Decadic where
>
> import Control.Monad.State

Now, we define u1 to be 5, and the function nextU takes $n$ and $u_n$ as input and produces $u_{n+1} = u_n^2 \bmod{10^{n+1}}$ .

> u1 :: Integer
> u1 = 5
>
> nextU :: Integer -> Integer -> Integer
> nextU n un = (un^2) `mod` (10^(n+1))

Now we can define the (infinite!) list of $u_n$ , which we’ll call us. generateUs takes the current values of $n$ and $u_n$ and generates a list with $u_n$ as its first element, and the remainder of the list computed by a recursive call to generateUs (with $n + 1$ and $u_{n+1}$ as its inputs). us just gets things started by calling generateUs with $n=1$ and $u_1 = 5$ .

> us :: [Integer]
> us = generateUs 1 u1
>   where generateUs n un = un : generateUs (n+1) (nextU n un)

Let’s see if it works by asking for the first ten elements of the list us:

*Decadic> take 10 us
[5,25,625,625,90625,890625,2890625,12890625,212890625,8212890625]

Seems good!

However, all is not roses. Let’s think about what our program had to do in order to calculate $u_{10}$ . Given $u_9 = 212890625$ , it first squared it to get $45322418212890625$ … but wait a minute! $45322418212890625$ ends in $212890625$ —of course it does, since that’s the defining property of the $u_n$ . But that means that half (more or less) of the work needed to square $212890625$ was completely wasted! There’s no point in computing the last half of $u_n^2$ , since we already know it is going to be $u_n$ . We only care about the rest of it.

Realizing that our program is doing too much work is one thing; turning this intuition into actual improvements is quite another! Initially it was far from obvious to me how to avoid the repeated work. But I finally figured it out.

The key is to remember at each step not just $u_n$ , but $u_n^2$ . However, since we know that $u_n^2$ ends with $u_n$ , we can remember $u_n$ along with the prefix (call it $p_n$ ) of $u_n^2$ that comes before the ending $u_n$ . For example, for $n = 9$ we’ll remember $p_9 = 45322418$ and $u_9 = 212890625$ (since $u_9^2 = 45322418212890625$ ). That is, we will always insist that

$u_n^2 = 10^n p_n + u_n$ .

The key question now is: given $p_n$ and $u_n$ , how do we compute $p_{n+1}$ and $u_{n+1}$ ? Well, computing $u_{n+1}$ is easy: it’s just $u_n$ with an extra digit tacked on, namely, the final digit of $p_n$ . In particular, we can define $d_{n+1} = p_n \bmod{10}$ to be the last digit of $p_n$ ; then $u_{n+1} = d_{n+1} 10^n + u_n$ .

It’s computing $p_{n+1}$ which is the tricky part. Since we want to have $u_{n+1}^2 = p_{n+1} 10^{n+1} + u_{n+1}$ , the idea is to expand out $u_{n+1}^2$ , manipulate it into the desired form, and see what we get for $p_{n+1}$ :

$\begin{array}{rcl} u_{n+1}^2 & = & (d_{n+1} 10^n + u_n)^2 \\ & = & d_{n+1}^2 10^{2n} + 2d_{n+1} 10^n u_n + u_n^2 \\ & = & d_{n+1}^2 10^{2n} + 2d_{n+1} 10^n u_n + p_n 10^n + u_n \end{array}$

Now, note that we can break $p_n$ into two parts: the last digit (which we’re calling $d_{n+1}$ ) and the rest. That is,

$p_n = 10 \lfloor p_n/10 \rfloor + d_{n+1}$

where $\lfloor - \rfloor$ denotes rounding down to the nearest integer, so $\lfloor p_n/10 \rfloor$ means to divide $p_n$ by 10 and throw away the remainder—that is, $p_n$ without its final digit. Continuing,

$\begin{array}{rcl} u_{n+1}^2 & = & d_{n+1}^2 10^{2n} + 2d_{n+1} 10^n u_n + p_n 10^n + u_n \\ & = & d_{n+1}^2 10^{2n} + 2d_{n+1} 10^n u_n + \lfloor p_n / 10 \rfloor 10^{n+1} + d_{n+1} 10^n + u_n \\ & = & d_{n+1}^2 10^{2n} + 2d_{n+1} 10^n u_n + \lfloor p_n / 10 \rfloor 10^{n+1} + u_{n+1} \\ & = & (d_{n+1}^2 10^{n-1} + 2 d_{n+1} u_n / 10 + \lfloor p_n / 10 \rfloor) 10^{n+1} + u_{n+1} \end{array}$

and there we have it! We have expressed $u_{n+1}^2$ in the form $p_{n+1} 10^{n+1} + u_{n+1}$ . In particular,

$p_{n+1} = d_{n+1}^2 10^{n-1} + 2 d_{n+1} u_n / 10 + \lfloor p_n / 10 \rfloor$

Notice that $2 d_{n+1} u_n / 10$ is actually guaranteed to be an integer, since $u_n$ is divisible by 5—this is the same phenomenon we ran across in proving that our simple algorithm for computing $u_n$ works correctly. So we can factor out a division by ten (division is generally slow so we’d rather only do it once):

$p_{n+1} = \lfloor (d_{n+1}^2 10^n + 2 d_{n+1} u_n + p_n) / 10 \rfloor$

And now for some code! As one final optimization, instead of keeping track of $n$ , we keep track of $10^n$ , so we don’t have to keep recomputing $10^n$ every iteration.

First, the data structure we use to keep track of the current values of $p_n$ , $u_n$ , and $10^n$ :

> -- State for incrementally constructing u_n.
> -- Invariant: curT = 10^n; un^2 = pn*curT + un
> data UState = UState { pn     :: Integer
>                      , un     :: Integer
>                      , curT   :: Integer
>                      }
>   deriving Show

The initial state for $n = 1$ has $p_1 = 2$ , $u_1 = 5$ , and $10^1 = 10$ .

> -- u_1 = 5;  5^2 = 25 = 2*10 + 5
> initUState = UState 2 5 10

The uStep function is the heart of the algorithm: it takes the current values of $p_n$ , $u_n$ , and $10^n$ and updates them to $p_{n+1}$ , $u_{n+1}$ , and $10^{n+1}$ . It also returns $d_{n+1}$ as its result, which we can use to build up $u$ digit-by-digit.

> uStep :: State UState Int
> uStep = do
>   u <- gets un
>   p <- gets pn
>   t <- gets curT
>
>   let d   = p `mod` 10      -- next digit
>       u'  = d * t + u       -- prepend the next digit to u
>       p'  = (p + 2*d*u + d*d*t) `div` 10   -- see above proof
>
>   put (UState p' u' (10*t)) -- record the new values
>
>   return $ fromInteger d    -- return the new digit

So, did we gain anything? As a concrete comparison, let’s see how long it takes to compute $u_{10001}$ . Using our first, simple code, it takes 7.2 seconds:

*Decadic> :set +s
*Decadic> length . show $ us !! 10000
10001
(7.23 secs, 208746872 bytes)

(I just had it print the number of digits of $u_{10001}$ to avoid wasting a ton of space printing out the entire number.) And using our new and hopefully improved code?

*Decadic> length . show . un . flip execState initUState
          $ replicateM_ 10000 uStep
10001
(1.55 secs, 225857080 bytes)

Only 1.5 seconds! Nice!

The other nice thing about uStep is that it spits out one digit of $u$ at a time, which we can use to define $u$ as an (infinite) list of digits—as if the digits were coming one at a time down a "tube", a $u$ -tube, one might say… get it, a $u$ … tube… heh… never mind.

> type TenAdic = [Int]
>
> u :: TenAdic
> u = 5 : evalState (sequence $ repeat uStep) initUState

*Decadic> reverse $ take 20 u
[9,2,2,5,6,2,5,9,9,1,8,2,1,2,8,9,0,6,2,5]

Nifty! Next time the real fun begins, as I show off some cool things we can do with our shiny new implementation of $u$ .

Posted in computation, convergence, infinity, iteration, modular arithmetic, number theory, programming | Tagged decadic, Haskell, idempotent, streaming, u | 1 Comment

The Math Less Traveled

Combinatorial proofs

Differences of powers of consecutive integers, part II

17×17 4-coloring with no monochromatic rectangles

Book review: Nine Algorithms that Changed the Future

Computing with decadic numbers

Differences of powers of consecutive integers

u-tube

Pages

Tags

Blogroll

Fun

Reference

Categories

Archives

Meta

Follow “The Math Less Traveled”