将小数提高到十进制的幂？

https://stackoverflow.com/questions/429165

06-07-2019
|

题

.net框架在Math类中提供了一个为double提供动力的方法。但是根据精度要求，我需要将小数提高到十进制幂[Pow（十进制a，十进制b）]。框架是否具有这样的功能？有没有人知道有这种功能的图书馆？

解决方案

为了解决我的问题，我找到了一些扩展系列，以及我把它们用来解决方程X ^ n = e ^（n * ln x）。

// Adjust this to modify the precision
public const int ITERATIONS = 27;

// power series
public static decimal DecimalExp(decimal power)
{
    int iteration = ITERATIONS;
    decimal result = 1; 
    while (iteration > 0)
    {
        fatorial = Factorial(iteration);
        result += Pow(power, iteration) / fatorial;
        iteration--;
    }
    return result;
}

// natural logarithm series
public static decimal LogN(decimal number)
{
    decimal aux = (number - 1);
    decimal result = 0;
    int iteration = ITERATIONS;
    while (iteration > 0)
    {
        result += Pow(aux, iteration) / iteration;
        iteration--;
    }
    return result;
}

// example
void main(string[] args)
{
    decimal baseValue = 1.75M;
    decimal expValue = 1/252M;
    decimal result = DecimalExp(expValue * LogN(baseValue));
}

Pow（）和Factorial（）函数很简单，因为幂总是一个int（在de power系列中）。

其他提示

对于正整数Exponent和小数基数，这应该是最快的：

// From http://www.daimi.au.dk/~ivan/FastExpproject.pdf
// Left to Right Binary Exponentiation
public static decimal Pow(decimal x, uint y){
    decimal A = 1m;
    BitArray e = new BitArray(BitConverter.GetBytes(y));
    int t = e.Count;

    for (int i = t-1; i >= 0; --i) {
        A *= A;
        if (e[i] == true) {
            A *= x;
        }
    }
    return A;
}

这是一个C＃程序，用于手动实现Math.Pow（），其精度高于.NET基于double的实现。剪切并粘贴到linqpad中立即运行，或将.Dump（）更改为Console.WriteLines。

我已经包含了对结果的测试。测试如下：

