Problem 104: Pandigital Fibonacci Ends

View on Project Euler

Project Euler Problem 104 Solution

EulerSolve provides an optimized solution for Project Euler Problem 104, Pandigital Fibonacci Ends, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for the first Fibonacci index \(k\) such that the first nine digits of \(F_k\) and the last nine digits of \(F_k\) are both 1-9 pandigital. In other words, each of those 9-digit blocks must contain every digit from 1 through 9 exactly once, with no zero and no repetition. The obstacle is scale. By the time the answer appears, the Fibonacci number itself has tens of thousands of digits, so constructing it exactly at every step would be wasteful. The implementations therefore split the question into two smaller numerical views of the same sequence: the tail is handled exactly with modular arithmetic, while the head is reconstructed from logarithms. Mathematical Approach Let the Fibonacci sequence be defined by \(F_1=1\), \(F_2=1\), and \(F_n=F_{n-1}+F_{n-2}\) for \(n \ge 3\). The search condition is naturally the conjunction of two tests on the same index \(n\): one test on \(F_n \bmod 10^9\), and one test on the leading decimal digits of \(F_n\). Splitting the Condition into Tail and Head If the last nine digits are all that matter, then the entire problem lives modulo \(10^9\). If the first nine digits are all that matter, then only the decimal order of magnitude and mantissa matter. Those two facts let us avoid big integers in the main loop....

Detailed mathematical approach

Problem Summary

The Fibonacci numbers are defined by \(F_1=F_2=1\) and \(F_n=F_{n-1}+F_{n-2}\). A block of nine digits is called \(1\)–\(9\) pandigital if it contains each of the digits \(1,2,\dots,9\) exactly once, in any order; such a block has no \(0\). The statement tells us that \(F_{541}\), with \(113\) digits, is the first Fibonacci number whose last nine digits are pandigital, and \(F_{2749}\), with \(575\) digits, the first whose first nine digits are. We must find the index \(k\) of the first Fibonacci number \(F_k\) whose first nine and last nine digits are both pandigital. [1]

The whole solution in one sentence: the last nine digits of a Fibonacci number depend only on the last nine digits of the two numbers before it, so they can be carried along exactly with small numbers, while the first nine digits can be read off with a logarithm from Binet's formula \(F_n\approx\varphi^n/\sqrt5\), where \(\varphi\approx1.618\) is the golden ratio; we walk through \(n=3,4,5,\dots\) and compute the first digits only at the rare indices whose last nine digits are already pandigital.

Mathematical Approach

1. Small Fibonacci numbers and their two ends

The first twenty Fibonacci numbers are \(1, 1, 2, 3, 5, 8, 13, 21, 34, 55,\)
\(89, 144, 233, 377, 610,\)
\(987, 1597, 2584, 4181, 6765.\)
From \(F_6\) on, each number is about \(1.6\) times the one before it, and since \(1.6^4\approx6.6\) is below \(10\) while \(1.6^5\approx10.5\) is above it, a new digit appears every four or five steps. To prove this, apply the recurrence four and five times: \(F_{n+4}=5F_n+3F_{n-1}\) and \(F_{n+5}=8F_n+5F_{n-1}\). For \(n\ge3\), \(F_n=F_{n-1}+F_{n-2}\) with \(1\le F_{n-2}\le F_{n-1}\) gives \(F_n/2\le F_{n-1}\lt F_n\), hence \(F_{n+4}\lt8F_n\) and \(F_{n+5}\ge10.5\,F_n\). If \(F_m\) is the first number with \(d\ge2\) digits, then \(F_{m-1}\lt10^{d-1}\), so \(F_{m+3}\lt8F_{m-1}\) is still below \(10^d\), while \(F_{m+5}\gt10F_m\ge10^d\) has more digits: for every \(d\gt1\) there are four or five Fibonacci numbers with \(d\) digits. [2]

We call the first nine digits of \(F_n\) its head and the last nine digits its tail; for example, \(F_{50}=12586269025\) has the head \(125862690\) and the tail \(586269025\). Nine digits appear for the first time in \(F_{40}=102334155\), whose head and tail are the same block; up to \(F_{83}\), which has \(17\) digits, the two blocks overlap, and from \(F_{84}=160500643816367088\) on they are separate. Since the first pandigital tail belongs to \(F_{541}\), every index that can solve the problem lies far beyond the overlap.

