# Project Euler 28: What is the sum of both diagonals in a 1001 by 1001 spiral?

Problem 28 of Project Euler reads

Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:

21 22 23 24 25
20  7  8  9 10
19  6  1  2 11
18  5  4  3 12
17 16 15 14 13

It can be verified that the sum of the numbers on the diagonals is 101.

What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral formed in the same way?

I am pretty sure it is possible to brute force a solution for this problem, but the thought of having to construct the spiral, so let instead of that let us attack it from a more analytical approach and search for a function of the ring number. This is done in two steps, first we will derive a formula where the calculation for the n’th ring depends on the n-1’th ring. And after that we will derive a formula which is independent on the previous ring.

## Deriving the formula for the corners of the ring

Lets start with the definitions as any good mathematician would do.

n is the ring number, starting from 0. Such that n=1 is the ring

7 8 9
6   2
5 4 3

and so on.

f(n) is the sum of diagonals for the the rings up to and including n.

We will need to find the corner values in the ring, and since n=0 is not a ring but a point I will treat it as a special case. Such that f(0) = 1.

The sides of a ring is 2n+1 wide, and thus the upper right corner in the ring is (2n+1)2, since it is equal to the area of the square.

Moving counter clockwise, in order to get to the second corner, we will take the value of the first corner and deduct the side length minus 1 from. since both numbers are included in the side. Such that the formula for the second corner is (2n+1)2 – 2n.

In order to get the third corner we once again deduct 2n, and the same thing for the fourth corner.

Summing all the corners we get

4(2n+1)2 – 12n

That means the sum of the diagonals can be written as

f(n) = 4(2n+1)2 – 12n + f(n-1)

Now it should be pretty easy to make a for loop for calculating f(500). However, I want to take a bit further.

## Deriving a non-iterative formula

Since we are working with a fairly nice formula, there is good chance that we can fit a polynomial to get an analytical solution for f(n). We used the same trick in the solution for problem 6 though we never proved it. The method uses differences and can also be seen at Ken Ward’s Mathematics Pages.

In order to figure out what order polynomial we need, we need to calculate the function value of f(n) and we need to calculate the differences. Δ1(n) represents the differences between f(n) and f(n-1). Δ2(n) = Δ1(n) – Δ1(n-1). and so on

 n 0 1 2 3 4 5 6 f(n) 1 25 101 261 537 961 1565 Δ1 24 76 160 276 424 604 Δ2 52 84 116 148 180 Δ3 32 32 32 32

Since a the third difference gives us a constant, we will need a third order polynomial function in order to fit f(n). We want to find the factors for the function

f(n) = ax3 + bx2 + cx + d

That means a total of 4 unknown factors, which can be found if we solve four equations. So lets solve the equations related to f(0) – f(3).

d = 1

a + b+ c + d = 25

8a + 4b + 2c + d = 101

27a + 9b + 3c + d = 261

This can be solved in a multitude of ways.  My preferred way would be to state it as a linear algebra problem, but we could also solve it by isolating and inserting into the equations.

I have the videos from a course taught by Gilbert Strang on Linear algebra, from which you will need the first two lectures in order to be able to solve this set of linear equations.  I wont show you how to solve the set of equations here, but recommend you the videos. Along with that I will recommend you two books. The one I used in university Linear Algebra and Its Applications by David C. Lay and the other is by Gilbert Strang Introduction to Linear Algebra.

Many modern calculators like TI 83 is fully capable of calculating the unknowns for you. However, I will strongly recommend you to do it by hand until you master it. Only then should you do it on a calculator or a program. It sounds old and harsh to say that, but over and over again, I keep realising that I understand the math building on top of the fundamentals only when I master the fundamentals.

After my little trip on the soap box, the result of solving the four equations is that the function

f(n) = 16/3x3 + 10x2 + 26/3x + 1

So now we have a complete analytical solution to any size spiral. The solution to the specific question is

f(500) =16/3*5003 + 10*5002 + 26/3*500 + 1 = 669 171 001

No computer code this time, just one single formula.

### 35 comments

Ronnie

I felt good when I saw my procedure on ur site….
My analysis..

starting from 1, going up-right, we get $f(x)=4x^2$
going up-left, $g(x)= (2x)^2$
going down right, $h(x)= (2x-1)^2+2x$
going down left, $i(x)= (2x)^2+1$

