## Monday, 31 March 2008

### Project Euler Problem 3

This post tackles Problem 3 in the Project Euler series.

Do you remember making Factor Trees in maths lessons? They are a neat way of showing how a composite number (a number that is not prime) is composed of prime factors. Suppose we wanted to factor the number 60. We might start by noticing that it is divisible by 2: 2 x 30 = 60. So we write

```        60
/  \
2   30
```

2 is prime, so we're finished with that branch of the tree. Look at the other branch: we've now got to solve the same problem, but with a smaller number. 30 is also divisible by 2, so we add another layer to the tree:

```        60
/  \
2   30
/  \
2    15
```

Applying the same procedure to 15 gives the final tree:

```        60
/  \
2   30
/  \
2    15
/  \
3    5```

The procedure stops when every branch of the tree has a prime number hanging from it. I'm not one for fortune telling, but reading these leaves will tell you something useful: the prime factorisation of the root number. In this case 60 = 2 x 2 x 3 x 5.

Problem 3 wants us to find the largest prime factor of a composite number. The Factor Trees methods suggests a nice recursive algorithm, and recursion has a very Functional feel to it. The algorithm is quite straight-forward.

We start with the number that we want to factorise, and feeling friendly, we call it n. We then try to find the first number bigger than 1 which divides n: call this number p. p must be prime: if it wasn't we'd have found another number first that divides n. We next cut n down to size by dividing it by p as may times as we can. Now for the recursive part: we start the procedure over again, this time with the smaller number that we found when dividing n by p. Any prime factor of this smaller number will be larger than p, and is also a factor of n so this trick will give us the answer we need.

There are just two special cases: first, there might not be a number p that divides n. In that case, n must be prime, so n is the answer we are looking for. Second, it might be that when we have divided n by p, we end up with 1, which won't factor any further: but in this case p is the largest prime that divides n, so p is the answer we need.

There's one other thing to note: when we are hunting for p we don't need to try all the numbers as far as n. We only need to go as far as √n. This is because n = √n x √n, so if it has a factor greater than √n, it will have a corresponding factor smaller than √n which we will find by searching only as far as √n.

We could be even cleverer and say that we only need to try dividing n by odd numbers (because if it is divisible by an even number, it will be divisible by 2, and we're chopping out 2 as a factor right at the beginning of the search). I want to keep the code simplish though, so for now I won't bother.

Now for some code:

```public void Solve()
{
long factor = FindGreatestPrimeFactor(1, 600851475143);
Console.WriteLine(factor);
}

public long FindGreatestPrimeFactor(long factorGreaterThan, long number)
{
long upperBound = (long)Math.Ceiling(Math.Sqrt(number));

// find next factor of number
long nextFactor = Range(factorGreaterThan + 1, upperBound)
.SkipWhile(x => number % x > 0).FirstOrDefault();

// if no other factor was found, then the number must be prime
if (nextFactor == 0)
{
return number;
}
else
{
// find the multiplicity of the factor
long multiplicity = Enumerable.Range(1, Int32.MaxValue)
.TakeWhile(x => number % (long)Math.Pow(nextFactor, x) == 0)
.Last();

long quotient = number / (long)Math.Pow(nextFactor, multiplicity);

if (quotient == 1)
{
return nextFactor;
}
else
{
return FindGreatestPrimeFactor(nextFactor, quotient);
}
}
}

private IEnumerable<long> Range(long first, long last)
{
for (long i = first; i <= last; i++)
{
yield return i;
}
}```

A few words of explanation might be called for.

First, you'll recognise the Range method as being an iterator. I created this one, rather than using Enumerable.Range primarily because I needed to be able to specify first and last as long rather than int. I also took the opportunity to make it more convenient to use by making it take an last parameter, rather than a count as Enumerable.Range does.

In order to find the next factor, we generate a sequence of all the numbers greater than any factor we found previously, up to √number. SkipWhile is another Extension method that works on a sequence (it kind of does the opposite of TakeWhile which we used last time). It ignores the first part of a sequence, and only includes the part after the condition stops being met. In this case, it skips all numbers that don't divide number, and then returns the rest of the input sequence. FirstOrDefault(), as you might expect from its name takes a sequence, and returns the first element of it; or if there are no elements in the sequence, it returns a default value - the default for the Type of element in the sequence, in this case 0.

Finding the multiplicity of the prime factor is similarly straightforward (multiplicity in this case means, how many times the factor divides into number). We generate the sequence of numbers, m, where nextFactor^m divides number, and then take the last one in that sequence.

And there we have it. Not as compact a solution as for previous problems, but interesting non-the-less. Anonymous said...

Here's a more compact (more naive) solution based on the LINQ-implemented prime number generator from: http://jacobcarpenter.wordpress.com/2008/03/26/linq-to-prime-numbers/

var odds =
from n in Enumerable.Range(0, int.MaxValue)
select 3 + (long) n * 2;

var primes = (new[] { 2L }).Concat(
from p in odds
where ! odds.TakeWhile(odd => odd * odd <= p).Any(odd => p % odd == 0)
select p);

var num = 600851475143;
Console.WriteLine(primes.TakeWhile(p => p * p <= num).Where(p => num % p == 0).Last()); Anonymous said...

Whoops! Apparently, I needed to manually tag those urls as hyperlinks:

- LINQ prime number generator

Unknown said...

Jacob,
I like it! I'd seen your post before when I was researching Euler problem 10.

I've got a slightly more efficient version of your prime generator which I'll post soon. Anonymous said...

Cool! Can't wait to see it.

Sævar said...

Hey - awesome site. I'm learning C# myself, and I had no idea about extensions and the lambda expressions, possibly due to looking at outdated material, but thanks to you I've seen the light.

Really enjoy your code; lots of nifty stuff which I didn't know C# could do.

My solution for Euler 3 was creating an extension to list all divisors and using Where(x => x.isPrime()).Max().. Of course, manually checking for primes using the loop until sqrt method isn't insanely efficient, but even so.

I'll keep on reading. Keep it coming.

Unknown said...

Sævar,
I'm glad you like my posts. Solving the problems in C#3.0 was a good learning experience for me too!

I've been pretty busy recently, but I will try to get more Project Euler posts up in the next month or so. Anonymous said...

Sieve's Algorithm worked best for me.
Make the code re-usable, seeing as problem 7 and onwards requires a lot of the work done in problem 3.

My hints and tips are here : http://www.ultralord.za.net/2010/05/14/project-euler-problem-3-solved/

Jacob said...

Cool! Can't wait to see it. Sujith487 said...

Same problem solved in C# language .

William Blackburn said...

I am having trouble printing the answer with Console.WriteLine.  I copied the code into VS and added the tried to call Solve from Main and I get an error: Error 1 An object reference is required for the non-static field, method, or property.
static void Main(string[] args)        {            Solve();        }        public void Solve()        {            long factor = FindGreatestPrimeFactor(1, 600851475143);            Console.WriteLine(factor);            Console.Read();        } Ujjwal Gulecha said...

sir i have solved this in 0.148 sec and the logic is very simple and here it is:
#include#includemain(){ long long j,k,count; long long n=600851475143LL; for(j=3;j<sqrt(n);j++) { if (n%j==0) { for(k=2;k<j;k++) { if (j%k==0) count=1; } if (count!=1) printf("%lld\n",j); } } } Jason said...

I have a javascript function that solved this in 0.016s :P