q_pi.txt

用tc来计算pi的精确计算

/*********************Nrutil.h*****************************/
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1)>(dmaxarg2)?\
		   (dmaxarg1):(dmaxarg2))

static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1)<(dminarg2)?\
		   (dminarg1):(dminarg2))

/*************************** 内存分配 **********************************/
double *dvector(long nl,long nh);
unsigned char *cvector(long nl,long nh);
void free_dvector(double *v,long nl,long nh);
void free_cvector(unsigned char *v,long nl,long nh);

/*********************Nrutil.c*****************************/
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1;
#define FREE_ARG char*

extern nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
  fprintf(stderr,"Numerical Recipes run-time error...");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"....now exitinf to system....");
  exit(1);
}
double *dvector(long nl,long nh)
/* allocate a double vector with subcript rang v[nl..nh]*/
{
  double *v;
  v=(double *) malloc((size_t) ( (nh-nl+1+1)*sizeof(double) ) );/* NR_END=1;*/
  if(!v) nrerror("alloction failure in dvector()");
  return v-nl+NR_END;
}

unsigned char *cvector(long nl,long nh)
/* allocate a unsigned char vector with subcript rang v[nl..nh]*/
{
  unsigned char *v;
  v=(unsigned char*) malloc((size_t) ( (nh-nl+1+1)*sizeof(unsigned char) ) );/* NR_END=1;*/
  if(!v) nrerror("alloction failure in cvector()");
  return v-nl+NR_END;
}

void free_dvector(double *v,long nl,long nh)
/* free a double vector allocated dvector() */
{
  free((FREE_ARG)(v+nl-1));   /* NR_END=1*/
}

void free_cvector(unsigned char *v,long nl,long nh)
/* free an unsigned char vector allocated cvector() */
{
  free((FREE_ARG)(v+nl-1));   /* NR_END=1*/
}
/*********************s_pi_chu.c*****************************/
#include "nrutil.h"
#define MF 4
#define BI (1.0/256)
#define MACC 3

/*字符串v[1..m]是以256为基数的数，其在v[1]后有基数的小数点;
  u[1..n]是设置为它的倒数的最高有效数字,在u[1]前有基数的小数点.*/
void Sinv(unsigned char u[],unsigned char v[],int n,int m)
{
  void Scpy(unsigned char u[],unsigned char v[],int n);
  void Smul(unsigned char w[],unsigned char u[],unsigned char v[],
	    int n,int m);
  void Sneg(unsigned char u[],int n);
  unsigned char *rr,*s;
  int i,j,maxmn,mm;
  float fu,fv;

  maxmn=IMAX(m,n);
  rr=cvector(1,1+(maxmn<<1));
  s=cvector(1,maxmn);
  mm=IMIN(MF,m);
  fv=(float) v[mm]; /* 采用原来的浮点运算以得到一个近似值*/
  for(j=mm-1;j>=1;j--){
      fv*=BI;
      fv+=v[j];
  }
  fu=1.0/fv;
  for(j=1;j<=n;j++){
      i=(int) fu;
      u[j]=(unsigned char) i;
      fu=256.0*(fu-i);
  }
  for(;;){                   /* 迭代牛顿法则到收敛*/
      Smul(rr,u,v,n,m);      /* 在s中构造2-u*v    */
      Scpy(s,&rr[1],n);
      Sneg(s,n);
      s[1]-=254;             /* 将s*u 的乘积赋于u */
      Smul(rr,s,u,n,n);
      Scpy(u,&rr[1],n);
      for(j=2;j<n;j++){      /* 若s的小数部分不是零,它就没有收敛到1 */
	  if(s[j]) break;
      }
      if(j==n){
	 free_cvector(s,1,maxmn);
	 free_cvector(rr,1,1+(maxmn<<1));
	 return;
      }
  }
}

Sdiv(unsigned char q[],unsigned char r[],unsigned char u[],
     unsigned char v[],int n,int m)
/* 以256为基数的无符号整数u[1..n]被v[1..m]除(要求m<n);
   产生2商为q[1..n-m+1]和余数r[1..m]  */
{
  int is;
  unsigned char *rr,*s;

  rr=cvector(1,(n+MACC)<<1);
   s=cvector(1,n+MACC);
  Sinv(s,v,n-m+MACC,m);            /* s=1/v */
  Smul(rr,s,u,n-m+MACC,n);         /* q=s*u */
  Scpy(q,&rr[1],n-m+1);
  Smul(rr,q,v,n-m+1,m);            /* 相乘相减得余数 */
  Ssub(&is,&rr[1],u,&rr[1],n);
  if(is) nrerror("MACC too small in Sdiv");
  Scpy(r,&rr[n-m+1],m);
  free_cvector(s,1,n+MACC);
  free_cvector(rr,1,(n+MACC)<<1);
}
/*********************s_pi_ren.c*****************************/
/* 将比255大的数划分成低位字节和高位进位,可对译成以256为基数的字符串
   进行多精度运算*/
