Euler Problem 97: Find last 10 digits of a big prime number

First, I'm glad that I've managed to solve 39 problems by now. Since the last problem, I don't find really interesting example to blog, but this problem 97 deserves a blog. Not really because Scala is good in answering this, but because the problem itself is interesting.


The problem is expressed in a little bit history of prime number, but finally the real problem is here: "What is the 10 last digits of the following prime number: 284332^7830457+ 1. Looks intimidating, right ? Try give a try calculating 28433 * (2: BigInt).pow(7830457) won't work, simply because its too large. Another thing is needed.

The other thing is the modular exponentiation, excellently explained here: http://en.wikipedia.org/wiki/Modular_exponentiation . Basically, the techniques are very useful to compute b^p mod m by observing that (a*b) mod m = ( (a mod m)  * (b mod m) ) mod m. 

Taking that into account,   b² mod m can be calculated as ( (b mod m) * (b mod m) )  mod m, and 

b³ mod m =( ( b mod m) * ( b² mod m) )mod m and so on. 

When this is translated to Scala it becomes:

 1 dev div(b: Long, exp: Long, m: BigInt): BigInt = {
 2   def div(b: Long, counter: Long, 
 3         m: BigInt, result: BigInt): BigInt = {
 4    if (counter == 0) 
 5       result
 6 
 7     else div(b, counter - 1, m, (result * b) % m)
 8   }
 9   div(b, exp, m, 1)
10 }

With this in our hand, we can then get the last 10 digit of 2^7830457 which nothing less than  2^7830457

mod 10000000000. Using div method above, we can get then the last 10 digits of the term. We can multiply the result with 28433 and then add 1 to get the answer to the problem.

OK. That was my 39th solved problem.

Problem: Euler Task 12: Number of factors of triangle number.

The following problem (euler 12) is a good example of the use of Stream collection in Scala. The use of Stream is explained in http://www.scala-lang.org/docu/files/collections-api/collections_14.html .
Basically, using Stream, we can create a list that computes its elements lazily. For example, to create a list of fibonacci number we can use 

def fibFrom(a: Int, b: Int): Stream[Int] = 

             a #:: fibFrom(b, a + b)

The problem is available here: http://projecteuler.net/index.php?section=problems&id=12 :

"The sequence of triangle numbers is generated by adding the natural numbers. So the 7^th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

 1: 1
 3: 1,3
 6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?"

The triangle number here is a good candidate for Stream collection.  We can define the triangle function as follow:

 1 def triangle(a: Int, b: Int) : Stream[Int] = 
 2    a #:: triangle(a + b + 1, b + 1)

Which means the first element of triangle number list is a and the next element is (a + b + 1), followed by ( a + b + 1 + b + 1 + 1), and so on .

triangle(1, 1).take(7).foreach(println) prints the first 7 triangle numbers.

The following function computes the number of factors of an integer:

1 def nfactor(n: Int): Int = 
2        (1 to csqrt(n)).foldLeft(0)( 
3         _ + incrfactor(n, _))

where incrfactor is a function that receives two parameters (n, x) and returns 0 if x does not divide n, 2 if x divides n and n/x != x , 1 otherwise (because when n/x = x ,we have only found one factor = x , otherwise if we have found 2: x and n/x).

1 def incrfactor(n: Int, x: Int) =
2     if (n % x == 0) (
3        if (x == n/x) 1 else 2

4     ) else 0    

With all these on our hands, we're ready to answer the problem. The answer to our problem is

triangle(1,1).find( nfactor(_) > 500)

which returns the first triangle number that has more than 500 factors. In my computer, the code above needs 6 s.

Problem: Euler Task 14: Which starting number, under one million, produces the longest chain?

I pick this problem 14 because it provides an interesting example of tuple used in combination with fold. Let's see how it works nicely.

First, the problem:

"The following iterative sequence is defined for the set of positive integers:

n n/2 (n is even)
n 3n + 1 (n is odd)

Using the rule above and starting with 13, we generate the following sequence:

13 40 20 10 5 16 8 4 2 1

It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.
Which starting number, under one million, produces the longest chain?"

Now, we're able to attack the problem. We will need a function that receives a term and returns the length of the chain. Let's say def chainLength(n : Int) : Int . Unfortunately, this will not work out since the term can be bigger than Integer.MAX_VALUE, which is the case when I first tested the problem. So, we'll need instead chainLength(n : Long).  We're ready for the code then.  We have basically four cases, when n <=0, n == 1,  n is odd, n is even. We'll have then something like :

def chainLength(n :Long): Int = {
   if (n <= 0) {
       throw new IllegalArgumentException("Negative term: " + n)
   } else if ( n == 1 ){
      1
   } else if (n% 2 == 0) {
      1 + chainLength(n/2)        
   } else {
     1 + chainLength(n *3  + 1)
   }
}

Almost OK. But the recursion is not tail recursion. We need to make it tail recursive by introducing an accumulator as the second parameter.

