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


\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

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
(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
(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

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

This entry was posted in computation, convergence, infinity, iteration, modular arithmetic, number theory, programming and tagged , , , , . Bookmark the permalink.

2 Responses to u-tube

  1. Pingback: Computing with decadic numbers | The Math Less Traveled

  2. Ray says:

    Superb analysis on how to work with decadic numbers, from a computing standpoint. I just passed the link to my colleagues. Nice work.

Comments are closed.