用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);
    }
}

 

马青公式计算圆周率程序

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

Wallis公式及其应用

今天我要讲的主要内容是什么是Wallis公式,以及它的推导过程。然后再讲述Wallis公式的两个重要应用,即推导 Stirling公式和求解Euler-Poisson积分。   Contens    ...
  • ACdreamers
  • ACdreamers
  • 2014年11月25日 16:33
  • 10693

泰勒公式的推导过程(日期推算星期)

星期制度是一种有古老传统的制度。据说因为《圣经·创世纪》中规定上帝用了六  天时间创世纪,第七天休息,所以人们也就以七天为一个周期来安排自己的工作和生  活,而星期日是休息日。从实际的角度来讲,以...
  • blackgooes
  • blackgooes
  • 2016年05月10日 16:06
  • 1962

谈谈JAVA如何计算字符串公式

因为最近忙着面试,也没太多时间分享自己的博客,但因为面试的过程中,碰到了几道类似的算法题,都是输入一个字符串数学公式,输出结果。自己整理了下,希望能和大家交流交流。我是利用谭浩强那本C语言上的算法写的...
  • qqylc0011
  • qqylc0011
  • 2015年04月11日 15:59
  • 2718

复化梯形公式

复化梯形公式重在理解梯形公式的概念,计算f(X)=#include //梯形公式 #include //#define f(x) (sin(x)/x) double f(double...
  • Mikeoperfect
  • Mikeoperfect
  • 2017年04月29日 10:38
  • 833

Excel个人所得税简洁计算公式

iamlaosong文 我们知道,最新的个人所得税起征点是 3500元,税率从3%到45%,分有7个等级,为提高计算速度,与不同税率相对应有7个速算扣除数,如下表所示: 级数 含税级距 不含...
  • iamlaosong
  • iamlaosong
  • 2016年12月19日 17:19
  • 6779

各种求圆周率π的算法(蒙特卡洛法的Java实现)

算法就是有穷规则构成的用于解决某一类问题的运算序列或执行步骤。要解决一个问题可能会有不同的方法,针对求圆周率π的近似值这个问题你能想到多少种算法呢?本文讨论三种计算方法,并在Java中实现其中基于蒙特...
  • baimafujinji
  • baimafujinji
  • 2017年03月02日 12:52
  • 2580

用割圆术求圆周率π

以半径为1的圆开始,依次在圆内画正六边形、正十二边形、正6n边形...。      只有n趋近于无穷大,便可以得到足够接近圆的多边形,计算的圆周率π也就越精确。 完整代码如下:...
  • sanqima
  • sanqima
  • 2014年06月12日 09:28
  • 1728

一起来算圆周率

转载自:http://www.guokr.com/blog/444081/ 自古计算圆周率是数学界一项流行运动,各大数学家争相破记录以名垂青史。想象有人为圆周率15年如一日地算,算的不是圆周率而...
  • memray
  • memray
  • 2015年09月18日 05:01
  • 2505

Java处理数学公式得出结果

支持含有字母的复杂表达式、不含字母的纯数学表达式先上测试程序,后面有源码: 1.纯数字的算式测试 public static void main(String[] args) thro...
  • YangRunkangBla
  • YangRunkangBla
  • 2015年12月08日 18:08
  • 2412
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:用Machin公式计算圆周率的源程序
举报原因:
原因补充:

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