Jason Chen的计算圆周率的程序

最新推荐文章于 2021-05-22 13:46:17 发布
最新推荐文章于 2021-05-22 13:46:17 发布
阅读量1k 收藏
点赞数
文章标签： fp file null div 语言 include


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include <time.h>
#include <conio.h>
#include <io.h>

int     *arctg5, *arctg239, *tmp;
int     DecLen, BinLen;
int     x, n, sign, NonZeroPtr;
int     Step;
char    fn_status[]   = "~pi_sts.___";
char    fn_arctg5[]   = "~atg5.___";
char    fn_arctg239[] = "~atg239.___";
char    fn_tmp[]      = "~tmp.___";
time_t TotalTime, time1, time2;

 

void __cdecl FirstDiv(int *arctg_array)
{
    __asm {
        mov esi, arctg_array
        mov dword ptr [esi], 1
        xor edx, edx
        mov ebx, x
        mov ecx, BinLen
fd1:    mov eax, [esi]
        div ebx
        mov [esi], eax
        add esi, 4
        loop fd1
        mov esi, arctg_array
        mov edi, tmp
        mov ecx, BinLen
        pushf
        cld
        rep movsd
        popf
    }
}

 

void __cdecl arctgx(int *arctg_array)
{
    int         NonZeroBytePtr;
    int         key;
    FILE        *fp;

    for (;NonZeroPtr<BinLen;) {
        NonZeroBytePtr = NonZeroPtr * 4;

        if (_kbhit()) {
            key=_getch();
            if (key==0 || key==0xe0) _getch();
            if (key=='p') {
                printf("Swap data to disk .../n");
                if (x==25)
                    Step = 1;
                else
                    Step = 2;
                time(&time2);
                TotalTime += time2 - time1;
                if ((fp=fopen(fn_status,"wt"))==NULL) {
                    printf("Create file %s error!!/n", fn_status);
                    exit(0);
                }
                fprintf(fp, "%d %d %d %d %d %d %d/n",
                    Step, DecLen, BinLen, n, sign, NonZeroPtr, TotalTime);
                fclose(fp);
                if ((fp=fopen(fn_arctg5,"wb"))==NULL) {
                    printf("Create file %s error!!/n", fn_arctg5);
                    exit(0);
                }
                if (fwrite(arctg5, 4, BinLen, fp) != (unsigned)BinLen) {
                    printf("Write file %s error!!/n", fn_arctg5);
                    exit(0);
                }
                fclose(fp);
                if ((fp=fopen(fn_arctg239,"wb"))==NULL) {
                    printf("Create file %s error!!/n", fn_arctg239);
                    exit(0);
                }
                if (fwrite(arctg239, 4, BinLen, fp) != (unsigned)BinLen) {
                    printf("Write file %s error!!/n", fn_arctg239);
                    exit(0);
                }
                fclose(fp);
                if ((fp=fopen(fn_tmp,"wb"))==NULL) {
                    printf("Create file %s error!!/n", fn_tmp);
                    exit(0);
                }
                if (fwrite(tmp, 4, BinLen, fp) != (unsigned)BinLen) {
                    printf("Write file %s error!!/n", fn_tmp);
                    exit(0);
                }
                fclose(fp);
                printf("Exit./n");
                exit(0);
            }
        }

        __asm {
            mov esi, tmp
            add esi, NonZeroBytePtr
            xor edx, edx
            mov ebx, x
            mov ecx, BinLen
            sub ecx, NonZeroPtr
arctg1:     mov eax, [esi]
            div ebx
            mov [esi], eax
            add esi, 4
            loop arctg1

            cmp sign, 1
            jne sub_

            mov esi, tmp
            add esi, NonZeroBytePtr
            mov edi, arctg_array
            add edi, NonZeroBytePtr
            xor edx, edx
            mov ebx, n
            mov ecx, BinLen
            sub ecx, NonZeroPtr
add_1:      mov eax, [esi]
            div ebx
            add [edi], eax
            adc dword ptr [edi-4], 0
            jnc add_3
            push edi
            sub edi, 4
add_2:      sub edi, 4
            add dword ptr [edi], 1
            jc add_2
            pop edi
add_3:      add esi, 4
            add edi, 4
            loop add_1
            jmp adj_var

sub_:       mov esi, tmp
            add esi, NonZeroBytePtr
            mov edi, arctg_array
            add edi, NonZeroBytePtr
            xor edx, edx
            mov ebx, n
            mov ecx, BinLen
            sub ecx, NonZeroPtr
sub_1:      mov eax, [esi]
            div ebx
            sub [edi], eax
            sbb dword ptr [edi-4], 0
            jnc sub_3
            push edi
            sub edi, 4
sub_2:      sub edi, 4
            sub dword ptr [edi], 1
            jc sub_2
            pop edi
sub_3:      add esi, 4
            add edi, 4
            loop sub_1

adj_var:    add n, 2
            neg sign
            mov esi, tmp
            add esi, NonZeroBytePtr
            cmp dword ptr [esi], 0
            jne adj_var_ok
            inc NonZeroPtr
adj_var_ok:
           }
      }
}