Simply let, $i=\lfloor k/2 \rfloor$
$\sum_{x=1}^{i} f(x) + \sum_{x=1}^{i} g(x) + \sum_{x=1}^{i} h(x) + \sum_{x=1}^{i} i(x) +1$
$\implies \sum_{x=1 }^{i} {4(4x^2+x+)} +1$
$\implies \frac{2} {3} {i(8i^2+15i+13)} +1$

So, putting k=1001 => i=500
we get the answer easily

I used this to calculate actually, since I was lazy to pick the calculator

from math import floor
lng=1001;
i= int(floor(lng/2))
print (2*i*(8*i*i+15*i+13))/3+1

I have done basically the same, in a different way…. derived the same formula you can see

Kristian

Nice way of deducing it. The only thing I can see is that rather than having the iterator bgin n, you should make it x, or change the first set of functions to be function of n.

Ronnie

I am extrmely sorry.. Since I used in my computer as f(n) rather than x.
I pasted and forgot to edit
The last two are
$h(x)=(2x-1)^2+2x$
$i(x)=4x^2+1$

Now, hope you have got, and if comfortable, change the last two as well..

And I must say, U r friendly…. and addicted to mathematics like I do.

P.S: I must compliment, ur eyes are like compiler, always looks for tokens

Kristian

I just updated your comment. My main goal of the blog, besides having fun is to encourage people to like math. And I don’t think I achieve that by being a jerk. In fact I don’t think you will accomplish much of anything in life by being a jerk.

Ronnie

May I know the “jerkish” thing I have done…………..
Yes, I do some experiment, And I have been correct everytime.. I pretty much follow the same way as you do..

Do you know, No one taught me Programming(its just my passion), And I do feel proud for that..
I said I am addicted to mathematics and programming, but I dont think that by being a “jerk”…

In some questions in PE, I used some blind shots, I shot them, because I knew even if I miss, it wont hit me back since the one I’m shooting is a target not a human with another gun. and in both the time, My shot was correct

And since I, in fact we dont know much about primes. We got to guess on that topic. In fact, I research much on primes alone.. thats why, by my experience I guessed them to be prime… and look, both were..

If It is the Jerkish thing.. I am “Jerk”..

Neither you have hurt me. nor I meant to do.. If I have done. Sorry,

Kristian

I think you misunderstand me. I don’t mean you are a jerk. I was just giving a bit of my philosophy to the world upon you complementing me for being nice.

Actually you seem like quite the opposite of a jerk to me. I think you are doing good and not least since you dare to ask questions and look for patterns.

Yes you might miss once in a while. No matter how accomplished a researcher you are, you will miss. Otherwise it wouldn’t me research.

I’m a very amateur self taught coder and I’m really bad at maths. Due to this, I find solving Project Euler problems pretty tricky. Any solutions I’ve managed to hack together so far are, from a mathematical standpoint, pretty bad I would think.

I solved this problem because I noticed a pattern when I stared at the example sprial for a while: the numbers at the corners of the spiral are related. Every time you find the next upper right corner number you can then derive the other three corners from it.

You can work out what each next right hand corner number will be with a very simple formula:

Call the centre number 1 the current number. Each movement diagonally up and right is one iteration. Starting with iteration = 1, the right corners can be found like this:

currentNumber + (iteration * 8) = next value of currentNumber

For each currentNumber you can then derive each of the other three corner numbers as follows:

lower right corner = currentNumber – (iteration * 6)
lower left corner = currentNumber – (iteration * 4)
upper left corner = currentNumber – (iteration * 2)

So, to find the sum of all the diagonals of a 1001 x 1001 spiral, in code:

int spiralSize = 1001;
int upperRightCornersToFind = (spiralSize - 1) / 2;
int currentNumber = 1;
int total = currentNumber;

