# Duncan's blog

## December 18, 2014

### Project Euler: problem 46 – Goldbach’s other conjecture

Filed under: PHP,Project Euler — duncan @ 11:57 pm
Tags: , , ,

I’m doing these Project Euler mathematical puzzles as a simple practical exercise for teaching myself PHP, and I’d appreciate any feedback on my code

It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a prime and twice a square.

9 = 7 + 2×12
15 = 7 + 2×22
21 = 3 + 2×32
25 = 7 + 2×32
27 = 19 + 2×22
33 = 31 + 2×12

It turns out that the conjecture was false.

What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?

This was another of those problems that I’d initially passed over, having decided it looked a bit tricky.  Then looking at it again realised it probably wasn’t too difficult.  And I turned out to be right!

A composite number is basically any positive integer that isn’t a prime number.

My logic is simply:

• Loop through odd numbers, from 3 upwards.
• For each number, if it’s prime, add it to an array of primes for later reference.
• If it’s not prime, work out if it matches the conjecture:
• Subtract the next largest prime, then examine the remainder
• Divide the remainder by 2.  if it’s a square number then we’re meeting the conjecture.  Move onto the next odd number.
• Otherwise keep subtracting the primes, checking the remainders.
• If we’ve looped through all the primes then we must have reached a number that doesn’t meet the conjecture.

```<?php
\$primes = [2];

\$i= 1;

while (true) {
\$i+= 2;

if (isPrime(\$i)) {
\$primes[] = \$i;
} else {
for(\$j = count(\$primes)-1; \$j > 0; \$j--) {
\$remainder  = \$i - \$primes[\$j];
\$remainder = \$remainder / 2;
\$root = sqrt(\$remainder);
if ((int) \$root == \$root) {
continue 2;
}
}

echo \$i;
break;
}
}

function isPrime(\$x)
{
\$root = sqrt(\$x);

for (\$i = 3; \$i <= \$root; \$i += 2) {
if (\$x % \$i == 0) {
return false;
}
}

return true;
}
```

I’m using a modified version of my original isPrime function, just because I know I don’t need to check for even numbers.

The loop backwards through the primes was maybe a bit of overkill; using a simple foreach loop through the primes in ascending order took more like 100ms.

What else… I check if a square root is the same as when it’s cast to an integer (not sure this is the best approach).  Then use continue 2; to get out of our inner-most loop and move onto the next value in our parent loop.

## December 14, 2014

### Project Euler: problem 33 (PHP) – Digit cancelling fractions

Filed under: PHP,Project Euler — duncan @ 10:54 pm
Tags: ,

Picture by jessicakelly

I’m doing these Project Euler mathematical puzzles as a simple practical exercise for teaching myself PHP, and I’d appreciate any feedback on my code

The fraction 49/98 is a curious fraction, as an inexperienced mathematician in attempting to simplify it may incorrectly believe that 49/98 = 4/8, which is correct, is obtained by cancelling the 9s.

We shall consider fractions like, 30/50 = 3/5, to be trivial examples.

There are exactly four non-trivial examples of this type of fraction, less than one in value, and containing two digits in the numerator and denominator.

If the product of these four fractions is given in its lowest common terms, find the value of the denominator.

It took me a couple of tries to get this, I think as the question doesn’t give enough details about what the rules are for anomalous fraction cancellation.

Supposing we say our fractions are of the form ab / cd.  Initially I was checking each of the following:

• a / c = ab / cd
• a / d = ab / cd
• b / c = ab / cd
• b / d = ab / cd

i.e. each possible combination of the first and second numerator and denominator digits.

Then I realised I could omit two of these, as we didn’t want to look at those cases where it’s the first numerator digit divided by the first denominator digit, or the second numerator digit divided by the second denominator digit.  Reducing what I was checking down to:

• a / d = ab / cd
• b / c = ab / cd

So I was now just looking at the first numerator over the second denominator, and the second numerator over the first denominator.  However this was still giving me too much possible fractions.  Further reading illustrated that the example given, 49 / 98, where the second numerator digit and the first denomator digit are cancelled out, is the rule for all cases.  Reducing what I was checking to just:

• a / d = ab / cd

The question also doesn’t really explain what else you can ignore.  Firstly if either digit is zero.  Secondly, the trivial examples include where a = b and c = d, e.g. 11 / 22.

And finally, I thought you could look at any values for ab and cd, e.g. I’d got 16 / 96 = 1 / 6.  But this meant I’d still got more than four ‘matching’ fractions.  What the question failed to mention is that the digits being cancelled out had to be identical, so really what I was checking was changed to:

