01 Oct 2009
Yet another problem involving primes. Guess I have to make a better prime generator soon…
The sum of the primes below 10 is 2 3 5 7 = 17.
Find the sum of all the primes below two million.
Soon, but not yet…
Using the same generator from problem 7, I go for an easy, but slow, solution for this one. There is one thing, though. Our generator generates up to the nth prime, but in this problem, we need to find every prime below 2000000. Already decided not to make a new prime generator, we cheat by asking it generate a lot of primes, and then jump out as soon as we hit a prime to big.
sum := 0;
for prime in PrimeGenerator(500000) do
begin
if prime >= 2000000 then
break;
sum := sum + prime;
end;
I choose to generate 500000 primes, pretty sure one prime out of four should be enough.
And that’s it – it’s all there at http://svn.vi-kan.net/euler.
30 Sep 2009
On our way to Lofoten last weekend, my wife and I had a small brainstorming on problem no. 9:
A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,
a2 b2 = c2
For example, 32 42 = 9 16 = 25 = 52.
There exists exactly one Pythagorean triplet for which a b c = 1000.
Find the product abc.
Both of us agreed there had to be some kind of mathematical formula for this kind of problem. The problem, though, was that non of us new this formula, and there wasn’t that many math books in the car that day… So instead of trying to guess the magic formula, we tried to optimize the search for the right numbers.
We could put three for loops inside each other, scanning the numbers from 1 to 1000 for each. Trying 1 1 1 as a possible solution isn’t that clever, though. Neither is the combination 1000 1000 1000.
Let’s take a look at the maximum possible values for a, b and c. The fact that a < b < c makes this kind of easy. Since all three numbers must be positive, c can’t be larger than 997 to make room for a = 1 and b = 2. The number b can’t be larger than 499, which makes room for a = 1 and c = 500. And finally, a can’t be larger than 332, which makes room for b = 333 and c = 335.
for a := 1 to 332 do
for b := a 1 to 498 do
for c := b 1 to 998 do
begin
...
end;
Now, the test it self should be easy:
if (a + b + c = 1000) and (a*a + b*b = c*c) then
...
We can do one more easy optimization. I our current c makes for a larger sum than 1000, we can break out of the innermost loop. This actually saves us a lot of cycles.
for c := b + 1 to 998 do
begin
sum := a + b + c;
if (sum) > 1000 then
break;
if (sum = 1000) and (a*a + b*b = c*c) then
begin
...
end;
end;
And that’s it. The complete solution is available as usual. Check it out at http://svn.vi-kan.net/euler.
30 Sep 2009
A couple of days ago, I solved eulers problem no. 8. I finally found some time to blog about my solution.
Find the greatest product of five consecutive digits in the 1000-digit number.
73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
I can’t see any smart solution for this problem. A simple brute force scan through the series of digits will have to do.
I start out with the number as a string. That makes it easy to loop through digit by digit. To make it even easier, I have made an array of five bytes that I place ‘over’ the current character. This gives me a ‘view’ of the five current digits that needs to be multiplied as bytes instead of characters.
const
SUBJECT = '73167176531330624919225119674426574742355349 .... 57530420752963450';
type
PFiveDigits = ^TFiveDigits;
TFiveDigits = array[0..4] of byte;
begin
for I := 1 to length(Subject) - 6 do
begin
FiveDigits := @SUBJECT[i];
end;
end;
A character in the range 0 to 9 has byte value in the range of 48 to 57. Before we multiply the values, we have to normalize them. We can do that by subtracting 48 from each byte value.
const
NULLCHAR = ord('0');
...
present := (FiveDigits[0]-NULLCHAR) * (FiveDigits[1]-NULLCHAR) * (FiveDigits[2]-NULLCHAR) * (FiveDigits[3]-NULLCHAR) * (FiveDigits[4]-NULLCHAR);
Now we just have to keep track of the highest produced product do find the answer.
largest := 0;
for I := 1 to length(Subject) - 6 do
begin
FiveDigits := @SUBJECT[i];
present := (FiveDigits[0]-NULLCHAR) * (FiveDigits[1]-NULLCHAR) * (FiveDigits[2]-NULLCHAR) * (FiveDigits[3]-NULLCHAR) * (FiveDigits[4]-NULLCHAR);
if present > largest then
begin
largest := present;
position := i;
end;
end;
As usual, you will find the code at http://svn.vi-kan.net/euler.
25 Sep 2009
It want be easy to keep up with my brother on this. Well, here’s my take on Project Eulers problem no 7:
By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10001st prime number?
Algorithms conserning primes is a large and complex field. I guess the are many different solutions to this problem, so the question will be how efficient the one that I find will be.
One thing is for sure, though. We need a way of finding primes. As far as I know, there is no formula to find the nth prime. We have to find the first 10 000 primes first before we can find our wanted one.
A Prime Generator
I will use a enumerator-pattern for our prime generator. For that, we will need two classes. First we will need a factory-class that implements the GetEnumerator function. Then we will need the class that handles the actual enumeration.
type
TPrimeGenerator = class
public
function GetCurrent: extended;
function MoveNext: boolean;
property Current: extended read GetCurrent;
end;
TPrimeGeneratorFactory = class
public
function GetEnumerator: TPrimeGenerator;
end;
With this implementation, we will be allowed to write code like this:
for prime in TPrimeGeneratorFactory.Create do
begin
...
end;
This code would generate an awful lot of primes, though. It would be nice to give it some kind limit. We can add that to the constructor, and check in the MoveNext method if the limit is reached. We would also need some live-time management on both the factory and the generator. The generator it self can be created and destroyed by the factory, but we can also make it easier by letting the factory implement an interface. If we do, delphi will manage it’s lifetime for us.
type
TPrimeGenerator = class
public
constructor Create(numberOfPrimes: integer);
...
end;
IPrimeGeneratorFactory = interface
function GetEnumerator: TPrimeGenerator;
end;
TPrimeGeneratorFactory = class(TInterfacedObject, IPrimeGeneratorFactory)
public
constructor Create(numberOfPrimes: integer);
destructor Destroy; override;
end;
function PrimeGenerator(numberOfPrimes: integer): IPrimeGeneratorFactory;
Now, we can rewrite our for … in – loop, and just get the number of primes that we want. We will know longer have any leaks either.
for prime in PrimeGenerator(10001) do
begin
...
end;
Finding Primes
Our generator doesn’t actually generate primes yet, though. Let’s see if we can make it do so. So how do we find primes? What is a prime? Wikipedia gives us the following definition of prime numbers:
In mathematics, a prime number (or a prime) is a natural number which has exactly two distinct natural number divisors: 1 and itself.
So if we loop through all numbers, and test if it’s possible to devide it by any other number, we should have our primes. The problem is, though, that this takes time. A lot of time. We have to find a way of shortening down the ammount of numbers we need to check. The first thing we should remove, is all even numbers. Every single even number is divisible by the number 2, so they can’t be primes. Except for the number 2 it self – it’s a prime, cause has only two divisors: 1 and itself.
Since 2 is the first prime number, it would be an easy exception to handle, and we have effectively cut shortened down possible primes to the half of what we started with.
Let’s write our FindNextPrime method:
function TPrimeGenerator.FindNextPrime(PreviousPrimeFound: extended): extended;
var
found: boolean;
currentNumber: extended;
begin
currentNumber := PreviousPrimeFound;
repeat
if (CurrentNumber = 1) or (currentNumber = 2) then
begin
currentNumber := currentNumber 1;
result := true;
end
else
begin
currentNumber := currentNumber 2;
found := CheckIfPrime(currentNumber);
end;
until found;
result := currentNumber;
end;
We here assume that PreviousPrimeFound equals 1 when no previous prime is found. So if the previous prime is 1 or 2, we just increment by one, knowing that both 2 and 3 is prime numbers, and signal a successfull found. If previous prime was not 1 or 2, it has to be 3 or any other odd number, so we increment by two to skip even numbers.
Let’s take a look at the CheckIfPrime( ) method. How do we check if a number is divisible by any other number? I guess we have to check:
function TPrimeGenerator.CheckIfPrime(number: integer): boolean;
var
i: integer;
begin
result := true;
for i := 2 to number - 1 do
begin
if number mod i = 0 then
begin
result := false;
break;
end;
end;
end;
Can you see a bottleneck here? For every number we test, there is one more number to devide. How can we optimize this?
First of all, every natural number can be expressed as a product of primes. 12 = 2 x 2 x 3, 20 = 2 x 2 x 5 and so on. So if a number is divisible by a non prime number, it should also be divisible by a prime number. So if we can find a way of just dividing by primes we would speed up things a lot.
So let’s add a list of previous found primes to the generator class, and loop through that list instead of all numbers.
function TPrimeGenerator.CheckIfPrime(number: extended): boolean;
var
itr: DIterator;
prevPrime: extended;
found: boolean;
begin
found := false;
itr := FPrimes.start;
while IterateOver(itr) do
begin
prevPrime := getExtended(itr);
if Frac(number / prevPrime) = 0 then
begin
found := true;
break;
end;
end;
result := not found;
end;
This makes things a lot faster, but I think we still can gain something. In the same wikipedia article that I refered to earlier, it is sufficient to test primes smaller than the sqaure root of the number we are currently testing. So let’s add that to our loop:
function TPrimeGenerator.CheckIfPrime(number: extended): boolean;
var
itr: DIterator;
prevPrime: extended;
found: boolean;
root: extended;
begin
root := sqrt(number);
...
while IterateOver(itr) do
begin
prevPrime := getExtended(itr);
if prevPRime > root then
break;
if ...
end;
...
end;
The Solution
Just like the solution for problem 6, you can find the complete code at http://svn.vi-kan.net/euler. It gives the answer in about 40–50 ms.
22 Sep 2009
Inspired by my brother over at geekality.net, I thought I should do an attempt on the various problems presented at projecteuler.net/. Since he already solved the first five, I jumped right in at problem 6:
The sum of the squares of the first ten natural numbers is,
12 22 … 102 = 385
The square of the sum of the first ten natural numbers is,
(1 2 … 10)2 = 552 = 3025
Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is
3025 385 = 2640.
Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.
This problem has two distinct parts. First, there is the sum of a sequence of the squared natural numbers, then there is the sum of a sequence of natural numbers, squared.
Sum of a series of squared natural numbers
The problem tells to summarize the square of the first hundred natural numbers. That gives us the following sequence:
12 22 … 992 1002
The simplest solution would be a brute force loop like this one:
function sumSquaredNumbers( ): longword;
var
i: integer;
begin
result := 0;
for i := 1 to 100 do
result := result (power(i, 2));
end;
That wouldn’t be much fun, though. There has to be a more general, optimized algorithm for this kind of work. First of all, we should make the function workable for any series length. That should be easy – just switch the hardcoded 100 with an given parameter:
function sumSquaredNumbers(serieslength: integer): longword;
var i: integer;
begin
for i := 1 to seriesLength do
result := …
end;
But the code is still the same, though. Nothing was optimized in any way. We still have the loop, making our function linear.
Google to the rescue!
math.com have a nice section on Series Expansions, which includes a table of power summations. This table states that
n
∑
k2
= 1 4 9 … n2
= (1/3)n3 (1/2)n2 (1/6)n
k=1
Why this should be true, I honestly can’t tell, but a fast comparison to our brute force solution shows that it gives us the right numbers. That’s good enough for now:
function sumSquaredNumbers(seriesLength: integer): longword;
begin
result := trunc(
((1/3) * power(seriesLength, 3))
((1/2) * power(seriesLength, 2))
((1/6) * seriesLength)
);
end;
I think this is a nice solution, and a constant one as well. It doesn’t matter if the series is of length one, hundred or thousand – the computation should take the same amount of time.
This part should be a lot easier. We could start with a brute force solution to this part too, but I think we could come up with something smarter right away.
I already know of a simple optimization. If you add the first and last number in the series, you will get the same as if you add the second and next to last number:
1
2
3
4
5
6
7
8
9
10
10
9
8
7
6
5
4
3
2
1
=
11
11
11
11
11
11
11
11
11
11
As you can see, every pair of numbers makes the same sum when added together. So for any given series, you can take the first and the last number and add them together, then you can multiply the sum by the length of the series, and finally divide by two.
n
∑
k
= 1 2 3 … n
= (1 n) * n / 2
k=1
I know there are other solutions to this, but this one will do.
function sumSquared(seriesLength: longword): longword;
begin
result := trunc(power(((1 seriesLength) * seriesLength / 2), 2));
end;
The solution
Now that we have solved each of the two sub parts of the problem, I guess we’re ready to solve the main thing.
Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.
Since we know the sum of the squares, and we now the square of the sum, it’s just a matter of substract the one from the other:
function diff(serieslength: longword): longword;
begin
result := sumSquared(serieslength) - sumSquaredNumbers(serieslength);
end;
I have uploaded my complete solution to http://svn.vi-kan.net/euler. The code also contains a simple stopwatch to time the code. This code is so simple, though, that it’s hard to get any measures.
And that’s it. Let’s see if we can solve the next one too.