void __cdecl mul_array(int *array, int multiplicator)
{
    __asm {
        mov esi, BinLen
        dec esi
        shl esi, 2
        add esi, array
        mov ecx, BinLen
        mov ebx, multiplicator
        xor edi, edi
mul1:   mov eax, [esi]
        mul ebx
        add eax, edi
        adc edx, 0
        mov [esi], eax
        mov edi, edx
        sub esi, 4
        loop mul1
        mov [esi], edx
    }
}

 

void __cdecl sub2array(int *array1, int *array2)
{
    __asm {
        mov esi, array1
        mov edi, array2
        mov ecx, BinLen
        dec ecx
sub1:   mov eax, [edi+ecx*4]
        sbb [esi+ecx*4], eax
        loop sub1
    }
}

 

void test(void)
{
    struct tm   *ts;
    FILE        *pi, *fp;
    int         i, tail, p10tail;

    printf("/nProgram to compute PI, by Jason Chen, May 1999./n");
    printf("        Dec Length               Time(h:m:s)/n");
    printf("            20000                 00:00:07/n");
    printf("           100000                 00:02:54/n");
    printf(" (Running on PII-233, 128MB, Win98 Dos mode)/n");
    printf("        Homepage: http://www.jason314.com/n");
    printf("           Email: jason@szonline.net/n/n");

    if ((fp=fopen(fn_status,"rt"))==NULL) {
        printf("Decimal length位数 = ");
        scanf("%d", &DecLen);
        if (DecLen < 100) DecLen = 100;
        BinLen = (int)(DecLen / log10(2) / 32) + 2;
        Step = 0;
        TotalTime = 0;
    }
    else {
        printf("Reading previous data .../n");
        fscanf(fp, "%d %d %d %d %d %d %d",
            &Step, &DecLen, &BinLen, &n, &sign, &NonZeroPtr, &TotalTime);
        fclose(fp);
        if (Step*DecLen*BinLen*n*sign*NonZeroPtr*TotalTime == 0) {
            printf("File %s error !!/nExit!/n", fn_status);
            exit(0);
        }
    }

    arctg5   =(int *) calloc(BinLen, 4);
    arctg239 =(int *) calloc(BinLen, 4);
    tmp      =(int *) calloc(BinLen, 4);
    if (arctg5==NULL || arctg239==NULL || tmp==NULL) {
        printf("Not enough memory !!/n");
        exit(0);
    }

    if (Step == 0) {
        memset(arctg5, 0, BinLen*4);
        memset(arctg239, 0, BinLen*4);
        memset(tmp, 0, BinLen*4);
    }
    else {
        if ((fp=fopen(fn_arctg5,"rb"))==NULL) {
            printf("Open file %s error!!/n", fn_arctg5);
            exit(0);
        }
        if (fread(arctg5, 4, BinLen, fp) != (unsigned)BinLen) {
            printf("File %s error!!/n", fn_arctg5);
            exit(0);
        }
        fclose(fp);
        if ((fp=fopen(fn_arctg239,"rb"))==NULL) {
            printf("Open file %s error!!/n", fn_arctg239);
            exit(0);
        }
        if (fread(arctg239, 4, BinLen, fp) != (unsigned)BinLen) {
            printf("File %s error!!/n", fn_arctg239);
            exit(0);
        }
        fclose(fp);
        if ((fp=fopen(fn_tmp,"rb"))==NULL) {
            printf("Open file %s error!!/n", fn_tmp);
            exit(0);
        }
        if (fread(tmp, 4, BinLen, fp) != (unsigned)BinLen) {
            printf("File %s error!!/n", fn_tmp);
            exit(0);
        }
        fclose(fp);
    }

    printf("Working ....../n");
    printf("稍候......");
    time(&time1);

    if (Step == 0) {
        x = 5;
        FirstDiv(arctg5);
        x = x * x;
        n = 3;
        sign = -1;
        NonZeroPtr = 1;
        arctgx(arctg5);
    }
    else if (Step == 1) {
        x = 5 * 5;
        arctgx(arctg5);
    }

    if (Step == 0 || Step == 1) {
        x = 239;
        FirstDiv(arctg239);
        x = x * x;
        n = 3;
        sign = -1;
        NonZeroPtr = 1;
        arctgx(arctg239);
    }
    else {
        x = 239 * 239;
        arctgx(arctg239);
    }

    mul_array(arctg5, 16);
    mul_array(arctg239, 4);
    sub2array(arctg5, arctg239);

    if ((pi=fopen("pi.txt","wt"))==NULL) {
        printf("Create file pi.txt error!!/n");
        exit(0);
    }

    printf("正在将结果写入文本文件: pi.txt .../n");
    fprintf(pi, "%d", arctg5[0]);
    for (i=1; i<=DecLen/9; i++) {
        arctg5[0] = 0;
        mul_array(arctg5, 1000000000);
        fprintf(pi, "%09d", arctg5[0]);
    }
    tail = DecLen % 9;
    p10tail = 1;
    for (i=1;i<=tail;i++) p10tail *= 10;
    arctg5[0] = 0;
    mul_array(arctg5, p10tail);
    fprintf(pi, "%0*d", tail, arctg5[0]);
    fclose(pi);

    time(&time2);
    TotalTime += time2 - time1;
    printf("已经完成了 !!/n");
    ts = gmtime(&TotalTime);
    ts->tm_mon --;
    ts->tm_mday = ts->tm_mday - 1 + ts->tm_mon * 31;
    printf("耗时 : ");
    if (ts->tm_mday > 0) printf("%d Day(s) ", ts->tm_mday);
    printf("%02d:%02d:%02d/n", ts->tm_hour, ts->tm_min, ts->tm_sec);

    if (_unlink(fn_status) == 0) {
        _unlink(fn_arctg5);
        _unlink(fn_arctg239);
        _unlink(fn_tmp);
    }
}

