Day: April 23, 2008

I thought I’d post a few code snippets in C# that have to do with basic number theory, since I use them in my programs from time to time. These snippets are by no means optimized, and the use of C# pretty much precludes their use in high-performance applications. Still, for relatively small arguments, these routines run surprisingly fast.

Prime numbers

To generate a list of prime numbers, we use the familiar Sieve of Eratosthenes. We can write one routine to generate the sieve:

bool[] GetPrimeSieve(long upTo)
{
    long sieveSize = upTo + 1;
    bool[] sieve = new bool[sieveSize];
    Array.Clear(sieve, 0, (int)sieveSize);
    sieve[0] = true;
    sieve[1] = true;
    long p, max = (long)Math.Sqrt(sieveSize) + 1;
    for (long i = 2; i <= max; i++)
    {
        if (sieve[i]) continue;
        p = i + i;
        while (p < sieveSize) { sieve[p] = true; p += i; }
    }
    return sieve;
}

The above function returns an array of booleans (the size of the given parameter), each of which is false if the array index is a prime number, or true if it's not a prime number. To get an actual list of prime numbers, we can use a function such as this:

long[] GetPrimesUpTo(long upTo)
{
    if (upTo < 2) return null;
    bool[] sieve = GetPrimeSieve(upTo);
    long[] primes = new long[upTo + 1];

    long index = 0;
    for (long i = 2; i <= upTo; i++) if (!sieve[i]) primes[index++] = i;

    Array.Resize(ref primes, (int)index);
    return primes;
}

The above function returns an actual list of primes less than or equal to the specified high limit. This means that we can easily compute the prime-counting function $$\pi(x)$$ for any integer x by counting the number of elements in the array returned by the function. We can write a similar function to generate a list of composites:

long[] GetCompositesUpTo(long upTo)
{
    if (upTo < 2) return null;
    bool[] sieve = GetPrimeSieve(upTo);
    long[] composites = new long[upTo + 1];

    long index = 0;
    for (long i = 2; i <= upTo; i++)
        if (sieve[i]) composites[index++] = i;

    Array.Resize(ref composites, (int)index);
    return composites;
}

As for determining if a single certain number is prime (without having to generate a giant sieve), we can simply use a function that attempts to factor the number. If the number happens to be divisible by an integer greater than 1 and less than or equal to its square root, then the number is not prime:

bool IsPrime(long n)
{
    long max = (long)Math.Sqrt(n);
    for (long i = 2; i <= max; i++)
        if (n % i == 0)
            return false;
    return true;
}

Greatest Common Divisor and Totient

To obtain the greatest common divisor (GCD) of two numbers, we use the usual Euclidean algorithm:

long gcd(long a, long b)
{
    long temp;
    while (b != 0)
    {
        temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

There is a slightly different binary algorithm for computing the GCD which is theoretically more efficient, but in practice (at least in C#) it's actually slightly less efficient than the algorithm above:

ulong gcdb(ulong a, ulong b)
{
    int shift;
    if (a == 0 || b == 0) return a | b;

    for (shift = 0; ((a | b) & 1) == 0; ++shift) { a >>= 1; b >>= 1; }
    while ((a & 1) == 0) a >>= 1;

    do {
        while ((b & 1) == 0) b >>= 1;

        if (a < b) {
            b -= a;
        } else {
            ulong temp = a - b;
            a = b;
            b = temp;
        }
        b >>= 1;
    } while (b != 0);
    return a << shift;
}

With the GCD readily available, determining whether two numbers are coprime is just a matter of telling whether or not their GCD is equal to 1. Also, calculating the LCM (least common multiple) of two numbers becomes trivial:

long lcm(long a, long b)
{
    long ret = 0, temp = gcd(a, b);
    if (temp != 0)
    {
        if (b > a) ret = (b / temp) * a;
        else ret = (a / temp) * b;
    }
    return ret;
}

Also using the GCD algorithm, it becomes easy to calculate Euler's totient function for a certain number, since the totient function is simply the number of integers less than or equal to n that are coprime to n:

long eulerTotient(long n)
{
    long sum = 0;
    for (long i = 1; i <= n; i++)
        if (gcd(i, n) == 1) sum++;
    return sum;
}

The above is a really naive algorithm. A much more efficient algorithm (one that is usually given in textbooks) is as follows:

int eulerTotient(int n)
{
    primes = GetPrimesUpTo(n+1);    //this can be precalculated beforehand
    int numPrimes = primes.Length;

    int totient = n;
    int currentNum = n, temp, p, prevP = 0;
    for (int i = 0; i < numPrimes; i++)
    {
        p = (int)primes[i];
        if (p > currentNum) break;
        temp = currentNum / p;
        if (temp * p == currentNum)
        {
            currentNum = temp;
            i--;
            if (prevP != p) { prevP = p; totient -= (totient / p); }
        }
    }
    return totient;
}