33,010
社区成员
发帖
与我相关
我的任务
分享自守数的计算探讨
何谓自守数？
在郭先强先生编写的Automorphic Number程序中进行了如下的定义：
【自守数】(在十进位制中,) 若一个 k 位正整数 N (可含前置 "0" ), 若满足如下性质: 任意两个或多个均以该字串 N 结尾的整数相乘, 其结果的最后 k 位数字一定还是 N, 那么, 则称 N 为 "k 位自守数".
前几日偶然访问到下面的页面
http://zhjyx.hfjy.net.cn/Resource/Book/Edu/SZJY/TS007047/0004_ts007047.htm
在这个页面中介绍了下面的方程的两个无限解（0和1是两个显然的解）
x^2=x
由此受到启发，编写的一段代码，经过一周左右的反复修改，实现了用计算机来快速求小于等于1,000,000,000的自守数。（事实上以现在的计算机速度和内存，计算10^8以下的数据比较现实一些）
我的代码是以VC编写的一个控制台程序，使用内联汇编，未采用MMX或SSE等指令集，只使用了普能的X86指令。速度与郭先生的程序相比，在位数较小时，比如10万以下，略慢一点(小于1s)，但在50万以上，速度略快，大约快2－3%左右。本人对MMX和SSE指令不熟悉，如果用其来优化核心代码，可能速度还能提高。
算法的核心就是如何求5的2＾n次方（但不需要求出全部数据，而只需要取低n位十进制数）。
如有兴趣者，可以查看上面的页面，也可与我讨论.需要代码请发邮件给我:
myh9999@vip.sina.com
待我整理好代码后，将贴出我的代码
在郭先强先生编写的Automorphic Number程序中进行了如下的定义：
【自守数】(在十进位制中,) 若一个 k 位正整数 N (可含前置 "0" ), 若满足如下性质: 任意两个或多个均以该字串 N 结尾的整数相乘, 其结果的最后 k 位数字一定还是 N, 那么, 则称 N 为 "k 位自守数".
前几日偶然访问到下面的页面
http://zhjyx.hfjy.net.cn/Resource/Book/Edu/SZJY/TS007047/0004_ts007047.htm
在这个页面中介绍了下面的方程的两个无限解（0和1是两个显然的解）
x^2=x
由此受到启发，编写的一段代码，经过一周左右的反复修改，实现了用计算机来快速求小于等于1,000,000,000的自守数。（事实上以现在的计算机速度和内存，计算10^8以下的数据比较现实一些）
我的代码是以VC编写的一个控制台程序，使用内联汇编，未采用MMX或SSE等指令集，只使用了普能的X86指令。速度与郭先生的程序相比，在位数较小时，比如10万以下，略慢一点(小于1s)，但在50万以上，速度略快，大约快2－3%左右。本人对MMX和SSE指令不熟悉，如果用其来优化核心代码，可能速度还能提高。
算法的核心就是如何求5的2＾n次方（但不需要求出全部数据，而只需要取低n位十进制数）。
如有兴趣者，可以查看上面的页面，也可与我讨论.需要代码请发邮件给我:
myh9999@vip.sina.com
待我整理好代码后，将贴出我的代码
...全文
请发表友善的回复…
发表回复
mathe 2004-11-09
- 打赏
- 举报
20s now. So 27% speedup.
gxqcn 2004-11-08
- 打赏
- 举报
我又改进了算法，在我的机器上计算 n=2^24 时速度提高了 30% 左右。
to myh9999 (虎子) & mathe():
我已将最新版发送给你们。
to mathe():
由于没有 ftp 权限，所以无法下载Intel C/C++编译器和VTune。
可否请你对新版再用 VTune 进行测试？谢谢！
to myh9999 (虎子) & mathe():
我已将最新版发送给你们。
to mathe():
由于没有 ftp 权限，所以无法下载Intel C/C++编译器和VTune。
可否请你对新版再用 VTune 进行测试？谢谢！
mathe 2004-11-05
- 打赏
- 举报
能不能提供一个没有图形界面的，像我的程序一样，只有命令行的，这样分析程序的性能比较简单。
比如我的程序，在n=10^7 （或2^24)时，估计有一半的时间花费在访问内存上。其中L2 Miss次数大概在100M
比如我的程序，在n=10^7 （或2^24)时，估计有一半的时间花费在访问内存上。其中L2 Miss次数大概在100M
mathe 2004-11-05
- 打赏
- 举报
呵呵，计算出10^8，花费2分46.516秒，但是不知道为什么，结果没有输出出来，估计可能输出列表框出问题，放不下了。计算过程中最多也就花700多M内存，不错:)
mathe 2004-11-05
- 打赏
- 举报
呵呵，不过你的程序很好，内存消耗很小，在我的计算机上计算n=10^8肯定没有问题，我的程序就Outofmemory了:)
mathe 2004-11-05
- 打赏
- 举报
呵呵，我知道可以那样计算，不过我偷懒，不愿意去修改。
在我的Windows 机器上做了试验，你的程序使用的内存明显比较少。
CPU 2.8G, L2 Cache 512M, Bus Clock 400M
在没有运行我们的程序时，内存使用450M
运行我的程序后，内存使用900M,时间15.436s (包含文件IO时间）
运行你的程序后，内存使用493M,时间17.141s
由于我的程序没有做位数优化，所以我有做了下n=2^24=16777216时候的比较
运行我的程序后，内存不变,时间是16.156s (同前面的唯一差别应该是在IO上）
运行你的程序后，内存使用534M,时间27.375s
在我的Windows 机器上做了试验，你的程序使用的内存明显比较少。
CPU 2.8G, L2 Cache 512M, Bus Clock 400M
在没有运行我们的程序时，内存使用450M
运行我的程序后，内存使用900M,时间15.436s (包含文件IO时间）
运行你的程序后，内存使用493M,时间17.141s
由于我的程序没有做位数优化，所以我有做了下n=2^24=16777216时候的比较
运行我的程序后，内存不变,时间是16.156s (同前面的唯一差别应该是在IO上）
运行你的程序后，内存使用534M,时间27.375s
xdspower 2004-11-05
- 打赏
- 举报
我觉得自守数表示为表达式时需要定义一些算术符号，比如位数运算{x},表示十进制数x的位数,高位可以有不定个数的0，但整个数是确定位数的，比如{0001}=4，{001}=3，则定义自守数可以表示为一个
对于非负整数N，有{N}=k,(N^2 mod 10^k）=N
一位的自守数有0,1,5,6
对于非负整数N，有{N}=k,(N^2 mod 10^k）=N
一位的自守数有0,1,5,6
gxqcn 2004-11-05
- 打赏
- 举报
其实，可以反推每次的有效位数。比如：计算 n = 10^7，
可以记录如下值：10000000，5000000，2500000，1250000，625000，312500，156250，156250，78125，39063，19532，9766，4883，2442，...，直到达到一个可直接初始化的值。
而后反算即可。
to mathe():
我已将最新版发到你的信箱，请注意查收。请你在你的系统上运行测试，并希望得到相应的测试结果，如计算速度，输入、输出耗时，内存消耗等，谢谢！
可以记录如下值：10000000，5000000，2500000，1250000，625000，312500，156250，156250，78125，39063，19532，9766，4883，2442，...，直到达到一个可直接初始化的值。
而后反算即可。
to mathe():
我已将最新版发到你的信箱，请注意查收。请你在你的系统上运行测试，并希望得到相应的测试结果，如计算速度，输入、输出耗时，内存消耗等，谢谢！
mathe 2004-11-05
- 打赏
- 举报
现在我用的算法还有一点比较浪费，就是计算10^7时，实际上计算的是
2^24=16777216位，所以很浪费，
比如如果只计算
2^23=8388608位，那么，内存就可以节省一半多的时间了。
其实如果计算时不是每次都增加一倍的位数，那么可以加快很多的。
2^24=16777216位，所以很浪费，
比如如果只计算
2^23=8388608位，那么，内存就可以节省一半多的时间了。
其实如果计算时不是每次都增加一倍的位数，那么可以加快很多的。
mathe 2004-11-05
- 打赏
- 举报
要看你如何进行乘法运算。在使用FFT进行乘法运算时，计算平方花费的时间是普通乘法的2/3,所以
用(x^2-1)^2 mod 10^(2k) 会比其他公式要快一些。
我上面的程序主要的缺点是内存消耗比较大，原因很简单，我使用了Intel的通用的FFT程序，程序没有专门为数据量很大的时候考虑过，所以在进行大数运算时，消耗内存过大。而这时，内存访问的速度成为主要制约因素了。其实这时，决定程序运行速度的主要因素在于内存访问的速度，而不是CPU的速度。不过我用的计算机是服务器 （Intel Xeon CPU)，L2 Cache的大小是512K,所以在内存访问的速度上肯定也比你的机器占优势的。
从运行时间的比较，在n=5,6,7,基本上是我的计算机上的数据比你的快了一倍，所以算法复杂度看来
是一样的，很可能运行速度差不多（假设我的计算机综合性能是你的一倍是比较合理的）。但是我现在的程序对内存的要求比较高，这个是明显的缺点。如果在你的计算机上运行，对于10^7肯定会比较慢，因为你的计算机的内存只有256K,如果有512K，那就会好很多。
我记得用DCT(离散余弦变换）也应该可以做卷积运算的，如果改成用DCT运算应该可以节省一倍的内存的，但是对于如何用DCT做卷积我不是很确定，所以没有做试验。
你应该是用Windows吧，可以从：
https://shale.intel.com/EvalCenter/EvalForm.aspx?ProductID=254
下载Intel的MKL for Windows, 要求有个Email地址，Intel会将一个试用的版权文件发到你的信箱的，里面有如何下载和使用MKL的说明。
然后还要到http://www.intel.com/software/products/noncom/
下载MKL的库文件。
用(x^2-1)^2 mod 10^(2k) 会比其他公式要快一些。
我上面的程序主要的缺点是内存消耗比较大，原因很简单，我使用了Intel的通用的FFT程序，程序没有专门为数据量很大的时候考虑过，所以在进行大数运算时，消耗内存过大。而这时，内存访问的速度成为主要制约因素了。其实这时，决定程序运行速度的主要因素在于内存访问的速度，而不是CPU的速度。不过我用的计算机是服务器 （Intel Xeon CPU)，L2 Cache的大小是512K,所以在内存访问的速度上肯定也比你的机器占优势的。
从运行时间的比较，在n=5,6,7,基本上是我的计算机上的数据比你的快了一倍，所以算法复杂度看来
是一样的，很可能运行速度差不多（假设我的计算机综合性能是你的一倍是比较合理的）。但是我现在的程序对内存的要求比较高，这个是明显的缺点。如果在你的计算机上运行，对于10^7肯定会比较慢，因为你的计算机的内存只有256K,如果有512K，那就会好很多。
我记得用DCT(离散余弦变换）也应该可以做卷积运算的，如果改成用DCT运算应该可以节省一倍的内存的，但是对于如何用DCT做卷积我不是很确定，所以没有做试验。
你应该是用Windows吧，可以从：
https://shale.intel.com/EvalCenter/EvalForm.aspx?ProductID=254
下载Intel的MKL for Windows, 要求有个Email地址，Intel会将一个试用的版权文件发到你的信箱的，里面有如何下载和使用MKL的说明。
然后还要到http://www.intel.com/software/products/noncom/
下载MKL的库文件。
gxqcn 2004-11-05
- 打赏
- 举报
受益非浅，谢谢！
mathe 2004-11-05
- 打赏
- 举报
一般来说，用汇编不见得会有很好的效果，除非你对计算机体系结构非常了解，即使这样，很多东西还是人很难理解的。用好的编译器得到的效果可能会比用自己写汇编程序会好。当然，在写c程序时如何写代码也是比较重要的，如果你能够对编译器比较了解，会有很大的帮助。
上面的数据是用Intel的VTune来分析的，通上面发给你的网址上相关链接也可以找到VTune.还有Intel C/C++编译器也不错，对优化程序肯定有帮助。上面还有一些现成的算法也非常有用(比如我直接使用了FFT）。
上面cycles是指整个程序运行的时钟周期数目，其实这个从程序运行的时间也可以估算出来，比如机器主频2.8G,运行45.8G个cycles,那么就运行了45.8G/2.8G秒了。
instructions是只总共执行的机器指令数目（同汇编指令一一对应），而循环里面被反复执行的指令要重复计算的。由于现在的机器一般都有好多个处理单元，如果内存访问比例不高，优化的好的程序的IPC (instructions/cycles)会大于1的。 L1,L2 miss是指计算机访问内存数据时，这个数据没有在对应的Cache里面，所以计算机要首先将数据从内存读入Cache,或者从L2 Cache读入L1 Cache.而计算机直接访问内存是非常慢的，反问L1要比L2快很多。上面两个程序的Cache Miss Ratio很高（将近100%)，也就是说，大部分需要访问的数据都不在Cache里面。DTLB是一个用来将线性地址（计算机中程序所以的内存地址）转化为物理地址（实际内存地址）的表格。同样这个表格也挺大的，不能全部放入Cache,所以也会有Miss,这时，访问内存就更慢了。
上面的数据是用Intel的VTune来分析的，通上面发给你的网址上相关链接也可以找到VTune.还有Intel C/C++编译器也不错，对优化程序肯定有帮助。上面还有一些现成的算法也非常有用(比如我直接使用了FFT）。
上面cycles是指整个程序运行的时钟周期数目，其实这个从程序运行的时间也可以估算出来，比如机器主频2.8G,运行45.8G个cycles,那么就运行了45.8G/2.8G秒了。
instructions是只总共执行的机器指令数目（同汇编指令一一对应），而循环里面被反复执行的指令要重复计算的。由于现在的机器一般都有好多个处理单元，如果内存访问比例不高，优化的好的程序的IPC (instructions/cycles)会大于1的。 L1,L2 miss是指计算机访问内存数据时，这个数据没有在对应的Cache里面，所以计算机要首先将数据从内存读入Cache,或者从L2 Cache读入L1 Cache.而计算机直接访问内存是非常慢的，反问L1要比L2快很多。上面两个程序的Cache Miss Ratio很高（将近100%)，也就是说，大部分需要访问的数据都不在Cache里面。DTLB是一个用来将线性地址（计算机中程序所以的内存地址）转化为物理地址（实际内存地址）的表格。同样这个表格也挺大的，不能全部放入Cache,所以也会有Miss,这时，访问内存就更慢了。
gxqcn 2004-11-05
- 打赏
- 举报
to mathe():
可否解释一下上面一排各项指标的含义？以及介绍一下如何使它们各自向好的方向优化的技术与经验？另外，这些数据是如何测得的？我想为以后的开发引入更科学的评判数据，谢谢！
可否解释一下上面一排各项指标的含义？以及介绍一下如何使它们各自向好的方向优化的技术与经验？另外，这些数据是如何测得的？我想为以后的开发引入更科学的评判数据，谢谢！
gxqcn 2004-11-05
- 打赏
- 举报
上述报告我未完全看明白，但大致了解了一个梗概，谢谢！
需要说明的是，我的程序虽然用汇编优化过(包括输出我都做过特别优化)，但却未用到MMX和SSE指令，与硬件的衔接不够，这是我的一个薄弱环节。
另外，高精度乘法会根据不同规模选用不同的算法，其临界点的设定以及算法中的某些常量只是满足在我的机器上接近最佳，比方说我未充分用到 L1/L2 Cache.
需要说明的是，我的程序虽然用汇编优化过(包括输出我都做过特别优化)，但却未用到MMX和SSE指令，与硬件的衔接不够，这是我的一个薄弱环节。
另外，高精度乘法会根据不同规模选用不同的算法，其临界点的设定以及算法中的某些常量只是满足在我的机器上接近最佳，比方说我未充分用到 L1/L2 Cache.
mathe 2004-11-05
- 打赏
- 举报
上面是n=2^24的数据：cycles/2.8G就是运行时间。
一般来说，L2 Miss的代价在100-200个cycles,L1 Miss可能10个cycles以内,DTLB Miss的代价比L2还要大，应该在200以上。CPU主频在2.8G,可以看出调用Intel FFT的程序时间基本上花费在内存访问上，所以它的Instruction/Cycles非常低，才0.2左右。所以这个代码在计算上已经没法优化了，如果要优化应该在内存访问上做优化。
gxqcn的代码花费在内存访问上的时间大概在20％吧（估计的，肯定不准）。其中DTLB miss ratio特别高，是不是说明访问数据比FFT版本更加随机。现在计算所占的比例还是比较高的，可以试一试换一个好一点的编译器看看，我觉得可以有所提高的。
一般来说，L2 Miss的代价在100-200个cycles,L1 Miss可能10个cycles以内,DTLB Miss的代价比L2还要大，应该在200以上。CPU主频在2.8G,可以看出调用Intel FFT的程序时间基本上花费在内存访问上，所以它的Instruction/Cycles非常低，才0.2左右。所以这个代码在计算上已经没法优化了，如果要优化应该在内存访问上做优化。
gxqcn的代码花费在内存访问上的时间大概在20％吧（估计的，肯定不准）。其中DTLB miss ratio特别高，是不是说明访问数据比FFT版本更加随机。现在计算所占的比例还是比较高的，可以试一试换一个好一点的编译器看看，我觉得可以有所提高的。
mathe 2004-11-05
- 打赏
- 举报
cycles instructions UseMMX UseSIMD L1 Cache Miss L2 Cache Miss L1 Miss Ratio L2 Miss Ratio DTLB Miss DTLB Miss Ratio
Intel FFT 45.8G 9.55G Yes NO 568M 177M 95.08% 95.89% 41.8M 82.10%
gxqcn's code 76.2G 46.22G NO NO 421M 115M 91.39% 92.58% 112M 88.78%
Intel FFT 45.8G 9.55G Yes NO 568M 177M 95.08% 95.89% 41.8M 82.10%
gxqcn's code 76.2G 46.22G NO NO 421M 115M 91.39% 92.58% 112M 88.78%
gxqcn 2004-11-05
- 打赏
- 举报
to mathe():
AutomorphicNumber_Console.exe 已发给你了，请分析其性能。
AutomorphicNumber_Console.exe 已发给你了，请分析其性能。
gxqcn 2004-11-04
- 打赏
- 举报
to mathe():
我曾用过 (x^2-1)^2 mod 10^2k 去计算，它的优点是代码简单，但效率欠佳，因为第二次平方将浪费太多。
可否与我的程序进行一下横向对比？我想借此机会检验一下我的算法库，包括计算时间、输入输出、内存消耗等各个方面。
我的 E-mail 是 gxqcn@163.com. 我可以将最新的 AutomorphicNumber.exe 及算法库 HugeCalc 提供给你使用；或者，你可以把你的程序及那个“Intel的Math Kernel Library”传给我，以便在相同的操作环境下对比。（从你提供的数据分析，我的 HugeCalc 似乎更优一点，不过这需要同台对比才可下定论。）
我曾用过 (x^2-1)^2 mod 10^2k 去计算，它的优点是代码简单，但效率欠佳，因为第二次平方将浪费太多。
可否与我的程序进行一下横向对比？我想借此机会检验一下我的算法库，包括计算时间、输入输出、内存消耗等各个方面。
我的 E-mail 是 gxqcn@163.com. 我可以将最新的 AutomorphicNumber.exe 及算法库 HugeCalc 提供给你使用；或者，你可以把你的程序及那个“Intel的Math Kernel Library”传给我，以便在相同的操作环境下对比。（从你提供的数据分析，我的 HugeCalc 似乎更优一点，不过这需要同台对比才可下定论。）
mathe 2004-11-04
- 打赏
- 举报
重新运行了一次n^7,时间14.77s,其中文件输出占1.1s.
mathe 2004-11-04
- 打赏
- 举报
试着下载了一个Intel的Math Kernel Library, 用里面的FFT来计算平方运算，使用算法：
(x^2-1)^2 (mod 10^2k)
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mkl_fft.h>
#define MOD (10000)
double *r;
double *tmp;
#ifdef _WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif
int main()
{
int i,n;
int k,k1,k2,c2;
scanf("%d",&k);
if(k<=1)return -1;
k2=k,c2=0;
while(k2){
c2++;k2>>=1;
}
k2=1<<(c2-1); //k2<=k and k2 is power of 2.
if(k2<k)k2<<=1;
r=(double *)malloc(sizeof(double)*(k2+2));
if(r==NULL){
printf("Out of memory\n");
return -1;
}
tmp=(double *)malloc(sizeof(double)*(2*k2+4));
if(tmp==NULL){
printf("Out of memory\n");
free(r);
return -1;
}
n=1;
r[0]=625.0;
while(4*n<k2){
//Step 1. calcuate r*r;
for(i=n;i<2*n;i++)r[ i ]=0.0;
dzfft1dc(r,2*n,0,tmp);
dzfft1dc(r,2*n,1,tmp);//The forward FFT
r[0]*=r[0];
for(i=1;i<=n;i++){
double re=r[ i ];
double im=r[i+n+1];
double nr=re*re-im*im;
double ni=re*im*2.0;
r[ i ]=nr;
r[i+n+1]=ni;
}
zdfft1dc(r,2*n,0,tmp);
zdfft1dc(r,2*n,1,tmp);
for(i=0;i<2*n-1;i++){
longlong w=(longlong)(r[ i ]+0.5);
r[ i ]=(double)(w%MOD);
w/=MOD;
r[i+1]+=(double)w;
}
//Step 2. r*r-1.0
r[0]-=1.0;
//Step 3. (r*r-1.0)Mod 10^(2n)
for(i=2*n;i<4*n;i++)r[ i ]=0.0;
dzfft1dc(r,4*n,0,tmp);
dzfft1dc(r,4*n,1,tmp);//The forward FFT
r[0]*=r[0];
for(i=1;i<=2*n;i++){//normalize
double re=r[ i ];
double im=r[i+2*n+1];
double nr=re*re-im*im;
double ni=re*im*2.0;
r[ i ]=nr;
r[i+2*n+1]=ni;
}
zdfft1dc(r,4*n,0,tmp);
zdfft1dc(r,4*n,1,tmp);
for(i=0;i<2*n;i++){//normalize
longlong w=(longlong)(r[ i ]+0.5);
r[ i ]=(double)(w%MOD);
w/=MOD;
r[i+1]+=(double)w;
}//discard digits after 2*n
n*=2;
}
k1=k/4;
printf("Result for k=%d is\n",k);
if(k%4){//Output first several digits:
int limit=1;
longlong val=(longlong)(r[k1]+0.5);
for(i=0;i<k%4;i++)limit*=10;
val%=limit;
printf("%0*d",k%4,(int)val);
}
for(i=k1-1;i>=0;i--){
int val=(int)(r[ i ]+0.5);
printf("%04d",val);
if(i%10==0)printf("\n");
}
printf("\n");
free(r);
free(tmp);
}
在我的计算机上 (CPU 2.8G RedHat8.1, Memory 256M)，对于
n=10^6
时间：1.032s (（包含输入输出时间）
但是在n=10^7时候很慢，需要将近1分钟，这是因为这个程序消耗内存比较厉害，
在n=10^7时，加上动态联接库，需要的内存超过300M,所以就要使用虚拟内存了。
上面的算法是nlog(n)的。
在另外一台主频3G,内存2G的机器上
计算n^6是0.765s,n^7是16s. (包含输入输出时间）
不过n^8就没法计算了，要out of memory.
(x^2-1)^2 (mod 10^2k)
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mkl_fft.h>
#define MOD (10000)
double *r;
double *tmp;
#ifdef _WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif
int main()
{
int i,n;
int k,k1,k2,c2;
scanf("%d",&k);
if(k<=1)return -1;
k2=k,c2=0;
while(k2){
c2++;k2>>=1;
}
k2=1<<(c2-1); //k2<=k and k2 is power of 2.
if(k2<k)k2<<=1;
r=(double *)malloc(sizeof(double)*(k2+2));
if(r==NULL){
printf("Out of memory\n");
return -1;
}
tmp=(double *)malloc(sizeof(double)*(2*k2+4));
if(tmp==NULL){
printf("Out of memory\n");
free(r);
return -1;
}
n=1;
r[0]=625.0;
while(4*n<k2){
//Step 1. calcuate r*r;
for(i=n;i<2*n;i++)r[ i ]=0.0;
dzfft1dc(r,2*n,0,tmp);
dzfft1dc(r,2*n,1,tmp);//The forward FFT
r[0]*=r[0];
for(i=1;i<=n;i++){
double re=r[ i ];
double im=r[i+n+1];
double nr=re*re-im*im;
double ni=re*im*2.0;
r[ i ]=nr;
r[i+n+1]=ni;
}
zdfft1dc(r,2*n,0,tmp);
zdfft1dc(r,2*n,1,tmp);
for(i=0;i<2*n-1;i++){
longlong w=(longlong)(r[ i ]+0.5);
r[ i ]=(double)(w%MOD);
w/=MOD;
r[i+1]+=(double)w;
}
//Step 2. r*r-1.0
r[0]-=1.0;
//Step 3. (r*r-1.0)Mod 10^(2n)
for(i=2*n;i<4*n;i++)r[ i ]=0.0;
dzfft1dc(r,4*n,0,tmp);
dzfft1dc(r,4*n,1,tmp);//The forward FFT
r[0]*=r[0];
for(i=1;i<=2*n;i++){//normalize
double re=r[ i ];
double im=r[i+2*n+1];
double nr=re*re-im*im;
double ni=re*im*2.0;
r[ i ]=nr;
r[i+2*n+1]=ni;
}
zdfft1dc(r,4*n,0,tmp);
zdfft1dc(r,4*n,1,tmp);
for(i=0;i<2*n;i++){//normalize
longlong w=(longlong)(r[ i ]+0.5);
r[ i ]=(double)(w%MOD);
w/=MOD;
r[i+1]+=(double)w;
}//discard digits after 2*n
n*=2;
}
k1=k/4;
printf("Result for k=%d is\n",k);
if(k%4){//Output first several digits:
int limit=1;
longlong val=(longlong)(r[k1]+0.5);
for(i=0;i<k%4;i++)limit*=10;
val%=limit;
printf("%0*d",k%4,(int)val);
}
for(i=k1-1;i>=0;i--){
int val=(int)(r[ i ]+0.5);
printf("%04d",val);
if(i%10==0)printf("\n");
}
printf("\n");
free(r);
free(tmp);
}
在我的计算机上 (CPU 2.8G RedHat8.1, Memory 256M)，对于
n=10^6
时间：1.032s (（包含输入输出时间）
但是在n=10^7时候很慢，需要将近1分钟，这是因为这个程序消耗内存比较厉害，
在n=10^7时，加上动态联接库，需要的内存超过300M,所以就要使用虚拟内存了。
上面的算法是nlog(n)的。
在另外一台主频3G,内存2G的机器上
计算n^6是0.765s,n^7是16s. (包含输入输出时间）
不过n^8就没法计算了，要out of memory.
加载更多回复（24）
