Project Euler 57: Investigate the expansion of the continued fraction for the square root of two.

Project Euler 57: Investigate the expansion of the continued fraction for the square root of two.

I found Problem 57 of Project Euler to be a rather interesting problem, with more than one solution. The problem description reads

It is possible to show that the square root of two can be expressed as an infinite continued fraction.

√ 2 = 1 + 1/(2 + 1/(2 + 1/(2 + … ))) = 1.414213…

By expanding this for the first four iterations, we get:

1 + 1/2 = 3/2 = 1.5
1 + 1/(2 + 1/2) = 7/5 = 1.4
1 + 1/(2 + 1/(2 + 1/2)) = 17/12 = 1.41666…
1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 = 1.41379…

The next three expansions are 99/70, 239/169, and 577/408, but the eighth expansion, 1393/985, is the first example where the number of digits in the numerator exceeds the number of digits in the denominator.

In the first one-thousand expansions, how many fractions contain a numerator with more digits than denominator?

First thing I will present is a brute force solution, which for all practical purposes are fast enough. The second solution is a closed form approximation to the problem, so it can be solved as fast as you can punch the calculator.

Brute Force Solution

The approach to the brute force solution is to do a bit of pattern recognition of the fractions, a bit like the once you get on an IQ test, please continue the this sequence….

We have

From this we can derive the following two formulas

 

 

where n means the numerator, and d means the denominator.  To calculate this in C# we would need two temporary variables to store the old values of the numerator and denominator, but it can be rewritten to

which means we only need two variables in order to store it. For all practical purposes this means nothing for the performance, but it makes the code a bit cleaner.
To check the number of digits we can use the log10 function and compare the integer parts. The total code for this looks like

const int limit = 1000;
int result = 0;

BigInteger den = 2;
BigInteger num = 3;

for (int i = 1; i < limit; i++) {
     num += 2 * den;
     den = num - den;
     if ((int)BigInteger.Log10(num) >(int)BigInteger.Log10(den))
         result++;
}

Unfortunately the expansions result in some rather huge numbers, and thus we have to use BigInteger if we want to get the result. This slows down the performance a bit which shows from the 8ms the solution took:

The number of iterations with more digits in numerator is 153
Solution took 8 ms

Approximate Solution

We know from the problem description that the first time we encounter a numerator with more digits than the denominator is after 8 expansions. For some reason this got me a bit curious to see when the next solution would arise. This happens after expansion 13 and then after 21 and 26. So after a few iterations a pattern started emerging that we would have 8,5,8,5,8,5 expansions between the solutions to the problem. I tried to iterate several more times and got the following

1234567x1234x1234567x1234x
1234567x1234x1234567x
1234567x1234x1234567x1234x
1234567x1234x12x1234x
1234567x1234x1234567x1234x
1234567x1234567x1234x
1234567x1234x1234567x1234x
1234567x1234567x1234x
1234567x1234x1234567x
1234567x1234x1234567x1234x
1234567x1234x1234567x

which shows that there is a trend of 13 solutions per 13 expansions, but the pattern is not completely regular. However, for fun and giggles lets try to assume that the pattern was regular and see how close we get. We could form a solution looking like

where the x is the limit and the means floor.

Plotting the error term of this approximation against the actual value of problem shows that the approximation is within +/- 1 of the actual solution.

In 68% of the cases between 1 and 1000 it hits the target, in 8,4% it shoots 1 above and in 23,6% it is just one short. For the value 1000 it hits spot on. So even if it is not a completely accurate approximation it is rather close.

Wrapping up

I presented you with two solutions, an accurate and rather fast brute force approach and an approximate but closed form solution. It might be possible to enhance the closed form solution, I haven’t investigated that in depth yet.

As usual you can find the source code, if you want a close peek at it.

How did you approach the problem?

Posted by Kristian

15 comments

I found the same formula for generating n/d using the last n/d after inspecting the problem for a bit.
And then I just looped 1000 times and incremented when log10(n) > log10(d).

Exactly the same as you did. 🙂

I found different (but equivalent) formulas for N and D: [latex]D_k=2D_{k-1}+D_{k-2}, \quad N_k= D_k+D_{k-1}[/latex]. Now for D this is a second order linear recurrence, which means we can get closed forms:

[latex]
D_k=\frac{(2-\sqrt2)(1-\sqrt2)^k+(2+\sqrt2)(1+\sqrt2)^k}{4}
[/latex]

[latex]
N_k=\frac{(1-\sqrt2)^{k+1}+(1+\sqrt2)^{k+1}}{2}
[/latex]

However I didn’t manage to find a clever way to use those closed forms to solve the problem at hand (even for brute forcing they are less efficient than the recurrence formulas).