2. Pandigital blocks

A pandigital block is an arrangement of the digits \(1\) to \(9\); read as numbers, such blocks are also called zeroless pandigital numbers. [3] The first place can hold any of the nine digits, the second any of the eight that are left, and so on, so there are \(9\cdot8\cdots1=362880\) such blocks, a number called nine factorial and written \(9!\). A tail can be any of the \(10^9\) blocks from \(000000000\) to \(999999999\), so if tails were random, about one in \(10^9/9!\approx2756\) would be pandigital. The tail \(839725641\) of \(F_{541}\) is pandigital; the tail \(002250249\) of \(F_{2749}\) is not.

3. The tail: arithmetic modulo \(10^9\)

The last two digits of a sum depend only on the last two digits of the numbers that are added: \(F_{11}=89\) and \(F_{12}=144\) end in \(89\) and \(44\), and \(89+44=133\) ends in \(33\), just like \(F_{13}=233\). In general, \(a\equiv b\pmod m\) means that \(a-b\) is a multiple of \(m\), and congruences can be added: from \(a\equiv a'\) and \(b\equiv b'\) it follows that \(a+b\equiv a'+b'\pmod m\), because \((a+b)-(a'+b')\) is a sum of two multiples of \(m\). [4] The remainder after division by \(10^9\) is the last nine digits, read as a number, so if \(t_n\) is the remainder of \(F_n\), then

$$t_n=(t_{n-1}+t_{n-2})\bmod 10^9.$$

A program keeps only the two latest tails; their sum stays below \(2\cdot10^9\), so it fits into an unsigned 32-bit integer, whose largest value is \(2^{32}-1=4294967295\). By induction on \(n\), each \(t_n\) is exactly the tail of \(F_n\), even when \(F_n\) has tens of thousands of digits.

The tails repeat. Only \(10^{18}\) different pairs of tails exist, so among the first \(10^{18}+1\) pairs \((t_n,t_{n+1})\) two are equal, by the pigeonhole principle: when more objects than boxes are put into boxes, some box receives two. [5] Take a repeat \((t_i,t_{i+1})=(t_j,t_{j+1})\) with \(i\lt j\) and \(i\) as small as possible. Every pair determines the pair before it, since \(t_{n-2}=(t_n-t_{n-1})\bmod 10^9\), so for \(i\gt1\) we would also have \((t_{i-1},t_i)=(t_{j-1},t_j)\), an earlier repeat. Hence \(i=1\): the tails return to their start \(1,1\) and repeat. The length of one round is the Pisano period \(\pi(10^9)\). [6]

Small moduli show how such periods combine. Modulo \(2\) the Fibonacci numbers are \(1,1,0,1,1,0,\dots\), so \(\pi(2)=3\). Modulo \(5\) they are \(1,1,2,3,0,3,3,1,4,0,\)
\(4,4,3,2,0,2,2,4,1,0\) and then \(1,1\) again, so \(\pi(5)=20\). Since \(10\) divides a number exactly when \(2\) and \(5\) both do, the last digits are back at \(1,1\) exactly when the remainders modulo \(2\) and modulo \(5\) are, that is, after a common multiple of \(3\) and \(20\) steps; the least common multiple gives \(\pi(10)=\mathrm{lcm}(3,20)=60\).

These are the cases \(e=1\) of the known periods \(3\cdot2^{e-1}\) of \(2^e\) and \(4\cdot5^e\) of \(5^e\), and the same argument combines the coprime factors \(2^9\) and \(5^9\) of \(10^9\): [6]

$$\begin{aligned}\pi(10^9)&=\mathrm{lcm}\bigl(3\cdot2^8,\,4\cdot5^9\bigr)\\&=1.5\cdot10^9.\end{aligned}$$

Running the recurrence through a whole period confirms this value. Within one period exactly \(567168\) indices have a pandigital tail, about \(4\%\) more than the \(1.5\cdot362880=544320\) that random tails would give.

4. The direct method and its cost

