计算欧拉数在 C# 上给出了有限的小数位数
Calculating the Euler's Number gives a limited number of decimals on C#
我正在尝试计算欧拉数 (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);
}
}
}
一个相当简单的改进是首先从最小的项开始执行加法,因此它们的大小相似,并且处理器可以在不损失太多精度的情况下执行加法。
注: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 位小数位
我很无聊所以我将 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 在后台使用
我正在尝试计算欧拉数 (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);
}
}
}
一个相当简单的改进是首先从最小的项开始执行加法,因此它们的大小相似,并且处理器可以在不损失太多精度的情况下执行加法。
注: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 位小数位
我很无聊所以我将 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 在后台使用