for (int i = 1; i <= rightCornersToFind; i++ {

currentNumber += 8 * i;
int bottomRight = currentNumber - (i * 6);
int bottomLeft = currentNumber - (i * 4);
int topLeft = currentNumber - (i * 2);

total += currentNumber + bottomRight + bottomLeft + topLeft;
}


I’m not sure if in mathematical terms this amounts to what you said already (I lack the mathematical vocabulary to understand the algebra) but it works. 🙂

Kristian

You are on the right track. What I do in the section “Deriving the formula for the corners of the ring” Is that instead of relaying on the previous number I derive an expression where I can calculate it based only on the current iteration.

In the next section I then use that expression to derive a formula where I can compute the answer by just inputting the iteration that we should sum over.

But as I said, finding the pattern you did is definitely on the right track. I can only encourage you to read the explanation a couple of times more to see if you can pick it up. The algebra is really pretty basic in the first section.

Joy

I think this is the sum of two simple series.
1+3+7+13+21+31+43+57+………and
1+5+9+17+25+37+49+65+……….
This is easier process. I wrote the code in python 3.3.

#1+3+7+13+21+31+43+57+……
a,sum1,dif=1,0,2
for i in range(1,1002):
sum1+=a
a+=dif
dif+=2

#1+5+9+17+25+37+49+65+…..
b,sum2,dif,k=1,0,4,0
for i in range(1,1002):
sum2+=b
b+=dif
if(i%2==0):
dif+=4
print(‘Here is the result:’,sum1+sum2-1)

Joy

Sorry. Indentations were missed in the for loops.

Kristian

Yep they work as series as well.

Gabriel

I did pretty much the same thing, but I got 668 670 001. What did I do wrong?
I found a quadratic equation for each diagonal, 4x^2-2x+1 and 4x^2-3+1. I added these together and got that the sum should be $\sum\limits_{x=-500}^500 {8x^2-5x+2}-1$
Using square pyramidal numbers, I found that equals $2*8*{2*500^3+3*500^2+500/over 6}+2002-1=668670001$

Kristian

It seems that I have a different formula for the quadratic equations of the diagonals. So maybe there is an error there?

One more code snippet which assignd values to every cell of the square accurately. Thus, one can compute not only the sum of diagnols but any patttern (say sume of all alternate rows/columns or sum of alternate cells etc.)

#define SQUARESIZE 1001
int main()
{
int SquareSize = SQUARESIZE;
int a[SQUARESIZE][SQUARESIZE];
int row=SquareSize/2,col=SquareSize/2,i,j,value=1,sum=0;
a[row][col] = value++;
for (i = 1; i < SquareSize; i+=2)
{
for (j=1;j<=i;j++)
a[row][++col] = value++;
for (j=1;j<=i;j++)
a[++row][col] = value++;
for (j=1;j<=i+1;j++)
a[row][–col] = value++;
for (j=1;j<=i+1;j++)
a[–row][col] = value++;
}
i–;
for(j=0;j<i;j++)
a[row][++col]=value++;
for (i=0;i0;i–)
sum+=a[SquareSize-i][i-1];
printf(“sum=%d\n”,sum-1);
return 0;
}
o/p:
sum=669171001

Chris Goldfrap

Like Gabriel I’ve been looking into equations for the straight line progressions, or spokes, to and from the Ulam centre. Not just the diagonals but also the vertical and horizontal rows, resulting in a figure similar to the British flag.

For a centre 0, they go 4x^2 (ie 4x^2 + 0x), 4x^2 + x, 4x^2 + 2x, 4x^2 + 3x, 4x^2 + 4x. For a centre 1, just add + 1 to make a standard quadratic, indeed for any centre n, add + n.

Five equations? But there are eight spokes, or if you will, four intersecting lines. That’s because there are in fact five lines: 4x^2 is just one straight line beginning and ending at the centre. It’s not bent at all. 4x^2 + x goes through the centre and is bent there through 1 right angle. 4x^2 + 2x through 2 right angles, so looks straight. 4x^2 + 3x through 3 right angles, so is bent the other way from 4x^2 + x, and 4x^2 + 4x is bent through 4 right angles, ie bent double so looks like one straight line.

For a hexagon number spiral on a chicken wire grid where the cells start off packed 6 round a central cell, the equations are 3x^2 etc. In an Ulam spiral on a squared sheet they’re packed 8 round the central cell. So maybe if P is the packing number, then the x^2 coefficient is P/2.

Madferret

I actually built the whole spiral in an matrix, noticing that in order to construct it, with the 1 already in place, the movements follow this pattern R D LL UU, RRR DDD LLLL UUUU, RRRRR DDDDD LLLLLL UUUUUU… and so on. I am not very good at analazing mathematical functions but I was able to code an algorithm.

while (counter <= movements - 1) {

 for (var a = 1; a < 5; a++) { if (a == 3) { rep++; } //log("Im going to write a " + map(a) + " " + rep + " times"); 

 for (var i = 0; i < rep; i++) { if (counter < movements) { makeMove(a); counter++; } } } rep++; } 

Kristian

I was considering that approach as well but decided against it after looking at the spiral. I guess I was too lazy to actually build it 🙂

Teun Pronk

I know this is an old post but I still wanted to post this and see if I can get a respond since im a bit puzzled here.

My method is different than your approach and my answer is different too.
Although if I run my method with a 5×5 grid the answer is also 101 like the example.

My approach.
There is a grid of NxN where N is 1001.
The center of the grid holds the value 1 and clockwise the value increases by 1.
Meaning the highest value is in the top right corner and has the value of N*N (Squareroot(N)).
The difference between the highest value and the next value is N-1.

You do this 5 times.
25-(5-1)=21; 21-(5-1)=17; 17-(5-1)=13 (Out ring complete).
To the next ring is the same interval.
13-(5-1)=9.
Since the ring is smaller you take your interval-2 to get to the next corner and you keep doing that till you reach 1.

In code:

function Problem28:string;
const
num=1001;
highNum=num*num;
var
sum:int64;
i,moves,step:integer;
begin
sum:=0;moves:=1;
i:=highNum;
step:=num-1;
repeat
sum:=sum+i;
i:=i-step;
inc(moves);
if moves=5 then
begin
step:=step-2;
moves:=0;
end;
until i<1;
writeln(sum);
readln;
end;


Yet I come to a total of 605700665.
How did this happen? Is there a flaw in my theory?
Any thoughts on this?

Thanks in advance 🙂

Kristian

I think the problem is that you check on moves=5, I think it should be equal to 4.

In your example you will have
25-4 = 21, moves = 1
21-4 = 17, moves = 2
17-4 = 13, moves = 3
13-4 = 9 , moves = 4
9-4 = 5 , moves = 5 (this is wrong)

Teun Pronk

Indeed, I noticed this but I couldnt edit my post anymore.
Stupid mistake lol…

hongtengyin

my approach:(python)

number = 1
sum_diagonals = -3

for n in range(0,501):
for x in range(0,4):
number += 2*n
sum_diagonals += number

print sum_diagonals

exalt

my solution of wich im proud of 🙂
still i feel like a noob compared to you guys…

#include &lt;iostream&gt;
#include &lt;cmath&gt;
using namespace std;
int main()
{
int limit = 1001, counter = 3, side = 1, sum = 1, square = 1;
while(side &lt; limit)
{
square = pow(counter, 3) / counter;
side = sqrt(square);
sum += square*4 – (side-1)*6;
counter += 2;
}
cout &lt;&lt; sum &lt;&lt; endl;
return 0;
}

exalt

my solution of wich im proud of 🙂
still i feel like a noob compared to you guys…

#include &lt;iostream&gt;
#include &lt;cmath&gt;
using namespace std;
int main()
{
int limit = 1001, counter = 3, side = 1, sum = 1, square = 1;
while(side &lt; limit)
{
square = pow(counter, 3) / counter;
side   = sqrt(square);
sum   += square*4 - (side-1)*6;
counter += 2;
}
cout &lt;&lt; sum &lt;&lt; endl;
return 0;
}


note: please remove my previous post 🙂

deepu

Int32 sum = 1, n = 2, p = 1, c = 4;
for (Int32 i = 2; i <= 501; i++)
{
c = 4;
while (c != 0)
{
p = p + n;
sum = sum + p;
c–;
}
n = n + 2;
}

Console.WriteLine("the sum is " + sum);
Console.ReadKey();

my solution :)!