The simplest program computes every \(F_n\) exactly with big integers, takes its tail as the remainder modulo \(10^9\) and converts the whole number to decimal for its head only when the tail is pandigital. \(F_n\) has about \(0.209\,n\) digits, as Section 6 explains, an addition or a remainder costs time proportional to the number of digits, and the digit counts of \(F_3,\dots,F_N\) add up to about \(0.1045\,N^2\): roughly \(10^9\) digit operations for \(N=10^5\) and a hundred times more for \(N=10^6\), besides the rare conversions. In Python such a direct scan took \(0.6\) seconds for \(N=10^5\) and \(93\) seconds for \(N=10^6\); the method of the next sections needs \(0.03\) and \(0.3\) seconds, because its work grows like \(N\) instead of \(N^2\).

5. Binet's formula

The two numbers

$$\varphi=\frac{1+\sqrt5}{2}\approx1.6180339887,$$

$$\psi=\frac{1-\sqrt5}{2}\approx-0.6180339887$$

are the two solutions of \(x^2=x+1\): \(\varphi\) is the golden ratio, an irrational number, and \(\psi\) is its conjugate. [7] Because \(\varphi^2=\varphi+1\), every power of \(\varphi\) can be reduced to a whole multiple of \(\varphi\) plus a whole number: \(\varphi^3=\varphi^2+\varphi\)\({}=2\varphi+1\), \(\varphi^4=3\varphi+2\), and in general, with \(F_0=0\),

$$\varphi^n=F_n\,\varphi+F_{n-1}\qquad(n\ge1).$$

Induction proves it. For \(n=1\) it reads \(\varphi=1\cdot\varphi+0\), and if it holds for \(n\), then

$$\begin{aligned}\varphi^{n+1}&=F_n\varphi^2+F_{n-1}\varphi\\&=F_n(\varphi+1)+F_{n-1}\varphi\\&=F_{n+1}\varphi+F_n.\end{aligned}$$

The same steps work for \(\psi\), since \(\psi^2=\psi+1\) as well, so \(\psi^n=F_n\psi+F_{n-1}\). [2] Subtracting the two identities removes \(F_{n-1}\), and since \(\varphi-\psi=\sqrt5\) we obtain Binet's formula

$$F_n=\frac{\varphi^n-\psi^n}{\sqrt5}.$$

For example, \(\varphi^{10}=55\varphi+34\approx122.99\), so \(\varphi^{10}/\sqrt5\approx55.0036\), and the small term \(\psi^{10}/\sqrt5\approx0.0036\) is exactly what turns it into \(F_{10}=55\). Because \(|\psi|\lt1\), this correction shrinks quickly: it is about \(-0.106\) for \(n=3\), \(0.0036\) for \(n=10\) and \(1.6\cdot10^{-11}\) for \(n=50\). Its size is below \(\tfrac12\) for every \(n\ge0\), so \(F_n\) is the whole number nearest to \(\varphi^n/\sqrt5\). [2]

6. The head from a logarithm

Logarithms separate the size of a number from its digits. Take \(F_{20}=6765=10^3\cdot6.765\); its common logarithm, the logarithm to base \(10\), is \(\log_{10}6765\approx3.83027\). The integer part \(3\), called the characteristic, is the exponent of the power of \(10\), so the number has \(3+1=4\) digits. The fractional part \(0.83027\), called the mantissa, is the logarithm of \(6.765\) and carries the digits: \(10^{0.83027}\approx6.765\), and \(\lfloor10^{0.83027+2}\rfloor=676\) gives the first three digits, where \(\lfloor y\rfloor\) is the largest whole number not above \(y\). Numbers that differ by a factor of a power of \(10\) share their mantissa. [8]

A large \(F_n\) is never computed in full; instead, Binet's formula gives its logarithm almost exactly. Write \(A_n=\varphi^n/\sqrt5\) for the main term of the formula and \(x=\log_{10}A_n\). Since logarithms turn powers into multiples and quotients into differences,

$$x=n\log_{10}\varphi-\log_{10}\sqrt5,$$

