Chudnovsky 二元分裂和因式分解
Chudnovsky binary splitting and factoring
在this article中,给出了使用二进制拆分的 Chudnovsky pi 公式的快速递归公式。在 python 中:
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);
}
}
我没有看完整的代码,但我快速浏览了一下,以便更好地理解您在问题中提供的摘录。
为了回答您问题中的一些要点,请先看一下这段代码:
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
与前面描述的类型的名称有关;下面两个词很容易理解:将两个fac
类型的整数除以它们的gcd.
C 代码迭代两个列表中的两个因子列表;由于因子是有序的(else if
和 else
语句周围的代码部分),因此很容易同步两个列表;每当检测到两个公因数(初始 if
语句)时,除法就是一个简单的减法问题:在该因数的两个幂列表中减去最小的幂(例如,a=2*5^3*17和b=3*5^5*19,在a和b的幂列表中相应的位置减去值3因子 5 导致 a=2*5^0*17 和 b=3*5^2*19).
在此操作期间,创建并调用了一个数字(相同 fac
类型)并调用 fmul
;这显然是两个数字的gcd。
在此步骤之后,名为 fmul
且属于 fac
类型的 gcd 被转换为 GMP 大整数,函数(也在程序代码中)称为 bs_mul
.这允许将 GCD 计算为一个大整数,以便同步两种形式的除数的新值:大整数和特殊 fac
类型。一旦 GCD 被计算为一个大整数,就很容易将两个初始数字除以 GCD。
因此函数作用于每个初始数字的两个版本。
希望对您有所帮助。
以下评论是关键:
/*
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)
*/
观察:
- 如果我们在最后的除法中将
p
和 q
都缩放 k
,这不会影响结果。
- 在计算对应于第二组评论的合并操作时,如果我们将
p(a,m)
、g(a,m)
和 q(a,m)
中的每一个缩放为 k
那么该因素将简单地传递到 p(a,b)
、g(a,b)
、q(a,b)
;类似地,如果我们要对 p(m,b)
、g(m,b)
和 q(m,b)
中的每一个进行缩放,那么该因子将继续存在。
行
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 中因式分解的方法。
- 使用不可变列表跟踪因式分解会极大地减慢速度,因为维护列表的工作太忙了。
- 使用 gmpy 的
gcd
不会加速
- 预先计算一个最多
6 * N
的素数列表(即可能整除 g
)并测试每个素数会使速度减慢 2 到 3 倍。
我的结论是,通过这种方法获得加速需要使用可变状态来跟踪因式分解,因此这对可维护性来说是一个很大的打击。
在this article中,给出了使用二进制拆分的 Chudnovsky pi 公式的快速递归公式。在 python 中:
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);
}
}
我没有看完整的代码,但我快速浏览了一下,以便更好地理解您在问题中提供的摘录。
为了回答您问题中的一些要点,请先看一下这段代码:
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
与前面描述的类型的名称有关;下面两个词很容易理解:将两个fac
类型的整数除以它们的gcd.
C 代码迭代两个列表中的两个因子列表;由于因子是有序的(else if
和 else
语句周围的代码部分),因此很容易同步两个列表;每当检测到两个公因数(初始 if
语句)时,除法就是一个简单的减法问题:在该因数的两个幂列表中减去最小的幂(例如,a=2*5^3*17和b=3*5^5*19,在a和b的幂列表中相应的位置减去值3因子 5 导致 a=2*5^0*17 和 b=3*5^2*19).
在此操作期间,创建并调用了一个数字(相同 fac
类型)并调用 fmul
;这显然是两个数字的gcd。
在此步骤之后,名为 fmul
且属于 fac
类型的 gcd 被转换为 GMP 大整数,函数(也在程序代码中)称为 bs_mul
.这允许将 GCD 计算为一个大整数,以便同步两种形式的除数的新值:大整数和特殊 fac
类型。一旦 GCD 被计算为一个大整数,就很容易将两个初始数字除以 GCD。
因此函数作用于每个初始数字的两个版本。
希望对您有所帮助。
以下评论是关键:
/*
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)
*/
观察:
- 如果我们在最后的除法中将
p
和q
都缩放k
,这不会影响结果。 - 在计算对应于第二组评论的合并操作时,如果我们将
p(a,m)
、g(a,m)
和q(a,m)
中的每一个缩放为k
那么该因素将简单地传递到p(a,b)
、g(a,b)
、q(a,b)
;类似地,如果我们要对p(m,b)
、g(m,b)
和q(m,b)
中的每一个进行缩放,那么该因子将继续存在。 行
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 中因式分解的方法。
- 使用不可变列表跟踪因式分解会极大地减慢速度,因为维护列表的工作太忙了。
- 使用 gmpy 的
gcd
不会加速 - 预先计算一个最多
6 * N
的素数列表(即可能整除g
)并测试每个素数会使速度减慢 2 到 3 倍。
我的结论是,通过这种方法获得加速需要使用可变状态来跟踪因式分解,因此这对可维护性来说是一个很大的打击。