How do I evaluate the division of BigIntegers as decimal rather than double in C#?

Question

TL;DR

F# PowerPack's BigRational type can be cast as a double to evaluate the value. However, after the numerator and denominator reach a certain size, the value returned is double.NaN.

Since BigRational keeps track of both the numerator and the denominator as System.Numerics.BigIntegers, you can use a workaround using properties of logarithms:

a / b = e^(ln(a) - ln(b))

With our BigInteger numerator and denominator, we can call

BigInteger num = myBigRational.Numerator;
BigInteger den = myBigRational.Denominator;
double value = Math.Exp(BigInteger.Log(num) - BigInteger.Log(den));

Due to the limitations of the structure of the double type for values close to 0, I would rather use a decimal. I just haven't figured out how.

Context

Just for fun, I'm writing a program that calculates pi using the Taylor series of arctan.

arctan(x) = x - x^3/3 + x^5/5 - x^7/7 +  x^9/9 - ...

If we evaluate the series at 1, we get

arctan(1) = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...

Since arctan(1) = pi/4, we can then multiply our series by 4 to calculate pi.

The goal of my program is to calculate how many terms of the series it takes to converge to pi accurate to n digits. For example, in order for the series to be accurate to one digit (3), it requires the first three terms:

1 term:  4 * (1) = 4
2 terms: 4 * (1 - 1/3) = 2.666666666666667
3 terms: 4 * (1 - 1/3 + 1/5) = 3.466666666666667

To be accurate to 2 digits (3.1), it requires the first 19 terms. Accurate to 3 digits (3.14) requires the first 119 terms, and so on.

I originally wrote my program using C#'s decimal type:

const int MaxDigits = 20;
private static void RunDecimalCalculation()
{
    decimal pi = 0m; // our current approximation of pi
    decimal denominator = 1m;
    decimal addSubtract = 1m;
    ulong terms = 0;

    for (int digits = 0; digits < MaxDigits; digits++)
    {
        decimal piToDigits, upperBound;
        GetBounds(digits, out piToDigits, out upperBound);

        while (pi >= upperBound | pi < piToDigits)
        {
            pi += addSubtract * 4m / denominator;
            denominator += 2m;
            addSubtract *= -1m;
            terms++;
        }

        PrintUpdate(terms, digits, pi);
    }
}

/// <summary>
/// Returns the convergence bounds for <paramref name="digits"/> digits of pi.
/// </summary>
/// <param name="digits">Number of accurate digits of pi.</param>
/// <param name="piToDigits">Pi to the first <paramref name="digits"/> digits of pi.</param>
/// <param name="upperBound">same as <paramref name="piToDigits"/>, but with the last digit + 1</param>
/// <example>
/// <code>GetBounds(1)</code>:
///     piToDigits = 3
///     upperBound = 4
/// 
/// <code>GetBounds(2)</code>:
///     piToDigits = 3.1
///     upperbound = 3.2
/// </example>
private static void GetBounds(int digits, out decimal piToDigits, out decimal upperBound)
{
    int pow = (int)Math.Pow(10, digits);
    piToDigits = (decimal)Math.Floor(Math.PI * pow) / pow;
    upperBound = piToDigits + 1m / pow;
}

I realized, though, that due to rounding errors in each iteration, after enough terms, the number of terms required might be off. So, I started looking into F# PowerPack's BigRational and rewrote the code:

// very minor optimization by caching common values
static readonly BigRational Minus1 = BigRational.FromInt(-1);
static readonly BigRational One = BigRational.FromInt(1);
static readonly BigRational Two = BigRational.FromInt(2);
static readonly BigRational Four = BigRational.FromInt(4);

private static void RunBigRationalCalculation()
{
    BigRational pi = BigRational.Zero;
    ulong terms = 0;
    var series = TaylorSeries().GetEnumerator();

    for (int digits = 0; digits < MaxDigits; digits++)
    {
        BigRational piToDigits, upperBound;
        GetBounds(digits, out piToDigits, out upperBound);

        while (pi >= upperBound | pi < piToDigits)
        {
            series.MoveNext();
            pi += series.Current;
            terms++;
        }

        double piDouble = Math.Exp(BigInteger.Log(pi.Numerator) - BigInteger.Log(pi.Denominator));
        PrintUpdate(terms, digits, (decimal)piDouble);
    }
}

// code adapted from http://tomasp.net/blog/powerpack-numeric.aspx
private static IEnumerable<BigRational> TaylorSeries()
{
    BigRational n = One;
    BigRational q = One;

    while (true)
    {
        yield return q * Four / n;

        n += Two;
        q *= Minus1;
    }
}

Unsurprisingly, this version runs incredibly slowly, which is fine. (The decimal version took 34 seconds to get to 9 accurate digits; the BigRational version took 17 seconds to get to 5 accurate digits, and has been running for about half an hour and hasn't reached 6 accurate digits). What's frustrating me, though, is that double is less accurate than decimal, so while the number of terms will always be correct, the value displayed from

double piDouble = Math.Exp(BigInteger.Log(pi.Numerator) - BigInteger.Log(pi.Denominator));

is inaccurate. Is there a way around this either through mathematical wizardry or some library that has decimal versions of Math.Exp() and BigInteger.Log()?

Answer 1

How do I evaluate the division of BigIntegers as decimal rather than double in C#?

You have BigIntegers N and D and wish to approximate the exact value N / D as a decimal. WOLOG assume both are positive.

That's easy enough. First solve this problem:

• Given N and D, find the BigInteger M that minimizes the absolute value of N / D - M / 10²⁸.

Second, solve this problem:

• Given M, find the BigInteger M' and the non-negative integer E such that M / 10 ²⁸ = M' / 10 ^E and E is minimized.

Third, solve this problem:

• Given M', find non-negative BigIntegers I0, I1, I2 such that M' = I0 + I1 * 2³² + I2 * 2⁶⁴ and each is strictly smaller than 2³². If there are no such I0, I1, I2, then your number cannot be represented as a decimal.

Convert I0, I1 and I2 to unsigned ints and then to signed ints.

Now you have everything you need to call

https://msdn.microsoft.com/en-us/library/bb1c1a6x(v=vs.110).aspx

And hey, you have a decimal in hand.

That said: you should not need to do this in the first place. You are doing your math in BigRationals; why ever would you want to drop out of them into decimals or doubles? Just approximate pi to whatever level you want in big rationals, and compare your slowly-converting series to that.

FYI, this series converges very slowly. Once you have empirically determined how slowly it converges, can you provide a proof of any bound on the rate of convergence?