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()?