#define LOBYTE(x) ( (unsigned char)((x)&0xff) )
#define HIBYTE(x) ( (unsigned char)((x)>>8&0xff) )


/* 加法*/
/* 无符号整数u[1..n]和v[1..n]相加,产生无符号整数w[1..n+1],以256为基数*/
void Sadd(unsigned char w[],unsigned char u[],unsigned char v[],int n)
{
 int j;
 unsigned short ireg=0;
 for(j=n;j>=1;j--){
     ireg=u[j]+v[j]+HIBYTE(ireg);
     w[j+1]=LOBYTE(ireg);
 }
 w[1]=HIBYTE(ireg);
}

/* 除法*/
/* 无符号整数u[1..n]中减去v[1..n],产生无符号整数w[1..n],以256为基数*/
/* 若相减结果为负,则is返回-1;反之,则返回0*/
void Ssub(int *is,unsigned char w[],unsigned char u[],unsigned char v[],int n)
{
 int j;
 unsigned short ireg=256;
 for(j=n;j>=1;j--){
     ireg=255+u[j]-v[j]+HIBYTE(ireg);
     w[j]=LOBYTE(ireg);
 }
 *is=HIBYTE(ireg)-1;
}
/* 捷加法(加一单个字节到一串数中) */
/* 将整数iv(0<=iv<=255)与以256为基数的无符号整数u[1..n]相加,产生w[1..n];*/
void Sjadd(unsigned char w[],unsigned char u[],int n,int iv)
{
 int j;
 unsigned short ireg;

 ireg=256*iv;
 for(j=n;j>=1;j--){
     ireg=u[j]+HIBYTE(ireg);
     w[j+1]=LOBYTE(ireg);
 }
 w[1]=HIBYTE(ireg);
}

/* 捷乘法*/
/* 以256为基数的无符号整数u[1..n]乘以整数iv(0<=iv<=255),产生w[1..n+1]*/
void Sjmul(unsigned char w[],unsigned char u[],int n,int iv)
{
 int j;
 unsigned short ireg=0;
 for(j=n;j>=1;j--){
     ireg=u[j]*iv+HIBYTE(ireg);
     w[j+1]=LOBYTE(ireg);
 }
 w[1]=HIBYTE(ireg);
}

/* 捷除法____以256为基数的无符号整数u[1..n]被整数iv除(而0<=iv<=255)*/
/* 产生商为w[1..n+1]和余数ir(0<=ir<=255)*/
void Sjdiv(unsigned char w[],unsigned char u[],int n,int iv,int *ir)
{
 int i,j;

 *ir=0;
 for(j=1;j<=n;j++){
     i=(*ir)*256+u[j];
     w[j]=(unsigned char)(i/iv);
     *ir=i%iv;
 }
}

/* 补码非操作____就是非操作*/
/* 以256为基数的无符号整数u[1..n]执行一次补码非操作*/
void Sneg(unsigned char u[],int n)
{
 int j;
 unsigned char ireg=256;

 for(j=n;j>=1;j--){
     ireg=255-u[j]+HIBYTE(ireg);
     u[j]=LOBYTE(ireg);
 }
}

/* 复制*/
/* 将v[1..n]复制到u[1..n]中*/
void Scpy(unsigned char u[],unsigned char v[],int n)
{
 int j;
 for(j=1;j<=n;j++)  u[j]=v[j];
}

/* 左移一串数*/
/* 将u[2..n+1]左移成u[1..n]*/
Szmov(unsigned char u[],int n)
{
 int j;
 for(j=1;j>=n;j++) u[j]=u[j+1];
}
/*********************s_pi_xia.c*****************************/
#include "nrutil.h"  /* 快速乘法*/
#define RX 256.0

