33,027
社区成员
发帖
与我相关
我的任务
分享全精度大整数阶乘快速算法[带源码]
算法的来源是：
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
但是作者有意的隐藏了一些关键的函数，本来打算写信给那个站长的，可是有碍于面子（写了人家不一定有空理你，有空也不一定会告诉你），只好自己研究了（好在不是很难）
算法源代码：
http://www.77studio.net/files/fac_src.rar
程序里很多很多的部分都是可以再优化的，由于速度已经接近HugeCalc，所以我都是懒得去深化，有兴趣者可以自行研究，有更好的idea请你告诉我
最新的报告，使用这个代码，把相应的"*"改为"*="，套用在gxqcn(GxQ)的CHugeInt上，速度已经和他的HugeCalc 3.0.0.1 beta(07-27版)几乎一摸一样了，当然他的算法还是没有公开的，真不知道两者有多少相同之处，我的有些关键代码都是胡诌出来的（就是通过交流，不断尝试得出的结果），呵呵
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
但是作者有意的隐藏了一些关键的函数，本来打算写信给那个站长的，可是有碍于面子（写了人家不一定有空理你，有空也不一定会告诉你），只好自己研究了（好在不是很难）
算法源代码：
http://www.77studio.net/files/fac_src.rar
程序里很多很多的部分都是可以再优化的，由于速度已经接近HugeCalc，所以我都是懒得去深化，有兴趣者可以自行研究，有更好的idea请你告诉我
最新的报告，使用这个代码，把相应的"*"改为"*="，套用在gxqcn(GxQ)的CHugeInt上，速度已经和他的HugeCalc 3.0.0.1 beta(07-27版)几乎一摸一样了，当然他的算法还是没有公开的，真不知道两者有多少相同之处，我的有些关键代码都是胡诌出来的（就是通过交流，不断尝试得出的结果），呵呵
...全文
请发表友善的回复…
发表回复
junmayang 2004-09-07
- 打赏
- 举报
学习！
yaos 2004-07-30
- 打赏
- 举报
宝宝：
考虑过FNT的具体实现的效率问题么？
a = 2可是最简单，不过对你们的10进制用处不打 ：）
另外，似乎10 ^ 8比10 ^ 9某些方面好
考虑过FNT的具体实现的效率问题么？
a = 2可是最简单，不过对你们的10进制用处不打 ：）
另外，似乎10 ^ 8比10 ^ 9某些方面好
liangbch 2004-07-30
- 打赏
- 举报
shines(郭子)过奖了,其实你的某些算法也很高明。
我在的阶乘计算上可以说是下足了功夫，并且不断改版，改动主要有2方面，一是更优化，速度更快，在敏感的程序段甚至使用汇编，而是算法更简单，条理更清晰。
阶乘计算器这场竞赛，还没有走完，我目前刚刚完成3.2版，准备继续出3个版本，3.5,4.2,4.5,这几个版本的关系如下,3.2为3.0优化版，程序架构比3。0更简单，更合理，但速度比3.0稍慢，3.5在3.2基础上修改而成，即先将n!分解质因数，将所有质因数相乘，尽可能使用平方，就是shine在本帖中提到的方法。4.2和3.2算法几乎相同，仅仅是当乘数大于某个数（超过4096个DWORD）时，采用FNT+CRT算法。5.5和3.5算法几乎相同，仅仅是当乘数大于某个数（超过4096个DWORD）时，采用FNT+CRT算法。
3.2---->3.5
| |
| |
4.2 4.5
我在的阶乘计算上可以说是下足了功夫，并且不断改版，改动主要有2方面，一是更优化，速度更快，在敏感的程序段甚至使用汇编，而是算法更简单，条理更清晰。
阶乘计算器这场竞赛，还没有走完，我目前刚刚完成3.2版，准备继续出3个版本，3.5,4.2,4.5,这几个版本的关系如下,3.2为3.0优化版，程序架构比3。0更简单，更合理，但速度比3.0稍慢，3.5在3.2基础上修改而成，即先将n!分解质因数，将所有质因数相乘，尽可能使用平方，就是shine在本帖中提到的方法。4.2和3.2算法几乎相同，仅仅是当乘数大于某个数（超过4096个DWORD）时，采用FNT+CRT算法。5.5和3.5算法几乎相同，仅仅是当乘数大于某个数（超过4096个DWORD）时，采用FNT+CRT算法。
3.2---->3.5
| |
| |
4.2 4.5
yaos 2004-07-30
- 打赏
- 举报
另外，你实现最小的两数字相乘似乎不必用这么复杂的数据结构啊
可以采用指针数组啊，保存每个数的指针，至于排序，因为都是整数，简单的箱排序足够了，可是O(n)的算法啊，每次取末尾的数字，然后，数组长度减2，插入更简单，应该是O(logn)算法吧
可以采用指针数组啊，保存每个数的指针，至于排序，因为都是整数，简单的箱排序足够了，可是O(n)的算法啊，每次取末尾的数字，然后，数组长度减2，插入更简单，应该是O(logn)算法吧
gxqcn 2004-07-30
- 打赏
- 举报
很遗憾，我没有看到liangbch(宝宝)3.0的代码。
如果你仅仅是为了计算阶乘时调用，“群组乘法”的某些特别算法确实提高不大；但如果参与计算的数非常不“齐整”，它的优势立即可以突显出来。
yaos(等待WoW前玩私服吧·无心人) ：
说实话，HugeCalc v2.x 在计算阶乘时用的算法就是“两两规约”法，而不是一般的分治法，在参与计算的数很“齐整”时其效率还很不错。
如果你仅仅是为了计算阶乘时调用，“群组乘法”的某些特别算法确实提高不大；但如果参与计算的数非常不“齐整”，它的优势立即可以突显出来。
yaos(等待WoW前玩私服吧·无心人) ：
说实话，HugeCalc v2.x 在计算阶乘时用的算法就是“两两规约”法，而不是一般的分治法，在参与计算的数很“齐整”时其效率还很不错。
xiang0123610 2004-07-30
- 打赏
- 举报
哪位能够列出一下这个算法计算几个数的大概时间吗,小弟我现在用的不是自己的电脑,不能测试.比如像1000,10000的阶乘所花的时间大概是多少
shines77 2004-07-30
- 打赏
- 举报
在这个PDF的222页，在“擂台：超大整数高精度快速算法(续)”那个帖子我也提到过
shines77 2004-07-30
- 打赏
- 举报
to yaos:
这个PDF资料不错，Matrix FFT有讲，但不多，但已经讲得很详细了，具体实现要自己来
http://www.jjj.de/fxt/
http://www.jjj.de/fxt/fxtbook.pdf
这个PDF资料不错，Matrix FFT有讲，但不多，但已经讲得很详细了，具体实现要自己来
http://www.jjj.de/fxt/
http://www.jjj.de/fxt/fxtbook.pdf
yaos 2004-07-30
- 打赏
- 举报
这么亲切，这不是浮点FFT流程么 ：）
liangbch 2004-07-30
- 打赏
- 举报
最朴素的fnt正变换和逆变换代码。
void fnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, z, a, b;
if (nn < 2) return;
if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);
mmax = nn >> 1;
while (mmax)
{
istep = mmax << 1;
// Optimize first step when z = 1
for (i = 0; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = a - b;
}
z = w;
for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = z * (a - b);
}
z *= w;
}
w *= w;
mmax >>= 1;
}
if (s) scramble (data, nn);
}
void ifnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, wt, wr, wtemp;
if (nn < 2) return;
if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);
if (s) scramble (data, nn);
mmax = 1;
while (nn > mmax)
{
istep = mmax << 1;
wr = wt = pow (w, nn / istep);
// Optimize first step when wr = 1
for (i = 0; i < nn; i += istep)
{
j = i + mmax;
wtemp = data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}
for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
wtemp = wr * data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}
wr *= wt;
}
mmax = istep;
}
}
下面是使用最普通的 fnt 计算 大数相乘的代码，代码中调用的函数请参照apfloat。
void convolution1 (const DWORD *buf1, const DWORD *buf2, DWORD *prod,size_t nn, int *i)
{
modint *rst1;
modint *rst2;
modint *rst3;
modint *tmp;
size_t n2 = (nn & -nn);
assert(nn == n2);
//----------------------------------------------
tmp = new modint[nn];
//----------------------------------------------
setmodulus (moduli[2]);
rst1 = new modint[nn];
memcpy(rst1,buf1,sizeof(modint)*nn);
fnt (rst1, primitiveroots[2], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[2], 1, nn);
multiplyinplace (rst1, tmp,nn);
fnt (rst1, primitiveroots[2], -1, nn);
fnt_DivN (rst1, nn);
//====================================
setmodulus (moduli[1]);
rst2 = new modint[nn];
memcpy(rst2,buf1,sizeof(modint)*nn);
fnt (rst2, primitiveroots[1], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[1], 1, nn);
multiplyinplace (rst2, tmp,nn);
fnt (rst2, primitiveroots[1], -1, nn);
fnt_DivN (rst2, nn);
//====================================
setmodulus (moduli[0]);
rst3 = new modint[nn];
memcpy(rst3,buf1,sizeof(modint)*nn);
fnt (rst3, primitiveroots[0], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[0], 1, nn);
multiplyinplace (rst3, tmp,nn);
fnt (rst3, primitiveroots[0], -1, nn);
fnt_DivN (rst3, nn);
//-------------------------------------------------
*i = carrycrt (rst3, rst2, rst1, nn);
memcpy(prod,rst3,sizeof(modint)*nn);
delete[] rst1;
delete[] rst2;
delete[] rst3;
delete[] tmp;
}
void fnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, z, a, b;
if (nn < 2) return;
if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);
mmax = nn >> 1;
while (mmax)
{
istep = mmax << 1;
// Optimize first step when z = 1
for (i = 0; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = a - b;
}
z = w;
for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = z * (a - b);
}
z *= w;
}
w *= w;
mmax >>= 1;
}
if (s) scramble (data, nn);
}
void ifnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, wt, wr, wtemp;
if (nn < 2) return;
if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);
if (s) scramble (data, nn);
mmax = 1;
while (nn > mmax)
{
istep = mmax << 1;
wr = wt = pow (w, nn / istep);
// Optimize first step when wr = 1
for (i = 0; i < nn; i += istep)
{
j = i + mmax;
wtemp = data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}
for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
wtemp = wr * data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}
wr *= wt;
}
mmax = istep;
}
}
下面是使用最普通的 fnt 计算 大数相乘的代码，代码中调用的函数请参照apfloat。
void convolution1 (const DWORD *buf1, const DWORD *buf2, DWORD *prod,size_t nn, int *i)
{
modint *rst1;
modint *rst2;
modint *rst3;
modint *tmp;
size_t n2 = (nn & -nn);
assert(nn == n2);
//----------------------------------------------
tmp = new modint[nn];
//----------------------------------------------
setmodulus (moduli[2]);
rst1 = new modint[nn];
memcpy(rst1,buf1,sizeof(modint)*nn);
fnt (rst1, primitiveroots[2], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[2], 1, nn);
multiplyinplace (rst1, tmp,nn);
fnt (rst1, primitiveroots[2], -1, nn);
fnt_DivN (rst1, nn);
//====================================
setmodulus (moduli[1]);
rst2 = new modint[nn];
memcpy(rst2,buf1,sizeof(modint)*nn);
fnt (rst2, primitiveroots[1], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[1], 1, nn);
multiplyinplace (rst2, tmp,nn);
fnt (rst2, primitiveroots[1], -1, nn);
fnt_DivN (rst2, nn);
//====================================
setmodulus (moduli[0]);
rst3 = new modint[nn];
memcpy(rst3,buf1,sizeof(modint)*nn);
fnt (rst3, primitiveroots[0], 1, nn);
memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[0], 1, nn);
multiplyinplace (rst3, tmp,nn);
fnt (rst3, primitiveroots[0], -1, nn);
fnt_DivN (rst3, nn);
//-------------------------------------------------
*i = carrycrt (rst3, rst2, rst1, nn);
memcpy(prod,rst3,sizeof(modint)*nn);
delete[] rst1;
delete[] rst2;
delete[] rst3;
delete[] tmp;
}
yaos 2004-07-30
- 打赏
- 举报
监测L2 Cache，一个cpuid足够了阿
yaos 2004-07-30
- 打赏
- 举报
在考虑短除法（n DWORD / DWORD），有点麻烦，你们有什么经验
普通的短除法汇编俺会，说的是倒数后相乘，发现不太容易
普通的短除法汇编俺会，说的是倒数后相乘，发现不太容易
yaos 2004-07-30
- 打赏
- 举报
:)
高德纳德书上有么？
那里有资料？
高德纳德书上有么？
那里有资料？
shines77 2004-07-30
- 打赏
- 举报
to yaos:
这个你就不懂了，FNT六步变换和FFT的6步矩阵变换是一样的（搬过来的），变换本身我了解，但是对于w的模变换我还不是很清楚，w即类似FFT中的w(变换为cos,sin的复数)
如果你看了FFT的文章，就知道6步变换是为了适应L1, L2 Cache，为了什么要置换矩阵（因为在处理连续的数据可以利用Cache Line），为什么要用矩阵的FFT，因为这样可以降低到纯粹的FFT的长度由N变为Sqrt(N)，矩阵FFT的实现和你上面说的FNT(其实也是FFT)的变换是
没有冲突的，矩阵FFT/FNT只是把FFT/FNT分解来计算而已，总体结构还是和你上面说的一样的
f(a) = FNT(a)
f(b) = FNT(b)
f(c) = f(a) * f(b)
c = IFNT(c)
则 c = a * b
题外话：apfloat的程序依赖机器Cache等参数，却没有自己检测硬件的能力，如果加上效率会更好一些，P4的L1 Cache是64字节对齐的，不是32
这个你就不懂了，FNT六步变换和FFT的6步矩阵变换是一样的（搬过来的），变换本身我了解，但是对于w的模变换我还不是很清楚，w即类似FFT中的w(变换为cos,sin的复数)
如果你看了FFT的文章，就知道6步变换是为了适应L1, L2 Cache，为了什么要置换矩阵（因为在处理连续的数据可以利用Cache Line），为什么要用矩阵的FFT，因为这样可以降低到纯粹的FFT的长度由N变为Sqrt(N)，矩阵FFT的实现和你上面说的FNT(其实也是FFT)的变换是
没有冲突的，矩阵FFT/FNT只是把FFT/FNT分解来计算而已，总体结构还是和你上面说的一样的
f(a) = FNT(a)
f(b) = FNT(b)
f(c) = f(a) * f(b)
c = IFNT(c)
则 c = a * b
题外话：apfloat的程序依赖机器Cache等参数，却没有自己检测硬件的能力，如果加上效率会更好一些，P4的L1 Cache是64字节对齐的，不是32
yaos 2004-07-30
- 打赏
- 举报
：）
看来apfloat的FNT效率不高啊
标准的FNT不用这么麻烦，就是编码可能多
f(a) = FNT(a)
f(b) = FNT(b)
f(c) = f(a) * f(b)
c = IFNT(c)
则 c = a * b
看来apfloat的FNT效率不高啊
标准的FNT不用这么麻烦，就是编码可能多
f(a) = FNT(a)
f(b) = FNT(b)
f(c) = f(a) * f(b)
c = IFNT(c)
则 c = a * b
liangbch 2004-07-30
- 打赏
- 举报
我的FNT是由apfloat改的，FNT六步变换我也没有看懂。等以后有时间可以考虑自己重写一个FNT算法。
FNT 乘法的基本过程是这样的，计算 计算数组a[i](i=0..N-1) *B[i](i=0..N-1)
1.将a 使用模m2,源根r2 做一个长度为2N的变换，在变换前，将后半部分补0
2.将b 使用模m2,源根r2 做一个长度为2N的变换，在变换前，将后半部分补0
3.计算 c2[i](0=0..2N-1)=a[i]*b[i]
4.将c2 进行逆变换，仍然叫做 c2
5.将a 使用模m1,源根r1 做一个长度为2N的变换，在变换前，将后半部分补0
6.将b 使用模m1,源根r1 做一个长度为2N的变换，在变换前，将后半部分补0
7.计算 c1[i](0=0..2N-1)=a[i]*b[i]
8.将c1 进行逆变换，仍然叫做 c1
9. 将a 使用模m0,源根r0 做一个长度为2N的变换，在变换前，将后半部分补0
10.将b 使用模m0,源根r0 做一个长度为2N的变换，在变换前，将后半部分补0
11.计算 c0[i](0=0..2N-1)=a[i]*b[i]
12.将c0 进行你变换，仍然叫做 c0
13 利用CRT，通过c0,c1,c2 求a*b的乘积 c.
在对序列做正变换或者逆变换时，使用6步变换算法。
使用10^8 造成序列增长，速度变慢，空间耗费增加。另外，在我的程序中，大数始终是一每个DWORD表示9位十进制数的方式表示的，如果此处使用10^8,必然导致和数的存储不兼容，所以需要修改现有的数的表示法和相关程序。固个人觉得使用10^8意义不大。
FNT 乘法的基本过程是这样的，计算 计算数组a[i](i=0..N-1) *B[i](i=0..N-1)
1.将a 使用模m2,源根r2 做一个长度为2N的变换，在变换前，将后半部分补0
2.将b 使用模m2,源根r2 做一个长度为2N的变换，在变换前，将后半部分补0
3.计算 c2[i](0=0..2N-1)=a[i]*b[i]
4.将c2 进行逆变换，仍然叫做 c2
5.将a 使用模m1,源根r1 做一个长度为2N的变换，在变换前，将后半部分补0
6.将b 使用模m1,源根r1 做一个长度为2N的变换，在变换前，将后半部分补0
7.计算 c1[i](0=0..2N-1)=a[i]*b[i]
8.将c1 进行逆变换，仍然叫做 c1
9. 将a 使用模m0,源根r0 做一个长度为2N的变换，在变换前，将后半部分补0
10.将b 使用模m0,源根r0 做一个长度为2N的变换，在变换前，将后半部分补0
11.计算 c0[i](0=0..2N-1)=a[i]*b[i]
12.将c0 进行你变换，仍然叫做 c0
13 利用CRT，通过c0,c1,c2 求a*b的乘积 c.
在对序列做正变换或者逆变换时，使用6步变换算法。
使用10^8 造成序列增长，速度变慢，空间耗费增加。另外，在我的程序中，大数始终是一每个DWORD表示9位十进制数的方式表示的，如果此处使用10^8,必然导致和数的存储不兼容，所以需要修改现有的数的表示法和相关程序。固个人觉得使用10^8意义不大。
shines77 2004-07-29
- 打赏
- 举报
to gxqcn(GxQ):
没啦，我不是那个意思，我只是觉得在你和liangbch(宝宝)的讨论中，我从他身上获得的灵感更多一些，有些时候你讲了我没有注意，他讲了我就比较容易领悟，他总是给人一种特别的feeling，尤其是他的3。0那段代码，虽然我想我知道是什么意思，但我一般是不会那样写的，有种很特别的感觉，清新优雅，我可能直来直去的比较多（思维简单一些）
至于“群组乘法”（计算许多数的连乘），我也考虑过很多的，但是基本上要么过于复杂，真要求一个最优解也许不难（程序可能很复杂），但是值不值得就值得考虑了，所以我也是想得多，实践中没有去做太多尝试，而且从实践看来，这个部分对最终的效率影响不是很大（粗略的看）
没啦，我不是那个意思，我只是觉得在你和liangbch(宝宝)的讨论中，我从他身上获得的灵感更多一些，有些时候你讲了我没有注意，他讲了我就比较容易领悟，他总是给人一种特别的feeling，尤其是他的3。0那段代码，虽然我想我知道是什么意思，但我一般是不会那样写的，有种很特别的感觉，清新优雅，我可能直来直去的比较多（思维简单一些）
至于“群组乘法”（计算许多数的连乘），我也考虑过很多的，但是基本上要么过于复杂，真要求一个最优解也许不难（程序可能很复杂），但是值不值得就值得考虑了，所以我也是想得多，实践中没有去做太多尝试，而且从实践看来，这个部分对最终的效率影响不是很大（粗略的看）
yaos 2004-07-29
- 打赏
- 举报
两两规约呢？
就是相邻的整数相乘，重复直到出结果
是不是不如你的快？
就是相邻的整数相乘，重复直到出结果
是不是不如你的快？
gxqcn 2004-07-29
- 打赏
- 举报
to shines(郭子)：
请别误会，我没说你从我那里获得了灵感！你的水平我是比较佩服的，英雄所见略同嘛！（把自己也夸奖了一下）：）毕竟大家来这里是为了交流，而不应计较太多别的事物，能达到“双赢”的目的那时你我共同的心愿。（我自认为，通过这些交流，已达到了相互促进的目的，对此，我心常怀感激。）
关于“群组乘法”（计算许多数的连乘）算法问题，我的宗旨和思路与你的不一样，不妨也抛出来让大家评论一下（申明：为原创）。
我的算法宗旨是：始终处理当前长度最小的两个数。
实现原理：
1、用一个结构保存该“群组”，并按长度排序；
2、从结构中pop最小的两个元素，将积再push进该结构中，并保持原递增/递减顺序；
3、重复2，直到该结构中仅剩一个元素，其值即为所求。
数据结构：
因上述过程需频繁地增/减元素及排序，所以我采用了非常复杂的红-黑树数据结构。
实验&讨论：
当参与计算的“群组”长度相近时，上述算法略逊于普通分治法；当参与计算的“群组”长度波动很大时，该算法比普通分治法性能提高很多。
所以，在即将发布的 HugeCalc v3.0.0.1 中，对内使用“群组乘法”对内部的阶乘实现采用普通分治法，而供外部调用的则采用上述算法（它具有很稳定的效能）。
请别误会，我没说你从我那里获得了灵感！你的水平我是比较佩服的，英雄所见略同嘛！（把自己也夸奖了一下）：）毕竟大家来这里是为了交流，而不应计较太多别的事物，能达到“双赢”的目的那时你我共同的心愿。（我自认为，通过这些交流，已达到了相互促进的目的，对此，我心常怀感激。）
关于“群组乘法”（计算许多数的连乘）算法问题，我的宗旨和思路与你的不一样，不妨也抛出来让大家评论一下（申明：为原创）。
我的算法宗旨是：始终处理当前长度最小的两个数。
实现原理：
1、用一个结构保存该“群组”，并按长度排序；
2、从结构中pop最小的两个元素，将积再push进该结构中，并保持原递增/递减顺序；
3、重复2，直到该结构中仅剩一个元素，其值即为所求。
数据结构：
因上述过程需频繁地增/减元素及排序，所以我采用了非常复杂的红-黑树数据结构。
实验&讨论：
当参与计算的“群组”长度相近时，上述算法略逊于普通分治法；当参与计算的“群组”长度波动很大时，该算法比普通分治法性能提高很多。
所以，在即将发布的 HugeCalc v3.0.0.1 中，对内使用“群组乘法”对内部的阶乘实现采用普通分治法，而供外部调用的则采用上述算法（它具有很稳定的效能）。
shines77 2004-07-29
- 打赏
- 举报
omegaSwingHighBound(n) 是每一层中multiList的元素个数，
如果上面的
第7层：2*11*13*17*19*29*31*53*59*61*67*71*73*79*83*89*97 '质数元素个数是17个
本来omegaSwingHighBound(n)这个函数只存在于FactorialPrimeSwingLuschny中，他的值等于FactorialPrimePartSchoenhage中的
unsigned long piN = 2 + (15 * n) / (8 * (lgN - 1));
在FactorialPrimePartSchoenhage中，原来的意思是比PI(X)（计算从0-X范围内的质数个数的函数，在讨论质数的帖子里有提过）稍大一个值，所以叫piN，因为这个算法中的primeList的长度至少等于PI(X)，计算过程中primeList还会变长，所以这个piN值比PI(X)稍大
在FactorialPrimeSwingLuschny和FactorialPrimeExpBinarySplit中这个omegaSwingHighBound(n)那个网站的作者并没有提供，但我们知道，在前者中，omegaSwingHighBound(n)最大值是比PI(X)稍大的值，而后者FactorialPrimeExpBinarySplit中肯定小于PI(X)，所以我们可以借用piN的值，稍微浪费一点内存而已，如果能够知道确切的omegaSwingHighBound(n)会更好，也许可以写个email去问问作者，相信他会乐意提供的
to gxqcn(GxQ)：
其实累乘的部分都是大家讨论中交流得到的，才有了上面的代码，我倒不是看了你那个《高精度快速阶乘算法》帖子，主要还是原来那个帖子的讨论，当然，思路都是一样的
我为了用简单的方式让相乘的数尽量的大小接近，还写了一个快速排序quickSort后，然后balanceList（即排序后，第一个和倒数第2个交换，第3个和倒数第4个交换。。。。），达到一个基本上（注意是基本上）前半部乘积和后半部乘积大小相仿的效果
但是效果非常不明显，快速排序和交换的时间损失基本上和平衡获得提升基本抵消（稍快了一小丁点），所以没用，感觉在这方面的设计优化不多，过于复杂反而难处理
如果上面的
第7层：2*11*13*17*19*29*31*53*59*61*67*71*73*79*83*89*97 '质数元素个数是17个
本来omegaSwingHighBound(n)这个函数只存在于FactorialPrimeSwingLuschny中，他的值等于FactorialPrimePartSchoenhage中的
unsigned long piN = 2 + (15 * n) / (8 * (lgN - 1));
在FactorialPrimePartSchoenhage中，原来的意思是比PI(X)（计算从0-X范围内的质数个数的函数，在讨论质数的帖子里有提过）稍大一个值，所以叫piN，因为这个算法中的primeList的长度至少等于PI(X)，计算过程中primeList还会变长，所以这个piN值比PI(X)稍大
在FactorialPrimeSwingLuschny和FactorialPrimeExpBinarySplit中这个omegaSwingHighBound(n)那个网站的作者并没有提供，但我们知道，在前者中，omegaSwingHighBound(n)最大值是比PI(X)稍大的值，而后者FactorialPrimeExpBinarySplit中肯定小于PI(X)，所以我们可以借用piN的值，稍微浪费一点内存而已，如果能够知道确切的omegaSwingHighBound(n)会更好，也许可以写个email去问问作者，相信他会乐意提供的
to gxqcn(GxQ)：
其实累乘的部分都是大家讨论中交流得到的，才有了上面的代码，我倒不是看了你那个《高精度快速阶乘算法》帖子，主要还是原来那个帖子的讨论，当然，思路都是一样的
我为了用简单的方式让相乘的数尽量的大小接近，还写了一个快速排序quickSort后，然后balanceList（即排序后，第一个和倒数第2个交换，第3个和倒数第4个交换。。。。），达到一个基本上（注意是基本上）前半部乘积和后半部乘积大小相仿的效果
但是效果非常不明显，快速排序和交换的时间损失基本上和平衡获得提升基本抵消（稍快了一小丁点），所以没用，感觉在这方面的设计优化不多，过于复杂反而难处理
加载更多回复（7）
