Lucas Lehmer optimization

Question

I've been working to optimize the Lucas-Lehmer primality test using C# code (yes I'm doing something with Mersenne primes to calculate perfect numbers. I was wondering it is possible with the current code to make further improvements in speed. I use the System.Numerics.BigInteger class to hold the numbers, perhaps it is not the wisest, we'll see it then.

This code is actually based on the intelligence found on: http://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test

This page (at the timestamp) section, some proof is given to optimize the division away.

The code for the LucasTest is:

public bool LucasLehmerTest(int num)
{
  if (num % 2 == 0)
     return num == 2;
  else
  {
     BigInteger ss = new BigInteger(4);
     for (int i = 3; i <= num; i++)
     {
        ss = KaratsubaSquare(ss) - 2;
        ss = LucasLehmerMod(ss, num);
     }
     return ss == BigInteger.Zero;
  }

}

Edit: Which is faster than using ModPow from the BigInteger class as suggested by Mare Infinitus below. That implementation is:

public bool LucasLehmerTest(int num)
{
  if (num % 2 == 0)
     return num == 2;
  else
  {
     BigInteger m = (BigInteger.One << num) - 1;
     BigInteger ss = new BigInteger(4);
     for (int i = 3; i <= num; i++)
       ss = (BigInteger.ModPow(ss, 2, m) - 2) % m;
     return ss == BigInteger.Zero;
  }

}

The LucasLehmerMod method is implemented as follows:

public BigInteger LucasLehmerMod(BigInteger divident, int divisor)
{
   BigInteger mask = (BigInteger.One << divisor) - 1;  //Mask
   BigInteger remainder = BigInteger.Zero;
   BigInteger temporaryResult = divident;

   do
   {
      remainder = temporaryResult & mask;
      temporaryResult >>= divisor;
      temporaryResult += remainder;
   } while ( (temporaryResult >> divisor ) != 0 );

   return (temporaryResult == mask ? BigInteger.Zero : temporaryResult);
}

What I am afraid of is that when using the BigInteger class from the .NET framework, I am bound to their calculations. Would it mean I have to create my own BigInteger class to improve it? Or can I sustain by using a KaratsubaSquare (derived from the Karatsuba algorithm) like this, what I found on Optimizing Karatsuba Implementation:

public BigInteger KaratsubaSquare(BigInteger x)
{
   int n = BitLength(x);                                 

   if (n <= LOW_DIGITS) return BigInteger.Pow(x,2);        //Standard square

   BigInteger b = x >> n;                                  //Higher half
   BigInteger a = x - (b << n);                            //Lower half
   BigInteger ac = KaratsubaSquare(a);                     // lower half * lower half
   BigInteger bd = KaratsubaSquare(b);                     // higher half * higher half
   BigInteger c = Karatsuba(a, b);                         // lower half * higher half

   return ac + (c << (n + 1)) + (bd << (2 * n));          
}

So basically, I want to look if it is possible to improve the Lucas-Lehmer test method by optimizing the for loop. However, I am a bit stuck there... Is it even possible?

Any thoughts are welcome of course.

Some extra thoughs:

I could use several threads to speed up the calculation on finding Perfect numbers. However, I have no experience (yet) with good partitioning. I'll try to explain my thoughts (no code yet):

First I'll be generating a primetable with use of the sieve of Erathostenes. It takes about 25 ms to find primes within the range of 2 - 1 million single threaded.

What C# offers is quite astonishing. Using PLINQ with the Parallel.For method, I could run several calculations almost simultaneously, however, it chunks the primeTable array into parts which are not respected to the search.

I already figured out that the automatic load balancing of the threads is not sufficient for this task. Hence I need to try a different approach by dividing the loadbalance depending on the mersenne numbers to find and use to calculate a perfect number. Has anyone some experience with this? This page seems to be a bit helpful: http://www.drdobbs.com/windows/custom-parallel-partitioning-with-net-4/224600406

I'll be looking into it further.

As for now, my results are as following. My current algorithm (using the standard BigInteger class from C#) can find the first 17 perfect numbers (see http://en.wikipedia.org/wiki/List_of_perfect_numbers) within 5 seconds on my laptop (an Intel I5 with 4 cores and 8GB of RAM). However, then it gets stuck and finds nothing within 10 minutes.

This is something I cannot match yet... My gut feeling (and common sense) tells me that I should look into the LucasLehmer test, since a for-loop calculating the 18th perfect number (using Mersenne Prime 3217) would run 3214 times. There is room for improvement I guess...

What Dinony posted below is a suggestion to rewrite it completely in C. I agree that would boost my performance, however I choose C# to find out it's limitations and benefits. Since it's widely used, and it's ability to rapidly develop applications, it seemed to me worthy of trying.

Could unsafe code provide benefits here as well?

Nightie

Do you really need the BigInteger here? And if you answer that with yes, perhaps ModPow is your friend. – Avigation

Mare Infinitus, that might be interesting. Yes I am in need of a BigInteger, since the perfect numbers have a tendency to grow quite big... The 17th perfect number consists of 1373 digits. I'll see what I can do with ModPow! – Nightie

I've written a test myself some ages ago, but in C# and with BigInteger. The ModPow definitly made the difference. Think it was for Project Euler. The question is: do you need that special algorithm or just one that does the job? – Avigation

Well I've just put it up to a test, it significantly increases performance! How could I have missed it... Thank you very much Mare Infinitus! I'll get back with some comparison results. – Nightie

Mare Infinitus, due to a small mistake in my code (it also found non mersenne primes), which I just fixed, I must admit that it works. It finds the first 20 perfect numbers in 138004 ms. However, my assumption (and what I found on Wiki) was right... Using the LucasLehmerMod function which I described above, finds the 20th perfect number in 116674 ms. That is 14 seconds faster... – Nightie

I meant almost 22 seconds... I had something else in mind while typing that I guess :-) – Nightie

Answer 1

One possible optimization is to use BigInteger ModPow

It really increases performance significantly.

Avigation

Mare Infinitus, I've updated your earlier suggestion in my original post. Implementing and using ModPow is actually slower than my optimization which removes divisions... However, plus one for the suggestion! – Nightie

I've accepted your answer thus far. It is the best suggestion when testing primes before entering the Lucas Lehmer test! – Nightie

Answer 2

Just a note for info... In python, this

ss = KaratsubaSquare(ss) - 2

has worse performance than this:

ss = ss*ss - 2

Snaggletooth

If any library implements large integer operations, then I would assume that it does it in the fastest possible way. So Python will use either the Karatsuba method itself, but highly optimised and tuned, or it will use faster methods, like the generalisation of the Karatsuba method, or FFT methods. – Amaleta

Answer 3

What about adapting the code to C? I have no idea about the algorithm, but it is not that much code.. so the biggest run-time improvement could be adapting to C.

Dotdotage

Yeah this doesn't seem like it would be too difficult to convert, either. Not much infrastructure that C# makes easier than C. – Yoakum

That is something I can agree with totally @dinony. However, I felt like programming it in C#. If that is the true optimization (which I do not doubt at all), I must look further. Somehow it must be possible to speed up the algorithm. I also noticed that the Karatsuba could be implemented a bit more efficient by using bitmasks. – Nightie

Yap I know, programming in C# is really tempting and fun, but when I hear optimization I think C/C++ :) .. I know it's not the answer you wanted, maybe someone who knows the algorithm can give you C# optimization hints – Dotdotage

Haha, Dinony, I know you are right. Perhaps that would be the next choice of course. Mare Infinitus put on a very helpful suggestion. I'll put that up for testing at the moment! – Nightie