梅森素数(Mersenne prime)判断, FFT 大数乘法 (非递归), O(n^2 l

原创代码，请勿转载！

梅森素数判定：
卢卡斯-莱默检验法：参考https://zh.wikipedia.org/wiki/%E5%8D%A2%E5%8D%A1%E6%96%AF-%E8%8E%B1%E9%BB%98%E6%A3%80%E9%AA%8C%E6%B3%95


卢卡斯-莱默检验法 是迭代算法，需要用到高精度乘法运算。
而现有的乘法运算算法最快的要数 快速傅里叶变换 算法，
可以将两位n-bit的大整数乘法在O(n lg n)计算出来。


FFT简介：快速傅里叶变换,参考 http://beige.ucs.indiana.edu/B673/node12.html


n = 8 为例子,计算两个4位数的乘法。
寻找根w使得 w^n = 1.
A[0..7][0..7]: A[i][j] = w^{i*j},其中 w^n = 1.
A^{-1} = w^{-1} ^ {i*j}.


对于大整数X和Y，首先做矩阵乘法运算
X' = A*X
Y' = A*Y
然后对X' Y' 对应bit位 对应相乘
Z = X' .* Y' (对应相乘)
最后还原
Z' = A^{-1}*Z


快速傅里叶变换的高效之处在于，A矩阵是具有很强的结构性，因此矩阵乘法运算过程可以优化到O(n lg n).


关于如何选择根w.
采用数论中的 root https://en.wikipedia.org/wiki/Root_of_unity_modulo_n
避免浮点数运算。?


对于n=2^k,寻找大素数MOD，使得 存在 w， w^{2^k}=1 (mod MOD).

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<ctime>
#include<cmath>

#define EXC(a,b) a ^ b ? a ^= b ^= a ^= b : 0
#define EXCB(a,b) a > b ? a ^= b ^= a ^= b : 0
#define FOR(i,a,b) for(int i=(a);i<(b);i++)
#define FORE(i,b) for(int i=(a);i<=(b);i++)
#define FORI(TYPE,it,s) for(TYPE::iterator it=s.begin();it!=s.end();it++) 
#define MST(a,b) memset(a,b,sizeof(a))
#define SCF(a) scanf("%d",&a)
#define SCFS(a) scanf("%s",a)
#define SCF2(a,b) scanf("%d%d",&a,&b)
#define SCF3(a,c) scanf("%d%d%d",&b,&c)
#define min2(a,b) ((a<b)?(a):(b))
#define max2(a,b) ((a>b)?(a):(b))
#define mk(a,b) make_pair(a,b)
using namespace std;
const int INF = 0x7FFFFFFF;
const double eps = 1e-8;
const int N = 1000000+5;//1<<27;
//#define SOUT "/dev/tty"
#define SOUT "CON"

typedef unsigned long long LL;
typedef bool BB;

int add(BB a[],int na,BB b[],int nb,BB c[],bool update);
int multi(BB a[],BB c[]);

//3221225473 = 2^30 * 3 + 1,11000000000000000000000000000001
const LL MOD = 3221225473LL;
// wa#Principal 2**k-th roots of unity: w^{2^(k-1)} = -1;
const LL wa[31]={ 1LL,3221225472LL,1013946479LL,1031213943LL,2526611335LL,66931328LL,651814490LL,1292405718LL,1958494276LL,764652596LL,1855261384LL,1168849724LL,211283056LL,1734477367LL,1445148504LL,1405588668LL,2519378301LL,1800384970LL,1358427440LL,511777184LL,3009749949LL,448192265LL,3182672916LL,140749519LL,538601489LL,1838654099LL,247425463LL,229807484LL,244140625LL,15625LL,125LL};
// wb#Inverses of principal 2**k-th roots of unity: wb = wa^{-1};
const LL wb[31]={ 1LL,2207278994LL,613406496LL,741264833LL,3098715598LL,704887665LL,764521865LL,2935727869LL,2450347685LL,532203874LL,2829850173LL,1588903488LL,919321183LL,1736059882LL,3000663231LL,2793181474LL,2702431710LL,615116962LL,2055537301LL,618423934LL,762996773LL,670776757LL,2925583959LL,1181743895LL,210645833LL,1294972776LL,2980752801LL,1869040087LL,121221157LL,2267742733LL};

#define IsComp(n)  (flag[n>>6]&(1<<((n>>1)&31)))
#define SetComp(n) flag[n>>6]|=(1<<((n>>1)&31))
#define MulAddMod(ans,p) ((ans  + a * b) % p)
#define MulAddModZ(a,p) ((a * b) % p)
#define bits(x) ceil(log((x)+0.5)/log(2.0))

using namespace std;

unsigned int flag[(N>>6)];
LL p,IW;			// prime p,2^IW = [p+p]
LL wps[N],wpr[N];	// w[i] = w^i;
LL arr[2][N];		// tmp array for calculation
BB D[N];			// tmp array for calculation
BB val[N],cont[N];	// *it:val = p*p-2,cont = -2;
// --- count prime ---
void count_Prime(int n){
	int i,j;
	MST(flag,0);
	for (i = 3; (j = i*i) < n; i += 2)
		if (!IsComp(i))
			for (; j < n; j += i<<1)
				SetComp(j);
}
// --- check prime ---
bool check_Prime(int x){
	if(x&1){
		return !IsComp(x);
	}
	else return (x==2);
}
// --- set w^i,n = 2^IW ---
void setw(int n,LL wt[],LL w){
	wt[0] = 1;
	FOR(i,1,n)	wt[i] = MulAddModZ(wt[i-1],w,MOD);
}
// --- set up enviroment for Mersene Prime test ---
void setup(int _p){
	p = _p;
	IW = bits(p<<1);
	setw(1<<IW,wps,wa[IW]);
	setw(1<<IW,wpr,wb[IW]);
	MST(val,0);
	MST(cont,0);
	//val = 4,(100),in the beginning.
	val[2] = 1;
	//cont = -2,(2222210),constant
	FOR(i,p) cont[i] = 1;
}
// --- big integer:add() ---
int add(BB a[],bool update){
//update = true,real calcuation in multi.
//update = false,handle -2 operation.
if(na<nb) return add(b,nb,na,c,update);
bool sum,tmp,flag = 0;
int n = na;//max2(na,nb);
FOR(i,nb){
	sum = a[i] ^ b[i] ^ flag;
	flag = (a[i] & b[i]) | (flag & (a[i] ^ b[i]));
	c[i] = sum; // in case c = a,or c = b
}
FOR(i,n){
	sum  = a[i] ^ flag;
	flag = a[i] & flag;
	c[i] = sum;
}
while(flag && update)
	FOR(i,n){
		sum = c[i] ^ flag;
		flag = c[i] & flag;
		c[i] = sum;
		if(!flag) break;
	}
while(n>0 && !c[n-1]) n--;
return n;
}
// --- DFT() X[] = from[] * matrix[][] ---
void DFT(int k,LL **from,LL **X,LL w[]){
	// k-bit,vector size is 2^k
	LL *a = *from,*b = *X;
	int d,n,mask;
	int u[3];//{i<<c,((i<<1)&mask)<<c,(((i<<1)|1)&mask)<<c}

	for(int c = k-1;c>=0;c--){
		d = 1<<(k-c);
		n = 1<<c;
		mask = (d-1)<<c;
		FOR(i,d){
			u[0] = i<<c;
			u[1] = (u[0]<<1)&(mask);
			u[2] = u[1] | n;
			FOR(j,n) b[u[0]|j] = MulAddMod(a[u[1]|j],a[u[2]|j],w[u[0]],MOD);//w^u[0]
		}
		swap(a,b);
	}
	if(!(k&1)){
		a = *from;
		*from = *X;
		*X = a;
	}
}
// --- big integer:multi() ---
int multi(BB a[],BB c[]){ //c = a*a
	int n = 1<<IW;// 2^IW >= p + p
	LL *A = arr[0];
	LL *C = arr[1];
	FOR(i,na) A[i] = a[i];
	FOR(i,n) A[i] = 0;

	DFT(IW,&A,&C,wps);
	FOR(i,n) A[i] = MulAddModZ(C[i],C[i],MOD);
	DFT(IW,wpr);

	int flag = 0;
	LL cur;
	FOR(i,n){
		cur = C[i]>>IW;
		c[i] = (cur+flag)&1;
		flag = (cur+flag)>>1;
	}
	while(n>p){
		MST(D,0);
		FOR(i,p,n){
			D[i-p] = c[i];
			c[i] = 0; //!important
		}
		n = add(c,D,n-p,true);
	}
	for(n = p;n>0 && !c[n-1];n--);
	return n;
}
// --- print array[],high bit first ---
void print(BB val[],int n){
LL ans = 0;
for(int i = n-1; i >= 0; i--){
	ans = (ans << 1) + val[i]; 
	printf("%d",val[i]);
}
printf(" -- %un",ans);
}
// --- check Mersene Prime ---
int check_Mersene(int x){
	if(!check_Prime(x)) return -1;
	int n = 3;//100
	setup(x);
	//print(val,n);
	FOR(i,2,p){
		n = multi(val,val);
		n = add(val,cont,val,false);
	//	print(val,n);
	}
	//printf("%d # %dn",x,n);
	return (!n);
}
// --- main ---
int main(int argc,char *argv[]){
	count_Prime(N-1);
	//check_Mersene(5);
	FOR(i,3,2000){
		long start = clock();
		int status = check_Mersene(i);
		if(status == -1){
			//printf("Not a Prime %d.",i);
			continue;
		}
		else if(status == 1){
			printf("tt++++Mersenne Prime p = %d.",i);
		}
		else{
			continue;
			//printf("%d.",i);
		}
		cout<< "t time:" << clock()-start<<endl;
	}
	return 0;
}
// -----------------------------------
//freopen("in","r",stdin);
//freopen("out2","w",stdout);

（编辑：邯郸站长网）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!