def chainLength(n :Long, len: Int): Int = {
   if (n <= 0) {
       throw new IllegalArgumentException("Negative term: " + n)
   } else if ( n == 1 ){
      len
   } else if (n% 2 == 0) {
           chainLength( n/2, len + 1)        
   } else {
           chainLength( n *3  + 1, len + 1)
   }
}

To get the chainLength for 13 we call then chainLength(13, 1)

Great. The only thing missing is now to find the lowest starting term between 2 to 1000000.
Here we're going to use combination of tuple (a, b) where a is the term and b is the length of sequence for the term. Here is the final code:

(2 to 1000000).foldLeft( (2,0) )( (x,y) => {
    val len = chainLength(y, 1)
    if (len > x._2 ) (y, len)
    else x
  })

The code took me less than 6 seconds in my machine.

Problem: Euler Task 25: What is the first term in the Fibonacci sequence to contain 1000 digits ?

Problem:
"What is the first term in the Fibonacci sequence to contain 1000 digits?"


Well, it should be only a variation of problem number 2. We need to modify the fibonacci function a little bit to stop when the number exceeds a value, for example 10⁹⁹.

The function becomes: 

1 def fibiter(max: BigInt) : BigInt = {
2   def fibiter(x: BigInt, acc: BigInt, index: Int) : BigInt = {
3       if (x > max)
4           index
5       else fibiter(x + acc,x, index + 1)
6    }
7    fibiter(1, 1, 2)
8 }

The line 3 checks if the last fibonacci number exceeds the maximum number. If it is, then its index is returned. Otherwise, the next fibonacci number is generated (x + acc) and the index is incremented.

The call starts with (1, 1, 2) to tell that the first two fibonacci numbers are already created. The solution to the problem is fibiter( (10: BigInt).pow 99)

Problem: Euler Task 16 - What is the sum of the digits of the number 2^1000?

"2¹⁵ = 32768 and the sum of its digits is 3 + 2 + 7 + 6 + 8 = 26.

What is the sum of the digits of the number 2¹⁰⁰⁰?"

One line:

 ((2 :BigInt) pow 1000).toString.foldLeft(0)(_ + _ - '0')

Problem: Euler Task 6 - Difference between Square of Sums and Sum of Squares

Well, task no 6 is pretty straight forward, so I write the solution quickly. Promise. I will stop coding this evening even when it is still raining  out there.

The problem is the following:

"The sum of the squares of the first ten natural numbers is,

1² + 2² + ... + 10² = 385

The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)² = 55² = 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."

First let's define square function. Easy:

def square(x : Int) = x * x

We define then the diff method that calculate the difference. The function receives a Range as its input. It is defined as follow:

def diff(r: Range) = 

   square( r.foldLeft(0)(_ + _)) - r.foldLeft(0)(_ + square(_))

The answer to the problem is then diff( 1 to 100).

While the exercise is easy, this is kind of problem where Scala, or functional programming in general, at its best: only three lines to solve this problem (I can even tweet the solution). That's why I find this exercise quite interesting.

Problem: Euler Task 3 - Largest Prime Factor