with \(\log_{10}\varphi\approx0.2089876402\) and \(\log_{10}\sqrt5\approx0.3494850022\); for \(n=20\) this gives \(x\approx3.83027\) again. Split \(x=j+f\) into its integer part \(j\) and its fractional part \(f\), with \(0\le f\lt1\). Then \(A_n=10^j\cdot10^f\) with \(1\le10^f\lt10\), so, as for \(6765\), \(A_n\) has \(j+1\) digits before the decimal point and the same digits as \(10^f\). Multiplying \(10^f\) by \(10^8\) brings eight more digits in front of the decimal point, so the first nine digits of \(A_n\) are

$$\bigl\lfloor10^{f+8}\bigr\rfloor.$$

For \(n=50\) we get \(x\approx10.0998970103\), so \(F_{50}\) has \(11\) digits, and \(f\approx0.0998970103\) gives \(10^{f+8}\approx125862690.25\) and the head \(125862690\), in agreement with \(F_{50}=12586269025\). For \(n=2749\) we get \(x\approx574.157538045\), so \(575\) digits, and \(10^{f+8}\approx143726895.53\) gives the head \(143726895\), which is pandigital, as the statement says. Since \(x\) grows by \(\log_{10}\varphi\approx0.209\) per step, the number of digits grows by \(0.209\) per step on average: by one at some steps and not at all at the others.

One complete step of the search looks like this. At \(n=541\) the tail \(839725641\) is pandigital, so the head is computed: \(x\approx112.71\) shows that \(F_{541}\) has \(113\) digits, and \(10^{f+8}\approx516212329.27\) gives the head \(516212329\). It repeats the digits \(1\) and \(2\), so the search moves on to \(n=542\).

7. Why dropping \(\psi^n\) is safe in the search

The formula gives the head of \(A_n\), not of \(F_n\), and these can differ. Let \(F_n\) have \(D\ge10\) digits and split it as \(F_n=H\cdot10^{D-9}+r\), where \(H\) is the true head and \(0\le r\lt10^{D-9}\) is the number formed by the other \(D-9\) digits. Then \(A_n=F_n+\varepsilon\) with \(\varepsilon=\psi^n/\sqrt5\), and \(|\varepsilon|\lt0.11\) for \(n\ge3\). Dividing by \(10^{D-9}\) gives

$$\frac{A_n}{10^{D-9}}=H+\frac{r+\varepsilon}{10^{D-9}}.$$

If \(r\ge1\), the numerator \(r+\varepsilon\) lies strictly between \(0\) and \(10^{D-9}\), so the last fraction lies strictly between \(0\) and \(1\) and the floor is exactly \(H\). The same bounds show that \(A_n\) lies strictly between \(10^{D-1}\) and \(10^D\), so the integer part \(j\) of Section 6 is \(D-1\) and \(10^{f+8}\) is exactly the left-hand side above. Hence the formula gives the true head whenever the digits after the head are not all zero.

If they are all zero and \(n\) is odd, then \(\varepsilon\lt0\) and the formula gives \(H-1\). This really happens: \(F_{45}=1134903170\) has ten digits and ends in \(0\), and \(10^{f+8}\approx113490316.99999999998\), so the formula claims the head \(113490316\) instead of \(113490317\). In the search this case cannot occur, because the head is computed only when the tail is pandigital, and a pandigital tail ends in one of the digits \(1\) to \(9\), so \(r\ge1\).

8. How many correct digits a double carries

Python and Java evaluate \(x\) with doubles, and so does the C++ program wherever long double is just a double, as on Apple silicon; most x86 compilers make that type the 80-bit extended format instead, and under Linux on 64-bit ARM it has 128 bits. [9] A double has a \(53\)-bit significand, good for \(15\) to \(17\) significant decimal digits; between \(2^e\) and \(2^{e+1}\) neighbouring doubles are \(2^{e-52}\) apart. [10] These digits must also hold the integer part of \(x\approx0.209\,n\), so the rounding error of \(f\) grows in proportion to \(n\), and an error \(\delta\) in \(f\) changes \(10^{f+8}\) by up to about \(10^9\cdot\ln10\cdot\delta\). In Python, compared with exact values, the computed \(10^{f+8}\) is never off by more than \(0.0065\) for \(n\le10^5\) and \(0.037\) for \(n\le10^6\).