//
//此程序代码可复制粘贴到vc++6.0调试运行。该程序为Jason Chen于1999年编的。我用过国内流行的
//Turbo c 2.0来运行它。有些标识符，它不认识，无法编译。而在vc++6.0中则顺利通过。
//要读懂此程序，需一定的功力。但照搬是没有问题的。程序用的是马青Machin公式:
// π=16arctg （1/5） - 4arctg（ 1/239 ）。

//arctgx=x-(1/3)(x^3)+(1/5)x^5-(1/7)x^7+(1/9)x^9-.....    Gregory's series 级数。

//用它来 计算π值到10万位在耗时上是没有什么问题的。快得很！当然还有更先进的方法。如日本的金田教授的方法等等。不管快慢如何，问题是这么庞大的位数如何存放在电脑内存中？而且还要保持其准确度！用过此程序后再仔细研究其中的c语言和汇编语言吧。显然,把数值变量转换为字符数组变量,然后逐位计算。并且还用到文件读写管理的方法来记录结果。大致就是这个思路。
// ,美国某大学的一位教授就是这么提示给他的学生的 。这是个关键。
//程序的开头说明部分和程序中的中文是我修改的。我已经运行多次。估计使用它时不会有大问题。
//在VC++中运行调试C语言程序的方法，在此都作了详细的说明。运行时出现命令提示符，
//在等号后面键入你想要的圆周率的位数，回车。等候。完成后电脑自动形成一个文本文件：Pi.txt.存
//放在所编出的程序文件的同一文件夹中。打开它即可看到π值文本。计算π值至10万位所耗时间约为25秒
//（最近出/的电脑）。调试成功后就会产生一个可脱离开发环境而独立执行的自己命名的.exe程序。

 

