计算欧拉数在 C# 上给出了有限的小数位数

Question

我正在尝试计算欧拉数 (e = 2.718281828459045235)，方法是执行阶乘函数，然后使用 while 循环调用它来计算常量。问题是我声明了一个 i 变量以使循环在 i<50 的情况下工作，如下所示，它给了我以下输出：

The number of euler is: 2.7182818284590455

但是，我想让程序输出更多的小数，但是无论我怎么增加循环中的i<50条件, 它总是给出相同数量的小数。

我尝试将 while 循环的条件更改为 i<150 并且我修改了 "long" 变量到一个 "ulong" 一个，但没有任何反应，它总是给出相同的输出。知道如何改变它吗？

如有任何帮助，我们将不胜感激，谢谢！

我的代码：

using System;

namespace eulerNumber
{
    class eulerCalculator
    {

        public static double factorial(long n)
        {
            if (n == 0)
            {
                return 1;
            }

            else if (n == 1)
            {
                return 1;
            }
            else
            {
                return n * factorial(n - 1);
            }
        }
        static void Main(string[] args)
        {
            double addition = 0, i = 0;

            while(i<50)
            {
                addition = addition + (1 / factorial((long)i)); 
                i++;
            }
            Console.WriteLine("The number of euler is: " + addition);
        }
    }
}

Answer 1

一个相当简单的改进是首先从最小的项开始执行加法，因此它们的大小相似，并且处理器可以在不损失太多精度的情况下执行加法。

注：factorial(50)为50 * factorial(49)，最后两项为1/factorial(49) + 1/factorial(50)。应用一些代数得到 1/factorial(49) + 1.0/50/factorial(49)，即 (1 + 1.0/50) / factorial(49)

假设你只计算分子，并在不计算的情况下跟踪分母中出现的阶乘。现在你有两个非常好的属性：

您永远不必计算溢出的数字（如阶乘 (i)）和
无论更新方程式中存在什么舍入误差，对最终答案都无关紧要，因为它是项中的误差，会被进一步划分并变得更小

这导致以下代码：

double accumulator = 1.0;

for( int i = 50; i > 0; --i )
{
    accumulator = 1.0 + accumulator / i;
}

演示：https://rextester.com/FEZB44588

使用 .NET 的 BigInteger class 的扩展允许我们拥有更多的精度

BigInteger scale = BigInteger.Pow(10, 60);
BigInteger accumulator = scale;

for( int i = 75; i > 0; --i )
{
    accumulator = scale + accumulator / i;
}

结果（插入小数点）：

2.718281828459045235360287471352662497757247093699959574966967

Wikipedia 的前 50 位小数：

2.71828182845904523536028747135266249775724709369995...

请注意，维基百科的措辞略有错误，这不是四舍五入到小数点后 50 位的值，这些是连续序列的前 50 位小数位

Answer 2

我很无聊所以我将 e=(1+1/x)^x approach 编码成简单的 C++（抱歉不是 C# 编码器），没有任何库或直接在字符串上计算的有趣的东西...

我还将算法从二进制移植到十进制基数，代码如下：

//---------------------------------------------------------------------------
// string numbers are in format "?.??????????????"
// where n is number of digits after decimal separator
void str_addmul(char *xy,char *x,char *y,int k,int n)   // xy = x+y*k, k = 0..9
    {
    int i,a,cy;
    for (cy=0,i=n+1;i>=0;i--)
        {
        if (i==1) i--;          // skip decimal separator
        a=(x[i]-'0')+((y[i]-'0')*k)+cy;
        cy   = a/10;
        xy[i]=(a%10)+'0';
        }
    xy[1]='.';
    xy[n+2]=0;
    }
//---------------------------------------------------------------------------
// string numbers are in format "?.??????????????"
// where n is number of digits after decimal separator
void str_mul(char *xy,char *_x,char *_y,int n)  // xy = x*y
    {
    int i,j;
    // preserve original _x,_y and use temp x,y instead
    char *x=new char[n+3]; for (i=0;i<n+3;i++) x[i]=_x[i];
    char *y=new char[n+3]; for (i=0;i<n+3;i++) y[i]=_y[i];
    // xy = 0.000...000
    i=0;
    xy[i]='0'; i++;
    xy[i]='.'; i++;
    for (;i<n+2;i++) xy[i]='0';
    xy[i]=0;
    // xy = x*y
    for (i=0;i<n+2;i++)
        {
        if (i==1) i++;          // skip decimal separator
        str_addmul(xy,xy,x,y[i]-'0',n);
        // x/=10
        for (j=n+1;j>2;j--) x[j]=x[j-1];
        x[2]=x[0]; x[0]='0';
        }
    delete[] x;
    delete[] y;
    }