dumindutraa

the required corner numbers are in an ordered sequense arent them???

1 (2) 3 (2) 5 (2) 7 (2) 9 (4) 13 (4) 17 (4) 21 (4) 25 (6) 31 (6) 37 (6) 43 (6) 49……

so jst a matter of taking the sum upto the number 1001*1001 Isnt it???

andrew

This can also be approached another way.

There are as many rings as it takes to move out of the center to the edge in any direction. Starting from the center, and considering width as one of the dimensions, there are a number of rings of
(width – 1) / 2. I accounted for the center value of one, but otherwise it’s just (width – 1) / 2 + 1.

Each ring is two points bigger in any direction than the ring before it. So the distance, and therefore the difference between diagonals is +2 more than it was the last time. It is a little like “drawing” the spiral while skipping all of the insignificant values between the four corners of each ring. You simply increase dx by two on each iteration and loop it four times. Here’s a sample code:

double width = 1001;
double sum = 1;
double currentValue = 1;
double dx = 0;
double squares = (width – 1) / 2;

for (int s = 0; s < squares; s++){
dx = dx + 2;
for (int i = 0; i < 4; i++){
currentValue = currentValue + dx;
sum = sum + currentValue;
}
}
System.out.println(sum);

I liked this approach because it was pretty straightforward and applied a very simple geometry, and in my case I was too lazy to do the math.