(Décidement, je suis en forme pour coder ce weekend. Après tout, avec le mauvais temps à Cote d'Azur ...)

Still from euler project. Now for task 3: largest prime factor. The integer factorization is an important subject and indeed it finds its application in cryptography. The fact that integer factorization is hard is used in many cryptography algorithms. But, to answer the question in euler project task 3, brute force approach should be enough since the integer to factorize is relatively small.

So, here we are. Here is the problem (http://projecteuler.net/index.php?section=problems&id=3)
"The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?"

First, the observation. The largest prime factor of 2 is 2. In general, I consider also that the largest prime factor of a prime number is the number itself. 

Then, we will try to check that the number is divided by a number. We check for every number from 2 until the square root of the checked number. Once found, the number is divided by the checked number and the result is again checked with the same algorithms, but with the checked number as the starting point (not 2). 

This yields the following implementation:

 1 
 2 import scala.collection.immutable.NumericRange
 3 
 4 def largestPrimeFactor(x: BigInt) = { 
 5   def largestPrimeFactor(x: BigInt, max: BigInt) : BigInt = {
 6      if (x == 2)  x 
 7      else {
 8         val nextFactor = 
 9           new NumericRange.Inclusive[BigInt](max, csqrt(x), 1).
10              find(x % _ == 0)
11         nextFactor match {
12           case Some(fct) => largestPrimeFactor( x / fct, fct)
13           case None => x
14         } 
15     }
16   }
17   largestPrimeFactor(x, 2)
18 }

Note the use of csqrt function that is  sqrt function implementation that works with BigInt. Then, line 8 - 10 is to get the first number that divides x. Line 11 - 13 acts accordingly, that is, if a factor is found, then it must be a prime number. Then x is divided by the found factor followed by recursive call to the same function.
If the find returns None, then it is the largest prime number.

As mentioned before, the brute force approach above is not practical for real problem with bigger number, especially for a number that is obtained by multiplying two large prime numbers. More sophisticated algorithms are required, and maybe I'll back to some of them one day in this blog. But for now, it's brute force.

Problem : Euler Task 2

Let's continue the effort of answering euler project challenge. Here, we answer the challenge no. 2.


Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 

89, ...

Find the sum of all the even-valued terms in the sequence which do not exceed four million.

First, let's investigate the implementation of fibonacci. A simple recursive program is the following

1 def fib(x: Int) : Long = {
2    if (x >= 2) {
3       fib(x -1) + fib(x -2)
4    } else {
5       1
6    }
7 }

The problem with the code above is the performance. fib(45) in my computer took 77 s. Fortunately, Scala can optimise tail recursion and indeed fibonacci can be implemented using tail recursive algorithms. The implementation is shown in the following code:

1 def fibiter(index: Int) : Long = {
2  def fibiter(x: Long, acc: Long, count: Int) : Long = {
3      if (count == index) 
4          acc
5      else fibiter(x + acc,x, count + 1)
6   }
7   fibiter(1, 0, 0)
8 }

To implement the algorithms, we need to define an auxiliary method, fibiter/3 that receives the "next" fibonacci number, the "current" fibonacci number, and the index of the count that needed to be compared against the requested index.
Scala optimises the code above so that it is treated just like a loop instead of a recursive call. Executing fibiter(45) does not take even 1ms. 
Having defined fibiter method, we're now ready to attack the problem, that is to sum of all even-valued terms.

First, just sum everything, not necessarily the even ones. The code is implemented as follow:

 1 def fibSum(max: Long) : Long = {
 2   def fibSum(x: Long, acc: Long, sum: Long) : Long = {
 3      if (x > max) {
 4          sum
 5      }
 6      else { 
 7        fibSum(x + acc, x, sum + x) 
 8      }
 9  }
10   fibSum(1, 0, 0)
11 }
12 
Like fibiter, we use here an auxilary method, fibSum with 3 parameters. acc is used as the accumulator of fibonacci number, just like in fibiter. In addition, sum is used as the accumulator of the sum of fibonacci number computed so far. The condition at line 3-5 is used to see if the next fibonacci number exceeds the max. If it is the case, the sum so far is returned, otherwise it recursively computes the next fibonacci number.
With the above algorithms in hand, we're now ready to attack the problem, that is, to sum only even fibonacci number:


 1 def fibSum(max: Long) : Long = {
 2   def fibSum(x: Long, acc: Long, sum: Long) : Long = {
 3     if (x > max) {
 4          sum
 5      }
 6      else if (x % 2 == 0) {
 7          fibSum(x + acc, x, sum + x) 
 8      } else {
 9          fibSum(x + acc, x, sum) 
10      }
11   }
12   fibSum(1, 0, 0)
13 }

The lines 6-10 are the key of the solution. The sum is updated only when the current fibonacci number is even. 


Finally, to further generalise the method, the fibSum is modified to accept a function Long => Boolean as follow:
 1 def fibSum(max: Long, filter: Long => Boolean) : Long = {
 2   def fibSum(x: Long, acc: Long, sum: Long) : Long = {
 3     if (x > max) {
 4          sum
 5      }
 6      else if (filter(x)) {
 7          fibSum(x + acc, x, sum + x) 
 8      } else {
 9          fibSum(x + acc, x, sum) 
10      }
11   }
12   fibSum(1, 0, 0)
13 }
To sum all even fibonacci numbers under 4000000, fibSum is called as follow: fibSum(4000000, x => x % 2 ==0).

Problem: Euler task 1

From Euler project, task 1.


"If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000."


Pretty Simple:

(1 to 1000).filter( x => (x %3 == 0) || ( x%5 == 0) ).

            foldLeft(0)(_ + _)

or:

(1 to 1000).foldLeft(0)

   ( (x,y) => x + (if (y % 3 == 0 || y % 5 ==0) y else 0))

The first looks better since the addition is only done when the number is divided by 3 or 5. The second one always executes the addition and the comparison.

However, the first one creates the list of number divided by 3 or 5 first before folding.

So, I think that the second version is better in term of number of created list.

Anyway, I think the following is the best one:

(1 to 1000 / 3).foldLeft(0)( _ + _ * 3) + 

(1 to 1000 / 5).foldLeft(0)( _ + _ * 5) -
(1 to 1000 / 15 ).foldLeft(0)( _ + _ *15)

The first two versions execute 1000 comparisons inside the (1 to ...) , while the second has no comparison. So, I believe the last version should be better than the first two.


In conclusion, my guess of the best performing implementation is the following:
1. Third implementation.
2.  Second implementation.
3.  First implementation.

Maybe the first implementation performs better than the second but the first consumes more memory.

Quick profiling on time of average execution , using 1000 000 as limit instead of 1000 in my machine: 
1. First implementation: 130 ms
2. Second implementation: 65 ms
3. Third implementation: 32 ms.

Finally, I found this one as even more better solution:

(1 to 1000 / 3).foldLeft(0)( _ + _ * 3) + (

 1 to 1000/ 5).foldLeft(0)( 

   (x,y) => x + (if (y % 3 ==0) 0 else (y * 5)))

Any thought ?