• a / c = ab / bc

Once I’d got all that cleared up, it was easy.

```<?php
\$product = 1;

for (\$numerator = 10; \$numerator <= 99; \$numerator++) {
for (\$denominator = \$numerator+1; \$denominator <= 99; \$denominator++) {
\$fraction = \$numerator / \$denominator;

\$numeratorAsString = (string)\$numerator;
\$denominatorAsString = (string)\$denominator;

if (
!isTrivial(\$numeratorAsString, \$denominatorAsString) &&
\$numeratorAsString[1] == \$denominatorAsString[0] &&
\$numeratorAsString[0] / \$denominatorAsString[1] == \$fraction
) {
\$product *= \$fraction;
}
}
}

function isTrivial(\$numerator, \$denominator) {
if (\$denominator[1] == 0) {
return true;
}

if (\$numerator[0] == \$numerator[1]) {
return true;
}

return false;
}

echo \$product;```

I cast my numerators and denominators to strings, enabling me to then reference the digits within each using array notation.

## December 10, 2014

### Project Euler: problem 27 (PHP) – Quadratic primes

Filed under: PHP,Project Euler — duncan @ 8:00 am
Tags: , ,

Photo by Ianqui Doodle

I’m doing these Project Euler mathematical puzzles as a simple practical exercise for teaching myself PHP, and I’d appreciate any feedback on my code.

Euler discovered the remarkable quadratic formula:

n² + n + 41

It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 39. However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41, and certainly when n = 41, 41² + 41 + 41 is clearly divisible by 41.

The incredible formula  n² − 79n + 1601 was discovered, which produces 80 primes for the consecutive values n = 0 to 79. The product of the coefficients, −79 and 1601, is −126479.

n² + an + b, where |a| < 1000 and |b| < 1000
where |n| is the modulus/absolute value of n
e.g. |11| = 11 and |−4| = 4

Find the product of the coefficients, a and b, for the quadratic expression that produces the maximum number of primes for consecutive values of n, starting with n = 0.

This was another puzzle I’d ignored previously, then looked at it again and realised it wasn’t that hard after all.  First time I’ve really done anything with quadratics since school probably!

```<?php
include 'isPrime.php';

\$maxPrimes = 0;
\$maxCoefficients = [];

calculatePrimes(1, 1);
calculatePrimes(1, -1);
calculatePrimes(-1, 1);
calculatePrimes(-1, -1);

function calculatePrimes(\$incrementA, \$incrementB) {
global \$maxPrimes;
\$a = 0;
while (abs(\$a) < 1000) {
\$b = 0;
while (abs(\$b) < 1000) {
\$n = 0;

while (true) {
\$q = (\$n * \$n) + (\$n * \$a) + \$b;

if (!isPrime(\$q)) {
if (\$n > \$maxPrimes) {
\$maxPrimes = \$n;
setMaxCoefficients(\$a, \$b);
}
break;
}

\$n++;
}
\$b+= \$incrementB;
}
\$a+= \$incrementA;
}
}

function setMaxCoefficients(\$a, \$b) {
global \$maxCoefficients;
\$maxCoefficients = [\$a, \$b];
}

echo \$maxCoefficients[0] * \$maxCoefficients[1];```

So I’ve wrapped up most of the code into one function, which I call four times.  Each time I’m either adding or subtracting 1 from each of the two coefficients.  Then there’s nested loops so we’re going from zero to +/- 999 for both coefficients.  Within that we loop again, incrementing \$n until our quadratic equation doesn’t return a prime number.  At that point, we check if the value of \$n is more than the previous maximum number of primes.  And if it is, I update a global array with those coefficients.

After we finish looping, I then calculate the production from those stored coefficients.  Simple, but not particularly fast.

## December 3, 2014

### Project Euler: problem 32 (PHP) – Pandigital products

Filed under: PHP,Project Euler — duncan @ 1:24 pm
Tags:

Photo by Lars Dahlin

I’m doing these Project Euler mathematical puzzles as a simple practical exercise for teaching myself PHP, and I’d appreciate any feedback on my code.

We shall say that an n-digit number is pandigital if it makes use of all the digits 1 to n exactly once; for example, the 5-digit number, 15234, is 1 through 5 pandigital.

The product 7254 is unusual, as the identity, 39 × 186 = 7254, containing multiplicand, multiplier, and product is 1 through 9 pandigital.

Find the sum of all products whose multiplicand/multiplier/product identity can be written as a 1 through 9 pandigital.

HINT: Some products can be obtained in more than one way so be sure to only include it once in your sum.