barry

 def euler28(): a = [] for i in range(0,1001): a += [[]] for j in range(0,1001): a[i] += [None]

 a[500][500] = 1 val = 2 orientation = [0,1] position = [500,501] def getCell(): return a[position[0]][position[1]] def setCell(v): a[position[0]][position[1]] = v def moveForward(): position[0] += orientation[0] position[1] += orientation[1] def rotateClockwise(): orientation[0], orientation[1] = orientation[1], -orientation[0] def outsideBoard(): return position[0] > 1000 or position[1] > 1000 while not outsideBoard(): setCell(val) rotateClockwise() moveForward() if getCell() is not None: rotateClockwise() rotateClockwise() moveForward() rotateClockwise() moveForward() val += 1 return sum(a[i][i] for i in range(0,1001)) + \ sum(a[i][1000-i] for i in range(0,1001)) - \ a[500][500] 

barry
def euler28():
a = []
for i in range(0,1001):
a += [[]]
for j in range(0,1001):
a[i] += [None]

a[500][500] = 1
val = 2
orientation = [0,1]
position = [500,501]

def getCell():
return a[position[0]][position[1]]

def setCell(v):
a[position[0]][position[1]] = v

def moveForward():
position[0] += orientation[0]
position[1] += orientation[1]

def rotateClockwise():
orientation[0], orientation[1] = orientation[1], -orientation[0]

def outsideBoard():
return position[0] &gt; 1000 or position[1] &gt; 1000

while not outsideBoard():
setCell(val)
rotateClockwise()
moveForward()
if getCell() is not None:
rotateClockwise()
rotateClockwise()
moveForward()
rotateClockwise()
moveForward()

val += 1

return sum(a[i][i] for i in range(0,1001)) + \
sum(a[i][1000-i] for i in range(0,1001)) - \
a[500][500]


David Th.

I found this one pretty easy, I solved it with:

sum=1
ult=1
for n in range(1,501):
sum+=20*n+4*ult
ult=ult+8*n

We get sum = 669171001

Each corner adds 2n from the last one, so in each ring we have 2n+4n+6n+8n=20n, and we have to add 4 times the last corner from the previous ring. That’s pretty much all.

[…] Although this problem is entirely solvable by using math alone, we chose to instead traverse the spiral and see if we notice a pattern. Let us start and 1 and […]

Miles

That’s brilliant!

I spent hours on a variant of this problem:

the N is not 1001, id could be up to 10^18; and the sum we wanted is too big, we just need sum % (1000000007);

I solved it in C++, however, I found the formula on some website, most of the time I spent is on Modular Arithmetic.

I could get the “Summing all the corners” formula myself, but don’t know a way to get the final formula.

You saved my life!

lioan

this is my solution, is not like yours cause i am almost a noob but works

def spiral_number_sum_deduction(n):
high_number = n ** 2 #this is the bigger number u can reach
number = 1000 #number of columns and rows – 1
result = 1
while number >= 1:
for i in range(1, 4):
result += high_number – number + 1 # plus 1 cause the number of rows and columns used was less 1
high_number -= number
number -= 2
result += high_number – number + 1
high_number -= number
return result

print(spiral_number_sum_deduction(1001))

meme

My approach was to derive the formula for the sum of each ring
which gave me the formula : n=1 is 1 otherwise 4n^2 – 6n + 6, where n is the nth ring
so just need to get the summation of the previous rings

javascript

var n = 1001;
var total = -3; //was set to -3 to get 1 for n=1
for (var i = 1; i <= n; i=i+2) {
total += (4*(Math.pow(i, 2))) - (6 * i) + 6;
}
console.log(total); //669171001


I was amazed on how you got a non-iterative formula (non really much of a math person :))
thanks for this post.

ksjv

your solution works if num is odd if n*n matrix is even means n is even then?

PostedBefore

Shouldn’t it be 16/3 * n^3 – 6*n^2 + 14/3*n – 3?
where you have 16/3*n^3 + 10*n^2 + 26/3*n + 1.