[BZOJ3512] DZY Loves Math IV -- 杜教筛,搜索

#include"bits/stdc++.h"
#define last last_

using namespace std;
typedef long long LL;

const int P=1000000007;

struct node{
	int n,m;
	inline LL hash(){return (LL)n<<30|m;}
	node(int _n=0,int _m=0){n=_n,m=_m;}
};

const int N=5000005,M=500005,hmod=1000003,L=5000000,X=100005;

struct HashMap{
	node w[M];
	int val[M],flagnull,h[M],nxt[M],ndc;
	HashMap(){flagnull=-1,ndc=0;}
	int&find(node x){
		for(int i=h[x.hash()%hmod];i;i=nxt[i])
			if(x.n==w[i].n&&x.m==w[i].m)return val[i];
		return flagnull;
	}
	int&operator[](node x){
		int&r=find(x);
		if(r!=-1)return r;
		w[++ndc]=x,nxt[ndc]=h[x.hash()%hmod],h[x.hash()%hmod]=ndc;
		return val[ndc];
	}
}q;

int n,m,mu[X],eu[N],seu[N],ed[X],pri[N>>3],pcnt,ans;
bool fl[N];
vector<int>di[X];

int G(int n){
	if(n<=L)return seu[n];
	node x(1,n);
	int r=q.find(x);
	if(r!=-1)return r;
	r=((LL)n*(n+1)>>1)%P;
	for(int i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		if(n/i==1)continue;
		r=(r-(LL)(n/i-1)*(G(j)-G(i-1)))%P;
	}
	return q[x]=(r+P)%P;
}

int F(int n,int m){
	if(n==1)return G(m);
	if(m==0)return 0;
	if(mu[n]==0)return (LL)F(ed[n],m)*(n/ed[n])%P;
	node x(n,m);
	int r=q.find(x);
	if(r!=-1)return r;
	r=0;
	for(int i=0,d;i<di[n].size();i++){
		d=di[n][i];
		r=(r+(LL)eu[n/d]*F(d,m/d))%P;
	}
	return q[x]=r;
}

int main(){
#ifndef ONLINE_JUDGE
	freopen("bzoj3512.in","r",stdin);
	freopen("bzoj3512.out","w",stdout);
#endif
	scanf("%d%d",&n,&m);
	mu[1]=ed[1]=eu[1]=seu[1]=1;
	for(int i=2;i<=L;i++){
		if(i<=n)di[i].push_back(1);
		if(!fl[i]){
			pri[++pcnt]=i;
			if(i<=n){
				ed[i]=i;
				mu[i]=-1;
			}
			eu[i]=i-1;
		}
		for(int j=1;i*pri[j]<=L;j++){
			fl[i*pri[j]]=1;
			if(i%pri[j]==0){
				if(i*pri[j]<=n){
					mu[i*pri[j]]=0;
					ed[i*pri[j]]=ed[i];
				}
				eu[i*pri[j]]=eu[i]*pri[j];
				break;
			}
			if(i*pri[j]<=n){
				mu[i*pri[j]]=-mu[i];
				ed[i*pri[j]]=ed[i]*pri[j];
			}
			eu[i*pri[j]]=eu[i]*(pri[j]-1);
		}
		for(int j=1;i*j<=n;j++)di[i*j].push_back(i);
		seu[i]=(seu[i-1]+eu[i])%P;
	}
	for(int i=1;i<=n;i++)ans=(ans+F(i,m))%P;
	cout<<ans<<endl;
	return 0;
}