This if the first one of these Project Euler puzzles I’ve completed in the last month.  I’d been sidelined by a much harder one, given up on it for now, and tried this instead for the first time.  It wasn’t too hard either; I’m not sure why I hadn’t attempted it previously.

```<?php
\$sums = [];

for (\$i = 1; \$i < 9; \$i++) {
for (\$j = 1234; \$j < 9876; \$j++) {
getPandigitalProduct(\$i, \$j, \$sums);
}
}

for (\$i = 12; \$i < 98; \$i++) {
for (\$j = 123; \$j < 987; \$j++) {
getPandigitalProduct(\$i, \$j, \$sums);
}
}

echo array_sum(\$sums);

function getPandigitalProduct(\$multiplicand, \$multiplier, &\$sums) {
\$product = \$multiplicand * \$multiplier;
\$number =  \$multiplicand . \$multiplier . \$product;

if (isPandigital(\$number)) {
\$sums[\$product] = \$product;
}
}

function isPandigital(\$number) {
\$length = strlen(\$number);

if (\$length != 9) {
return false;
}

for (\$i = 1; \$i <= \$length; \$i++) {
if (strpos(\$number, (string)\$i) === false) {
return false;
}
}

return true;
}
```

Not the cleanest code, but it’s reasonably fast.  In the isPandigital function I just loop from 1 to 9, checking that each of those digits is present in the number (after checking that it’s 9 digits long of course).

The important thing to understand is why I’m doing two separate nested loops.  The first one is for getting 9-digit identities which look like:

1 * 2345 = 6789

The second one is for getting identities which look like:

12 * 345 = 6789

Initially I had it in one set of nested loops, like:

```for (\$i = 1; \$i < 99; \$i++) {
for (\$j = 1; \$j < 9999; \$j++) {
...
}
}```

However that made 979804 iterations!  Doing it this way there’s a total of 143440 iterations (still way too many of course).

Also note I’m explicitly passing the \$sums array by reference by specifying the ampersand before the argument in the function declaration.

What else… I had to cast the value of \$i to a string before being able to use it succesfully with the strpos function.

## September 19, 2014

### Project Euler redux

Filed under: PHP,Project Euler — duncan @ 8:00 am
Tags: ,

A few years ago I did several of the Project Euler puzzles, a series of mathematical problems that are good programming exercises.  I started out doing them in ColdFusion, but ended up doing a a few in Python due to its better abilities at dealing with massive numbers and calculations.

I’ve now decided to convert at least some of them from CFML / Python code into PHP, as a little tutorial for myself to increase my familiarity with some of its basic syntax.  I’m doing this for basically the same reasons my colleague Adam Cameron is shifting his blog’s focus from ColdFusion to PHP.

I’m not going to try and solve the problems that I’ve already done by looking at them completely from scratch.  Instead I’ll assume my logic or approach from before was basically sound, and try to simply convert the code into PHP, perhaps with minor improvements.

What I’m  hoping is I’ll get used to common PHP functions and perhaps some gotchas that a newcomer to the language wouldn’t automatically be aware of.  If you happen to be reading any of these, have some experience in the language and spot anything daft I’m doing or know better approaches and so on, I’d appreciate any feedback.

So without further ado, here’s the first Project Euler puzzle, a load more will hopefully follow.

## December 13, 2009

### Project Euler problem 58

Filed under: Coldfusion,Project Euler — duncan @ 9:00 am
Tags: , , ,

Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.

```37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49
```

It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?

This problem has some similarities to problem 28. In that case the spiral went clockwise; this one it goes anti-clockwise. For our purposes that’s a red herring.

We want to loop round, calculating the value for the corners of our spiral.
For each value, work out if it is prime (re-using the IsPrime function I first used back in problem 3).
Keep count of how many primes there are.
Also keep track of how long the side of our spiral is.
At each loop round the spiral, check if the ratio of primes / length is less than 10%.

```<cfset numPrimes = 0>
<cfset length = 1>
<cfset number = 1>
<cfset increment = 2>
<cfset ratio = 100>

<cfloop condition="ratio GT 10">
<cfset length += 2>

<!--- add our increment 4 times, calculating each time if it's prime --->
<cfloop index="i" from="1" to="4">
<cfset number += increment>

<!--- the 4th corner is never prime --->
<cfif i LT 4 AND IsPrime(number)>
<cfset numPrimes++>
</cfif>
</cfloop>

<cfset ratio = numPrimes / ((2*length)-1) * 100>

<cfset increment += 2>
</cfloop>

<cfoutput>#length#</cfoutput>

<cfscript>
function isPrime(x)
{
var isPrime = true;
var i = 0;

if (x LT 2)
{
return false;
}

if ((NOT x MOD 2) AND (x GT 2))
{	// a multiple of 2, but not 2 itself
return false;
}

for (i = 3; i LTE SQR(x); i = i + 2)
{
if (NOT x MOD i)
{	// found a factor of x
return false;
}
}

return isPrime;
}
</cfscript>```

