22

在本文中，给出了使用二元分裂的 Chudnovsky pi 公式的快速递归公式。在蟒蛇中： 

C = 640320
C3_OVER_24 = C**3 // 24
def bs(a, b):
    if b - a == 1:
        if a == 0:
            Pab = Qab = 1
        else:
            Pab = (6*a-5)*(2*a-1)*(6*a-1)
            Qab = a*a*a*C3_OVER_24
        Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
        if a & 1:
            Tab = -Tab
    else:
        m = (a + b) // 2
        Pam, Qam, Tam = bs(a, m)
        Pmb, Qmb, Tmb = bs(m, b)

        Pab = Pam * Pmb
        Qab = Qam * Qmb
        Tab = Qmb * Tam + Pam * Tmb
    return Pab, Qab, Tab

N = int(digits/DIGITS_PER_TERM + 1)
# Calclate P(0,N) and Q(0,N)
P, Q, T = bs(0, N)
one = 10**digits
sqrtC = sqrt(10005*one, one)
return (Q*426880*sqrtC) // T

这种方法已经非常快了，但是提到了 GMP 库网站gmp-chudnovsky.c上的一个实现，在二进制拆分中也考虑了分子和分母。由于代码已经过优化并且我很难理解，所以这背后的总体思路是什么？我不知道分数是否被简化，数字是否保持因式分解而不是完全相乘，或两者兼而有之。

这是二进制拆分的代码示例： 

/* binary splitting */
void
bs(unsigned long a, unsigned long b, unsigned gflag, long int level)
{
  unsigned long i, mid;
  int ccc;

  if (b-a==1) {
    /*
      g(b-1,b) = (6b-5)(2b-1)(6b-1)
      p(b-1,b) = b^3 * C^3 / 24
      q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
    */
    mpz_set_ui(p1, b);
    mpz_mul_ui(p1, p1, b);
    mpz_mul_ui(p1, p1, b);
    mpz_mul_ui(p1, p1, (C/24)*(C/24));
    mpz_mul_ui(p1, p1, C*24);

    mpz_set_ui(g1, 2*b-1);
    mpz_mul_ui(g1, g1, 6*b-1);
    mpz_mul_ui(g1, g1, 6*b-5);

    mpz_set_ui(q1, b);
    mpz_mul_ui(q1, q1, B);
    mpz_add_ui(q1, q1, A);
    mpz_mul   (q1, q1, g1);
    if (b%2)
      mpz_neg(q1, q1);

    i=b;
    while ((i&1)==0) i>>=1;
    fac_set_bp(fp1, i, 3);  /*  b^3 */
    fac_mul_bp(fp1, 3*5*23*29, 3);
    fp1[0].pow[0]--;

    fac_set_bp(fg1, 2*b-1, 1);  /* 2b-1 */
    fac_mul_bp(fg1, 6*b-1, 1);  /* 6b-1 */
    fac_mul_bp(fg1, 6*b-5, 1);  /* 6b-5 */

    if (b>(int)(progress)) {
      printf("."); fflush(stdout);
      progress += percent*2;
    }

  } else {
    /*
      p(a,b) = p(a,m) * p(m,b)
      g(a,b) = g(a,m) * g(m,b)
      q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
    */
    mid = a+(b-a)*0.5224;     /* tuning parameter */
    bs(a, mid, 1, level+1);

    top++;
    bs(mid, b, gflag, level+1);
    top--;

    if (level == 0)
      puts ("");

    ccc = level == 0;

    if (ccc) CHECK_MEMUSAGE;

    if (level>=4) {           /* tuning parameter */
#if 0
      long t = cputime();
#endif
      fac_remove_gcd(p2, fp2, g1, fg1);
#if 0
      gcd_time += cputime()-t;
#endif
    }

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(p1, p1, p2);

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(q1, q1, p2);

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(q2, q2, g1);

    if (ccc) CHECK_MEMUSAGE;
    mpz_add(q1, q1, q2);

    if (ccc) CHECK_MEMUSAGE;
    fac_mul(fp1, fp2);

    if (gflag) {
      mpz_mul(g1, g1, g2);
      fac_mul(fg1, fg2);
    }
  }

  if (out&2) {
    printf("p(%ld,%ld)=",a,b); fac_show(fp1);
    if (gflag)
      printf("g(%ld,%ld)=",a,b); fac_show(fg1);
  }
}
4

