相关链接
题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4635
数据生成器：http://paste.ubuntu.com/24360378/
神犇题解：http://www.cnblogs.com/clrs97/p/5625508.html
解题报告
对于第一问，显然可以用莫比乌斯反演做到$O(Tm)$
对于第二问，考虑枚举$k$，然后$10^3$以内的数最多有$4$种不同的质因数
于是搞一个状压$DP$，用矩阵快速幂优化
单词询问时间复杂度：$O(32m^2 + 32^3 log (n))$
看起来蛮有希望过的，但卡了一下午常也没有卡过 QwQ
正解的话，我们可以学习Claris用容斥
具体来讲枚举$k$，然后枚举$k$的每一个质因数是否满足条件
然后配合一点预处理之类的，就可以做到$O(m^2 + m \log m)$了
Code
这份代码在BZOJ被卡常了 QwQ
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int MOD = 1000000007;
const int N = 10000009;
const int M = 32;
int n,m,SZ,s1[N],hc[N],to[1001];
int tot,pri[N>>3],mu[N]; bool vis[N];
inline int read() {
char c=getchar(); int f=1,ret=0;
while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
return ret * f;
}
inline int Pow(int w, int t) {
int ret = 1;
for (;t;t>>=1,w=(LL)w*w%MOD)
if (t&1) ret=(LL)ret*w%MOD;
return ret;
}
inline void GetMu() {
mu[1] = 1;
for (int i=2;i<N;i++) {
if (!vis[i]) mu[i] = -1, pri[++tot] = i;
for (int j=1;j<=tot&&pri[j]*i<N;j++) {
vis[i*pri[j]] = 1;
if (i%pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
mu[i] = mu[i-1] + mu[i];
}
}
inline int cal(int mx) {
int ret = 0;
for (int l=1,r,tmp;l<=mx;l=r+1) {
r = mx / (mx / l);
tmp = (LL)(mu[r] - mu[l-1]) * hc[mx/l] % MOD;
ret = (ret + tmp) % MOD;
}
return (ret + MOD) % MOD;
}
struct Matrix{
int a[M][M];
inline Matrix() {memset(a,0,sizeof(a));}
inline Matrix(const Matrix *A) {
for (int i=0;i<M;i++) for (int j=0;j<M;j++) a[i][j] = A->a[i][j];
}
inline Matrix(int x) {
memset(a,0,sizeof(a));
for (int i=0;i<SZ;i++) a[i][i] = x;
}
inline void operator *= (const Matrix &A) {
Matrix ret;
for (int i=0,*t1,*t3;i<SZ;i++) {
t1 = ret.a[i]; const int *t2 = A.a[i];
for (int j=0;j<SZ;j++) {
t3 = a[j];
for (int k=0;k<SZ;k+=4) {
t1[k] = (t1[k] + (LL)t3[k] * t2[j]) % MOD;
t1[k+1] = (t1[k+1] + (LL)t3[k+1] * t2[j]) % MOD;
t1[k+2] = (t1[k+2] + (LL)t3[k+2] * t2[j]) % MOD;
t1[k+3] = (t1[k+3] + (LL)t3[k+3] * t2[j]) % MOD;
}
}
}
for (int i=0;i<SZ;i++) for (int j=0;j<SZ;j++) a[i][j] = ret.a[i][j];
}
inline Matrix operator ^ (int t) {
Matrix ret(1),w(this);
for (;t;t>>=1,w*=w) if (t&1) ret*=w;
return ret;
}
}ans,trans;
inline int solve(int w) {
tot = 0;
for (int i=2,tmp=w;i<=tmp;i++) {
if (tmp % i == 0) {
pri[++tot] = i; tmp /= i;
while (tmp % i == 0) pri[tot] *= i, tmp /= i;
}
}
ans = Matrix(); trans = Matrix();
ans.a[0][0] = 1; SZ = 1 << tot;
memset(to,0,sizeof(to));
for (int i=1,t=1,ww;ww=pri[i],i<=tot;i++,t<<=1)
for (int j=ww;j<=m;j+=ww) to[j] |= t;
for (int i=1,tt;tt=to[i],i<=m;i++) {
for (int p=0,ww;ww=p,p<SZ;p++)
++trans.a[p|tt][p];
}
trans = trans ^ n;
ans *= trans;
return ans.a[SZ-1][0];
}
int main() {
int T = read(), type = read();
if (type == 1) {
GetMu();
for (int t=1,vout;vout=0,t<=T;t++) {
n = read(); m = read();
int L = read(), R = read();
for (int l=1,r;l<=m;l=r+1) {
r = m / (m / l);
hc[m / l] = Pow(m / l, n);
}
for (int l=L,r,tmp;l<=R;l=r+1) {
r = m / (m / l); tmp = cal(m / l);
vout = (vout + (min(r,R)-max(L,l)+1ll) * tmp) % MOD;
}
printf("%d\n",vout);
}
} else if (type == 2) {
for (int t=1,vout,l,r;vout=0,t<=T;t++) {
n = read(); m = read(); l = read(), r = read();
for (int w=l;w<=r;w++) vout = (vout + solve(w)) % MOD;
printf("%d\n",vout);
}
}
return 0;
}