/* 运用快速傅里叶变换计算*/
/* 以256为基数的无符号整数u[1..n]与v[1..m]相乘,产生积为w[1..n+m]*/
void Smul(unsigned char w[],unsigned char u[],unsigned char v[],int n,int m)
{
 int j,mn,nn=1;
 double cy,t,*a,*b;     /*宣告fft程序,为双精度版本*/
 void drealft(double data[],unsigned long n,int isign);
 mn=IMAX(m,n);
 while(nn<mn)  nn<<=1; /* 为变换找出最小可用的2的幂次 */
 nn<<=1;
 a=dvector(1,nn);
 b=dvector(1,nn);
 for(j=1;j<=n;j++)       a[j]=(double) u[j];  /* 将u移入双精度浮点型数组*/
 for(j=n+1;j<=nn;j++)    a[j]=0.0;
 for(j=1;j<=m;j++)       b[j]=(double) v[j];  /* 将v移入双精度浮点型数组*/
 for(j=m+1;j<=nn;j++)    b[j]=0.0;
 /* 执行卷积:首先求出二个傅里叶变换*/
 drealft(a,nn,1);
 drealft(b,nn,1);
 /*复数相乘的结果(实部和虚部)
 b[1]*=a[1];  b[2]*=a[2];
 for(j=3;j<nn;j++){
     b[j]=  (t=b[j])*a[j]-b[j+1]*a[j+1];
     b[j+1]=t*a[j+1]+b[j+1]*a[j];
 }
 drealft(b,nn,-1);          /* 进行傅里叶反变换*/

 /* 执行最后完成所有进位的过程*/
 cy=0.0;
 for(j=nn;j>=1;j--){
     t=b[j]/(nn>>1)+cy+0.5; /* 这0.5是为了舍入误差而允许的*/
     b[j]=t;
     while(b[j]>=RX) b[j]-=RX;
     cy=(int)(t/RX);
 }
 if(cy>=RX) nrerror("cannot happen in fftmul");
 w[1]=(unsigned char) cy;
 for(j=2;j<=n+m;j++)
   w[j]=(unsigned char) b[j-1];
 free_dvector(b,1,nn);
 free_dvector(a,1,nn);
}
/*********************s_pi_fft.c*****************************/
#include <math.h>
#define SWAP(a,b) tempr=(a);(b)=(a);(b)=tempr;
void dfour1(double data[],unsigned long nn,int isign)
/* 假如isign=1,则制换后的data[1..2*nn]是输入的离散傅里变换;
  假如isign=-1,则制换后的data[1..2*nn]是输入的nn个离散傅里逆变换;
*/
{
  unsigned long n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  double tempr,tempi;
  n=nn<<1;
  j=1;
  for(i=1;i<n;i+=2){
    if(j>i){
      SWAP(data[j],data[i]);
      SWAP(data[j+1],data[i+1]);
    }
    m=n>>1;
    while(m>=2&&j>m){
      j-=m;
      m>>=1;
    }
    j+=m;
  }
  mmax=2;
  while(n>mmax){
    istep=mmax<<1;
    theta=isign*(6.28318530717959/mmax);
    wtemp=sin(0.5*theta);
    wpr=-2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for(m=1;m<mmax;m+=2){
      for(i=m;i<=n;i+=istep){
	j=i+mmax;
	tempr=wr*data[j]-wi*data[j+1];
	tempi=wr*data[j+1]+wi*data[j];
	data[j]=data[i]-tempr;
	data[j+1]=data[i+1]-tempi;
	data[i]+=tempr;
	data[i+1]+=tempi;
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
  }
}

void drealft(double data[],unsigned long n,int isign)
/*计算一组n个实值数据点傅里叶变换;用傅里叶变换的正半频率替换这些数据;
  复变换的第一个和最后一个分量的实数值分别返回单元data[1]和data[2]; */
{
  void dfour1(double data[],unsigned long nn,int isign);
  unsigned long i,i1,i2,i3,i4,np3;
  double c1=0.5,c2,h1r,h1i,h2r,h2i;
  double wr,wi,wpr,wpi,wtemp,theta;

  theta=3.141592653589793/(double)(n>>1);
  if(isign==1){
    c2=0.5;
    dfour1(data,n>>1,1);
  }else{
    c2=0.5;
    theta=-theta;
  }
  wtemp=sin(0.5*theta);
  wpr=-2.0*wtemp*wtemp;
  wpi=sin(theta);
  wr=1.0+wpr;
  wi=wpi;
  np3=n+3;
  for(i=2;i<=(n>>2);i++){
    i4=1+( i3=np3-( i2=1+(i1=i+i-1) ) );
    h1r=c1*(data[i1]+data[i3]);
    h1i=c1*(data[i2]+data[i4]);
    h2r=-c2*(data[i2]+data[i4]);
    h2i= c2*(data[i1]+data[i3]);
    data[i1]=h1r+wr*h2r-wi*h2i;
    data[i2]=h1i+wr*h2i+wi*h2r;
    data[i3]=h1r-wr*h2r+wi*h2i;
    data[i2]=-h1i+wr*h2i+wi*h2r;
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
  }
  if(isign==1){
    data[1]=(h1r=data[1])+data[2];
    data[2]=h1r-data[2];
  }else{
    data[1]=c1*((h1r=data[1])+data[2]);
    data[2]=c1*(h1r-data[2]);
    dfour1(data,n>>1,-1);
  }
}
/*********************s_pi_sqr.c*****************************/
#include <math.h>
#include "nrutil.h"
/* #define MF 3*/
#define BI (1.0/256)

#define IAZ 28

void Ssqrt(unsigned char w[],unsigned char u[],unsigned char v[],int n,int m)
/* 字符串v[1..m]是以256为基数,其在v[1]后有小数点;
   w[1..n]设为v的平方根(基数小数点在w[1]后)
*/
{
  int i,ir,j,mm;
  float fu,fv;
  unsigned char *r,*s;

  r=cvector(1,n<<1);
  s=cvector(1,n<<1);
  mm=IMIN(m,MF);
  fv=(float)v[mm];
  for(j=mm-1;j>=1;j--){
    fv*=BI;
    fv+=v[j];
  }
  fu=1.0/sqrt(fv);
  for(j=1;j<=n;j++){
    i=(int) fu;
    u[j]=(unsigned char) i;
    fu=256.0*(fu-1);
  }
  for(;;){                   /* 迭代牛顿法则收敛,构造s=(3-v*u*u)/2 */
    Smul(r,u,u,n,n);
    Szmov(r,n);
    Smul(s,r,v,n,m);
    Szmov(s,n);
    Sneg(s,n);
    s[1]=-253;
    Sjdiv(s,s,n,2,&ir);
    for(j=2;j<n;j++){
      if(s[j]){            /* 若s的小数部分不为零,它不能收敛到1 */
	Smul(r,s,u,n,n);   /* 用su代替u */
	Scpy(u,&r[1],n);
	break;
      }
    }
    if(j<n)   continue;
    Smul(r,u,v,n,m);
    Scpy(w,&r[1],n);
    free_cvector(s,1,n<<1);
    free_cvector(r,1,n<<1);
    return;
  }
}


void S2dfr(unsigned char a[],unsigned char s[],int n,int *m)
/* 将以256为基数的小数a[1..n],(小数点在a[1]前);转换成十进制,用s[1..m]表示;
   其中m是返回值;
*/
{
  int j;
  *m=(int)(2.408*n);
  for(j=1;j<=(*m);j++){
    Sjmul(a,a,n,10);
    s[j]=a[1]+IAZ;
    Szmov(a,n);
  }
}
/*********************s_pi_zpi.c*****************************/
#include <stdio.h>
#include "nrutil.h"
#define IAOFF 48
extern void Ssqrt(unsigned char w[],unsigned char u[],unsigned char v[],int n,int m);
extern void S2dfr(unsigned char a[],unsigned char s[],int n,int *m);


void Szpi(int n)
/* 计算和打印pi的前n位字节的多精度程序 */
{
  int ir,j,m;
  unsigned char mm,*x,*y,*sx,*sxi,*t,*s,*pi;

  x=cvector(1,n+1);
  y=cvector(1,n<<1);
  sx=cvector(1,n);
  sxi=cvector(1,n);
  t=cvector(1,n<<1);
  s=cvector(1,3*n);
  pi=cvector(1,n+1);
  t[1]=2;                      /*设 T=2   */
  for(j=2;j<=n;j++)  t[j]=0;
  Ssqrt(x,x,t,n,n);            /* x[1]=2的平方根 */
  Sadd(pi,t,x,n);            /* pi[1]=2+2的平方根 */;
  Szmov(pi,n);
  Ssqrt(sx,sxi,x,n,n);         /* 2的平方根的平方根 */
  Scpy(y,sx,n);
  for(;;){
    Sadd(x,sx,sxi,n);
    Sjdiv(x,&x[1],n,2,&ir);
    Ssqrt(sx,sxi,x,n,n);
    Smul(t,y,sx,n,n);
    Sadd(&t[1],&t[1],sxi,n);
    x[1]++;
    y[1]++;
    Sinv(s,y,n,n);
    Smul(y,&t[2],s,n,n);
    mm=t[2]-1;
    j=t[n+1]-mm;
    if(j>1||j<-1){
      for(j=3;j<=n;j++){
	if(t[j]!=mm){
	  Smul(s,pi,&t[1],n,n);
	  Scpy(pi,&s[1],n);
	  break;
	}
      }
      if(j<=n) continue;
    }
    printf("pi=\n");
    s[1]=pi[1]+IAOFF;
    s[2]='.';
    m=mm;
    S2dfr(&pi[1],&s[2],n-1,&m);
    s[m+3]=0;
    printf("%64s\n",&s[1]);
    free_cvector(pi,1,n+1);
    free_cvector(s,1,3*n);
    free_cvector(t,1,n<<1);
    free_cvector(sxi,1,n);
    free_cvector(sx,1,n);
    free_cvector(y,1,n<<1);
    free_cvector(x,1,n+1);
    return;
  }
}
/*********************s_pimain.c*****************************/
#include "s_pi_ren.c"
#include "s_pi_xia.c"
#include "s_pi_chu.c"
#include "s_pi_fft.c"
#include "s_pi_zpi.c"
#include "s_pi_sqr.c"
#include "nrutil.c"

main()
{
 int i=10;
 Szpi(i);
}