BigInteger/BigRational 转换为 double 和返回的问题

BigInteger/BigRational problems with converting to double and back

也许我不明白如何使用这两个类型:BigInteger/BigRational,但总的来说我想实现这两个方程式:

这是我的数据:n=235,K = 40,这个小 p(实际上称为 rho)是 5。一开始问题出在 Power 函数上:结果非常非常大 -所以这就是我使用 BigInteger 库的原因。但后来我意识到会有一个除法,结果将是一个双精度类型的数字 - 所以我改用 BigRational 库。

这是我的代码:

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));
    }

主要问题是它没有 return 任何“正常”值,例如 0.3 或 0.03。打印到控制台的结果是一个很长的数字(其中有 1200 位数字)...

有人可以看看我的代码并帮助我解决问题并能够通过代码解决这个方程式。谢谢

Console.WriteLine(pNResult); 在后台调用 BigRational.ToString(),它以 numerator/denominator.

的形式打印数字

考虑到本例中分子和分母都很大,很容易在输出中遗漏 /

BigRational 支持转换为 decimaldouble。不过,在这种情况下,结果太小而无法放入 decimal 中。转换为 double,得到结果 7.89682541396914E-177.

如果您需要更高的精度,您需要自定义转换为十进制格式的字符串,例如 this Whosebug answer.

使用该自定义转换例程将结果精确到 1000 位小数 -

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

- 给出的结果为:

0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000078968254139691306770128897137459492828971170349380336740935269651539684650525033676003134593283361305530675112470528408219177025044254116462798561450442318290046626248451723040397770263675109107145461310779641705093156106311143727608208629473359566457461384474633112850335950017209558136575135801388668687571284492241030561019606955986265585636660304889792027894460104216176719717671500843399685686146432982358441225578366059001576682388503227237202077881334695352338638383337717103303153521108812750644260562351186866587629456292506971252525125976755540274041651740194108430555751648707933592643410475214924394223640168857340953563111097979394441303100701008120008166339365089771585037880235325673143152814510586536335380671360865230428857049658368242543653234599817430185879648427434216378356518036776477170130227628307039


要检查您的计算代码是否正确运行,您可以为不同的函数添加单元测试(FactorialSeriesP 本身的计算)。

此处实用的方法是手动计算 knp 的某些小值的结果,并检查您的函数是否计算出相同的结果。

如果您使用 Visual Studio,则可以使用 this MSDN page 作为创建单元测试项目的起点。请注意,被测函数必须对单元测试项目可见,并且您的单元测试项目需要将引用添加到您正在进行计算的现有项目中,如 link 中所述.

从最容易检查的 Factorial 开始,您可以添加这样的测试:

[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));
    }
}

您问题中的代码通过了该测试。

然后您可以为 Series 函数添加一个测试方法:

[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));
}

这表明您的代码存在一些问题:

  • n 实际上不应该是这个函数的参数,因为在这种情况下 n 不是函数 P 的参数,而实际上只是 "index-of-summation". n 在此函数中的局部值由您的 i 变量表示。
  • 这意味着您的 Factorial(n) 呼叫需要更改为 Factorial(i)
  • 循环也是差一,因为求和的 Sigma 符号包含 Sigma 顶部的数字,所以你应该有 <= fin (或者你也可以写成这只是 < k).

这是更新后的 Series 函数:

// 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;
}

要测试 P(n) 计算,您可以将其移出到它自己的函数中进行测试(我在这里将其称为 ComputeP):

[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));
}

这在计算 P(n) 时出现了问题 - 您在其中对 double 进行了强制转换,但不应该出现(结果不准确 - 您需要保留所有BigRational 中的中间结果)。不需要转换,所以只需删除它就可以解决这个问题。

这是更新后的 ComputeP 函数:

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;
}

为了避免混淆,这里是整个更新的计算程序:

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 
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();
    }
}

下面是整个单元测试程序:

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));
    }
}

顺便说一下,P(n)P(n) 结果现在是 npk 输入值的更新程序:

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000593109980769066916025972569398424267669807629726200017375290861590898269902277869938365969961320969473356001666906480007119114830921839913623591124192047955091318951831902550404167336054683697071654765071519020060437129398945035521954738463786221029427589397688847246112810536958194364039693387170592425527136243952416704526069736811587380688876091926255908361275575249492845970903676492429684929779402600032481018886875698972533534890841796034626337674846620462046294537488580901129338625628349474358946962065227890599744775562637784553656488649841148591533557896418988044457914999854241038974478576578909626765823565817758792682480009619613438867365912697996527957775248350987801430141776875171808382272960426476953742528769626555642957093028553993908356226007570404005591174451216846471710162760343


注意:您应该将更多手工检查的结果添加到单元测试中,并检查我在这里将代数解释为代码的任何工作,以确保这是正确的。