Project Euler 80: Calculating the digital sum of the decimal digits of irrational square roots.

Project Euler 80: Calculating the digital sum of the decimal digits of irrational square roots.

Based on the headline Problem 80 of Project Euler sounds easy. It seems we just have to calculate the square root to a certain precision. So lets look at the problem description

It is well known that if the square root of a natural number is not an integer, then it is irrational. The decimal expansion of such square roots is infinite without any repeating pattern at all.

The square root of two is 1.41421356237309504880…, and the digital sum of the first one hundred decimal digits is 475.

For the first one hundred natural numbers, find the total of the digital sums of the first one hundred decimal digits for all the irrational square roots.

Indeed this was the core of the problem, calculating the square root of some numbers to a certain precision. The first thing we should note is that all square roots are irrational unless they are perfect squares. So except from 1,2,4,9,25,26,47,64,81 and 100 we need to calculate all the numbers.

Calculating the square root

The second thing is that we need a method to calculate the square root. Some programming language such as Java and Python have support for arbitrary precision floating point numbers. I know that SuprDewd made a library for BigDecimal  that as well but that felt kind of cheating if I used that. So I was looking at the wikipedia page for calculating square roots. And there are many good ones which are probably very fast. However, a foot note linking to a paper called square root by subtraction by Frazer Jarvis caught me eye.

The difference between this and other methods is that it keeps all calculations in integers. That means we can easily implement an algorithm without floating point, something I love since I find it much easier to handle precision in non-floating point numbers.

I will briefly go through the algorithm and suggest that you read the paper (which is an easy read) in order to dig into the algorithm.

If we want to calculate the square root of n, we need two numbers a and b which starts such that a = 5n and b = 5.

The we have the repeated steps which can be one of two

if a>= b: a = a-b and b = b+10.

if a< b: a = a*100 and add a zero just before the last digit of b  (which will always be ‘5’).

By repeating these steps b will approach the square root of n.  Apart from the last two digits of b this is the correct square root down to that precision. I don’t think it is the fastest algorithm out there, but it is one that can be easily implemented and one that can be done by hand if necessary.

My implementation in C# looks like this

private BigInteger Squareroot(int n, int digits) {
    BigInteger limit = BigInteger.Pow(10, digits + 1);
    BigInteger a = 5 * n;
    BigInteger b = 5;

    while (b < limit) {
        if (a >= b) {
            a -= b;
            b += 10;
        } else {
            a *= 100;
            b = (b/10) * 100 + 5;
        }
    }

    return b/100;
}

It takes two input the number we want to calculate the square root of and the precision. It uses the BigInteger class in order to handle the many digits, if it could operate on normal integers I don’t know.

Putting it all together

With the square root algorithm home we have solved the core of the problem, now me miss a digit sum calculation and a main loop. I have shamelessly stole the digit sum algorithm from Problem 56. The main loop of the program combines these and avoid the perfect squares. My implementation looks like

int result = 0;
int j = 1;

for (int i = 1; i <= 100; i++) {
    if (j * j == i) {
        j++;
        continue;
    }
    result += DigitSum(Squareroot(i, 100));
}

Executing this gives us

Total of the digital sums is 40886
Solution took 26,5812 ms

So not the fastest implementation around, but the ease of implementing more than makes up for the execution time I think.

Wrapping up

This was problem 80 we have now solved using an algorithm to calculate square roots to an arbitrary precision. The execution time is not that good, but the algorithm was rather easy to implement. I have only posted parts of the source code. You can find the full source code here.

Since algorithms for calculating square roots are numerous, how did you solve it?

Posted by Kristian

16 comments