I don’t mean to censor you, but I deleted your second comment and updated the second formula so it parses.

I have the exact same issues with the latex formulas as you. I just have the ability to edit 🙂

Nicely found. This would be a good thing if you had to answer for one specific number. But you are right, since we need to check every number anyway, closed form or recurrence does not make much of a difference.

Hello there i found this amazing pattern.
The first 8 sequence is S[3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408, and 1393/985].
S2 is formed from S1, S1 numerator + S1 denominator is S2’s denominator.
S2’s denominator plus S1’s denominator is S2’s numerator.

Very nice approximate !!!
Thanks

I have written the code in C++. The answer it is giving me is 123 which is wrong. Can anyone help me with it?

#include
#include
using namespace std;

int main()
{
long double d1=2,n1=3;
int count=0;
for (int i=0;i(long)log10(d1))
{
count++;
cout<<"num is "<<n1<<" and den is "<<d1<<"\n";
}
}
cout<<count;
}

Pretty sure your brute force solution will get a significant speed up if you do this instead of performing a log on 2000 numbers (some of which are ridiculously big):


proc solve {{n 1000}} {
puts "57. In the first $n expansions of the square root of 2, how many fractions contain a numerator with more digits than denominator?"

set answer 0
set num 3
set den 2
set power10 10

for {set expansion 1} {$expansion=$power10} {
# if the denominator is still less than 10/100/1000/etc. then increase the answer
if {$den<$power10} { incr answer }
set power10 [expr {$power10*10}]
}

# next expansion
set num [expr {$num+2*$den}]
set den [expr {$num-$den}]
}

return $answer
}

So most iterations, the only thing I do is check to see if the numerator has gone over the next power of 10, which is a simple comparison. On the occasions when it has, check if the denominator made it over too. If not, add to the answer.

That code’s Tcl btw

Jean-Marie Hachey

Table 1
Number of fractions containing a numerator with more digits
than the denominator (NFND)
Re: Project Euler Problem 57

http://img11.hostingpics.net/pics/838263pe57tab1.jpg

Table 1 shows that in the first 1000 expansions, 153 fractions contain
a numerator with more digits than the denominator; this quantity reaches a value of 1508 fractions after 10 000 expansions.
Slopes of consecutive points also show that the points are not collinear.

___

Table 2
Distribution of digits in the numerator and denominator.
Re: Project Euler Problem 57

http://img11.hostingpics.net/pics/138849pe57tabletwo.jpg

In Table 2, as predicted with the algo, we find 6 fractions (bold) where the number of digits in the numerator is greater.

___

Sources:

1) Project Euler 57
Kristian’s algorithm
http://www.mathblog.dk/files/euler/Problem57.cs
2) Microsoft Visual C# 2010 Express
(Reference added : System.Numerics)
3) Microsoft Office Excel (2007)
4) General interest:
The Square Root of Two to 1 Million Digits
http://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil

Hey Kristian,
The pattern you are looking for is Parabola. The repeating sequence 8,5,8,5… will make a parabola (Or hyperbola) like shape.

Graph

Last update was my bad
This is the sequence (tested upto 10K)

(8,5,8,5,8,5,8,8,5,8,5,8,5,3)
(5,8,5,8,5,8,8,5,8,5,8,5,8,8,5,8,5,8,8,5,8,5,8,5,8,8,5,8,5,8,5,3)*x

Graph

It follows a rather odd pattern, I played around with other square root convergents. Here’s a graph for the square root of 2 (I just plotted the difference on the y axis, n = 15000): http://imgur.com/a/9xr1e (it’s rather difficult to see the exact pattern, but it gives you an idea of what’s going on). Here’s one with n = 1500: http://imgur.com/a/CU1PA. And just for fun sqrt(10), n = 1500: http://imgur.com/a/5YWjr. It’s quite interesting.

Hi Kristian,

I tried to avoid log10(x) or len(str(x)) calculations by keeping extra variables that are powers of 10, to make comparing the number of digits cheap:

# Crazy fast, avoid log10(x) or len(str(x)): 0ms :D
def solution3(N=1000, verbose=False):
n_gt_d = 0
(num, den) = (1, 1)
# to avoid len(str(x)) calculations, we keep tenpowers for num, den:
(tpn, tpd) = (10, 10)
# so: num < tpn, den = tpn:
tpn *= 10
if den >= tpd:
tpd *= 10
if tpn != tpd:
n_gt_d += 1
if verbose:
print "%4d %d/%d %f" % (i+1, num, den, float(num)/den)
return n_gt_d

Indentation is essential in python, but got lost here. Should’ve been
<code language=”python”>

Leave a Reply