wanchuanhua
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
计算圆周率的经典程序
04-29
只用两行代码就可将圆周率精确到小数点后几千位 至于到精确到多少,由你自己设定,很经典的小程序
高级语言程序设计__第四章例子 用C语言实现计算π的值（近似值）
华工的专栏
10-21 8716
用下列公式计算π的值：π/4= 1- 1/3 + 1/5 -1/7+... ±1/n , 精度要求为 |1/n|&lt; 10的-4次方 分析 ： 设计变量s存放每项的正负符号， n 存放每项的分母， t 存放每项 在do-while 循环体中首先计算出当前一项 t=s/n, 将t 存放到结果pi中 =pi+t， 然后计算下一项符号s= s*-...
c语言编程计算π的值
li_641233642的博客
07-14 9260
根据下面的公式计算π的值 π=4-4/3+4/5-4/7+4/9-4/11+… 请打印出一个表格来显示：计算出来π的近似值。 #include <stdio.h> #include <math.h> double get_pi(double n) //返回值为double类型,不能写int...
计算圆周率 .exe
03-06
可以快速的计算指定精度的圆周率 输入100 输出： 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
Windows圆周率计算工具
06-18
基于OpenCL的圆周率计算工具，可用于GPU性能测试，支持nVidia、AMD、Intel等厂商GPU、CPU。GPi可充分利用GPU强大的并行计算能力。 操作简单，一键开始，使用主流显卡可在3秒内算得Pi小数点后100万位。
Java小程序计算圆周率代码
08-29
Java小程序计算圆周率代码 Java小程序计算圆周率代码是指使用Java语言编写的小程序来计算圆周率π的近似值。圆周率π是数学中一个非常重要的常数，它是圆的周长与其直径的比值。计算圆周率π的方法有很多种，本文...
分别用MPI和串行程序计算圆周率
07-27
C++计算圆周率，分别用穿行计算和通过MPI实现的并行计算来进行。并行计算课程实验代码，分别用MPI和串行程序实现圆周率的计算并输出时间
计算圆周率
05-30
描述中反复提到了“计算圆周率c代码”，这可能是指一个C语言程序，该程序实现了某种特定的算法来估算π的值。C语言是一种强大的、低级别的编程语言，适用于执行高效的数值计算任务。 在提供的标签“计算圆周率”下...
python计算圆周率pi的方法
09-21
主要介绍了python计算圆周率pi的方法,涉及Python针对数学运算的相关技巧,需要的朋友可以参考下
蒙特卡罗方法 计算圆周率程序的python实现
09-04
蒙特卡罗方法 计算圆周率程序的python实现。 蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法。当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，...
圆周率π小数点后一亿位数.txt
02-14
圆周率一亿位
c语言编程快速收敛的圆周率计算,[原创]圆周率PI的计算(精确到几十万位)
weixin_39653717的博客
05-18 854
[原创]圆周率PI的计算(精确到几十万位)//环境：VC6.0,Console Application//原理：π=2+1/3*(2+2/5*(2+3/7*(2+...//特点：内嵌汇编提速并扩大了计算范围//限制：位数ws原则上没有限制但因为本//算法的时间正比于ws的平方所以将位数//控制在二、三十万以内较好。本人曾用//奔Ⅳ2.6GHz算20万位，耗时10分钟左右#include#inclu...
c语言计算10亿位圆周率,C语言：圆周率的计算
weixin_39807691的博客
05-20 1124
Description根据下面关系式，求圆周率pi的值，直到最后一项的值小于epsilon。epsilon的值由键盘输入。pi/2=1+1/3+2!/(3*5)+3!/(3*5*7)+...+n!/(3*537*...*(2n+1))Input输入一个浮点数表示epsilon的值Output输出圆周率Pi的值，保留6位小数。Sample Input0.000001Sample Output3.14...
C语言求圆周率π
liudaoqun的博客
03-27 3972
C语言求圆周率π #include <stdio.h> #include <stdlib.h> #include <math.h> int main() { float s = 1; float pi = 0; float i = 1.0; float n = 1.0; while (fabs(i) >= 1e-6) { pi += i; ...
关于计算圆周率PI的经典程序
weixin_30445169的博客
07-27 1676
短短几行代码，却也可圈可点。如把变量s放在PI表达式中，还有正负值的处理，都堪称经典。尤其是处处考虑执行效率的思想令人敬佩。 /* pi/4=1-1/3+1/5-1/7+1/9-…… */ #include <stdio.h> int main(){ int s=1; float pi=0.,n=1.,t=1.; ...
C语言 迭代求圆周率Pi
incredibleimpact的博客
12-09 1983
求圆周率Pi 效果图 代码如下 #include<stdio.h> #include<math.h> int main(int argc,char const*argv[]) { const int a=2; double pi=1,k=1+1e-10,ret,item; int j=1; do{ ret=1; double go=2; for(int i=1;i<j;i++)
【C语言程序设计】C语言求圆周率π（三种方法）
热门推荐
HUYA69的博客
03-05 2万+
题目一： 利用公式①计求π的近似值，要求累加到最后一项小于10^(-6)为止。 程序代码： #include <stdio.h> #include <stdlib.h> #include <math.h> int main(){ float s=1; float pi=0; float i=1.0; float n=1.0; while(fabs(i)>=1e-6){ ...
【c语言】计算圆周率的近似值
小飞的博客
03-31 2万+
用公式π/4=1-1/3+1/5-1/7+1/9-…,求圆周率π的近似值，直到上述级数最后一项的绝对值小于10^(−4) 为止 #include<stdio.h> int main(void) { float pi=0,t=1.0,n=1.0; while(1/n>=1e-4)//精度可以缩小，但计算时间会增加 //建议不小于1e-6 { pi=pi+t...
C语言使用define定义圆周率,C语言中的宏处理
weixin_27442001的博客
05-22 5911
在C语言中使用宏，我们经常这么做，但是为什么使用宏，他可以做什么，我们或许只是一知半解，下面简单介绍如何在C语言中使用宏预处理器。#define语句对于define语句，我们已经非常熟悉，一个宏定义语句，主要用途是定义符号常量，如：#define PI 3.14159定义了圆周率PI 为3.14159，在以后用到3.14159的地方直接用PI代替即可，也方便我们编写代码，更大的好处是，如果我不想将...
matlab计算圆周率程序
最新发布
06-11
在MATLAB中，你可以使用数值积分方法来近似计算圆周率π。其中一种常见的方法是使用著名的马鞍点公式（也称为蒙特卡洛...运行这个程序后，会计算指定数量的随机点来估计π的值。增加`n_points`的值，结果会更加准确。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值