用Machin公式计算圆周率的源程序

原创 2004年09月08日 18:46:00

用Machin公式计算圆周率的源程序

 
/*  Program to compute PI, by Jason Chen, May 1999
**
**  Open VC++ IDE, new a win32 console application project,
**  add this file to the project and build it.
**
**  Homepage : http://www.jason314.com
**  Email    : jason@szonline.net
*/
#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 main(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   = calloc(BinLen, 4);
    arctg239 = calloc(BinLen, 4);
    tmp      = 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("Press 'p' to pause & exit/n");
    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("Writing result to file: 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("Done !!/n");
    ts = gmtime(&TotalTime);
    ts->tm_mon --;
    ts->tm_mday = ts->tm_mday - 1 + ts->tm_mon * 31;
    printf("Time : ");
    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);
    }
}

 

解读传说中计算π的超牛的C程序

  • xiao69
  • xiao69
  • 2009-11-23 21:46:00
  • 5318

PI值计算

HDU 2179 pi值计算  先发上大数版本的程序(java水的,不想写高精度了。。) import java.math.BigDecimal; import java.math.B...
  • swm8023
  • swm8023
  • 2011-10-08 01:54:40
  • 4676

圆周率PI的源程序(暂计算到100000位并且可以精确到无限位!)

对超范围的问题很感兴趣!常规的算法是无法处理这类问题的,在一位学长的帮助下终于解决的圆周率的问题!下面的源代码供大家学习交流!!!#include #include #include int main...
  • aa569658851
  • aa569658851
  • 2010-06-19 09:19:00
  • 217

计算圆周率的C程序

这是很早以前用C语言写的一个计算圆周率的程序,  算法是用泰勒公式计算反正切值。在命令行不跟参数执行该程序则使用Gauss公式计算前1000位圆周率的值,如果带一个命令行参数,则该值为要计算的位数。如...
  • u013948187
  • u013948187
  • 2015-08-04 07:10:46
  • 2982

马青公式计算圆周率程序

马青公式(梅钦公式、Machin)计算圆周率是这样的: 用这公式每增加计算一项,就可以增加约1.39位圆周率准确值,计算百万以下级别精度的圆周率,速度还算满意。 马青公式级数展开后,可以看...
  • MHL_1208980380
  • MHL_1208980380
  • 2017-01-09 05:30:52
  • 2038

<em>计算</em>圆周率 <em>源代码</em>

<em>计算</em>圆周率 <em>源代码</em> 2014-04-12 上传大小:41KB <em>计算</em>圆周率 能够快速<em>计算</em>圆周率,...C语言<em>计算</em>圆周率<em>Pi</em>的程序例子-CSDN下载 时间:2018-4-9 [免费]用canvas<em>计算</em>圆周率...
  • 2018年04月17日 00:00

Pi的计算程序

一直对PI的计算程序感到好奇,到网上找了找相关文章。照抄写了一个程序。用了第三个式子。ARRAYSIZE代表位数。:)#include #include #include #define ARRAYS...
  • skyofdreams
  • skyofdreams
  • 2006-08-15 12:18:00
  • 1551

利用公式计算π=4(1-1/3+1/5-1/7+1/9-...)的近似值,直到括号中最后一项的绝对值小于0.000001为止

#include #include using namespace std; int main( ) {int a=1;  double b=1,c=1,pi=0;  while(fabs...
  • w92a01n19g
  • w92a01n19g
  • 2012-09-15 18:24:40
  • 1500

并行<em>计算</em>中求<em>pi</em>的源程序

并行<em>计算</em>中求<em>pi</em>的源程序,比较简单,容易实现和理解。是比较经典的MPI求<em>pi</em>程序。综合评分:5 收藏评论(3)举报 所需: 3积分/C币 下载个数: 79 开通VIP 立即下载...
  • 2018年04月14日 00:00

转载:计算几何算法(含源代码)

㈠ 点的基本运算 1. 平面上两点之间距离 1 2. 判断两点是否重合 1 3. 矢量叉乘 1 4. 矢量点乘 2 5. 判断点是否在线段上 2 6. 求一点饶某点旋转后的坐标 2 7. 求矢量夹角 ...
  • lxmhuendan
  • lxmhuendan
  • 2008-01-11 14:58:00
  • 993
收藏助手
不良信息举报
您举报文章:用Machin公式计算圆周率的源程序
举报原因:
原因补充:

(最多只允许输入30个字)