3

Maybe I didn't understood how to use those two types: BigInteger/BigRational, but generally speaking I want to implement those two equations: enter image description here

This is my data: n=235, K = 40 and this small p (which actually is called rho) is 5. In the beginning the problem was with the Power function: the results were very very very big - so that is why I used the BigInteger library. But then I realize that there will be a division made and the result will be a number of type double - so I changed to BigRational library.

Here is my code:

static void Main(string[] args)
    {
        var k = 40;
        var n = 235;
        var p = 5;

        // the P(n) equation
        BigRational pnNumerator = BigRational.Pow(p, n);
        BigRational pnDenominator = BigRational.Pow(k, (n - k)) * Factorial(k);


        // the P(0) equation

        //---the right side of "+" sign in Denominator
        BigRational pk = BigRational.Pow(p, k);
        BigRational factorialK = Factorial(k);
        BigRational lastPart = (BigRational.Subtract(1, (double)BigRational.Divide(p, k)));
        BigRational factorialAndLastPart = BigRational.Multiply(factorialK, lastPart);
        BigRational fullRightSide = BigRational.Divide(pk, factorialAndLastPart);
        //---the left side of "+" sign in Denominator
        BigRational series = Series(k, p, n);


        BigRational p0Denominator = series + fullRightSide;
        BigRational p0Result = BigRational.Divide(1, p0Denominator);

        BigRational pNResult = BigRational.Divide((pnNumerator * p0Result), pnDenominator);
        Console.WriteLine(pNResult);
        Console.ReadKey();
    }

    static BigRational Series(int k, int p, int n)
    {
        BigRational series = new BigRational(0.0);
        var fin = k - 1;
        for (int i = 0; i < fin; i++)
        {
            var power = BigRational.Pow(p, i);
            var factorialN = Factorial(n);
            var sum = BigRational.Divide(power, factorialN);
            series += sum;
        }
        return series;
    }

    static BigRational Factorial(int k)
    {
        if (k <= 1)
            return 1;
        else return BigRational.Multiply(k, Factorial(k - 1));
    }

The main problem is that it does not return any "normal" value like for example 0.3 or 0.03. The result which is printed to the console is a very long number (like 1200 digits in it)...

Can someone please take a look at my code and help me fix the problem and be able to solve this equations by the code. Thank you

Eru
  • 332
  • 3
  • 17

1 Answers1

6

Console.WriteLine(pNResult); calls BigRational.ToString() under-the-hood, which prints the number in the form numerator/denominator.

It's easy to miss the / in the output given how large the numerator and denominator both are in this case.

BigRational supports conversions to decimal and to double. The result is too small to fit in a decimal in this case though. Converting to a double, gives the result 7.89682541396914E-177.

If you need better precision, you'll need a custom conversion to a decimal-format string, like the one in this Stackoverflow answer.

Using that custom conversion routine to get the result to 1000 decimal places -

Console.WriteLine(pNResult.ToDecimalString(1000));

- gives the result as:

0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000078968254139691306770128897137459492828971170349380336740935269651539684650525033676003134593283361305530675112470528408219177025044254116462798561450442318290046626248451723040397770263675109107145461310779641705093156106311143727608208629473359566457461384474633112850335950017209558136575135801388668687571284492241030561019606955986265585636660304889792027894460104216176719717671500843399685686146432982358441225578366059001576682388503227237202077881334695352338638383337717103303153521108812750644260562351186866587629456292506971252525125976755540274041651740194108430555751648707933592643410475214924394223640168857340953563111097979394441303100701008120008166339365089771585037880235325673143152814510586536335380671360865230428857049658368242543653234599817430185879648427434216378356518036776477170130227628307039

To check that your calculation code is operating correctly, you can add unit-tests for the different functions (Factorial, Series and the computation of P itself).

An approach that is practical here is to calculate the results by hand for certain small values of k, n and p and check that your functions compute the same results.

If you're using Visual Studio, you can use this MSDN page as a starting point for creating a unit-test project. Note that the functions under test must be visible to the unit-test project, and your unit-test project will need to have a reference added to your existing project where you're doing the computation, as explained in the link.