When I first read the description for this problem I had just discovered Wolfram Mathematica (similar to MATLAB, which you are most likely familiar with), so I solved this problem using it:
Map[Total[(RealDigits[Sqrt[N[#, 200]]] // First)[[1 ;; 100]]] &, Select[Range[2, 100], Sqrt[#] != Floor[Sqrt[#]] &]] // Total

But since you mentioned my BigDecimal class, I thought I should try using it for this problem:

int curSq = 1, curN = 1, sum = 0;
for (int i = 1; i <= 100; i++)
{
    if (i == curSq) // My perfect-square filter is way cooler ;D
    {
        curN += 2;
        curSq += curN;
        continue;
    }

    // I turn it into a string, since working with the mantissa directly would probably be pain
    string sqrt = BigDecimal.Sqrt(new BigDecimal(i, 100)).ToString();
    for (int j = 0, k = 0; j < sqrt.Length && k < 100; j++)
    {
        if (sqrt[j] != '.')
        {
            sum += sqrt[j] - '0';
            k++;
        }
    }
}

Console.WriteLine(sum);

This runs in about 200 ms.
Now I’ve done two pretty cheat-ish solutions. Not cool 🙁
But I’m going to read that paper. Thanks.

Yes I am familiar with Mathematica. However, I haven’t used it that much. I have used Maple a bit. I do believe that Mathematica (and Maple) are quite different from Matlab. The latter is for numerical calculations. The former are much more symbolic minded. So different purposes, but all of them useful tools.

Neither of your solutions is cheatish in my opinion. The first one uses the best tool for the problem, and the second one is made with code you did all by your self. Something which is kind of cool.

Mind helping me out? This is my program output, showing the number and its square root’s digital sum, with the answer on the last line. IDK what I’m doing wrong!

2 475
3 441
4 0
5 472
6 470
7 397
8 464
9 0
10 456
11 482
12 404
13 416
14 483
15 498
16 0
17 447
18 397
19 469
20 476
21 481
22 404
23 450
24 472
25 0
26 451
27 393
28 396
29 455
30 463
31 469
32 467
33 468
34 455
35 435
36 0
37 452
38 456
39 459
40 481
41 445
42 440
43 497
44 478
45 435
46 445
47 418
48 393
49 0
50 464
51 461
52 426
53 470
54 455
55 426
56 460
57 447
58 403
59 492
60 473
61 429
62 435
63 433
64 0
65 433
66 457
67 491
68 444
69 431
70 457
71 409
72 417
73 411
74 420
75 428
76 452
77 477
78 449
79 453
80 475
81 0
82 447
83 466
84 448
85 418
86 450
87 450
88 401
89 400
90 469
91 403
92 493
93 443
94 494
95 451
96 466
97 432
98 468
99 438
233 Milliseconds elapsed.
40320

I can’t tell you exactly. My guess was that you were making a rounding error, but some of the numbers are several places off. I have a list of the first 10 digit sums, that should help you debugging the problem:
2 475
3 441
4 0
5 473
6 471
7 398
8 465
9 0
10 459
11 484

What solution are you pursuing, and can you share it with us once you iron out the last few kinks?

/Kristian

I was trying to do this problem using an old java BigInteger square root (Heron’s Method) class that I had used before. You can look at it here: http://andrewkoroluk.com/euler/code/BigSquareRoot.java

Hey man thanks for your help. I figured out that I needed to set my Heron’s Method digits to 102, multiply the answer by 10^99, and add those digits. Also, I just figured out that I needed to discard all rational square roots from the final sum. That’s 102 problems down! Thanks!

Here’s my final Prob 80 source code: http://andrewkoroluk.com/euler/code/euler80.java

Well you mostly if not figured it out yourself. I just provided a few samples for you to test your code on. But thanks for the source code, always fun to see what other people do. Well done.

/Kristian

Abdullah

Here is my APL (NARS2000) short and fast implementation that took 228 ms

+/+⌿(100⍴10)⊤⍎,’ ‘,((90 1↑z),0 2↓z←90 101↑⍕⍪√0x+(⍳99)∼(⍳9)*2),’x’
40886

Just saw this thread and couldn’t help but respond. Here’s my APL (NARS2000) code (5 ms):

⎕PP←110 ⋄ ⎕FPC←2⍟10*⎕PP ⋄ ⎕IO←0
N←1..100v ⋄ +/’0123456789’⍳∊100⍴¨(⍕¨√N~N*2)~¨⊂’. ‘
40886

This solution takes advantage of NARS2000’s support for MPFR (variable precision floating point — VFP) numbers. The assignment to ⎕PP specifies how many digits of printing precision we want (a few more than we need so as to avoid round-off error). The assignment to ⎕FPC specifies the number of bits of precision MPFR is to use. The assignment to ⎕IO specifies that origin-dependent code is done in origin-zero (in particular the table lookup). The assignment to N specifies the numbers from 1 to 100 as VFP numbers.

N~N*2 is the set difference between N and its squares, thus leaving only those numbers whose square root is irrational.
√ is the square root primitive function.
⍕¨ formats as a character vector each value in its argument.
~¨⊂’. ‘ is the set difference of each formatted number and ‘. ‘ — that is it removes all dots and blanks, leaving digits only.
100⍴¨ reshapes each number to the first 100 digits.
∊ flattens the nested array into one long stream of digits.
‘0123456789’⍳ performs a table lookup of each digit in the long stream in the string ‘0123456789’ which in effect converts each digit to a number.
+/ sums the vector of numbers to produce the final result.

I did things a little differently. As people pointed out, there’s no BigDecimal class in C#. So rather than writing my own, I used BigInteger to take the square root of each number multiplied by 10^200, which will give 100 digits of precision after taking the square root. I started with an initial guess by taking the square root using doubles and multiplying that by 10^100; then finally polishing the root with Newton-Raphson.
Takes 7 ms. Here’s the code:

class Sum
{
public static int OfFirst100DigitsOfTheSquareRootOf(int x)
{
int sqrt = (int)Math.Sqrt(x);
if (sqrt*sqrt == x)
return 0;

// there is no BigDecimal class, so take square root of x*1E200,
// which will give us an answer with 100 digits of precision.
BigInteger n = x;
for(int i=0; i<100; ++i)
n *= 100;

// take double square root and convert to our first guess, using string/parse
var input = Math.Sqrt(1E26*(double)x).ToString();
int period = input.IndexOf('.');
if (period != -1) input = input.Substring(0, period);
while (input.Length <= 100) input = input + "0";

BigInteger n0 = BigInteger.Parse(input);
BigInteger n1 = 0;
while (true)
{
n1 = (n0 + n/n0)/2; // Newton-Raphson
if (n1 == n0)
break;
n0 = n1;
}

// sum of first 100 decimal
var s = n0.ToString();
int sum = 0;
for (int i=0; i<100; ++i)
sum += s[i] - '0';
return sum;
}
}

[TestClass]
public class UnitTest1
{
[TestMethod]
public void First100IrrationalSquareRoots()
{
int sum = 0;
for(int i=0; i<=100; ++i)
sum += Sum.OfFirst100DigitsOfTheSquareRootOf(i);
Assert.AreEqual(40886, sum);
}
}

Jean-Marie Hachey

Excerpt from the statement of Problem Euler 80:
“The square root of two is 1.41421356237309504880…, and the digital sum of the first one hundred decimal digits is 475.”
http://projecteuler.net/problem=80

Comment :
475 is the sum of the first 100 digits including integer part and decimal part of number 1.4142…
___

Table 1 shows the first 10 non-perfect-square numbers (2-13) and the sums (specific and cumulative) of their two-first-digit-square roots. One can see that the sums include both integer part and decimal part of the square roots.
Cumulative sums calculated are confirmed by the algorithm (2).
(Reference System.IO not used in this algo.)

http://img15.hostingpics.net/pics/377198PE80sqrttable1.jpg
___

Sources:

1) “1 followed by numbers that are not squares.”
http://oeis.org/A164514
2) Project Euler 80
Kristian’s algorithm.
http://www.mathblog.dk/files/euler/Problem80.cs
3) Microsoft Visual C# 2010 Express
(Reference added: System.Numerics)
4) Microsoft Office Excel (2007)

returns the Answer in “0.0337419509888 seconds”

from math import *
def sqroot(N) : #For irrational square roots.
r = int(sqrt(N))
k = N – (int(sqrt(N))**2)
Q = k*100
R = 20*r
i = 1
s = “”
while i <= 100 – len(str(r)) : ##100 specifies total digits in it.
j = 1
while (R+j)*j < Q :
j += 1
s = s + str(j-1)
Q = (Q – ((R+j-1)*(j-1)))*100
p = ((R+j-1))
R = (p + p%10)*10
i = i + 1
#return str(r) + "." + s ##Returns the square root as a "100"(can be changed) digit number.
u = 0
for m in xrange(len(s)) : ##Returns sum of the digits for the first k digits
u = u + int(s[m])
if u == 0 :
return 0
else :
return u + r

k = 0
for i in xrange(1,101) :
if i not in [1,4,9,16,25,36,49,64,81,100] :
k = k + sqroot(i)
print k

I have used the methodology of Method 2 in http://www.wikihow.com/Calculate-a-Square-Root-by-Hand.

returns the Answer in “0.0337419509888 seconds”

from math import *
def sqroot(N) : #For irrational square roots.
r = int(sqrt(N))
k = N – (int(sqrt(N))**2)
Q = k*100
R = 20*r
i = 1
s = “”
while i <= 100 – len(str(r)) : ##100 specifies total digits in it.
j = 1
while (R+j)*j < Q :
j += 1
s = s + str(j-1)
Q = (Q – ((R+j-1)*(j-1)))*100
p = ((R+j-1))
R = (p + p%10)*10
i = i + 1
#return str(r) + "." + s ##Returns the square root as a "100"(can be changed) digit number.
u = 0
for m in xrange(len(s)) : ##Returns sum of the digits for the first k digits
u = u + int(s[m])
if u == 0 :
return 0
else :
return u + r

k = 0
for i in xrange(1,101) :
if i not in [1,4,9,16,25,36,49,64,81,100] :
k = k + sqroot(i)
print k

Michael Press

That square root algorithm is an example of a spigot algorithm.
Since the algorithm computes the decimal digits directly it
looks like you can speed up your program by writing


int sum = 0;

while (b = b) {
a -= b;
b += 10;
sum++;
} else {
a *= 100;
b = (b/10) * 100 + 5;
}
}

return sum

Tnen


result += DigitSum(Squareroot(i, 100));

becomes


result += Squareroot(i, 100);

Also, why not


int sum = 0;

for(int i = 0; i = b) {
a -= b;
b += 10;
sum++;
} else {
a *= 100;
b = (b/10) * 100 + 5;
}
}

return sum

and eliminate the value limit
and all those big integer comparisons?

I implemented this scheme in C with the GNU Multiple Precision Libary
and it runs in 6 millisecond on an eight year old home desktop machine.

A simple optimization is to replace


b = (b/10) * 100 + 5;

with


b = 10 * b - 45;

and eliminate around 9000 divisions.

Stefan Gruenwald


from math import sqrt
from decimal import *
getcontext().prec = 105

result=0
for i in xrange(1,101):
temp=sum(map(int,list(str(int(Decimal(i).sqrt()*10**99)))))
if temp>=10:
result+=temp

print result

So I used the Wikipedia digit by digit calculation with some optimisations.

Firstly, you can completely ignore the decimal point. The algorithm will work using integers since we only need to know what the digits are, not where they are in relation to the decimal point. Secondly, since we are only calculating square roots of n where 2 <= n <= 99, only the first digit pair is non zero and it is equal to n.

My solution took under 0.5 seconds in Swift even though I had to roll my own Big Integer implementation and I found x simply by counting down from 9 instead of using a more clever method (I wanted to avoid having to implement division).

Leave a Reply