【自守数】(在十进位制中,) 若一个 k 位正整数 N (可含前置 "0" ), 若满足如下性质: 任意两个或多个均以该字串 N 结尾的整数相乘, 其结果的最后 k 位数字一定还是 N, 那么, 则称 N 为 "k 位自守数".
我们有结论:
如果x是k位自守数，那么
(x^2-1)^2 mod 10^(2k) 是2k位自守数。
关于整数的乘法，
假设N位L进制数
A= a0+a1*L+a2*L^2+...+a(N-1)*L^(N-1)
B=b0+b1*L+b2*L^2+...+b(N-1)*L^(N-1)
C=A*B=(a0*b0)+(a0*b1+a1*b0)*L+...+a(N-1)*a(N-1)*L^(2*N-2)
我们可以看到C的每个系数同卷积非常象
如果我们在A同B后面都添加N项0,
也就是A=a0+a1*L+...+a(2N-1)*L^(2N-1)
B=b0+b1*L+...+b(2N-1)*L^(2N-1)
其中a(N),a(N+1),...,a(2N-1);b(N),b(N+1),...,b(2N-1)都是0
那么
C=(a0*b0+a1*b(2N-1)+...+a(2N-1)*b1)+
(a0*b1+a1*b0+a2*b(2N-1)+...+a(2N-1)*b2)*L+
...+
(a0*b(2N-1)+a1*b(2N-2)+...+a(2N-1)*b(0))*L^(2N-1)
其中各个系数正好是A和B的卷积，只是计算结果有可能每项超出[0,L)之间的范围要求，需要最后做调整。
而卷积可以通过离散傅立叶变换来计算，如果计算两个数的平方，那么还可以少一次傅立叶变换。
下面是计算自守数的代码，使用了Intel MKL中傅立叶变换函数:
#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);
}
发表于 @ 2007年10月31日 10:10:00|评论(loading...)|编辑