Like in problem 28, we can see the corner values go:
3, 5, 7, 9, [increment: 2]
13, 17, 21, 25, [increment: 4]
31, 37, 43, 49 [increment: 6]
etc.

So we can work out the corner values by starting at 1, adding 2 to our increment each time round the spiral, then adding the increment four times to calculate the corners.

We don’t need to work out if the 4th corner is prime; it will always be a square number. e.g. on spiral with length 3, the 4th corner is 32 = 9. On length 5, the 4th corner is 52 = 25, and so on. So the if statement <cfif i LT 4 AND IsPrime(number)> will only execute the IsPrime() function if i LT 4 evaluates to true.

## December 12, 2009

### Project Euler problem 44

Filed under: Project Euler,Python — duncan @ 3:08 pm
Tags: , ,

Pentagonal numbers are generated by the formula, Pn=n(3n−1)/2. The first ten pentagonal numbers are:

1, 5, 12, 22, 35, 51, 70, 92, 117, 145, …

It can be seen that P4 + P7 = 22 + 70 = 92 = P8. However, their difference, 70 − 22 = 48, is not pentagonal.

Find the pair of pentagonal numbers, Pj and Pk, for which their sum and difference is pentagonal and D = |Pk − Pj| is minimised; what is the value of D?

I’d originally made a start on this in Python, then gave up on it. Came back to it much later and tried doing it in ColdFusion. My solution there took too long to run, so I went back to my Python version and finished it off.

Here’s the code:

```import math

def pentagonal(x):
return x * (3 * x - 1) / 2

def isPentagonal(x):
n = (math.sqrt((24 * x) + 1) + 1) / 6

if n == int(n) and n > 0:
return 1
else:
return 0

i = 0
pentagonals = []
solution = 0

while 1:
i += 1
p = pentagonal(i)

for j in pentagonals:
diff = abs(p - j)

if isPentagonal(diff) == 1:
intSum = p + j

if isPentagonal(intSum) == 1:
print("solution:", int(diff))
solution = 1
break

if solution:
break

pentagonals.append(p)
```

The formula for calculating the pentagonal numbers is given to us. I already had a formula for the reverse, checking if a number is pentagonal, from problem 45.
To calculate the square root, I needed to import the math library.

Basically I loop until I find a solution. I calculate the pentagonal value of i. Then I subtract each pentagonal number smaller from it (but use the abs() function to get the absolute value).
If that difference is pentagonal, then I do a sum as well.
If that is also pentagonal, then our solution is the difference.
As I go along I append each pentagonal number into the array.
This solution runs in about 12 seconds for me.

## December 9, 2009

### Project Euler problem 55

Filed under: Project Euler,Python — duncan @ 12:00 am
Tags: , , , ,

Problem 55:

If we take 47, reverse and add, 47 + 74 = 121, which is palindromic.

Not all numbers produce palindromes so quickly. For example,

349 + 943 = 1292,
1292 + 2921 = 4213
4213 + 3124 = 7337

That is, 349 took three iterations to arrive at a palindrome.

Although no one has proved it yet, it is thought that some numbers, like 196, never produce a palindrome. A number that never forms a palindrome through the reverse and add process is called a Lychrel number. Due to the theoretical nature of these numbers, and for the purpose of this problem, we shall assume that a number is Lychrel until proven otherwise. In addition you are given that for every number below ten-thousand, it will either (i) become a palindrome in less than fifty iterations, or, (ii) no one, with all the computing power that exists, has managed so far to map it to a palindrome. In fact, 10677 is the first number to be shown to require over fifty iterations before producing a palindrome: 4668731596684224866951378664 (53 iterations, 28-digits).

Surprisingly, there are palindromic numbers that are themselves Lychrel numbers; the first example is 4994.

How many Lychrel numbers are there below ten-thousand?

Yet another problem which was easy to code in ColdFusion, but it couldn’t deal with the large numbers generated. So rewrote it in Python.

Just for interest, here’s what I did in CFML. NB: there’s a logic error here, which I fixed once I moved to Python:

