2

I wrote a simple single-threaded application that checks whether a given number is prime. It is meant to support arbitrarily large integers, so I made use of the BigInteger class in C#. Here is the code:

using System;
using System.Diagnostics;
using System.Numerics;

namespace PrimeChecker {
    public static class Program {
        /// <summary>
        /// Gets the number of threads that parallel loops will use from the thread pool when called.
        /// </summary>
        private static BigInteger MinusOne = BigInteger.MinusOne,
            Zero = BigInteger.Zero,
            One = BigInteger.One,
            Two = new BigInteger(2),
            Three = new BigInteger(3),
            Five = new BigInteger(5),
            Six = new BigInteger(6);

        public static void Main(string[] args) {
            string value = args.Length == 0 ? null : args[0];
            BigInteger num, factor;
            Stopwatch elapsed = new Stopwatch();
            do {
                if (value == null) {
                    Console.Write("\nEnter integer to check if prime (or blank to exit): ");
                    value = Console.ReadLine();
                }
                if (value.Length == 0)
                    return;
                value.Trim();
                if (BigInteger.TryParse(value, out num)) {
                    Console.Write("Checking...");
                    elapsed.Restart();
                    if (num.IsZero)
                        Console.WriteLine("nope, divisible by 0. :/ (elapsed: 0ms)");
                    else if (num.IsOne)
                        Console.WriteLine("nope, divisible by 1. :/ (elapsed: 0ms)");
                    else {
                        factor = IsPrime(num);
                        elapsed.Stop();
                        if (factor.IsOne)
                            Console.WriteLine("prime! :D (elapsed: " + elapsed.Elapsed.ToString() + ")");
                        else
                            Console.WriteLine("nope, divisible by at least " + factor + " and " + (num / factor) + ". :/ (elapsed: " + elapsed.Elapsed.ToString() + ")");
                    }
                } else
                    Console.WriteLine("Not a valid integer... :*");
                value = null;
            } while (true);
        }

        private static BigInteger LogBase2(BigInteger num) {
            if (num <= Zero)
                return MinusOne; //does not support negative values
            BigInteger i = Zero;
            while (!(num >>= 1).IsZero)
                i++;
            return i;
        }

        public static BigInteger Sqrt(this BigInteger num) {
            if (num.IsZero)
                return Zero;
            else if (num.IsOne)
                return One;
            else if (num.Sign > 0) {
                BigInteger root = LogBase2(num);
                BigInteger lowerBound;
                while (num < (lowerBound = root * root) || num > lowerBound + root + root) {
                    root += num / root;
                    root >>= 1;
                }
                return root;
            } else
                return MinusOne;
        }

        /// <summary>
        /// Returns 0 if num is 0, 1 if the value is prime or 1, -1 if the value is negative,
        /// else returns the smallest factor of the num.
        /// </summary>
        /// <param name="num">The integer to check if prime.</param>
        public static BigInteger IsPrime(BigInteger num) {
            if (num.IsZero)
                return Zero;
            else if (num.Sign < 0)
                return MinusOne;
            else if (num <= Three || num == Five)
                return One;
            else if (num.IsEven)
                return Two;
            else if ((num % Three).IsZero)
                return Three;
            BigInteger root = Sqrt(num);
            if (!root.IsEven)
                root++;
            BigInteger i = Five;
            while (i < root && !((num % i).IsZero || (num % (i += Two)).IsZero)) //main loop here
                i += Four;
            return i < root ? i : One;
        }
    }
}

However, I wanted to investigate writing the same application in C++. I managed to accomplish this by making use of the MPIR library on Windows. The full code code and Visual Studio project can be found on GitHub: https://github.com/mathusummut/PrimeCheckerCpp, but here is the main loop:

while (mpz_cmp(i, root) == -1) {  //loop is here
    mpz_mod(temp, num, i);
    if (mpz_sgn(temp) == 0)
        break;
    mpz_add_ui(i, i, 2u);
    mpz_mod(temp, num, i);
    if (mpz_sgn(temp) == 0)
        break;
    mpz_add_ui(i, i, 4u);
}

When I run the C# executable to check whether 2305843009213693951 is prime, I get this output:

Enter integer to check if prime (or blank to exit): 2305843009213693951
Checking...prime! :D (elapsed: 00:00:52.2468554)

whereas when I run the C++ executable, I get this output:

Enter integer to check if prime (or blank to exit): 2305843009213693951
Checking...prime! :D (elapsed: 60.429640s)

Why is the C++ implementation slower? Is the MPIR library at fault here, or am I doing something wrong?

NB: Both are compiled and tested in Release mode 64-bit with full optimizations enabled with no debugger attached.

NB: @hatchet was right, when I tried an integer that is not representable with 64 bits, the C++ implementation was faster. For example, trying 108446744073709551647, the C# version took 13.5 minutes whereas the C++ version took only 6.5 minutes to execute!

MathuSum Mut
  • 2,765
  • 3
  • 27
  • 59
  • 1
    To me it seems in C# it seems you step by 6 (`i += Six;`) but in C++ you step by 4 (`mpz_add_ui(i, i, 4u);`)? – Thomas Weller Mar 06 '17 at 23:05
  • 1
    That's not exactly the case, as in the C++ version the increment by 2 is performed in place, whereas in C# it's not, that's why I increment by 6 not by 4 in C#. The code is functionally equivalent. – MathuSum Mut Mar 06 '17 at 23:06
  • 2
    I think this boils down to what optimizations have been cooked into the library implementations. I can't speak to MPIR, but .Net's BigInteger mod operation looks for cases where the values can fit in a uint or ulong, and if so, just do regular mod with those values, thus skipping the more expensive mod operation that would happen for larger numbers. Also, once you get into numbers greater than max ulong value, I think you'll see your isPrime times degrade pretty badly the larger the numbers get. – hatchet - done with SOverflow Mar 06 '17 at 23:08
  • 2
    By comparison, a miller-rabin prime check for the test value in your question (which can be done deterministically for all 64 bit values) runs about 4 orders of magnitude faster than the isPrime method you're using. – hatchet - done with SOverflow Mar 06 '17 at 23:11
  • I suspect that hatchet is correct. Try a number that's about twice the magnitude of the one you are using. Please post if you find this is the case (I'm curious now :-)) – Dweeberly Mar 06 '17 at 23:13
  • I am aware that the current implementation is simplistic and slow, and I seek to improve upon it, but truthfully the question is about whether using MPIR as a BigInteger library is worth the time. – MathuSum Mut Mar 06 '17 at 23:16
  • 1
    My point is that your test number is too small to run the code that matters in BigInteger, and because of that, it's not a useful test of those libraries. – hatchet - done with SOverflow Mar 06 '17 at 23:18
  • Looking at the Miller-Rabin code, I see that it uses random numbers...is it guaranteed to give the correct answer? – MathuSum Mut Mar 06 '17 at 23:20
  • @MathuSumMut - up to 64 bit values, witness values have been determined that guarantee deterministic answers. http://blog.softwx.net/2013/05/miller-rabin-primality-test-in-c.html .Yes, beyond that, it is probabilistic, but given enough witness values, has a high level of reliability. – hatchet - done with SOverflow Mar 06 '17 at 23:23
  • 1
    @MathuSumMut - try your comparison using 108446744073709551647. It will take about a half hour for BigInteger, but this is out of the ulong range and may be better for comparing to BigInteger to MPIR for your use. – hatchet - done with SOverflow Mar 07 '17 at 00:25
  • Your hypothesis was correct! I tried `108446744073709551647`, and it took 13.5 minutes for the C# version, and the C++ version only took 6.5 minutes. :) – MathuSum Mut Mar 07 '17 at 09:29
  • @MathuSumMut - I'm glad you tried the larger example. It's a minor point, but there is a minor bug in your code. In your isPrime test, there are two mod tests in your loop condition. If `num % i` is not zero, but `num % i+1` is zero, you return `i` as the factor, but actually it was `i+1` that was the factor. – hatchet - done with SOverflow Mar 08 '17 at 18:01
  • Ah, right, in the C# version, yes? I'll increment i in place then, thanks. (fixed) – MathuSum Mut Mar 08 '17 at 18:20

0 Answers0