So the ninth digit is in danger whenever \(10^{f+8}\) lies within a few hundredths of a whole number. At \(n=179054\), for example, the tail \(142563897\) is pandigital and \(F_{179054}\), with \(37420\) digits, begins with \(528995913994\), but Python computes \(10^{f+8}\approx528995914.00007\), giving the head \(528995914\). Which heads come out wrong even depends on the last binary digit of the constants and on the order of the roundings: on an Apple-silicon Mac, Python, Java and a GCC build of the C++ program get the head wrong at two pandigital tails before the search stops, this one included, and a Clang build at neither. Neither error changes a decision, since neither the true nor the computed head is pandigital there. A safe program would rely not on observed errors like these but on a proven bound for the rounding error, and would recompute the head with more digits, or exactly, whenever \(10^{f+8}\) lies within that bound of a whole number; none of the three programs does this.

9. Why there are infinitely many solutions

The C++ and Java programs have no upper limit on \(n\), and an exact search of this kind stops only if a solution exists. We now prove that there are infinitely many. That the programs, which compute heads in floating point, really stop at the first one is a separate fact, checked under Complexity and Verification: up to it no rounding error of Section 8 changes a decision, and it lies far below \(2^{31}-1\), the largest value of their int counters. Let \(P=1.5\cdot10^9\) be the period of Section 3. All indices \(n=541+iP\) with \(i=0,1,2,\dots\) have the tail \(839725641\) of \(F_{541}\), so by Section 7 their heads are \(\lfloor10^{f+8}\rfloor\), where \(x\) grows by \(\alpha=P\log_{10}\varphi\) from one \(i\) to the next. The number \(\log_{10}\varphi\) is irrational: if it were a fraction \(a/b\) with \(b\ge1\), then \(\varphi^b=10^a\) would be rational, but \(\varphi^b=F_b\varphi+F_{b-1}\) with \(F_b\ge1\) is irrational because \(\varphi\) is. So \(\alpha\) is irrational too, and by the equidistribution theorem of Weyl, Sierpiński and Bohl the fractional parts of \(\alpha,2\alpha,3\alpha,\dots\) spread evenly over the interval from \(0\) to \(1\): in the long run every subinterval receives a share of the points equal to its length. [11] Shifting all points by the same amount does not change this. The head is the pandigital block \(123456789\) exactly when \(f\) is at least \(\log_{10}1.23456789\) and less than \(\log_{10}1.2345679\), an interval of positive length, so infinitely many \(i\) give an index whose two ends are both pandigital.

10. How rare the solutions are

The same argument measures how often both ends are pandigital. For \(n\ge40\) let \(g\) be the fractional part of \(\log_{10}F_n\) itself. A block \(h\) is the head of \(F_n\) exactly when \(h\le10^{g+8}\lt h+1\), that is, when \(g\) lies in an interval of length \(\log_{10}(1+1/h)\). The values of \(f\) spread evenly, by the theorem of Section 9, and \(\log_{10}F_n\) differs from \(x\) only by \(\log_{10}(1+\varepsilon/F_n)\), which tends to \(0\). So one of \(f\) and \(g\) can lie in the interval while the other does not only when \(f\) is within a shrinking distance of one of its ends, and because \(f\) spreads evenly, such indices have the long-run share \(0\). Hence in the long run the length \(\log_{10}(1+1/h)\) is the share of indices with head \(h\); this is Benford's law for leading digits, which the Fibonacci numbers are known to obey. [12] Adding these lengths over the \(9!\) pandigital blocks gives the share \(\beta\) of pandigital heads,

$$\beta=\sum_h\log_{10}\Bigl(1+\frac1h\Bigr),$$

which is about \(3.873\cdot10^{-4}\), one index in \(2582\); pandigital tails have the share \(567168/P\), one in \(2645\). Applying the argument of Section 9 separately to each class of indices modulo \(P\) shows that both ends are pandigital for the share \(\beta\cdot567168/P\) of all indices, about one in \(6.8\) million. Small ranges come close to these shares: for \(n\le10^6\) there are \(387\) pandigital heads, just as many as expected, and \(345\) pandigital tails against \(378\) expected. A count with Python's doubles finds only \(385\), because the rounding error of Section 8 hides the heads of \(F_{613888}\) and \(F_{868622}\).