```<cfscript>
function isPalindrome(x)
{
if (x EQ Reverse(x))
return true;
else
return false;
}
</cfscript>

<cfset Lychrel = ArrayNew(1)>

<cfloop index="i" from="1" to="9999">
<cfset k = i>

<cfloop index="j" from="1" to="50">

<cfset temp = k + Reverse(k)>

<cfif IsPalindrome(temp)>
<cfoutput>#k# + #Reverse(k)# = #temp#</cfoutput><br>
<cfset ArrayAppend(Lychrel, i)>
<cfbreak>
</cfif>

<cfset k = temp>
</cfloop>
</cfloop>

<cfoutput><strong>#ArrayLen(Lychrel)#</strong></cfoutput>
```

Running this it quickly throws an error:
‘The value 210+E11200002108.1 cannot be converted to a number.’

Python doesn’t seem to have a built-in Reverse function like ColdFusion, so you have to write your own.

```def Reverse(x):
num = str(x)
newnum = ''

for i in range(len(num)-1, 0, -1):
newnum += num[i]

newnum += num[0]

return newnum

Lychrel = []
limit = 10000

for i in range(1, limit):
k = i

for j in range(50):
temp = k + int(Reverse(k))

if str(temp) == str(Reverse(temp)):
break

k = temp

if j == 49:
Lychrel.append(i)

print(len(Lychrel))
```

This code is pretty fast, but it could be better. For instance I’m casting everything to a string or an int as required, but that could maybe be tidied up somehow.
Also in my Reverse function, I somehow couldn’t work out how to properly loop backwards through the string, so ended up having to add on the first character outside of the loop.

Initially I was appending to my Lychrel array at the point where I break out of the loop. However on reading the question again I realised the Lychrel numbers are the numbers which aren’t solved within 50 iterations, not the other way around. I don’t even need an array, I could just have a counter that I increment. But having it in the array gives us the advantage of being able to output the array contents to make sure we’re on the right track.

## February 6, 2009

### Project Euler problem 56

Filed under: Project Euler,Python — duncan @ 12:00 am
Tags: , , ,

Problem 56:

A googol (10100) is a massive number: one followed by one-hundred zeros; 100100 is almost unimaginably large: one followed by two-hundred zeros. Despite their size, the sum of the digits in each number is only 1.

Considering natural numbers of the form, ab, where a, b < 100, what is the maximum digital sum?

So, let’s do nested loops, both going 1 to 100. Inside the inner loop, calculate a^b (or i^j in this case). Then loop through the digits of that value, adding them up. Keep track of which is the largest value this produces.

Running time about 3 seconds in Python:

```import time
tStart = time.time()

maxsum = 0

for i in range(100):
for j in range(100):
num = i ** j
total = 0

for k in str(num):
total += int(k)

if total > maxsum:
maxsum = total

print(maxsum)
print("time:" + str(time.time() - tStart))
```

## February 4, 2009

### Project Euler problem 53

Filed under: Project Euler,Python — duncan @ 12:00 am
Tags: , ,

Problem 53:

There are exactly ten ways of selecting three from five, 12345:

123, 124, 125, 134, 135, 145, 234, 235, 245, and 345

In combinatorics, we use the notation, 5C3 = 10.

In general,
nCr = n! / (r!(n-r)!)
where r <= n, n! = n*(n-1)*…*3*2*1, and 0! = 1.

It is not until n = 23, that a value exceeds one-million: 23C10 = 1144066.

How many, not necessarily distinct, values of nCr, for 1 <= n <= 100, are greater than one-million?

At first I looked at this and thought it was all about iterating through the combinations of digits, just like how it lists 123, 124, etc. Which sounds tricky. But it’s much simpler than that. In fact they give you the hardest part of the algorithm, the equation needed:
n! / (r!(n-r)!)

In Python this solution takes about 0.36 seconds. I’m going to start using the time module to calculate running times.

```import time
tStart = time.time()

values = 0

def factorial(x):
y = 1
for i in range(1,x+1):
y = y * i

return y

def calc(n, r):
return  factorial(n) / (factorial(r) * factorial(n-r))

for i in range(1,101):
for j in range(1, i):
if calc(i, j) > 1000000:
values += 1

print(values)
print("time:" + str(time.time() - tStart))
```

So, two functions, one calculates the factorial of x, the other just calculates our equation. The second function could instead have just been expressed as one line in our procedural code.

Then loop from 1 to 100. Within that, do an inner loop from 1 to whatever our current outer loop is on, passing the two loop counters into our function. We don’t really care which numbers produce values greater than a million; we only care how many numbers there are. So value is just a simple counter.

Next Page »

Create a free website or blog at WordPress.com.