Calculating a range of arbitrarily large prime numbers in F# 2.0

The library System.Numerics.dll, who was introduced in .Net 4.0, contains a System.Numerics.BigInteger structure. BigInteger represents a whole number of arbitrary size (or precision). Before .Net 4.0, the largest number that could be represented out of the box was System.Double.MaxValue. Written in decimal notation, this would be a whole number with 309 digits (more than 179 thousand centillion). However, using Double for whole number calculations is error-prone. Double is internally stored with floating point logic; the larger the number, the less exact it becomes. For instance, you cannot exactly represent the number ten quadrillion (i.e., ten million billion):

let x = 1.0e+16
let y = x - 1.0
printf "%b" (x = y)
// Prints true! 

Thanks to the introduction of System.Numerics.BigInteger, it is now possible for any .Net Framework language to exactly represent any whole number regardless of its size. In the following example, we will use this type to implement a function that calculates a range of arbitrarily large prime numbers. We will design it as a static library method which can be called from any .Net 4.0 client application. The signature, seen from C#, will be as follows:

using System.Collections.Generic;
using System.Numerics;

namespace MyCompany.MyProduct{
    public static class BigMath{
        public static IEnumerable<BigInteger> 
            PrimesBetween(BigInteger i, BigInteger j){
            // etc.

To implement the library method, proceed as follows:

• (If not done yet: Install the F# Power Pack)
• Create a new F# class library project
• Set .Net Framework 4 client profile as the compilation target
• Add references to System.Numerics.dll and FSharp.PowerPack.dll
• Add a new source file BigMath.fs with the following code:

/// Provides static methods for mathematical functions involving 
/// arbitrarily large numbers. 
module MyCompany.MyProduct.BigMath

/// <summary>Gets a lazily evaluated sequence of prime numbers.</summary>
/// <param name="i">The beginning of the interval for the prime numbers.</param>
/// <param name="j">The end of the interval for the prime numbers.</param>
/// <returns>A sequence with the prime numbers between <paramref name="i"/> and <paramref name="j"/>.</returns>
let PrimesBetween(i, j) =
    /// Approximates the square root of a BigInteger with a 
    /// BigRational, according to the "babylonian method". 
    let sqrtI i =
        let iN = BigRational.FromBigInt i
        let half n = n / 2N
        let rec approachWith n =
            let diff = BigRational.PowN(n, 2) - iN
            if diff >= 0N && diff < 1N then n
            else approachWith(half(n + iN / n))
        approachWith(half iN)

    /// Gets true if i can be divided by a number >= 2,
    /// where i is a non-negative BigInteger.
    let hasDivisor i =
        let maxDiv = 
            let sr = sqrtI i
            min (sr.Numerator / sr.Denominator + 1I) (i - 1I)
        let rec hasDiv div =
            if div > maxDiv then false
            elif i % div = 0I then true
            else hasDiv(div + 1I)
        hasDiv 2I   

    /// Gets true if i is a prime number. 
    let isPrime i =
        match i with
        | _ when i < 2I -> false
        | _ when i = 2I -> true
        | _ -> not (hasDivisor i)

    /// The sequence of numbers between i and j.
    let range =
        let step = i |> compare j
        match step with
        | 0 -> Seq.singleton i
        | _ -> seq {i .. bigint step .. j}

    /// The resulting sequence of prime numbers between i and j.
    range |> Seq.filter isPrime  

The compiled F# library can be used in C# as follows:

// using System;
// using MyCompany.MyProduct;
var primes = BigMath.PrimesBetween(10000200, 10000150);
foreach (var p in primes)
    Console.WriteLine(p.ToString("#,#"));
// Produces:
// 10,000,189
// 10,000,169

The parameters of the PrimesBetween function are declared in tupled (not curried) form, which conforms to the F# Component Design Guidelines' subsection Avoid the use of currying of parameters.

All recursive functions are implemented with tail recursion, using accumulator arguments. As a consequence, the compiled IL code can never produce a stack overflow.

The result is a strongly typed generic enumerable of BigInteger. At any given time, no more than one prime number is held in memory, regardless of the size of the interval between i and j.

The square root calculation in the sqrtI function is an optimization. To know whether a number has a divisor, it is sufficient to try division by all numbers between 2 and the square root. For the calculation of the square root, we use the type Microsoft.FSharp.Math.BigRational from FSharp.PowerPack.dll.

To calculate extremely large prime numbers in a productive situation, we would have to further speed up the function using a cache.