Starting with Factorial, which is the easiest to check, you could add a test like this:

[TestClass]
public class UnitTestComputation
{
    [TestMethod]
    public void TestFactorial()
    {
        Assert.AreEqual(1, Program.Factorial(0));
        Assert.AreEqual(1, Program.Factorial(1));
        Assert.AreEqual(2, Program.Factorial(2));
        Assert.AreEqual(6, Program.Factorial(3));
        Assert.AreEqual(24, Program.Factorial(4));
    }
}

The code in your question passes that test.

You can then add a test method for your Series function:

[TestMethod]
public void TestSeries()
{
    int k = 1;
    int p = 1;
    BigRational expected = 1;
    Assert.AreEqual(expected, Program.Series(k, p));

    k = 2;
    p = 1;
    expected += 1;
    Assert.AreEqual(expected, Program.Series(k, p));

    k = 3;
    p = 1;
    expected += (BigRational)1 / (BigRational)2;
    Assert.AreEqual(expected, Program.Series(k, p));

    k = 1;
    p = 2;
    expected = 1;
    Assert.AreEqual(expected, Program.Series(k, p));

    k = 2;
    p = 2;
    expected += 2;
    Assert.AreEqual(expected, Program.Series(k, p));
}

This showed up some problems in your code:

  • n shouldn't actually be a parameter to this function, because in this context n isn't the parameter to function P, but actually just the "index-of-summation". n's local value in this function is represented by your i variable.
  • This then means that your Factorial(n) call needs to change to Factorial(i)
  • The loop is also off-by-one, because the Sigma notation for the summation is inclusive of the number at the top of the Sigma, so you should have <= fin (or you could also have written this simply as < k).

This is the updated Series function:

// CHANGED: Removed n as parameter (n just the index of summation here)
public static BigRational Series(int k, int p)
{
    BigRational series = new BigRational(0.0);
    var fin = k - 1;
    // CHANGED: Should be <= fin (i.e. <= k-1, or < k) because it's inclusive counting
    for (int i = 0; i <= fin; i++)
    {
        var power = BigRational.Pow(p, i);
        // CHANGED: was Factorial(n), should be factorial of n value in this part of the sequence being summed.
        var factorialN = Factorial(i);
        var sum = BigRational.Divide(power, factorialN);
        series += sum;
    }
    return series;
}

