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

by Marc Sigrist 15. November 2010 20:30

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! 

Through the introduction of System.Numerics.BigInteger, it is now possible for any .Net Framework language to exactly represent any whole number whatsoever. 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 by any .Net 4.0 client application. The signature, seen from C#, should be as follows:

Desired Signature, as Seen From C#

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

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

Implementation in F#

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 the .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  

Usage from C#

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

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

Concluding Remarks

The parameters of the PrimesBetween function are declared in tupled (not curried) form, which corresponds to par. 4.2 "Object and Member Design" of the F# Component Design Guidelines, which says "Do not use currying of parameters in vanilla .NET APIs".

All recursive functions are implemented tail recursively, using accumulator arguments, so that no stack overflow will ever occur.

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 speed up the function using a cache. Perhaps unexpectedly, this cache would not necessarily store the resulting prime numbers. Instead, it would be gradually built with all prime numbers from the bottom up, starting at 2, whenever the recursive hasDiv function is called. The hasDiv function would then be modified to no more try out all whole numbers starting from 2 as potential divisors, but only all (recursively cached) primary numbers starting from 2, which is sufficient. This would enormously reduce the amount of divisions. The next step would then be to persist the cache in a database. The resulting library could be run for many years, returning larger and larger prime numbers, until (presumably) the currently calculated BigInteger number would be so large as to overstrain the system memory. I leave it up to you to speculate approximately when this might happen...