By Guy Jacobson and Kiran Kedlaya

This is a numerical steganography puzzle. Solvers need to have access to a tool that can manipulate very large integers; I will use Python here, since it has baked-in support for arbritrary precision arithmetic. The fraction shown is just what it looks like: a rational number. Solvers need to convert the number into its decimal expansion and then go on a brief treasure hunt through the digits found there, delving deep. The final step will require using a bit of mathematical finesse and not mere brute force.

Looking at the initial segment of the expansion just beyond the decimal point

.1911091620152008 . . .

solvers should notice that the digits taken in pairs are the alphabetic indices (A=01, B=02, . . .) of the beginning of a message:

 19 11 09 16 20 15 20 08 . . . S K I P T O T H . . .

Looking at more digit pairs, solvers can recover a complete message, terminated by 00:

 19 11 09 16 20 15 20 08 15 21 19 01 14 04 20 08 04 09 01 09 20 00 S K I P T O T H O U S A N D T H D I G I T

So it looks like we will need to venture considerably deeper into the decimal expansion. It is not difficult to compute a few thousand digits of this fraction; in Python, for example, you can just write

`10**2000 * a // b`

to get the first 2000 digits of the fraction a/b, expressed as a large integer.

The digits starting at position 1000 in the decimal expansion encode another alphabetic-index message, again terminated by 00:

 20 08 05 14 20 15 07 15 15 07 15 12 20 08 04 09 01 09 20 00 T H E N T O G O O G O L T H D I G I T

At this point, solvers need to find a way of computing specific digits of a/b without generating the entire decimal expansion starting from the decimal point, since there is no practical way of computing or storing a googol (10100) digits. You might think, “Hey, since this is a rational number, it’s going to have a repeating decimal, so we can just compute the repeating part and use that to find the googolth digit.” But the fraction was intentionally constructed to have a very large repeating sequence in its expansion (with a length of 16,859,147,591,006,830,465,663,878 digits, if you must know), so . . . good luck with that approach!

There are probably several ways to find digits deep inside the expansion of a rational number quickly, but here’s a fairly straightforward method that relies on a few modular arithmetic tricks. We want to calculate the kth digit of a/b. We noted earlier that we can shift the decimal point by simply multiplying by 10k to get 10k·a/b. This will be a number whose decimal value will have the desired digit in the ones place (just to the left of the decimal point). But it will also have a fractional part to the right of the decimal point. We can get rid of this fractional part by first subtracting off 10k·a mod b from 10k·a before dividing by b. Then

[10k·a - (10k·a mod b)] / b

is an integer whose least-significant digit is the kth digit of a/b. If we take this value mod 10, we will get the digit we seek:

([10k·a - (10k·a mod b)] / b) mod 10

Assuming that the original denominator b is relatively prime to 10, the multiplicative inverse of b mod 10 exists; let’s call this value b-1. (For our fraction, where the denominator ends in 1, b-1 = 1.) We can replace the division in our formula with a multiplication by b-1 to get:

[b-1·10k·a - b-1·(10k·a mod b)] mod 10

Since the first term inside the brackets has a factor of 10k, it will vanish when we mod by 10 at the end (assuming k > 0), leaving just

[-b-1·(10k·a mod b)] mod 10

Now the only tricky bit remaining is to compute 10k mod b efficiently for very large values of k. We cannot do the obvious thing and multiply 10 by itself a googol times, or we would be still be computing at the heat death of the universe. Fortunately, there is the well-known trick of exponentiation by squaring that allows us to compute powers quickly. In fact, this technique is used by Python’s three-argument `pow()` function which does exponentiation modulo a given base, so we don’t even need to code up the trick ourselves if we program in Python. (See the code at the end of this solution for an example implementation.)

Using this method, we can find the sequence of digits in a/b starting at position 10100 and recover the final message as we did previously:

 12 09 14 04 05 14 12 01 02 00 L I N D E N L A B

yielding the solution to this puzzle, LINDEN LAB.

Here is a bit of code in Python to solve the problem, using the digit-extraction technique described above:

```#!/usr/bin/python

def gcd(b, n):
x0, x1, y0, y1 = 1, 0, 0, 1
while n != 0:
q, b, n = b // n, n, b % n
x0, x1 = x1, x0 - q * x1
y0, y1 = y1, y0 - q * y1
return b, x0, y0

def mulinv(b, n):
g, x, _ = gcd(b, n)
if g == 1:
return x % n

# generate digits of num/denom starting at d in a given base
def digit (num, denom, d = 1, base = 10): # assumes denom is relatively prime to base
denominv = mulinv (denom, base)
n = (num * pow (base, d - 1, denom)) % denom
while True:
n = (n * base) % denom
yield (-n * denominv) % base

# decode ASCII from fraction
def decodefrac(num, denom, start = 1):
d = digit(num, denom, start)
s = ''
while True:
c = 10 * d.next() + d.next();
if (c == 0):    # terminated by 00
return s
s += chr(c + ord('A') - 1)

a = 155619367777177344298579121640577930358466048997080741489502115154736052423438639656...
b = 814295694336200301695257325143637445297122644140195688542610298982742658367669506161...

print decodefrac (a, b, 1)
print decodefrac (a, b, 1000)
print decodefrac (a, b, 10**100)```

## Construction Explanation

This is an explanation of how this puzzle was constructed. The knowledge used here is not needed in order to solve the puzzle.

In order to realize the intended solution path, we needed to find a rational number r/s with the following substrings in its decimal expansion.

• digits 1 to 44: 19110916201520081521190114042008040901092000
• digits 1000 to 1039: 2008051420150715150715122008040901092000
• digits 10100 to 10100+19: 12091404051412010200

Furthermore, we wanted r and s not to be too large.

To simplify the process, we decided to construct r/s to have the form a/b + c/d where the summands were subject to the following conditions.

• The number a/b has a long period and has a substring matching the first step of the path.
• The number c/d has a short period, has zeroes at the location of the first step of the path, and has substrings correctly complementing with a/b to give the second and third steps (plus a trailing zero in each case to prevent an unwanted carry into the final digit).

The second condition on a/b states that this fraction is at least

19110916201520081521190114042008040901092000 / 1044

with the difference being no more than 1/1044. We got some initial candidates from the convergents of the continued fraction of this quantity, with a and b each having 22 digits or slightly more, but these did not have as large a period as we wanted. (Reminder: the period is the multiplicative order of 10 modulo b, that is, the smallest positive integer n such that b divides 10n-1). However, standard heuristics (such as the prime number theorem and Artin’s conjecture on primitive roots) suggest that trying replacing the original fraction with nearby values would, in a reasonable time, yield a convergent with a suitable denominator.

To construct c/d, we took d to be a factor of 10p-1 for some relatively small prime p, so that (as long as d > 9) the period of c/d would equal p. In order to execute the construction, we needed to ensure that the three prescribed substrings correspond to distinct locations in the period; that is, the quantities

1 to 44, 1000 to 1040, 10100 to 10100+20

had to be pairwise distinct modulo p. This of course forces p >= 44+41+21 = 106. We took p=131, for which the listed quantities correspond to

1 to 44, 83 to 123, 52 to 72

For convenience, we also imposed zeroes in positions 124 to 131. Now the decimal expansion of (108·c)/d, ignoring the digits before the decimal point, has zeroes in position 1 to 52 and prescribed digits in positions 60 to 80.

If we had chosen the remaining digits randomly, we are likely to have a denominator of 10131-1. However, since we had flexibility to choose them, we looked for a choice that would make the numerator have a common factor with 10131-1. The prime factorization of 10131-1 can be computed using any standard mathematical software package (we used Magma and Sage, but Mathematica should also work). Note that numbers of the form mn - 1 are easier to factor than random numbers using current techniques, but this is a long story outside the scope of this explanation. In any case, the factorization is

32 · 80173 · 109517 · 141811693 · 446790173 · 7370364319027 · 15594845538029429933 · 7317723970031057677693 · 131758351065116151205213 · 180222062287834025451247081

We opted to search for a multiple of 80173·109517. This amounted to solving a linear Diophantine equation of the form

Ax + By = C (mod 80173·109517)

for some A,B,C, where x is an integer from 0 to 107-1 (representing the digits from 45 to 51) and y is an integer from 0 to 109-1 (representing the digits from 73 to 82). Given the small sizes of the numbers involved, we did this by exhausting over the choices for x, and for each of these solving for y (using standard modular arithmetic) to see if it fell into the correct range.

In this way, we arrived at the values

• a = 122974716351446828004968
• b = 643478915687283605559691
• c = 130232918073045339638509705978224505361995762494579108317537000000000
• d = 1265458237223626203402473149639710975904321641786205068866013979603779880318998437005817381711485428451586671804079361871

yielding the values of r and s in the statement of the puzzle:

• r = 155619367777177344298579121640577930358466048997080741489502115154736052423438639656322566701951512888971130859442675019505459646950828457775128
• s = 814295694336200301695257325143637445297122644140195688542610298982742658367669506161728178242191616818890832959609127053599963426497407079941861