目标=。4％pa，每日复合10 000
答案=应为10 040
如何=十进制b = 10000; for（int i = 0; i <！> lt; 365; i ++）{b * = rate;其中rate =（1.004）^（1/365）

我测试了3种速率的实现：（1）手动计算（2）Excel（3）Math.Pow

手动计算具有最高的准确度。结果是：

Manually calculated rate:   1.0000109371043837652682334292
Excel rate:                 1.000010937104383712500000M [see formula =(1.004)^(1/365)]
Math.Pow rate:              1.00001093710438

Manual - .4%pa on R10,000:  10040.000000000000000000000131 
Excel - .4%pa on R10,000:   10039.999999999806627646709094 
Math.Pow - .4%pa on R10,000:10039.999999986201948942509648

我还在那里留下了一些额外的工作 - 我曾经用它来确定最高因子可以适应ulong（= 22）。

Linqpad代码：

/*
a^b = exp(b * ln(a))
    ln(a) = log(1-x) = - x - x^2/2 - x^3/3 - ...   (where |x| < 1)
        x: a = 1-x    =>   x = 1-a = 1 - 1.004 = -.004
    y = b * ln(a)
    exp(y) = 1 + y + y^2/2 + x^3/3! + y^4/4! + y^5/5! + ...
        n! = 1 * 2 * ... * n        
*/

/*
//
// Example: .4%pa on R10,000 with daily compounding
//

Manually calculated rate:   1.0000109371043837652682334292
Excel rate:                 1.000010937104383712500000M =(1.004)^(1/365)
Math.Pow rate:              1.00001093710438

Manual - .4%pa on R10,000:  10040.000000000000000000000131 
Excel - .4%pa on R10,000:   10039.999999999806627646709094 
Math.Pow - .4%pa on R10,000:10039.999999986201948942509648 

*/

static uint _LOOPS = 10;    // Max = 22, no improvement in accuracy after 10 in this example scenario
//  8: 1.0000109371043837652682333497
//  9: 1.0000109371043837652682334295
// 10: 1.0000109371043837652682334292
// ...
// 21: 1.0000109371043837652682334292
// 22: 1.0000109371043837652682334292

// http://www.daimi.au.dk/~ivan/FastExpproject.pdf
// Left to Right Binary Exponentiation
public static decimal Pow(decimal x, uint y)
{
    if (y == 1)
        return x;

    decimal A = 1m;
    BitArray e = new BitArray(BitConverter.GetBytes(y));
    int t = e.Count;

    for (int i = t-1; i >= 0; --i) {
        A *= A;
        if (e[i] == true) {
            A *= x;
        }
    }
    return A;
}

// http://stackoverflow.com/questions/429165/raising-a-decimal-to-a-power-of-decimal
// natural logarithm series
public static decimal ln(decimal a)
{
    /*
    ln(a) = log(1-x) = - x - x^2/2 - x^3/3 - ...   (where |x| < 1)
        x: a = 1-x    =>   x = 1-a = 1 - 1.004 = -.004
    */
    decimal x = 1 - a;
    if (Math.Abs(x) >= 1)
        throw new Exception("must be 0 < a < 2");

    decimal result = 0;
    uint iteration = _LOOPS;
    while (iteration > 0)
    {
        result -= Pow(x, iteration) / iteration;
        iteration--;
    }
    return result;
}

public static ulong[] Fact = new ulong[] {
    1L,
    1L * 2,
    1L * 2 * 3,
    1L * 2 * 3 * 4,
    1L * 2 * 3 * 4 * 5,
    1L * 2 * 3 * 4 * 5 * 6,
    1L * 2 * 3 * 4 * 5 * 6 * 7,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20,
    14197454024290336768L, //1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 * 21,        // NOTE: Overflow during compilation
    17196083355034583040L, //1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 * 21 * 22    // NOTE: Overflow during compilation
};

// http://stackoverflow.com/questions/429165/raising-a-decimal-to-a-power-of-decimal
// power series
public static decimal exp(decimal y)
{
    /*
    exp(y) = 1 + y + y^2/2 + x^3/3! + y^4/4! + y^5/5! + ...
    */

    uint iteration = _LOOPS;
    decimal result = 1; 
    while (iteration > 0)
    {
        //uint fatorial = Factorial(iteration);
        ulong fatorial = Fact[iteration-1];
        result += (Pow(y, iteration) / fatorial);
        iteration--;
    }
    return result;
}

void Main()
{   
    decimal a = 1.004M;
    decimal b = 1/365M;

    decimal _ln = ln(a);
    decimal y = b * _ln;
    decimal result = exp(y);
    result.Dump("Manual rate");

    decimal excel = 1.000010937104383712500000M;    // =(1.004)^(1/365)
    excel.Dump("Excel rate");


    decimal m = (decimal)Math.Pow((double)a,(double)b);
    m.Dump("Math.Pow rate");

    //(result - excel).Dump("Diff: Manual - Excel");
    //(m - excel).Dump("Diff: Math.Pow - Excel");

    var f = new DateTime(2013,1,1);
    var t = new DateTime(2014,1,1);
    Test(f, t, 10000, result, "Manual - .4%pa on R10,000");
    Test(f, t, 10000, excel, "Excel - .4%pa on R10,000");
    Test(f, t, 10000, m, "Math.Pow - .4%pa on R10,000");
}

decimal Test(DateTime f, DateTime t, decimal balance, decimal rate, string whichRate)
{
    int numInterveningDays = (t.Date - f.Date).Days;
    var value = balance;
    for (int i = 0; i < numInterveningDays; ++i)
    {
        value *= rate;
    }
    value.Dump(whichRate);
    return value - balance;
}

/*

// Other workings:

//
// Determine maximum Factorial for use in ln(a)
//

ulong max    =  9,223,372,036,854,775,807 * 2   // see http://msdn.microsoft.com/en-us/library/ctetwysk.aspx
Factorial 21 = 14,197,454,024,290,336,768
Factorial 22 = 17,196,083,355,034,583,040
Factorial 23 = 8,128,291,617,894,825,984 (Overflow)

public static uint Factorial_uint(uint i)
{
    // n! = 1 * 2 * ... * n
    uint n = i;
    while (--i > 1)
    {
        n *= i;
    }
    return n;
}

public static ulong Factorial_ulong(uint i)
{
    // n! = 1 * 2 * ... * n
    ulong n = i;
    while (--i > 1)
    {
        n *= i;
    }
    return n;
}

void Main()
{
    // Check max ulong Factorial
    ulong prev = 0;
    for (uint i = 1; i < 24; ++i)
    {
        ulong cur = Factorial_ulong(i);
        cur.Dump(i.ToString());
        if (cur < prev)
        {
            throw new Exception("Overflow");
        }
        prev = cur;
    }
}
*/

我认为这很大程度上取决于您计划插入的数量。如果'a'和'b'不是'nice'数字，那么您可能会获得一个无法终止的值，这是无法存储的如果C＃BigDecimal像Java BigDecimal一样行为，它可能会在这种情况下引发异常。

你确定你真的想这样做吗？ decimal乘法比double慢大约40倍，所以我希望小数Math.Pow()几乎无法使用。

但是，如果您只期望整数幂，我建议您使用已在SO上讨论过的基于整数的幂算法。

许可以下： CC-BY-SA 和归因

不隶属于 StackOverflow