Not Only Algorithm,不仅仅是算法,关注数学、算法、数据结构、程序员笔试面试以及一切涉及计算机编程之美的内容 。。
你的位置:NoAlGo博客 » 数学 » 

计算自然常数e

e是数学中最重要的常数之一,作为自然对数的底数,应用非常广泛。e的值约为2.718281828,但e的精确值究竟是多少呢?本文将使用一种最简单的方法求出e的前若干位的值。

泰勒展开

e的计算方法有很多种,其中最简单的是泰勒展开。泰勒公式是一个用函数在某点的信息描述其附近取值的公式,其可以表示成

f(x) = f(a)/0! + f'(a)/1!*(x-a) + f''(a)/2!*(x-a)^2 + ... + f(n)(a)/n!*(x-a)^n + ...

令f(x) = e^x,且令x = 1, a = 0,即可得到

e = 1 + 1/1! + 1/2! + 1/3! + ... + 1/n! + ...

于是可以利用这个公式求取e的值,但这里随着n的值越来越大,要提高计算的精度,只能够自己实现高精度计算过程。

本文介绍的方法参考Pascal Sebah在1999年7月实现的版本,关于更多的常数计算方法参考这里:Numbers, constants and computation

代码实现

根据泰勒展开,实现的过程非常直观。

使用数组表示高精度小数,根据e的特性,这里的高精度小数中第一个数字表示整数部分,其余部分表示小数部分,这样使得高精度的加法和除法都非常简单。注意,这里的高精度仅仅满足了计算e所需要的过程,并不对所有的高精度运算成立。
高精度小数使用10000进制,可以使用整数进行初始化,并进行两个高精度小数的加法和高精度小数除以低精度整数的运算。

然后就是计算的主函数,直接循环直到新添加的项为0即可。因为这里高精度小数的精度有限,于是随着n的增大,后面的项会越来越小,最终为0。

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>

const int BASE = 10000;      //高精度数字进制
const int NUM_DIGITS = 10000;//高精度数字长度

class BigReal {
public:
	BigReal (int val = 0) { //整数初始化
		n = 1 + static_cast<int>(NUM_DIGITS / log10(BASE));
		data = new int[n];
		memset(data, 0, n * sizeof(int));
		data[0] = val;
	}
	~BigReal() {
		delete(data);
	}

	bool IsZero() { //是否为0
		for (int i = 0; i < n; i++)
			if (data[i]) return false;
		return true;
	}

	void Add(const BigReal &x) { //加另一个大数
		for (int i = n - 1, carry = 0; i >= 0; i--) {
			data[i] += x.data[i] + carry;
			if (data[i] < BASE) {
				carry = 0;
			} else {
				data[i] -= BASE;
				carry = 1;
			}
		}
	}

	void Div(int d) { //除小整数
		for (int i = 0, carry = 0, q; i < n; i++) {
			data[i] = data[i] + carry * BASE;
			q = data[i] / d;
			carry = data[i] - q * d;
			data[i] = q;
		}
	}

	void Print() { //打印
		printf("%d.", data[0]);
		for (int i = 1; i < n; i++) {
			printf("%.4d", data[i]);
			if (i % 15 == 0) {
				printf("%8d\n", i * 4);
			}
		}
	}

private:
	int n, *data; //长度,指针
};

int main() {
	BigReal e(1), k(1);
	for (int i = 1; !k.IsZero(); i++) {
		// e = 1 + 1/1! + 1/2! + 1/3! + ...
		k.Div(i);
		e.Add(k);
	}
	e.Print();
}

计算结果

计算的结果以长度60为单位进行显示,下面截取其前面360位的计算结果进行显示:

2.718281828459045235360287471352662497757247093699959574966967      60
627724076630353547594571382178525166427427466391932003059921     120
817413596629043572900334295260595630738132328627943490763233     180
829880753195251019011573834187930702154089149934884167509244     240
761460668082264800168477411853742345442437107539077744992069     300
551702761838606261331384583000752044933826560297606737113200     360

另一种方法

这里介绍Xavier Gourdon在1999年7月实现的另一种方法,代码非常简短,可以打印出e的前面9000位,并且可以通过把9009改成更大的数字来计算更多的位。

#include <stdio.h>
int main() {
	int N=9009,n=N,a[9009],x;
	while(--n)a[n]=1+1/n;
	for(;N>9;printf("%d",x))
		for(n=N--;--n;a[n]=x%n,x=10*a[n-1]+x/n);
}

下面截取其在gcc编译后运行结果的前一部分进行展示:

2718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200

上一篇: 下一篇:
  1. 你个逗逼,a 是int型的,a 始终是1.你自己都没弄懂这个程序,胡乱贴个屁.

  2. x是使用gcc默认初值0。数组a并非全1,下标1的位置值为2。自己去跑一遍就知道了。

  3. 男生寝室忽然停电,同学们大喊:来电,来电!不一会,电果然来了。 神了男生们一起欢呼起来:来女人,来女人!管理寝室的阿姨来了,大吼一声:都给老娘闭嘴。

  4. 学好IT好就业选硅谷IT,学技能拿文凭事半功倍,紧跟专业老师一起冲浪IT行业。我们有建设学习型专业师资团队,教师领跑学生紧跟随其后。(QQ:800015777 电话0754—88989555)

  5. 从百度点进来的,支持一下,希望站长您多出一些好文章。

我的博客

NoAlGo头像编程这件小事牵扯到太多的知识,很容易知其然而不知其所以然,但真正了不起的程序员对自己程序的每一个字节都了如指掌,要立足基础理论,努力提升自我的专业修养。

站内搜索

最新评论