[BZOJ4635]数论小测验

Description
有一个长度为N的数组Ai，每个元素可以取1~M中的一个正整数。那么一共有M^N种可能的数组。因为 SHUXK 对数
论有特殊的爱好，所以他立刻想到了下面两个问题：

1. 对于给定的正整数K，有多少个数组Ai满足GCD(A1,A2...An) = K
2. 对于给定的正整数K，有多少个数组Ai满足K|Lcm(A1,A2...An)
   （我相信机智的你在看到这道题的英文名称时就已经猜得八九不离十了）当然，对于 SHUXK 来说，只询问一个K的
   情况很简单。于是他加强了一下：给定一个范围L~R，对L~R中所有的正整数K都求出答案，并且把它们加起来作为
   最终的结果。然后， SHUXK 就不会了……他需要你来帮助他计算出答案。 由于答案可能很大，我们只要求输出答
   案对1,000,000,007取模的结果。

Input
第一行包含两个正整数T和Type，分别表示测试数据的组数和数据类型
（ Type = 1时只询问 gcd 的问题， Type = 2时只询问 lcm 的问题）。
以下T行，每行包含四个正整数N,M,L,R（ 1 ≤ L ≤ R≤ M），含义如题面所述。
T<=500，Type=1时，M<=10^7，Type=2时M<=1000

Output
输出文件应包含T行，每行输出一个整数，表示一组测试数据的答案对
1,000,000,007取模的结果。答案的顺序应与输入数据的顺序保持一致。

Sample Input
5 1
5 5 1 5
5 5 2 5
5 5 3 5
5 5 4 5
5 5 5 5

Sample Output
3125
34
3
2
1

---

我们设(f[k]=sumlimits_{A_i=1}^ksumlimits_{A_2=1}^k...sumlimits_{A_n=1}^k[gcdlimits_{i=1}^n(A_i)=1])

反演一下得(f[k]=sumlimits_{x=1}^kmu(x)(lfloordfrac{k}{x} floor)^n)，那么我们能够用数论分块在(O(sqrt{k}))的时间内求出

则type=1时要求(sumlimits_{i=l}^rf[lfloordfrac{m}{i} floor])，这部分我们可以用分块在(O(sqrt{m}))的时间内求出

因此这部分时间复杂度为(O(m))

---

第二问考虑容斥，由于(k)很小，我们考虑将(k)分解质因数，lcm的意义是指数最大值，于是我们容斥，枚举哪些质因数的指数没有达到要求

于是问题转化为求出(m)以内有多少个数不满足以下约束，每一种约束形如一个质因子的指数不应该超过一个数

这个新的问题同样可以容斥，假定只有一个限制，那么答案显然为(m-dfrac{m}{p^a})，因为不符合条件的数只要包含了(p^a)，其余包含什么都行

考虑两个限制，那么答案即为(m-dfrac{m}{p_1^{a_1}}-dfrac{m}{p_2^{a_2}}+dfrac{m}{p_1^{a_1}p_2^{a_2}})，这个就是最基本的容斥方法，推广到多个约束亦然如此

得到(m)以内的方案数后(n)次方即可，因为每个位置都有这么多方案

一定记得预处理(n)次方快速幂，否则会TLE……

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
	static char buf[1000000],*p1=buf,*p2=buf;
	return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
inline int frd(){
	int x=0,f=1; char ch=gc();
	for (;ch<'0'||ch>'9';ch=gc())	if (ch=='-')	f=-1;
	for (;ch>='0'&&ch<='9';ch=gc())	x=(x<<3)+(x<<1)+ch-'0';
	return x*f;
}
inline int read(){
	int x=0,f=1; char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')	f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<3)+(x<<1)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x<0)	putchar('-'),x=-x;
	if (x>9)	print(x/10);
	putchar(x%10+'0');
}
const int N=1e7,M=1e3,p=1e9+7;
namespace Mobius{
	int prime[(N/10)+10],mu[N+10],sm[N+10];
	bool inprime[N+10];
	void prepare(){
		int tot=0; mu[1]=sm[1]=1;
		for (int i=2;i<=N;i++){
			if (!inprime[i])	mu[prime[++tot]=i]=-1;
			for (int j=1;j<=tot&&i*prime[j]<=N;j++){
				inprime[i*prime[j]]=1;
				if (i%prime[j]==0){
					mu[i*prime[j]]=0;
					break;
				}mu[i*prime[j]]=-mu[i];
			}
			sm[i]=sm[i-1]+mu[i];
		}
	}
	int mlt(int a,int b){
		int res=1;
		for (;b;b>>=1,a=1ll*a*a%p)	if (b&1)	res=1ll*res*a%p;
		return res;
	}
	int work(int m,int n){
		int res=0;
		for (int i=1,pos;i<=m;i=pos+1){
			pos=m/(m/i);
			res=(1ll*(sm[pos]-sm[i-1])*mlt(m/i,n)+res)%p;
		}
		return res;
	}
	void main(int T){
		prepare();
		for (;T;T--){
			int n=frd(),m=frd(),L=frd(),R=frd(),res=0;
			for (int i=L,pos;i<=R;i=pos+1){
				pos=min(m/(m/i),R);
				res=(1ll*(pos-i+1)*work(m/i,n)+res)%p;
			}
			printf("%d
",res);
		}
	}
};
const int P=1e9+7;
namespace IEP{//Inclusion-Exclusion Principle
	int prime[M+10],from[M+10];
	int A[10],mlu[M+10][10],g[M+10];
	bool inprime[M+10];
	int n,m,L,R,tot,res,Ans;//top;
	void prepare(){
		int tot=0;
		for (int i=2;i<=M;i++){
			if (!inprime[i])	from[prime[++tot]=i]=i;
			for (int j=1;j<=tot&&i*prime[j]<=M;j++){
				inprime[i*prime[j]]=1;
				if (i%prime[j]==0){
					from[i*prime[j]]=from[i];
					break;
				}from[i*prime[j]]=prime[j];
			}
		}
		prime[0]=tot;
	}
	int mlt(int a,int b){
		int res=1;
		for (;b;b>>=1,a=1ll*a*a%p)	if (b&1)	res=1ll*res*a%p;
		return res;
	}
	void solve(int k){
		int x=k,top=0;
		while (x!=1){
			int y=from[x],tmp=1;
			while (from[x]==y)	tmp*=y,x/=y;
			mlu[k][++top]=tmp;
		}
		mlu[k][0]=top;
	}
	void main(int T){
		prepare(); g[0]=1;
		for (int i=1;i<=M;i++)	solve(i);
		for (;T;T--){
			n=frd(),m=frd(),L=frd(),R=frd(),Ans=0;
			for (int i=1;i<=m;i++)	g[i]=mlt(i,n);
			for (int k=L;k<=R;k++){
				for (int sta=0;sta<1<<mlu[k][0];sta++){
					int num=0,tot=0,res=0;
					for (int i=0;i<mlu[k][0];i++)	if ((sta>>i)&1)	num++,A[tot++]=mlu[k][i+1];
					for (int s=0;s<1<<tot;s++){
						int cnt=0,tmp=m;
						for (int j=0;j<tot;j++)	if ((s>>j)&1)	cnt++,tmp/=A[j];
						res=((cnt&1?-1:1)*tmp+res)%p;
					}
					Ans=((num&1?-1:1)*g[res]+Ans)%p;
				}
			}
			printf("%d
",(Ans+p)%p);
		}
	}
};
int main(){
	int T=frd(),type=frd();
	if (type==1){
		Mobius::main(T);
		return 0;
	}
	if (type==2){
		IEP::main(T);
		return 0;
	}
	return 0;
}