To test the P(n) calculation you can move that out into its own function to test (I've called it ComputeP here):

[TestMethod]
public void TestP()
{
    int n = 1;
    int k = 2;
    int p = 1;
    // P(0) = 1 / (2 + 1/(2*(1 - 1/2))) = 1/3
    // P(1) = (1/(1/2 * 2)) * P(0) = P(0) = 1/3 
    BigRational expected = 1;
    expected /= 3;
    Assert.AreEqual(expected, Program.ComputeP(k, n, p));

    n = 2;
    k = 2;
    p = 1;
    // P(2) = (1/(1*2)) * P(0) = 1/6
    expected = 1;
    expected /= 6;
    Assert.AreEqual(expected, Program.ComputeP(k, n, p));
}

This showed up a problem with calculating P(n) - you had a cast to double in there which shouldn't have been present (the result is inaccurate then - you need to keep all the intermediate results in BigRational). There's no need for the cast, so just removing it fixes this problem.

Here is the updated ComputeP function:

public static BigRational ComputeP(int k, int n, int p)
{
    // the P(n) equation
    BigRational pnNumerator = BigRational.Pow(p, n);
    BigRational pnDenominator = BigRational.Pow(k, (n - k)) * Factorial(k);


    // the P(0) equation

    //---the right side of "+" sign in Denominator
    BigRational pk = BigRational.Pow(p, k);
    BigRational factorialK = Factorial(k);
    // CHANGED: Don't cast to double here (loses precision)
    BigRational lastPart = (BigRational.Subtract(1, BigRational.Divide(p, k)));
    BigRational factorialAndLastPart = BigRational.Multiply(factorialK, lastPart);
    BigRational fullRightSide = BigRational.Divide(pk, factorialAndLastPart);
    //---the left side of "+" sign in Denominator
    BigRational series = Series(k, p);


    BigRational p0Denominator = series + fullRightSide;
    BigRational p0Result = BigRational.Divide(1, p0Denominator);

    BigRational pNResult = BigRational.Divide((pnNumerator * p0Result), pnDenominator);
    return pNResult;
}

For avoidance of confusion, here is the whole updated calculation program:

using System;
using System.Numerics;
using System.Text;
using Numerics;

public class Program
{
    public static BigRational ComputeP(int k, int n, int p)
    {
        // the P(n) equation
        BigRational pnNumerator = BigRational.Pow(p, n);
        BigRational pnDenominator = BigRational.Pow(k, (n - k)) * Factorial(k);


        // the P(0) equation

        //---the right side of "+" sign in Denominator
        BigRational pk = BigRational.Pow(p, k);
        BigRational factorialK = Factorial(k);
        // CHANGED: Don't cast to double here (loses precision)
        BigRational lastPart = (BigRational.Subtract(1, BigRational.Divide(p, k)));
        BigRational factorialAndLastPart = BigRational.Multiply(factorialK, lastPart);
        BigRational fullRightSide = BigRational.Divide(pk, factorialAndLastPart);
        //---the left side of "+" sign in Denominator
        BigRational series = Series(k, p);


        BigRational p0Denominator = series + fullRightSide;
        BigRational p0Result = BigRational.Divide(1, p0Denominator);

        BigRational pNResult = BigRational.Divide((pnNumerator * p0Result), pnDenominator);
        return pNResult;
    }

    // CHANGED: Removed n as parameter (n just the index of summation here)
    public static BigRational Series(int k, int p)
    {
        BigRational series = new BigRational(0.0);
        var fin = k - 1;
        // CHANGED: Should be <= fin (i.e. <= k-1, or < k) because it's inclusive counting
        for (int i = 0; i <= fin; i++)
        {
            var power = BigRational.Pow(p, i);
            // CHANGED: was Factorial(n), should be factorial of n value in this part of the sequence being summed.
            var factorialN = Factorial(i);
            var sum = BigRational.Divide(power, factorialN);
            series += sum;
        }
        return series;
    }

    public static BigRational Factorial(int k)
    {
        if (k <= 1)
            return 1;
        else return BigRational.Multiply(k, Factorial(k - 1));
    }

    static void Main(string[] args)
    {
        var k = 40;
        var n = 235;
        var p = 5;
        var result = ComputeP(k, n, p);
        Console.WriteLine(result.ToDecimalString(1000));
        Console.ReadKey();
    }
}

// From https://stackoverflow.com/a/10359412/4486839
public static class BigRationalExtensions
{
    public static string ToDecimalString(this BigRational r, int precision)
    {
        var fraction = r.GetFractionPart();

        // Case where the rational number is a whole number
        if (fraction.Numerator == 0 && fraction.Denominator == 1)
        {
            return r.GetWholePart() + ".0";
        }

        var adjustedNumerator = (fraction.Numerator
                                           * BigInteger.Pow(10, precision));
        var decimalPlaces = adjustedNumerator / fraction.Denominator;

        // Case where precision wasn't large enough.
        if (decimalPlaces == 0)
        {
            return "0.0";
        }

        // Give it the capacity for around what we should need for 
        // the whole part and total precision
        // (this is kinda sloppy, but does the trick)
        var sb = new StringBuilder(precision + r.ToString().Length);

        bool noMoreTrailingZeros = false;
        for (int i = precision; i > 0; i--)
        {
            if (!noMoreTrailingZeros)
            {
                if ((decimalPlaces % 10) == 0)
                {
                    decimalPlaces = decimalPlaces / 10;
                    continue;
                }

                noMoreTrailingZeros = true;
            }

            // Add the right most decimal to the string
            sb.Insert(0, decimalPlaces % 10);
            decimalPlaces = decimalPlaces / 10;
        }

        // Insert the whole part and decimal
        sb.Insert(0, ".");
        sb.Insert(0, r.GetWholePart());

        return sb.ToString();
    }
}

And here is the whole unit-test program:

using System;
using Microsoft.VisualStudio.TestTools.UnitTesting;
using Numerics;


[TestClass]
public class UnitTestComputation
{
    [TestMethod]
    public void TestFactorial()
    {
        Assert.AreEqual(1, Program.Factorial(0));
        Assert.AreEqual(1, Program.Factorial(1));
        Assert.AreEqual(2, Program.Factorial(2));
        Assert.AreEqual(6, Program.Factorial(3));
        Assert.AreEqual(24, Program.Factorial(4));
    }

    [TestMethod]
    public void TestSeries()
    {
        int k = 1;
        int p = 1;
        BigRational expected = 1;
        Assert.AreEqual(expected, Program.Series(k, p));

        k = 2;
        p = 1;
        expected += 1;
        Assert.AreEqual(expected, Program.Series(k, p));

        k = 3;
        p = 1;
        expected += (BigRational)1 / (BigRational)2;
        Assert.AreEqual(expected, Program.Series(k, p));

        k = 1;
        p = 2;
        expected = 1;
        Assert.AreEqual(expected, Program.Series(k, p));

        k = 2;
        p = 2;
        expected += 2;
        Assert.AreEqual(expected, Program.Series(k, p));
    }

    [TestMethod]
    public void TestP()
    {
        int n = 1;
        int k = 2;
        int p = 1;
        // P(0) = 1 / (2 + 1/(2*(1 - 1/2))) = 1/3
        // P(1) = (1/(1/2 * 2)) * P(0) = P(0) = 1/3 
        BigRational expected = 1;
        expected /= 3;
        Assert.AreEqual(expected, Program.ComputeP(k, n, p));

        n = 2;
        k = 2;
        p = 1;
        // P(2) = (1/(1*2)) * P(0) = 1/6
        expected = 1;
        expected /= 6;
        Assert.AreEqual(expected, Program.ComputeP(k, n, p));
    }
}

Incidentally, the P(n) result with the updated program for your input values for n, p and k is now:

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000593109980769066916025972569398424267669807629726200017375290861590898269902277869938365969961320969473356001666906480007119114830921839913623591124192047955091318951831902550404167336054683697071654765071519020060437129398945035521954738463786221029427589397688847246112810536958194364039693387170592425527136243952416704526069736811587380688876091926255908361275575249492845970903676492429684929779402600032481018886875698972533534890841796034626337674846620462046294537488580901129338625628349474358946962065227890599744775562637784553656488649841148591533557896418988044457914999854241038974478576578909626765823565817758792682480009619613438867365912697996527957775248350987801430141776875171808382272960426476953742528769626555642957093028553993908356226007570404005591174451216846471710162760343

NOTE: You should add to the unit-tests with more results you've checked by hand, and also check any of my working here in interpreting the algebra as code to ensure this is correct.

Community
  • 1
  • 1
softwariness
  • 4,022
  • 4
  • 33
  • 41
  • Okay, thank you for the answer - can you tell me if my code and this result is counted correctly? For example does my functions Series and Factorial works correct? I checked them and also I'm getting the same result as you are, but still I thought that my result isn't correct... :( – Eru Mar 29 '15 at 07:11
  • I've updated my answer to show how you how you might go about testing code like this, and also made some corrections to your program code. You should ensure that my corrections are correct though with further checking and testing. – softwariness Mar 29 '15 at 14:23
  • 1
    I've tested this in a symbolic math program and I get 0.0067 for P(0). Hopefully that helps with double checking. – Asad Saeeduddin Mar 29 '15 at 15:09
  • 1
    @Eru Also, the final result is the same: 5.9311e-199 (didn't bother checking to high precision). – Asad Saeeduddin Mar 29 '15 at 15:16
  • @Asad, thanks, that tallies for `P(0)` (with `k=40`, `p=5`): `0.0067379469` to 10 decimal places. – softwariness Mar 29 '15 at 15:17
  • 1
    @softwariness to be honest I don't know how to thank you for the whole help! It was very clear what you wrote and helped me a lot. Thank you again for your help and time :) – Eru Mar 30 '15 at 09:36