2 回答 2

6

以下评论是关键：

/*
    g(b-1,b) = (6b-5)(2b-1)(6b-1)
    p(b-1,b) = b^3 * C^3 / 24
    q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
*/

/*
    p(a,b) = p(a,m) * p(m,b)
    g(a,b) = g(a,m) * g(m,b)
    q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
*/

/*
          p*(C/D)*sqrt(C)
    pi = -----------------
             (q+A*p)
*/

观察：

  1. 如果我们在最终除法中同时缩放p和缩放q因子k，则不会影响结果。
  2. 在计算对应于第二组注释的合并操作时，如果我们要缩放 , 中的每一个，p(a,m)那么到那时该因子将简单地传递到, , ; 类似地，如果我们要对 、 、 中的每一个进行缩放，那么该因素就会发挥作用。g(a,m)q(a,m)kp(a,b)g(a,b)q(a,b)p(m,b)g(m,b)q(m,b)

  3. 线

    fac_remove_gcd(p2, fp2, g1, fg1);

    有效地

    k = gcd(p2, g1);
p2 /= k; // p(m,b) /= k
g1 /= k; // g(a,m) /= k

    这具有缩小p(a,b),g(a,b)和q(a,b)的净效应gcd。通过前两次观察，这种缩小一直干净地进行到最终结果。

后记

我尝试了三种在 Python 中实现因式分解的方法。

  • 使用不可变列表跟踪分解会大大减慢速度，因为维护列表的工作太多了。
  • 使用 gmpygcd并不能提高速度
  • 预先计算一个素数列表6 * N（即可能除以g）并测试每个素数会使事情减慢 2 到 3 倍。

我得出的结论是，使用这种方法来加速需要使用可变状态来跟踪分解，所以它对可维护性来说是一个很大的打击。

于 2015-11-28T13:41:35.140 回答
3

我没有查看完整的代码，但我快速浏览了一下，以便更好地理解您在问题中提供的摘录。

为了回答您的问题中的某些观点，请先查看这段代码：

typedef struct {
  unsigned long max_facs;
  unsigned long num_facs;
  unsigned long *fac;
  unsigned long *pow;
} fac_t[1];

它将新类型定义为结构（如果您根本不懂 C，假设它就像一个基本的 Pyhton 对象嵌入变量但没有方法）。这种类型允许将整数处理为两个整数值和两个数组（例如：两个列表）：

  • 最大的因素
  • 因素数
  • 因素清单
  • 这些因素的（相应）权力清单

同时，代码保持与 libgmp 类型中的大整数相同的数字（这是函数参数的含义mpz_t p和含义mpz_t g）。

现在，您所展示的功能如何。它被称为fac_remove_gcd；首字母fac与前面描述的类型的名称有关；以下两个词很容易理解：将 type 的两个整数除以fac它们的 gcd。

C 代码遍历两个列表中的两个因子列表；else if由于因子是有序的（ andelse语句周围的代码部分），因此很容易同步两个列表；每当检测到两个公因子（初始if语句）时，除法就是一个简单的减法问题：在这个因子的两个幂列表中减去最小的幂（例如 a=2*5^3*17 和 b=3* 5^5*19，值 3 将在a和b的幂列表中减去对应于因子 5 的位置，导致 a=2*5^0*17 和 b=3*5^2*19） .

fac在此操作期间，会创建并调用一个数字（相同类型） fmul；这显然是两个数字的gcd。

在这一步之后，调用的 gcdfmul和fac类型被转换为 GMP 大整数，并调用函数（也在程序的代码中）bs_mul。这允许将 GCD 计算为一个大整数，以便以两种形式同步除数的新值：大整数和特殊fac类型。一旦 GCD 被计算为一个大整数，很容易将两个初始数字除以 GCD。

因此，这些函数作用于每个初始数字的两个版本。

希望它可以提供帮助。

于 2015-11-15T18:26:24.140 回答