Smaller versions of the question show how irregular the first solution is. The first index whose first \(\ell\) and last \(\ell\) digits are both arrangements of \(1,\dots,\ell\) is \(1\), \(8\), \(9842\), \(18479\), \(503214\), \(294699\), \(2988291\) and \(4736109\) for \(\ell=1,\dots,8\). For \(\ell=3\), \(F_{9842}\) has \(2057\) digits and begins and ends with \(321\); and the index for \(\ell=6\) is smaller than the one for \(\ell=5\).

How the Code Works

The C++ program [13] starts in main. parse_arguments accepts only the option --skip-checkpoints and rejects anything else with exit code \(1\). Unless the checkpoints are skipped, run_checkpoints checks the pandigital part of the two facts of the statement, not that \(541\) and \(2749\) come first: the tail from trailing_nine_digits_at_index(541), which runs the recurrence of Section 3, and the head from leading_nine_digits(2749) must be pandigital, or the program stops with exit code \(2\). Then solve searches.

is_pandigital_1_to_9 strips the digits of a tail or head, an unsigned 32-bit integer, with % 10 and / 10, sets bit \(d\) of a mask for each digit \(d\), returns false at once for a \(0\) or a repeated digit, and at the end requires nine digits and the mask 0x3FE, binary \(1111111110\), with exactly the bits \(1\) to \(9\) set. A tail with leading zeros has fewer than nine digits as a number and is rightly rejected, since its block contains a \(0\).

solve keeps \(t_{n-2}\) and \(t_{n-1}\) in the unsigned 32-bit variables a and b, starting from \(F_1=F_2=1\), and for \(n=3,4,\dots\) adds them in 64 bits and reduces the sum modulo \(10^9\), as in Section 3. Only when the new tail is pandigital does it call leading_nine_digits(n), the formula of Section 6, with the two logarithms as static constants and all arithmetic in long double. Before the conversion to an integer it adds \(10^{-10}\): in the 80-bit format, whose neighbouring values near \(10^8\) are about \(7\cdot10^{-12}\) apart, this lifts values lying within \(10^{-10}\) below a whole number, though it cannot repair errors as large as those of Section 8. For a double, whose neighbours near \(10^8\) are about \(1.5\cdot10^{-8}\) apart, it changes nothing. A ten-digit result would be divided by \(10\). The first index at which both tests pass is returned and printed.

The Python and Java programs follow the same plan with ordinary doubles, without checkpoints, and with \(\log_{10}\sqrt5\) computed as half of \(\log_{10}5\). Python pads the tail with zfill(9), compares the set of its characters with the digits \(1\) to \(9\), takes the first nine characters of str(int(10**(frac + 8))) as the head and scans \(3\le n\lt10^6\), printing nothing if no index there qualified. Java keeps the tails in long variables, pads with String.format("%09d", b), marks digits in a boolean array and cuts the head from the string of (long) Math.pow(10, frac + 8). String.format writes the digits of the default locale: under an Arabic or Persian locale they are Eastern Arabic digits, for which c - '0' gives an index far outside that array, and the program stops with an ArrayIndexOutOfBoundsException. It should therefore be run with a locale that writes ASCII digits, for example with the option -Duser.language=en.

Complexity and Verification

If the search stops at index \(k\), it makes \(k-2\) steps of the tail recurrence, each an addition, a remainder and a test of at most nine digits, so the time grows linearly with \(k\), against \(k^2\) for the direct method of Section 4. The head costs one multiplication, one floor and one power of \(10\) and is needed only at pandigital tails: about one index in \(2645\) over a whole period, and \(345\) times up to \(10^6\). The memory is constant. On an Apple-silicon Mac the C++ program takes about \(0.01\) seconds, Python \(0.13\) seconds and Java \(0.21\) seconds.