//---------------------------------------------------------------------------
char* compute_e(int n)      //  e = (1+1/x)^x where x is big power of 10
    {
    int i,x10,m=n+n+4;
    char* c=new char[m+3];  // ~double precision
    char* a=new char[m+3];  // ~double precision
    char* e=new char[n+3];  // target precision
    // x = 10^m
    // c = 1 + 1/x = 1.000...0001;
    i=0;
    c[i]='1'; i++;
    c[i]='.'; i++;
    for (;i<m+1;i++) c[i]='0';
    c[i]='1'; i++;
    c[i]=0;
    // c = c^x
    for (x10=0;x10<m;x10++)
        {
        str_mul(a,c,c,m);   // c^2
        str_mul(c,a,a,m);   // c^4
        str_mul(c,c,c,m);   // c^8
        str_mul(c,c,a,m);   // c^10
        }
    // e = c
    for (i=0;i<n+2;i++) e[i]=c[i]; e[i]=0;
    delete[] a;
    delete[] c;
    return e;
    }
//---------------------------------------------------------------------------

用法：

char *e=compute_e(100); // 100 is number of digits after decimal point
cout << e;  // just output the string somewhere
delete[] e; // release memory after

compute_e(100) 与参考的结果：

e(100) = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
e(ref) = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

代码很慢，因为它是在字符串和 base10 中计算的，而不是在 base2 中的整数数组上计算的，并且只使用了朴素的数学运算实现...但是在 335 ms 和 200 中仍然完成了 100 位数字数字 2612.525 ms 在我的旧 PC 上应该 比具有相同精度的迭代更快 ...

要获得 base10 算法，方程式为：

x = 10^digits
e = (1 + 1/x) ^ x

所以当你在 decadic 中写 x 和术语 (1 + 1/x) 你会得到：

x             =   1000000000000 ... 000000
c = (1 + 1/x) = 1.0000000000000 ... 000001

现在我刚刚将修改为 10ing 的电源？而不是 c = c*c 我迭代 c = c*c*c*c*c*c*c*c*c*c 就是这样......

多亏了这些东西是在以 10 为基数的字符串上计算的，不需要 ...

最后，为了获得 n 十进制数字的精度，我们必须使用 m = 2*n + 4 数字精度进行计算，并将最终结果剪切为 n 数字 ...

所以只需将这些东西移植到 C# 就可以使用字符串而不是 char* 这样你就可以摆脱 new[]/delete[] 其余的在 C# 中应该是相同的 ...

有点好奇所以我测量了一下：

[   0.640 ms] e( 10) = 2.7182818284
[   3.756 ms] e( 20) = 2.71828182845904523536
[  11.172 ms] e( 30) = 2.718281828459045235360287471352
[  25.234 ms] e( 40) = 2.7182818284590452353602874713526624977572
[  46.053 ms] e( 50) = 2.71828182845904523536028747135266249775724709369995
[  77.368 ms] e( 60) = 2.718281828459045235360287471352662497757247093699959574966967
[ 121.756 ms] e( 70) = 2.7182818284590452353602874713526624977572470936999595749669676277240766
[ 178.508 ms] e( 80) = 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759
[ 251.713 ms] e( 90) = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178
[ 347.418 ms] e(100) = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
              e(ref) = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

并使用 this 显示 ~O(n^2.8) 复杂性

这里 100 位的 e 在我的 arbnum（任意精度浮点数）上实现，二进制 (1+1/x)^x 和 Ben Voigt 方法进行比较：

//---------------------------------------------------------------------------
const int n=100;    // digits
//---------------------------------------------------------------------------
arbnum compute_e0()
    {
    arbnum c,x;
    int bit;
    const int bits =int(float(n)/0.30102999566398119521373889472449)+1;
    const int bits2=(bits<<1);
    // e=(1+1/x)^x  ... x -> +inf
    c.one(); c>>=bits; c++;  // c = 1.000...0001b = (1+1/x)          = 2^-bits + 1
    for (bit=bits;bit;bit--)        // c = c^x = c^(2^bits) = e
        {
        c*=c;
        c._normalize(bits2); // this just cut off the result to specific number of fractional bits only to speed up the computation instead you should cut of only last zeros !!!
        }
    c._normalize(bits);
    return c;
    }
//---------------------------------------------------------------------------
arbnum compute_e1()
    {
    // e = 1 + e/i ... i = { n, n-1, n-2, .... 1 }
    int i;
    arbnum e;
    const int bits =int(float(n)/0.30102999566398119521373889472449)+1;
    for (e=1.0,i=70;i;i--)
        {
        e/=i;
        e._normalize(bits);
        e++;
        }
    return e;
    }
//---------------------------------------------------------------------------

结果如下：

    reference e  = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
[   2.799 ms] e0 = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
[  10.918 ms] e1 = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

Karatsuba and approximation divider 在后台使用

计算欧拉数在 C# 上给出了有限的小数位数

Calculating the Euler's Number gives a limited number of decimals on C#

c#

math

eulers-number