用马青公式计算π的c程序
/*一个计算π的c程序
======================================
由于大多数计算机内置算法有一定精度限制,你想计算那么多位就会产生问题。
这里有一个c程序,允许计算要多少有多少。但马青公式在试图计算亿位时不理想。
下面就是这个程序。超过万位时此程序就不能胜任了。而Jason Chen的程序可算到10万位!
目前最快的方法是Chudnovsky、Ramanujan和金田康正的高斯-勒让德算法。
*/
///
/****************************************/
/* Compute pi to arbitrary precision */
/* Author Roy Williams February 1994 */
/* Uses Machin's formula... */
/* pi/4 = 4 arctan(1/5) - arctan(1/239) */
/****************************************/
/* compile with cc -O -o pi pi.c */
/* run as "pi 1000" */
/****************************************/
/* The last few digits may be wrong.......... */
#include <stdio.h>
#include <malloc.h>
#define BASE 10000
int nblock;
int *tot;
int *t1;
int *t2;
int *t3;
void copy(int *result, int *from)
{
int i;
for(i=0; i<nblock; i++)
result[i] = from[i];
}
int zero(int *result)
{
int i;
for(i=0; i<nblock; i++)
if(result[i])
return 0;
return 1;
}
void add(int *result, int *increm)
{
int i;
for(i=nblock-1; i>=0; i--){
result[i] += increm[i];
if(result[i] >= BASE){
result[i] -= BASE;
result[i-1]++;
}
}
}
void sub(int *result, int *decrem)
{
int i;
for(i=nblock-1; i>=0; i--){
result[i] -= decrem[i];
if(result[i] < 0){
result[i] += BASE;
result[i-1]--;
}
}
}
void mult(int *result, int factor)
{
int i, carry = 0;
for(i=nblock-1; i>=0; i--){
result[i] *= factor;
result[i] += carry;
carry = result[i]/BASE;
result[i] %= BASE;
}
}
void div(int *result, int denom)
{
int i, carry = 0;
for(i=0; i<nblock; i++){
result[i] += carry*BASE;
carry = result[i] % denom;
result[i] /= denom;
}
}
void set(int *result, int rhs)
{
int i;
for(i=0; i<nblock; i++)
result[i] = 0;
result[0] = rhs;
}
void print(int *result)
{
int i, k;
char s[10];
printf("%1d./n", result[0]);
for(i=1; i<nblock; i++){
sprintf(s, "%4d ", result[i]);
for(k=0; k<5; k++)
if(s[k] == ' ') s[k] = '0';
printf("%c%c%c%c", s[0], s[1], s[2], s[3]);
if(i%15 == 0) printf("/n");
}
printf("/n");
}
void arctan(int *result, int *w1, int *w2, int denom, int onestep)
{
int denom2 = denom*denom;
int k = 1;
set(result, 1);
div(result, denom);
copy(w1, result);
do{
if(onestep)
div(w1, denom2);
else {
div(w1, denom);
div(w1, denom);
}
copy(w2, w1);
div(w2, 2*k+1);
if(k%2)
sub(result, w2);
else
add(result, w2);
k++;
} while(!zero(w2));
}
//void main(int argc, char *argv[])
void test()
{
int ndigit = 10000; /*位数可以自己试改,如:1000....*/
/*
if(argc == 2)
ndigit = atoi(argv[1]);
else {
fprintf(stderr, "Usage: %s ndigit/n", argv[0]);
fprintf(stderr, "(Assuming 10000 digits)/n");
}
*/
if(ndigit < 20) ndigit = 20;
nblock = ndigit/4;
tot = (int *)malloc(nblock*sizeof(int));
t1 = (int *)malloc(nblock*sizeof(int));
t2 = (int *)malloc(nblock*sizeof(int));
t3 = (int *)malloc(nblock*sizeof(int));
if(!tot || !t1 || !t2 || !t3){
fprintf(stderr, "Not enough memory/n");
exit(1);
}
arctan(tot, t1, t2, 5, 1);
mult(tot, 4);
arctan(t3, t1, t2, 239, 2);
sub(tot, t3);
mult(tot, 4);
print(tot);
}
/
用马青公式计算π的c程序相关推荐
- 用马青公式计算圆周率,Python语言
马青公式: 简单实现 import times_time = time.time() #定义计算的位数 num = 100000 #多计算10位,以防出错 num1 = num + 10 #定义计算小 ...
- Python利用马青公式计算圆周率Π并写入文件
一.什么是马青公式 马青公式由英国天文学教授约翰·马青(John Machin ,1686 –1751)于1706年发现,他利用这个公式计算到了100位的圆周率. 马青公式每计算一项可 ...
- c语言马青公式计算圆周率,数学圆周率计算马青公式π/4=4arctan1/5-arctan1/239如何得出的?...
共回答了16个问题采纳率:87.5% 设 x = arctan A tan x = A tan 2x = (2 tan x) / (1 - tan^2 x) tan 2x = (2A) / (1 - ...
- Machin(梅钦/马青)公式计算圆周率π
作为一名计算机的初学者,因为老师作业的要求去完成圆周率的计算,因此突然产生了兴趣,想尝试自己用梅钦公式来完成这个任务.上网找了一些资料,也看见了短短几行就完成任务的代码,实在佩服,不过那样的代码实在对 ...
- 利用Java的BigDecimal与马青公式精确计算π后10000位,
首先给出公式如下: π=16arctan1/5−4arctan1/239: 即是 π=16×(1/(1×5)−1/(3×5的3次方)+1/(5×5的5次方)-)−4×(1/(1×239)−1/(3×2 ...
- 【Python、数学】计算任意位数的圆周率π(马青公式)
1. 公式准备 计算准确圆周率的马青公式: 对反正切进行级数展开: 就可以得到 π = 16(1/5 - 1/3/5^3 + 1/5/5^5 - ...) - 4(1/239 - 1/3/239^3 ...
- 利用马青公式输出π的后任意位数字
马青公式 π=16arctan15−4arctan1239π=16arctan15−4arctan1239\pi = 16arctan \frac{1}{5} - 4arctan \frac{1}{2 ...
- python求圆周率马青公式_Python 实现丘德诺夫斯基(Chudnovsky)法計算高精度圓周率...
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 源码:(需要将 @ 替换成 ASCII空格) # -*- coding: UTF-8 -*- # 丘德诺夫斯基法計算高精度圓周率程序 # Calculat ...
- matlab 流程计算方法,吸波材料LLG公式计算复磁导率的过程及matlab程序
看到一篇paper,利用Landau-Lifshitz-Gilbert 公式计算片状颗粒的复磁导率.(JAP 107,033913, 2010) http://scitation.aip.org/co ...
最新文章
- 5s的app显示无法连接服务器,苹果5s无法连接app store解决方法汇总
- 频率分布直方图组距如何确定_吃透教材理解教参,《直方图》教学反思
- mysql 去重求总数_Mysql获取去重后的总数
- HTML基础 --- HTML简介
- python和java一样吗-三分钟看懂Python和Java的区别
- 高性能Javascript:高效的数据访问
- java form表单传参_表单Form传参数
- 【博客话题】感谢您,我的老师
- windows CMD生成文件夹树状图(tree)命令(以图形显示驱动器或路径的文件夹结构)
- 黑色幽默:“新知青”电影《走着瞧》首映
- Beta冲刺——星期三
- js中构造函数与普通函数的区别
- python stack使用_python inspect.stack() 的简单使用
- 【学堂在线数据挖掘:理论方法笔记】第八天(4.2)
- python爬去微博签到数据_GitHub - fs6/weiboSpider: 新浪微博爬虫,用python爬取新浪微博数据...
- 论现场跟客户演示软件产品
- 移动应用开发者的阶级状况:多数是无产阶级
- View 的各种知识
- docker日志显示时间时区错误,时区UST问题/群晖docker日志时间不正确 寻找解答过程
- 第三届“网鼎杯”官方资格赛圆满结束,问鼎之战即将开启!