Besides the two checkpoints, separate code checked the following for this article. Exact big integers confirm the digit counts \(113\) and \(575\), that \(F_{541}\) and \(F_{2749}\) are the first numbers with a pandigital tail and a pandigital head, and the tail \(t_n\) for every \(n\le10^6\); at each of the \(345\) pandigital tails a \(40\)-digit evaluation of \(\lfloor10^{f+8}\rfloor\) gives the true head, as Section 7 predicts. Exactly one index up to \(10^6\) has both ends pandigital, and all three programs print it. On the Apple-silicon Mac the double-precision heads are wrong at \(1\) to \(8\) of these tails, depending on the language and the compiler, but never where it changes a decision. The \(387\) pandigital heads of Section 10 were counted without logarithms, from a lower and an upper bound for the leading \(60\) digits of every \(F_n\), carried as integers and rounded outwards whenever a digit is dropped; at every index from \(40\) to \(10^6\) the two bounds agree on the first nine digits. The period and the count \(567168\) come from running the recurrence through a whole period, and the smaller versions of Section 10 were found with \(50\)-digit heads and confirmed with exact integers.

Footnotes and References

The references below supply the background facts; the arguments about tails, heads, rounding and the existence of solutions are derived here.

  1. Project Euler 104 — Pandigital Fibonacci Ends. The official statement defines Fₙ = Fₙ₋₁ + Fₙ₋₂ with F₁ = F₂ = 1, notes that F₅₄₁ (113 digits) is the first Fibonacci number whose last nine digits are 1–9 pandigital and F₂₇₄₉ (575 digits) the first whose first nine digits are, and asks for the index k of the first Fibonacci number with both properties.
  2. Wikipedia — Fibonacci sequence. Binet's formula Fₙ = (φⁿ − ψⁿ)/√5; the identity φⁿ = Fₙφ + Fₙ₋₁, proved by induction; Fₙ is the integer closest to φⁿ/√5, since |ψⁿ/√5| < 1/2 for n ≥ 0; and for every d > 1 there are 4 or 5 Fibonacci numbers with d decimal digits.
  3. Wikipedia — Pandigital number. A pandigital number contains every digit of its base; a base-10 number without the digit 0, such as 923456781, may also be called pandigital, or zeroless pandigital.
  4. Wikipedia — Modular arithmetic. Integers are congruent modulo m when their difference is a multiple of m, and congruence is compatible with addition: a₁ ≡ b₁ and a₂ ≡ b₂ give a₁ + a₂ ≡ b₁ + b₂ (mod m).
  5. Wikipedia — Pigeonhole principle. If n items are put into m containers with n > m, at least one container holds more than one item; the principle is also known as Dirichlet's box principle.
  6. Wikipedia — Pisano period. The Fibonacci numbers modulo n repeat with the Pisano period π(n); π(2ᵏ) = 3·2ᵏ⁻¹, π(5ᵏ) = 4·5ᵏ, π(mn) is the least common multiple of π(m) and π(n) for coprime m and n, and the table lists π(10) = 60.
  7. Wikipedia — Golden ratio. φ ≈ 1.618 satisfies φ² = φ + 1, the other root of x² − x − 1 = 0 is (1 − √5)/2 ≈ −0.618, and φ is irrational.
  8. Wikipedia — Common logarithm. The characteristic, the integer part, counts how far the decimal point must be moved, and numbers that differ by a power-of-10 factor share the mantissa, the fractional part: log₁₀ 120 = 2 + log₁₀ 1.2 ≈ 2.07918.
  9. Wikipedia — Long double. Most x86 C compilers implement long double as the 80-bit extended precision type; on 64-bit ARM it is a synonym for double under Windows and macOS, and the quadruple-precision format under Linux.
  10. Wikipedia — Double-precision floating-point format. The 53-bit significand (52 bits stored) gives 15 to 17 significant decimal digits; from 2ⁿ to 2ⁿ⁺¹ the representable numbers are spaced 2ⁿ⁻⁵² apart.
  11. Wikipedia — Equidistribution theorem. For irrational a, the sequence a, 2a, 3a, … modulo 1 is uniformly distributed; this was proved in 1909 and 1910 separately by Hermann Weyl, Wacław Sierpiński and Piers Bohl.
  12. Wikipedia — Benford's law. A number begins with the digit string n with probability log₁₀(1 + 1/n); the Fibonacci numbers provably obey the law, and for kⁿ it follows from the equidistribution theorem when log₁₀ k is irrational.
  13. C++ — Euler104.cpp. The linked immutable revision contains the functions described above and the option --skip-checkpoints; it is the implementation record for this article.

Problem 104 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

