[quote]感谢楼上的，我正在写大数乘法，那个AGM算法是几何算术平均值法吗~~[/quote]
是的

[quote][quote]感谢楼上的，我正在写大数乘法，那个AGM算法是几何算术平均值法吗~~[/quote]
是的[/quote]

你做过?[em14]

以前查过资料啊，AGM就是几何算术平均值法，一样的

你仔细读一下陈先生的主页 www.jason314.com,文中对AGM算法描述的很详细，运用这些知识完全可以编程的。
  主要的难点是如何实现
    1。大数乘法（2，3，4的基础）
    2。大数的倒数（3.基础）
    3。大数除法
    4。大数的平方根
 要实现一个效率不算太高的大数乘法并不难，最容易硬乘法，其次是分治法，最难的FFT/FNT算法。这些在《数与编程》有所提及，但不够详细。
你也可以http://numbers.computation.free.fr/Constants/constants.html 学到一些大数计算的方法。
  原代码方面，关于使用 分治法的 计算乘法的代码 可以 从 NTL(http://www.shoup.net)，GMP(http://swox.com/gmp/) 中看到。
  关于各种FFT算法和相关算法可参照 fxtbook.pdf(http://www.jjj.de/fxt/)
  关于使用 FFT 计算乘法的代码 可以 从 GMP 和　ooura fft (http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html,内有完整的计算大数乘法，除法，平方根，ＡＧＭ算法）
)中看到。
  关于使用FNT(快速数论变换）计算乘法的代码 可以 从 apfloat(http://www.apfloat.org/apfloat/) 中看到。
  我曾看到一个较简单的使用分治法计算PI的源代码，这是我第一次看到如何使用分治法计算乘法的代码，获益非浅，可惜我不能在我的电脑上找到．

　ＡＧＭ中计算ＰＩ的源文件为　const_pi.c, ooura fft 中计算ＰＩ的源文件为 pi_fft.c, apfloat 中计算ＰＩ的源文件为aptest.cpp

下面给出 一个一个老外 使用 martin 公式计算 PI 的源代码，思想和你的程序是一样的，只是书写风格不一样，你可以学习一下别人的代码，这样进步会快一些。
  

/*
** Pascal Sebah : September 1999
** 
** Subject:
**
**    A very easy program to compute Pi with many digits.
**    No optimisations, no tricks, just a basic program to learn how
**    to compute in multiprecision.  
**
** Formulae:
**
**    Pi/4 =    arctan(1/2)+arctan(1/3)                     (Hutton 1)
**    Pi/4 =  2*arctan(1/3)+arctan(1/7)                     (Hutton 2)
**    Pi/4 =  4*arctan(1/5)-arctan(1/239)                   (Machin)
**    Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
**
**      with arctan(x) =  x - x^3/3 + x^5/5 - ...
**
**    The Lehmer's measure is the sum of the inverse of the decimal
**    logarithm of the pk in the arctan(1/pk). The more the measure
**    is small, the more the formula is efficient.
**    For example, with Machin's formula:
**
**      E = 1/log10(5)+1/log10(239) = 1.852
** 
** Data:
**
**    A big real (or multiprecision real) is defined in base B as:
**      X = x(0) + x(1)/B^1 + ... + x(n-1)/B^(n-1)
**      where 0<=x(i)<B
**
** Results: (PentiumII, 450Mhz)
**    
**   Formula      :    Hutton 1  Hutton 2   Machin   Gauss
**   Lehmer's measure:   5.418     3.280      1.852    1.786
**
**  1000   decimals:     0.2s      0.1s       0.06s    0.06s
**  10000  decimals:    19.0s     11.4s       6.7s     6.4s
**  100000 decimals:  1891.0s   1144.0s     785.0s   622.0s
**
** With a little work it's possible to reduce those computation
** times by a factor 3 and more:
**  
**     => Work with double instead of long and the base B can
**        be choosen as 10^8
**     => During the iterations the numbers you add are smaller
**        and smaller, take this in account in the +, *, /
**     => In the division of y=x/d, you may precompute 1/d and
**        avoid multiplications in the loop (only with doubles)
**     => MaxDiv may be increased to more than 3000 with doubles
**     => ...
*/

#include <time.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
long B=10000; /* Working base */
long LB=4;    /* Log10(base)  */
long MaxDiv=450;  /* about sqrt(2^31/B) */
/*
** Set the big real x to the small integer Integer 
*/
void SetToInteger (long n, long *x, long Integer) {
  long i;
  for (i=1; i<n; i++) x[i] = 0;
  x[0] = Integer;
}
/*
** Is the big real x equal to zero ?
*/
long IsZero (long n, long *x) {
  long i;
  for (i=0; i<n; i++)  
    if (x[i])    return 0;
    return 1;
}
/*
** Addition of big reals : x += y
**  Like school addition with carry management
*/
void Add (long n, long *x, long *y) {
  long carry=0, i;
  for (i=n-1; i>=0; i--) {
    x[i] += y[i]+carry;
    if (x[i]<B) carry = 0;
    else {
      carry = 1;
      x[i] -= B;
    }
  }  
}
/*
** Substraction of big reals : x -= y
**  Like school substraction with carry management
**  x must be greater than y
*/
void Sub (long n, long *x, long *y) {
  long i;
  for (i=n-1; i>=0; i--) {
    x[i] -= y[i];
        if (x[i]<0) {
          if (i) {    
        x[i] += B;
        x[i-1]--;
      }
        }
  }  
}
/*
** Multiplication of the big real x by the integer q 
** x = x*q.
**  Like school multiplication with carry management
*/
void Mul (long n, long *x, long q) {
  long carry=0, xi, i;
  for (i=n-1; i>=0; i--) {
    xi  = x[i]*q;        
    xi += carry;        
    if (xi>=B) {
      carry = xi/B;
      xi -= (carry*B);
    }
    else 
      carry = 0;
    x[i] = xi;
    }  
}
/*
** Division of the big real x by the integer d 
** The result is y=x/d.
**  Like school division with carry management
**  d is limited to MaxDiv*MaxDiv.
*/
void Div (long n, long *x, long d, long *y) {
  long carry=0, xi, q, i;
  for (i=0; i<n; i++) {
    xi    = x[i]+carry*B;
    q     = xi/d;
    carry = xi-q*d;   
    y[i]  = q;        
  }  
}
/*
** Find the arc cotangent of the integer p = arctan (1/p)
**  Result in the big real x (size n)
**  buf1 and buf2 are two buffers of size n
*/
void arccot (long p, long n, long *x, long *buf1, long *buf2) {
  long p2=p*p, k=3, sign=0;
  long *uk=buf1, *vk=buf2;
  SetToInteger (n, x, 0);
  SetToInteger (n, uk, 1);    /* uk = 1/p */
  Div (n, uk, p, uk);
  Add (n, x, uk);              /* x  = uk */

  while (!IsZero(n, uk)) {
    if (p<MaxDiv)
      Div (n, uk, p2, uk);  /* One step for small p */
    else {
      Div (n, uk, p, uk);   /* Two steps for large p (see division) */
      Div (n, uk, p, uk);  
    }
    /* uk = u(k-1)/(p^2) */
    Div (n, uk, k, vk);       /* vk = uk/k  */
    if (sign) Add (n, x, vk); /* x = x+vk   */
    else Sub (n, x, vk);      /* x = x-vk   */
    k+=2;
    sign = 1-sign;
  }
}
/*
** Print the big real x
*/
void Print (long n, long *x) {
  long i; 
  printf ("%d.", x[0]);
  for (i=1; i<n; i++) {
    printf ("%.4d", x[i]);
    if (i%25==0) printf ("%8d\n", i*4);
  }
  printf ("\n");
}
/*
** Computation of the constant Pi with arctan relations
*/
void main () {  
  clock_t endclock, startclock; 
  long NbDigits=10000, NbArctan;
  long p[10], m[10];
  long size=1+NbDigits/LB, i;
  long *Pi      = (long *)malloc(size*sizeof(long));
  long *arctan  = (long *)malloc(size*sizeof(long));
  long *buffer1 = (long *)malloc(size*sizeof(long));
  long *buffer2 = (long *)malloc(size*sizeof(long)); 
  startclock = clock();    
  /*
  ** Formula used: 
  **   
  **   Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
  */
  NbArctan = 3;
  m[0] = 12; m[1] = 8;  m[2] = -5;
  p[0] = 18; p[1] = 57; p[2] = 239; 
  SetToInteger (size, Pi, 0);
  /*
  ** Computation of Pi/4 = Sum(i) [m[i]*arctan(1/p[i])] 
  */
  for (i=0; i<NbArctan; i++) {
    arccot (p[i], size, arctan, buffer1, buffer2);
    Mul (size, arctan, abs(m[i]));
    if (m[i]>0) Add (size, Pi, arctan);  
    else        Sub (size, Pi, arctan);  
  }
  Mul (size, Pi, 4);
  endclock = clock ();
  Print (size, Pi);  /* Print out of Pi */
  printf ("Computation time is : %9.2f seconds\n",  
         (float)(endclock-startclock)/(float)CLOCKS_PER_SEC ); 
  free (Pi);
  free (arctan);
     free (buffer1);
     free (buffer2);
}

感谢感谢!

我的FFT已经测试完成了,那个站我原来看过;都没余项啊,难到都不是搞数学的,或者是藏了一手

对于Machin's formula，或者AGM算法，其余项很容易推出来的，非得人家给余项吗？别人99.9%都做了,那0.1%只好麻烦楼主了。

呵呵，说反了吧，余项是误差分析的基础，其它的只要给几个公式就行了，写那么多有个屁用啊~

[quote]对于Machin's formula，或者AGM算法，其余项很容易推出来的，非得人家给余项吗？别人99.9%都做了,那0.1%只好麻烦楼主了。[/quote]
梁兄真是好同志啊, 把所有相关资料都详细地提供了,呵呵. 你说的那个分治法的源码我手头有,需要的人可以找我,呵呵.