u-tube

[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$ .

u-tube

About Brent

2 Responses to u-tube

Follow Blog via Email

Meta

Categories

Archives

u-tube

Share this:

Like this:

Related

About Brent

2 Responses to u-tube

Follow Blog via Email

Meta

Categories

Archives