bool is_pandigital_1_to_9(std::uint32_t value) {
    int mask = 0;
    int count = 0;
    while (value > 0) {
        const int d = static_cast<int>(value % 10U);
        if (d == 0) {
            return false;
        }
        const int bit = 1 << d;
        if ((mask & bit) != 0) {
            return false;
        }
        mask |= bit;
        value /= 10U;
        ++count;
    }
    return count == 9 && mask == 0x3FE;
}

std::uint32_t leading_nine_digits(const int n) {
    static const long double log10_phi = std::log10((1.0L + std::sqrt(5.0L)) / 2.0L);
    static const long double log10_sqrt5 = std::log10(std::sqrt(5.0L));

    const long double x = static_cast<long double>(n) * log10_phi - log10_sqrt5;
    const long double fractional = x - std::floor(x);
    long double v = std::pow(10.0L, fractional + 8.0L);

    std::uint32_t lead = static_cast<std::uint32_t>(v + 1e-10L);
    while (lead >= 1000000000U) {
        lead /= 10U;
    }
    return lead;
}

std::uint32_t trailing_nine_digits_at_index(const int target_index) {
    if (target_index <= 2) {
        return 1;
    }

    std::uint32_t a = 1;
    std::uint32_t b = 1;
    constexpr std::uint32_t mod = 1000000000U;

    for (int n = 3; n <= target_index; ++n) {
        const std::uint32_t c = static_cast<std::uint32_t>((static_cast<std::uint64_t>(a) + b) % mod);
        a = b;
        b = c;
    }
    return b;
}

int solve() {
    constexpr std::uint32_t mod = 1000000000U;
    std::uint32_t a = 1;
    std::uint32_t b = 1;

    for (int n = 3;; ++n) {
        const std::uint32_t c = static_cast<std::uint32_t>((static_cast<std::uint64_t>(a) + b) % mod);
        a = b;
        b = c;

        if (!is_pandigital_1_to_9(b)) {
            continue;
        }

        const std::uint32_t lead = leading_nine_digits(n);
        if (is_pandigital_1_to_9(lead)) {
            return n;
        }
    }
}

bool run_checkpoints() {
    if (!is_pandigital_1_to_9(trailing_nine_digits_at_index(541))) {
        std::cerr << "Checkpoint failed for trailing pandigital index 541" << '\n';
        return false;
    }
    if (!is_pandigital_1_to_9(leading_nine_digits(2749))) {
        std::cerr << "Checkpoint failed for leading pandigital index 2749" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

# Problem 104: Pandigital Fibonacci ends
# Find the first Fibonacci number where both the first and last 9 digits are 1-9 pandigital.

import math

def is_pandigital(s):
    return len(s) == 9 and set(s) == set('123456789')

def solve():
    MOD = 10**9
    a, b = 1, 1
    log_phi = math.log10((1 + math.sqrt(5)) / 2)
    log_sqrt5 = math.log10(5) / 2
    
    for n in range(3, 1000000):
        a, b = b, (a + b) % MOD
        if not is_pandigital(str(b).zfill(9)):
            continue
        # Check leading digits using logarithms
        log_fn = n * log_phi - log_sqrt5
        frac = log_fn - int(log_fn)
        leading = str(int(10**(frac + 8)))[:9]
        if is_pandigital(leading):
            print(n)
            return

solve()

Java

public class Euler104 {
    static boolean isPan(String s) {
        if (s.length() != 9)
            return false;
        boolean[] seen = new boolean[10];
        for (char c : s.toCharArray()) {
            int d = c - '0';
            if (d == 0 || seen[d])
                return false;
            seen[d] = true;
        }
        return true;
    }

    public static void main(String[] args) {
        long MOD = 1000000000L, a = 1, b = 1;
        double logPhi = Math.log10((1 + Math.sqrt(5)) / 2), logSqrt5 = Math.log10(5) / 2;
        for (int n = 3;; n++) {
            long c = (a + b) % MOD;
            a = b;
            b = c;
            String tail = String.format("%09d", b);
            if (!isPan(tail))
                continue;
            double logFn = n * logPhi - logSqrt5;
            double frac = logFn - Math.floor(logFn);
            String lead = Long.toString((long) Math.pow(10, frac + 8)).substring(0, 9);
            if (isPan(lead)) {
                System.out.println(n);
                return;
            }
        }
    }
}