定義(a$b)i = Σaj * bi^j ，那麼題目相當於求A^N的第0項，其中Ai=[i為質數且i<=m]。

       顯然我們可以用快速冪+分治乘做到N^1.59logN，光榮TLE。

       UPD：後來發現這是fwt是Nlog^2N的好像能過。

       和普通多項式乘法一樣，我們尋找一種變換trans(a)，使得trans(a$b)=trans(a)*trans(b)，這裡有(X*Y)i=Xi*Yi。只要能O(NlogN)轉化和復原就能解決這個問題，先來看a,b含有兩個數的情況。

      a(x,y)，b(z,w)。trans(a)=(x-y,x+y)，trans(b)=(z-w,z+w)，那麼：

      trans(a$b)=trans((xz+yw,xw+yz))=(xz+yw-xw-yz,xz+yw+xw+yz)，trans(a)*trans(b)=((x-y)(z-w),(x+y)(z+w))，顯然兩者相同。

      對於一般情況，定義a=(a1,a2)，表示ai=a1i(i<n/2)或者a2(i-n/2)(i>n/2)，有trans(a)=(trans(a1)-trans(a2),trans(a1)+trans(a2))。這個和上面一樣用a,b展開即可證明。

      復原就是每一步都反著來。然後遞迴實現即可。不過這道題目非遞迴還是很簡單的應該會快不少。

      UPD：本質就是fwt。

AC程式碼如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
#define mod 1000000007
#define inv 500000004
#define N 70005
using namespace std;

int n,m,cnt,a[N],c[N]; bool vis[N];
int ksm(int x,int y){
	int t=1; for (; y; y>>=1,x=(ll)x*x%mod) if (y&1) t=(ll)t*x%mod;
	return t;
}
void trs(int l,int r){
	if (l==r) return; int mid=(l+r)>>1,i,j,x;
	trs(l,mid); trs(mid+1,r);
	for (i=l,j=mid+1; i<=mid; i++,j++){
		x=(a[i]+a[j])%mod;
		a[i]=(a[i]-a[j]+mod)%mod; a[j]=x;
	}
}
void rsto(int l,int r){
	if (l==r) return; int mid=(l+r)>>1,i,j,x;
	for (i=l,j=mid+1; i<=mid; i++,j++){
		x=(ll)(a[j]-a[i]+mod)*inv%mod;
		a[i]=(ll)(a[i]+a[j])*inv%mod; a[j]=x;
	}
	rsto(l,mid); rsto(mid+1,r);
}
int main(){
	int i,j;
	for (i=2; i<=50000; i++){
		if (!vis[i]) c[++cnt]=i;
		for (j=1; j<=cnt && i*c[j]<=50000; j++){
			vis[i*c[j]]=1;
			if (!(i%c[j])) break;
		}
	}
	c[cnt+1]=50001;
	while (~scanf("%d%d",&n,&m)){
		memset(a,0,sizeof(a));
		for (i=1; c[i]<=m; i++) a[c[i]]=1;
		j=1; while (j<=m) j<<=1; j--;
		trs(0,j);
		for (i=0; i<=j; i++) a[i]=ksm(a[i],n);
		rsto(0,j);
		printf("%d\n",a[0]);
	}
	return 0;
}
 

by lych

2016.5.17