数论习题（二）

（二）.筛子

- P4213 【模板】杜教筛

给定 $1\le n<2^{31}$，求下面两个东西：

$$S_{\phi }(n)=\sum _{i=1}^n\phi (i),S_{\mu }(n)=\sum _{i=1}^n\phi (i)$$

显然这个东西不能直接 $O(n)$ 求，那么我们就用杜教筛！

1.筛 $\phi$

我们知道 $\phi \ast 1=id$，于是我们取 $g=1$，有：

$$S_{\phi }(n)=\sum _{i=1}^n(\phi \ast 1)(i)-\sum _{i=2}^{n}1(i)\cdot S_{\phi }(\lfloor \frac{n}{i}\rfloor )$$

也就是：

$$\sum _{i=1}^{n}i-\sum _{i=2}^n{S}_{\phi }(\lfloor \frac{n}{i}\rfloor )=\frac{n(n+1)}{2}-\sum _{i=2}^n{S}_{\phi }(\lfloor \frac{n}{i}\rfloor )$$

就OK了，直接筛就行。

2.筛 $\mu$

我们知道 $\mu \ast 1=\epsilon$，于是我们取 $g=1$，有：

$$S_{\mu }(n)=\sum _{i=1}^n(\mu \ast 1)(i)-\sum _{i=2}^{n}1(i)\cdot S_{\mu }(\lfloor \frac{n}{i}\rfloor )$$

也就是：

$$\sum _{i=1}^n\epsilon (i)-\sum _{i=2}^n{S}_{\mu }(\lfloor \frac{n}{i}\rfloor )=1-\sum _{i=2}^n{S}_{\mu }(\lfloor \frac{n}{i}\rfloor )$$

就OK了，直接筛就行。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1e6+7e5+10;
ll T=1;
ll n;
ll prime[N],v[N],mu[N],phi[N],tot;
ll smu[N],sphi[N];
map <ll,ll> mp_mu,mp_phi;
void shai(ll n){
	phi[1]=1,mu[1]=1,v[1]=1;
	for(int i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i;prime[++tot]=i;mu[i]=-1;phi[i]=i-1;
		}
		for(int j=1;j<=tot;j++){
			if(prime[j]>v[i] || prime[j]>n/i) break;
			v[i*prime[j]] = prime[j];
			if(i%prime[j]==0) mu[i*prime[j]]=0,phi[i*prime[j]]=phi[i] * prime[j];
			else mu[i*prime[j]] = -mu[i],phi[i*prime[j]] = phi[i] * (prime[j]-1);
		}
	}
	for(int i=1;i<=n;i++){
		smu[i]=smu[i-1] + mu[i],sphi[i]=sphi[i-1]+phi[i];
	}
}
ll sum_mu(ll n){
	if(n<=N-10) return smu[n];
	if(mp_mu.find(n)!=mp_mu.end()) return mp_mu[n];
	ll res=1;
	for(ll l=2,r;l<=n;l=r+1){
		r= n/(n/l);
		res-=(r-l+1) * sum_mu(n/l);
	}return mp_mu[n]=res;
}
ll sum_phi(ll n){
	if(n<=N-10) return sphi[n];
	if(mp_phi.find(n)!=mp_phi.end()) return mp_phi[n];
	ll res=n*(n+1)/2;
	for(ll l=2,r;l<=n;l=r+1){
		r= n/(n/l);
		res-=(r-l+1) * sum_phi(n/l);
	}return mp_phi[n]=res;
}
int main(){
	shai(N-10);
	scanf("%lld",&T);
	for(int Case=1;Case<=T;Case++){
		scanf("%lld",&n);
		printf("%lld %lld\n",sum_phi(n),sum_mu(n));
	}
	return 0;
}

---

-BZOJ4176. Lucas的数论

给定 $1\le n\le {10}^9$，求：

$$\sum _{i=1}^n\sum _{j=1}^{n}d(ij)$$

这玩意看着和 P3327 [SDOI2015] 约数个数和 很像，故我们直接照搬结论：

$$\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=1]\lfloor \frac{n}{i}\rfloor \lfloor \frac{n}{j}\rfloor$$

莫反拆开：

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{d\mid gcd(i,j)}\mu (d)\lfloor \frac{n}{i}\rfloor \lfloor \frac{n}{j}\rfloor$$

枚举 $d$：

$$\sum _{d=1}^n\mu (d)(\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{n}{id}\rfloor )$$

然后两层都数论分块，$\mu$ 的前缀和用杜教筛筛除来就行了。

复杂度不会推，但时限 3s，能过。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=5e6+10,mod=1000000007;
ll n;
ll prime[N],tot,v[N],mu[N];
map <ll,ll> mp;
void shai(ll n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(v[i]==0){
			v[i]=i,prime[++tot]=i;mu[i]=-1;
		}
		for(int j=1;j<=tot;j++){
			if(prime[j] > v[i] || prime[j] > n/i) break;
			v[i*prime[j]] = prime[j];
			if(i%prime[j]==0) mu[i*prime[j]]=0;
			else mu[i*prime[j]] = -mu[i];
		}
	}
	for(int i=1;i<=n;i++) (mu[i] += mu[i-1])%=mod;
}
ll sum_mu(ll n){
	if(n<=N-10) return mu[n];
	if(mp.find(n)!=mp.end()) return mp[n];
	ll res=1;
	for(ll l=2,r;l<=n;l=r+1){
		r= n/(n/l);
		res-=(r-l+1) * sum_mu(n/l);
	}return mp[n]=res;
}
int main(){
	shai(N-10);
	scanf("%lld",&n);
	ll ans=0;
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ll cnt=0,m=n/l;
		for(int lj=1,rj;lj<=m;lj=rj+1){
			rj=m/(m/lj);
			(cnt+=(rj-lj+1) * (m/lj) %mod) %=mod;
		}
		(ans += cnt*cnt%mod*(sum_mu(r)-sum_mu(l-1))%mod)%=mod;
	}
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}

---

-BZOJ3512. DZY Loves Math IV

给定 $1\le n\le {10}^5$，$1\le m\le {10}^9$，求：

$$\sum _{i=1}^n\sum _{i=1}^m\phi (ij)$$

首先我们知道：

$$\phi (ij)=\frac{\phi (i)\phi (j)gcd(i,j)}{\phi (gcd(i,j))}$$

很可惜，这个题没法用这个东西推（我试了，会推出来一个很长且没法算的东西），我们要用这个：

若 $p\in \mathbb{P}$，则：

$$\phi (n)=\{\begin{array}{ll}(p-1)\phi (\frac{n}{p})& p\mid n\text{ 且 }{p}^2\nmid n\\ p\phi (\frac{n}{p})& p^2\mid n\end{array}$$

我们重新看这个式子，发现 $n$ 这一维很小，所以可以考虑先枚举 $i\in [1,n]$，分别求和。

我们记：

$$f(a,n)=\sum _{i=1}^n\phi (ai)$$

我们随便找一个 $p\in \mathbb{P}$ 且 $p\mid a$，尝试把 $a$ 拆成 ${\displaystyle \frac{a}{p}}\times a$，分情况讨论：

1.$p^2\nmid a$

此时我们无法保证是否存在 $p^2\nmid ai$，即有一部分 $p\mid i$ 的情况下，有 $p^2\mid ai$。

于是原式变为：

$$\sum _{i=1}^n[p\nmid i]\phi (p\cdot \frac{ai}{p})+\sum _{i=1}^n[p\mid i]\phi (p\cdot \frac{ai}{p})$$

带一下上面的式子，也就是：

$$\sum _{i=1}^n[p\nmid i](p-1)\cdot \phi (\frac{ai}{p})+\sum _{i=1}^n[p\mid i]p\cdot \phi (\frac{ai}{p})$$

我们可以把后半段的 $p-1$ 强行 -1，然后再加上：

$$\sum _{i=1}^n(p-1)\cdot \phi (\frac{ai}{p})+\sum _{i=1}^n[p\mid i]\phi (\frac{ai}{p})$$

也就是：

$$(p-1)\cdot \sum _{i=1}^n\phi (\frac{ai}{p})+\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\phi (ai)$$

依据 $f$ 的定义，也就是：

$$(p-1)\cdot f(\frac{a}{p},n)+f(a,\lfloor \frac{n}{p}\rfloor )$$

2.$p^2\mid a$

我们此时能保证 $p^2\mid ai$，于是可以直接套公式：

$$\sum _{i=1}^{n}p\cdot \phi (\frac{ai}{p})$$

也就是：

$$p\cdot f(\frac{a}{p},n)$$

---

于是我们把 $f$ 写成记忆化搜索，一直搜就行了

边界条件：

1.$a=1$

此时原式就是 $\sum _{i=1}^n\phi (i)$，跑杜教筛即可；

而这个 $m 又是之前的 $m$ 除以一个数下取整得来的，和杜教筛内部的数论分块分界点相同，所以杜教筛只会跑一遍。

2.$n=0$

这时候直接返回 $0$ 就行了，我因为没加这个而 MLE 了好久。

---

我们这个题要求的是：$\sum _{i=1}^{n}f(i,m)$，枚举就行了。

至于为什么复杂度是对的，我也不会证awa

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef int ll;
const int N=1e6+10,mod=1e9+7,inv2=500000004;
ll n,m;
ll tot,v[N],phi[N],maxp[100010];
unordered_map <ll,ll> mpphi,mpf[100010];
vector <ll> prime;
void shai(ll n){
	phi[1]=1,v[1]=1;
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,prime.push_back(i);tot++;
			phi[i]=i-1;
		}
		for(ll j=0;j<tot;j++){
			if(prime[j] > v[i] || prime[j] > n/i) break;
			v[i*prime[j]] = prime[j];
			if(i%prime[j]==0) phi[i*prime[j]]=prime[j]*phi[i];
			else phi[i*prime[j]]=(prime[j]-1)*phi[i];
		}
	}
	for(int i=1;i<=100000;i++) maxp[i]=i;
	for(ll i=1;i<=n;i++) (phi[i]+=phi[i-1])%=mod;
	for(ll i=0;i<tot;i++){
		for(ll j=prime[i];j<=100000;j+=prime[i]){
			while(maxp[j] % prime[i]==0 && maxp[j]!=prime[i]) maxp[j]/=prime[i];
		}
	}
}
ll S_phi(ll n){
	if(n<=N-10) return phi[n];
	if(mpphi.find(n)!=mpphi.end()) return mpphi[n];
	ll res=1ll*n*(n+1)%mod*inv2%mod;
	for(ll l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		(res-=1ll*(r-l+1)*S_phi(n/l)%mod)%=mod;
	}
	return mpphi[n]=(res+mod)%mod;
}
ll f(ll a,ll m){
	if(mpf[a].find(m)!=mpf[a].end()) return mpf[a][m];
	if(!m) return 0;
	if(a==1){
		return S_phi(m);
	}
	ll p=maxp[a],pa=a/p;
	if(pa % p){
		mpf[a][m]=(1ll*(p-1) * f(pa,m) %mod + f(a,m/p) ) %mod;
	}else{
		mpf[a][m]=(1ll*p * f(pa,m)) % mod;
	}
	return mpf[a][m];
}
int main(){
	shai(N-10);
	scanf("%d%d",&n,&m); 
	ll ans=0;
	for(ll i=1;i<=n;i++){
		(ans+=f(i,m))%=mod;
	}
	printf("%d\n",(ans+mod)%mod);
	return 0;
}

---

-P3768 简单的数学题

给定 $n$，求：

$$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$$

直接上套路推式子，枚举 $gcd$：

$$\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)=1]ij{d}^2$$

莫反拆开，$d^2$ 提出去：

$$\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{k\mid gcd(i,j)}\mu (k)ij$$

枚举 $k$：

$$\sum _{d=1}^n{d}^3\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)k^2\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{dk}\rfloor }ij$$

$i,j$ 上界有 $dk$，于是老套路令 $T=dk$，枚举 $T$：

$$\sum _{T=1}^n(\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }{)}^2\sum _{d\mid T}{d}^3\mu (\frac{T}{d})\cdot (\frac{T}{d}{)}^2$$

后面那个关于 $d$ 的式子，把 $T^2$ 提出来，就是：

$$T^2\sum _{d\mid T}d\cdot \mu (\frac{T}{d})$$

其中有 $id\ast \mu$，卷出来就是 $\phi$

所以原式：

$$\sum _{T=1}^n(\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }{)}^2{T}^2\phi (T)$$

典中典数论分块，我们只需要尝试求出 $f(T)=T^2\phi (T)$ 的前缀和就行了，杜教筛嘛。

我们令 $g=i{d}^2$，则：

$$(f\ast g)(T)=\sum _{d\mid T}{d}^2\phi (d)\cdot (\frac{T}{d}{)}^2=T^2\sum _{d\mid T}\phi (d)=T^3=i{d}^3(T)$$

也就是 $f\ast g=i{d}^3$，即 $i{d}^2\ast f=i{d}^3$。

直接跑杜教筛就行了。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1e7+10;
ll inv2,inv6,mod,n;
ll prime[N],tot,v[N],phi[N];
ll sf[N];
map <ll,ll> mp;
void shai(ll n){
	phi[1]=1;v[1]=1;
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,prime[++tot]=i;
			phi[i]=i-1;
		}for(ll j=1;j<=tot;j++){
			if(prime[j]>v[i]||prime[j]>n/i) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0) phi[i*prime[j]]=prime[j] *phi[i];
			else phi[i*prime[j]]=(prime[j]-1)*phi[i];
		}
	}
	for(ll i=1;i<=n;i++) sf[i]=(sf[i-1] + i%mod*i%mod*phi[i]%mod)%mod;
}
ll qpow(ll a,ll n){
	ll res=1;
	while(n){
		if(n&1) res=(res*a)%mod;
		a=(a*a)%mod,n>>=1;
	}return res;
}
ll S_2(ll n){
	n%=mod;
	return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
ll S(ll n){
	if(n<=N-10) return sf[n];
	if(mp.find(n)!=mp.end())return mp[n];
	ll res=n%mod*((n+1)%mod)%mod*inv2%mod;
	res=(res*res)%mod;
	for(ll l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		(res-=(S_2(r)-S_2(l-1))%mod*S(n/l)%mod)%=mod;
	}
	return mp[n]=(res+mod)%mod;
}
int main(){
	scanf("%lld%lld",&mod,&n);
	inv2=qpow(2,mod-2)%mod;
	inv6=qpow(6,mod-2)%mod;
	shai(N-10);
	ll ans=0;
	for(ll l=1,r;l<=n;l=r+1){
		ll nl=n/l;r=n/(n/l);
		ll tmp=nl%mod*((nl+1)%mod)%mod*inv2%mod,ttt=(S(r)-S(l-1))%mod;
		tmp=(tmp*tmp)%mod,ttt=(ttt+mod)%mod;
		(ans+=tmp* ttt %mod)%=mod;
	}printf("%lld\n",(ans%mod+mod)%mod);
	return 0;
}

---

- P5325 【模板】Min_25 筛

min_25 筛即可。

我们令 $f_1^{\prime }(x)=x,f_2^{\prime }(x)=x^2$，对俩函数筛一遍，算就行了。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=2e5+10,mod=1e9+7;
const int inv2=500000004,inv6=166666668;
ll n,qn;
ll p[N],tot,v[N];
ll s1[N],s2[N];//质数前缀和&质数平方前缀和 
ll m,id1[N],id2[N],w[N];//离散化 
ll g1[N],g2[N];
void shai(ll n){
	for(int i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,p[++tot]=i;
		}for(int j=1;j<=tot;j++){
			if(p[j] > v[i] || p[j] > n/i) break;
			v[i*p[j]]=p[j];
		}
	}
	for(int i=1;i<=tot;i++){
//		(s1[i] += p[i]) %=mod;
//		(s2[i] += p[i] * p[i] % mod) %=mod;
		s1[i] = (s1[i-1] + p[i]) %mod;
		s2[i] = (s2[i-1] + p[i] * p[i] % mod) %mod;
	}return;
}
ll id(ll x){
	if(x <= qn) return id1[x];
	else return id2[n/x];
}
ll S(ll n,ll j){
	if(p[j] > n) return 0;
	ll res = (g2[id(n)] - g1[id(n)])%mod - (s2[j] - s1[j])%mod;
	res=(res%mod+mod)%mod;
	for(int k=j+1;k<=tot && p[k] * p[k] <= n;k++){
		for(ll e=1,now=p[k];now<=n;now*=p[k],e++){
			(res += now%mod * (now%mod-1) %mod * (S(n/now,k) + (e>1)) %mod)%=mod;
		}
	}return res;
}
int main(){
	scanf("%lld",&n);qn=sqrt(n);shai(qn);
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l),w[++m]=n/l;
//		printf("l:%lld r:%lld\n",l,r);
		g1[m] = w[m]%mod * (w[m]%mod + 1) %mod * inv2 %mod -1;
		g2[m] = w[m]%mod * (w[m]%mod + 1) %mod * (2 * w[m]%mod + 1) %mod *inv6 %mod -1;
		if(w[m] <= qn) id1[w[m]]=m;
		else id2[n/w[m]] = m;
	}//离散化 
	//这里通过滚动数组优化。 
	for(ll j=1;j<=tot;j++){//j是阶段，我们先枚举j，当然直接枚到tot就行，因为筛的质数最大到根号n 
		for(ll i=1;i<=m && p[j] * p[j] <= w[i];i++){
			//根据转移术式，当pj*pj>wi时，根本就不用变！ 
			(g1[i] -= p[j] *              (g1[id(w[i]/p[j])] - s1[j-1] )%mod )%=mod;
			(g2[i] -= p[j] * p[j] % mod * (g2[id(w[i]/p[j])] - s2[j-1] )%mod )%=mod;
			(g1[i] += mod )%=mod;(g2[i] += mod) %=mod;
		}
	}
	ll ans=S(n,0) + 1;
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}
//10000000000

---

- loj6053. 简单的函数

给定积性函数 $f$，满足：

$$\begin{array}{rl}& f(1)=1\\ & f(p^c)=p\mathrm{xor}c\end{array}$$

给定 $1\le n\le$，求 $\sum _{i=1}^{n}f(i)$。

我们使用 min_25 筛。

对于质数 $p$，显然有：

$$f(p)=p\mathrm{xor}1=p-1+[p=2]2$$

可以先忽略 $[p=2]2$ 部分，对 $f_1^{\prime }(p)=1$ 和 $f_2^{\prime }(p)=p$ 的 $g$ 筛出来，然后拟合 $v$ 函数，做就行了。

（函数名参见此博客）

我们发现，$f(2)$ 处的值我们少算了 $2$，于是在最后把答案额外加二就行了。

但我突然好奇为什么？

比如一个例子：令 $n=2$。

因为假如说 $f(2)=1,f(3)=2$，那么就会算出错误的 $f(6)=1\times 2=2$，这显然是错误的，即为什么 $f(2)$ 少算的值不会影响其它的计算？

我们把 $S$ 函数的递推式和求解 $S$ 的函数代码搬过来：

$$S(n,j)=q(n)-\sum _{i=1}^{j}f(P_i)+\sum _{k>j}^{P_k^2\le n}\sum _{e=1}^{P_k^e\le n}f(P_k^e)\times (S(\frac{n}{P_k^e},k)+[e>1])$$

ll S(ll n,ll j){
	if(p[j] > n) return 0;
	ll res = (g2[id(n)] - g1[id(n)])%mod - (s2[j] - s1[j])%mod;
	res=(res%mod+mod)%mod;
	for(int k=j+1;k<=tot && p[k] * p[k] <= n;k++){
		for(ll e=1,now=p[k];now<=n;now*=p[k],e++){
			(res += (p[k] ^ e) %mod * (S(n/now,k) + (e>1)) %mod)%=mod;
		}
	}
	return res;
}

我们有 $q(n)=\sum _{i=1}^{n}f(i)[i\in \mathbb{P}]$。

注意到，我们只在第一次调用 $S(n,0)$ 时，$\sum _{i=1}^{j}f(P_i)$ 中不包含错误的 $f(2)$，而 $q(n)$ 中始终包含 $f(2)$。

这意味着，在我们调用 $S(n,j)$，其中 $j>0$ 时，$q$ 与 $\sum _{i=1}^{j}f(P_i)$ 中的 $f(2)$ 会抵消掉，且后面的循环从 $k=j+1$ 开始枚举，不会产生额外的影响。

在 $j=0$ 时，我们 $k$ 会从 $1$ 开始枚举，而 $e$ 也会从 $1$ 开始枚举，注意到代码里写的是 p[k]^e，即 $P_k\mathrm{xor}1$，即 $2\mathrm{xor}1$，并没有任何错误的 $f(2)$ 带来影响！！！

唯一出现影响的 $f(2)$，只在计算 (g2[id(n)] - g1[id(n)]) 中出现了一次！

所以我们最后只需要加上 $2$ 就能得到正确的答案了。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=2e5+10,mod=1e9+7;
const int inv2=500000004,inv6=166666668;
ll n,qn;
ll p[N],tot,v[N];
ll s1[N],s2[N];//质数前缀和&质数平方前缀和 
ll m,id1[N],id2[N],w[N];//离散化 
ll g1[N],g2[N];
void shai(ll n){
	for(int i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,p[++tot]=i;
		}for(int j=1;j<=tot;j++){
			if(p[j] > v[i] || p[j] > n/i) break;
			v[i*p[j]]=p[j];
		}
	}
	for(int i=1;i<=tot;i++){
		s1[i] = (i)%mod;
		s2[i] = (s2[i-1] + p[i] % mod) %mod;
	}return;
}
ll id(ll x){
	if(x <= qn) return id1[x];
	else return id2[n/x];
}
ll S(ll n,ll j){
	if(p[j] > n) return 0;
	ll res = (g2[id(n)] - g1[id(n)])%mod - (s2[j] - s1[j])%mod;
	res=(res%mod+mod)%mod;
	for(int k=j+1;k<=tot && p[k] * p[k] <= n;k++){
		for(ll e=1,now=p[k];now<=n;now*=p[k],e++){
			(res += (p[k] ^ e) %mod * (S(n/now,k) + (e>1)) %mod)%=mod;
		}
	}
	return res;
}
ll f1(ll x){
	return (x-1) %mod ;
}
ll f2(ll x){
	x%=mod;
	return x * (x+1) %mod *inv2 %mod -1;
}
int main(){
	scanf("%lld",&n);qn=sqrt(n);shai(qn);
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l),w[++m]=n/l;
		g1[m] = f1(w[m]);
		g2[m] = f2(w[m]);
		if(w[m] <= qn) id1[w[m]]=m;
		else id2[n/w[m]] = m;
	}//离散化 
	for(ll j=1;j<=tot;j++){//j是阶段，我们先枚举j，当然直接枚到tot就行，因为筛的质数最大到根号n 
		for(ll i=1;i<=m && p[j] * p[j] <= w[i];i++){
			(g1[i] -= ( 1) %mod * (g1[id(w[i]/p[j])] - s1[j-1] )%mod )%=mod;
			(g2[i] -= p[j] %mod * (g2[id(w[i]/p[j])] - s2[j-1] )%mod )%=mod;
			(g1[i] += mod )%=mod;(g2[i] += mod) %=mod;
		}
	}
	ll ans=S(n,0) + 1 ;
	if(n>=2) ans+=2;
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}

---

- 51nod 1847 奇怪的数学题

抄题解：小蒟蒻yyb，lcw

我们记 $\mathrm{sgcd}(n,m)$ 为 $n,m$ 的次大公约数。

给定 $1\le n\le {10}^9,1\le \alpha \le 50$，求：

$$\sum _{i=1}^n\sum _{j=1}^n\mathrm{sgcd}(i,j{)}^{\alpha }$$

对 $2^{32}$ 取模的值。

（为了防止变量名重复，我们把原题中的 $k$ 均替换为 $\alpha$）

这个题融合了杜教筛，min_25筛，莫反推式子，甚至从未听说过的第二类斯特林数求自然数幂和！

为了方便表示，我们令 $\mathrm{lcp}(a)=\underset{p|a}{min}p$

显然，我们有：

$$\mathrm{sgcd}(a,b)=\frac{gcd(a,b)}{\mathrm{lcp}(gcd(a,b))}$$

于是我们推式子，枚举 $gcd$：

$$\sum _{d=1}^n{\left(\frac{d}{\mathrm{lcp}(d)}\right)}^{\alpha }\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)=1]$$

如果 $i$ 和 $j$ 的上界不同，我们只能莫反拆开，但由于这里两者上界相同，我们直接考虑利用 $\phi$ 的性质将其转化成 $\phi$：

$$\sum _{d=1}^n{\left(\frac{d}{\mathrm{lcp}(d)}\right)}^{\alpha }(2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (i)-1)$$

（注意这个“$-1$”在关于 $i$ 的求和外面！）

于是我们对 $d$ 进行数论分块，后者可以杜教筛筛 $\phi$，我们重点考虑前者的前缀和怎么求。

我们设 $f(n)={\left({\displaystyle \frac{n}{\mathrm{lpf}(n)}}\right)}^{\alpha }$，求 $\sum _{i=1}^{n}f(i)$。

我们讨论 $n$ 为质数和合数，分别考虑：

$n$ 为质数

此时显然有 $f(n)=1$，故我们索求的其实就是 $n$ 以内的素数个数。

我们使用不完整的 min_25 筛：令 $f^{\prime }(p)=1$，我们令 $g$ 函数：

$$g(n,j)=\sum _{i=2}^n{f}^{\prime }(i)[i\in \mathbb{P}\mid \underset{p|i}{min}p>P_j]$$

对这个函数求：

$$g(n,\underset{P_k^2\le n}{max}k)$$

即 $\sum _{i=2}^n{f}^{\prime }(i)[i\in \mathbb{P}]$ 即可，显然这两者时等价的。

当然，我们最终求出来的 $g$ 数组，包含了 $n$ 取遍 $\lfloor {\displaystyle \frac{n}{1}}\rfloor ,\lfloor {\displaystyle \frac{n}{2}}\rfloor \dots \lfloor {\displaystyle \frac{n}{n}}\rfloor$ 这些值，所有的 $g$ 之值。这些数会在统计答案是有用。

（关于为什么把 $g$ 的定义写出来，是害怕自己看不懂）

$n$ 为合数

我们令 $f^{\prime }(n)=n^{\alpha }$，对 $f^{\prime }$ 进行不完整的 min_25 筛。

我们把所有的合数按照其最小质因数进行分类。

然后我们关注 min_25 筛中 $g$ 函数的递推过程，发现有：

当 $P_j^2\le n$ 时，$g(n,j-1)$ 比 $g(n,j)$ 多出来的贡献为：

$$f^{\prime }(P_j)\times (g(\frac{n}{P_j},j-1)-\sum _{i=1}^{j-1}{f}^{\prime }(P_i)))$$

我们这里其实找到了一类合数，它们有共同特征：最小质因子为 $P_j$。

我们要求这一系列合数，它们的 ${\left({\displaystyle \frac{n}{\mathrm{lpf}(n)}}\right)}^{\alpha }$ 之和，即 $f$ 之和（注意是 $f$ 不是 $f^{\prime }$！不要搞混！）。

其实我们发现上述贡献式中后面一部分，即 $g(\frac{n}{P_j},j-1)-\sum _{i=1}^{j-1}{f}^{\prime }(P_i))$，它的实际含义其实就是这些合数，除掉它们的最小质因子 $P_k$ 后，剩余东西的 $\alpha$ 次方和！也就是 $f$ 之和！！！

故我们在求 $g(n,0\sim \underset{P_k^2\le n}{max}k)$，时，把所有的 $g(\frac{n}{P_j},j-1)-\sum _{i=1}^{j-1}{f}^{\prime }(P_i))$加到一个 $h_n$ 里面，就能够得到对于 $n$ 以内所有的合数，它们的 $f$ 之和。

注意这里的 $n$ 会取遍 $\lfloor {\displaystyle \frac{n}{1}}\rfloor ,\lfloor {\displaystyle \frac{n}{2}}\rfloor \dots \lfloor {\displaystyle \frac{n}{n}}\rfloor$ 这些值，所以我们也会得到一个 $h$ 数组，存了 $h_{\lfloor \frac{n}{1}\rfloor },h_{\lfloor \frac{n}{2}\rfloor }\dots h_{\lfloor \frac{n}{n}\rfloor }$ 这些值，这在我们最终统计答案时有用。

---

还有一个东西差点忘了说了，在我们求合数部分时，我们需要为 $g$ 赋初始值，即：

$$g(n,0)=\sum _{i=2}^n{i}^{\alpha }$$

显然不能枚举，我们需要拉插，但模数不支持逆元，所以我们只能用别的方式。

这里就用了一个神奇的东西：第二类斯特林数求自然数幂和。

即：

$$\sum _{i=1}^n{i}^{\alpha }=\sum _{j=0}^{\alpha }\left\{\begin{array}{c}k\\ j\end{array}\right\}\frac{(n+1{)}^{\underset{―}{j+1}}}{j+1}$$

其中下降幂为 $x^{\underset{―}{n}}=x(x-1)(x-2)\dots (x-n+1)$。

$\left\{\begin{array}{c}k\\ j\end{array}\right\}$ 为第二类斯特林数，有个递推式：

$\left\{\begin{array}{c}n\\ m\end{array}\right\}=\left\{\begin{array}{c}n-1\\ m-1\end{array}\right\}+m\left\{\begin{array}{c}n-1\\ m\end{array}\right\}$

初始值：$\left\{\begin{array}{c}0\\ 0\end{array}\right\}=1$

这个东西可以 $O({\alpha }^2)$ 预处理处斯特林数，每次 $O(k^2)$ 查询自然数幂和，一共查询 $\sqrt{n}$ 次，复杂度可以接受。

---

最后我们可以数论分块原式子，通过手玩，可以发现数论分块的右端点 $r$ 组成的集合，与数论分块中 $\lfloor \frac{n}{l}\rfloor$ 组成的集合相同！于是我们只需要进行一次上述过程，就能把前半部分所有取值的前缀和都求出来！

而后面的 $\phi$ 通过杜教筛算就行。

---

说一点后话：啊啊啊啊啊啊啊啊啊啊啊啊啊啊逆天啊啊啊啊啊啊啊啊啊啊啊

我调了好久，原因是：

for(int j=1;j<=tot;j++){
	for(int i=1;i<=m &&1ll * prime[j] * prime[j] <= w[i];i++){
		int idx=id(w[i]/prime[j]);
		g1[i] -= (g1[idx] - s1[j-1] ) ;
		h[i] += g2[idx] - s2[j-1];
		g2[i] -= pk[j]  * (g2[idx] - s2[j-1] ) ;
	}
}

这一段代码，由于我筛出了 ${10}^6$ 以内的所有质数，所以 prime[j] * prime[j] 必须要乘以 1ll ！！！！！不然就会炸掉，然后莫名进入这个内层循环跑一遍！！！！！！！！！！

太离谱了！！！！！！！！！！！！

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef unsigned int ll;
const int N=1e9+10,M=1e6+10,K=1e6+10,ALP=55;
ll n,alp;
ll prime[K],tot,v[K],phi[K];
ll s1[M],s2[M],pk[M];
ll qpow(ll a,ll n){
	ll res=1;
	while(n){
		if(n&1) res *= a;
		a*=a;n>>=1;
	}return res;
}
void shai(ll n){
	phi[1]=1;
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,prime[++tot]=i;
			phi[i]=i-1;
		}for(ll j=1;j<=tot;j++){
			if(prime[j] > v[i] || prime[j] > n/i) break;
			v[i*prime[j]] = prime[j];
			if(i%prime[j] == 0) phi[i*prime[j]] = phi[i] * prime[j];
			else phi[i*prime[j]] = phi[i] * (prime[j]-1);
		}
	}
	for(ll i=1;i<=n;i++) phi[i] += phi[i-1];
	for(ll i=1;i<=tot;i++){
		s1[i] = s1[i-1] + 1;
		pk[i] = qpow(prime[i],alp);
		s2[i] = s2[i-1] + pk[i];
	}
}
ll stl[ALP][ALP];
ll zrsmh(ll n){
	ll res=0;
	if(n<=alp){
		for(ll i=1;i<=n;i++){
			ll tmp=1;
			for(ll j=1;j<=alp;j++) tmp *= i;
			res += tmp;
		}return res;
	}
	for(ll j=0;j<=alp;j++){
		ll tmp=1;
		for(ll l=0;l<j+1;l++){
			if((n+1-l) % (j+1) )tmp *= (n+1-l);
			else tmp *= (n+1-l) / (j+1);
		}res += tmp * stl[alp][j];
	}return res;
}
ll m,qn,w[M],id1[M],id2[M],g1[M],g2[M],h[M];
ll f1(ll x){
	return x-1;
}
ll f2(ll x){
	return zrsmh(x) - 1;
}
int id(int x){
	if(x <= qn) return id1[x];
	else return id2[n/x];
}
ll tph[M];
bool vis[M];
ll Sphi(ll x){
	if(x<=K-10) return phi[x];
	if(vis[id(x)]) return tph[id(x)];
	ll res=1ll * x*(x+1)/2;
	vis[id(x)]=1;
	for(ll l=2,r;l<=x;l=r+1){
		r= x/(x/l);
		res-=(r-l+1) * Sphi(x/l);
	}return tph[id(x)]=res;
}
void Init(){
	qn=sqrt(n);
	stl[0][0]=1;
	for(ll i=1;i<=alp;i++){
		for(ll j=1;j<=alp;j++){
			stl[i][j]=stl[i-1][j-1] + j * stl[i-1][j];
		}
	}
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l),w[++m]=n/l;
		g1[m] = f1(w[m]);
		g2[m] = f2(w[m]);
		if(w[m] <= qn) id1[w[m]]=m;
		else id2[n/w[m]] = m;
	}
	for(int j=1;j<=tot;j++){
		for(int i=1;i<=m &&1ll * prime[j] * prime[j] <= w[i];i++){
			int idx=id(w[i]/prime[j]);
			g1[i] -= (g1[idx] - s1[j-1] ) ;
			h[i] += g2[idx] - s2[j-1];
			g2[i] -= pk[j]  * (g2[idx] - s2[j-1] ) ;
		}
	}
}
ll ask(ll x){
	return g1[id(x)] + h[id(x)] ;
}
int main(){
	scanf("%u%u",&n,&alp);
	shai(K-10);Init();
	ll ans=0;
	Sphi(n);
	ll lst=0,now;
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		now=ask(r);
		ans += (2*Sphi(n/l)-1)*(now-lst);
		lst=now;
	}
	cout<<ans<<endl;
	return 0;
}

---

（三）.自然数幂和与拉格朗日插值

-BZOJ2137. submultiple

给定 ${\displaystyle m=\prod _{i=1}^n{p}_i^{a_i}}$，求：

$$\sum _{j\mid m}d(j{)}^k$$

其中测试点 $\mathbf{\text{A}}$ 中 $1\le a_i\le {10}^5,1\le k\le 2^{31}-1$，测试点 $\mathbf{\text{B}}$ 中 $1\le a_i\le 2^{63}-1,1\le k\le 12$

既然题目给好了 $m$ 的质因数分解，我们可以按照每个质因数选多少个，来枚举 $m$ 的因数。

而约数个数函数 $d$ 刚好有一个和质因数次数相关的公式，于是原式：

$$\sum _{i_1=0}^{a_1}\sum _{i_2=0}^{a_2}\dots \sum _{i_n=0}^{a_n}((i_1+1)(i_2+1)\dots (i_n+1){)}^k$$

也就是：

$$(\sum _{i_1=0}^{a_1}(i_1+1{)}^k)(\sum _{i_2=0}^{a_2}(i_2+1{)}^k)\dots (\sum _{i_n=0}^{a_n}(i_n+1{)}^k)$$

整个写成 $\prod$ 记号：

$$\prod _{j=1}^n\sum _{i=1}^{a_j+1}{i}^k$$

这其实就是算 $n$ 次自然数幂和，我们分别解决不同范围的测试点。

1.测试点 $\mathbf{\text{A}}$

这个点中 $a$ 小，$k$ 大，直接预处理+前缀和就行了。

复杂度 $O(k\log k\sum a)$

2.测试点 $\mathbf{\text{B}}$

这个点中 $a$ 大，$k$ 小，直接拉格朗日插值就行了。

复杂度 $O(n{k}^2)$

具体做法在这篇博客里。

这样这个题就做完了

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1e5+10,mod=1000000007;
ll n,k,a[N],maxn;
ll sum[N],inv[N];
ll qpow(ll a,ll n){
	ll res=1;
	while(n){
		if(n&1) res=(res*a)%mod;
		a=(a*a)%mod,n>>=1;
	}return res;
}
void getS(ll n){
	for(int i=1;i<=n;i++)sum[i]=(sum[i-1] + qpow(i,k))%mod;
}
void solve1(){
	ll ans=1;getS(maxn+1);
	for(int i=1;i<=n;i++){
		(ans*=sum[a[i]+1])%=mod;
	}printf("%lld\n",(ans+mod)%mod);
}
void solve2(){
	getS(k+2);ll ans=1;
	for(int i=-(k+2);i<=k+2;i++){
		inv[i+k+2]=qpow((i+mod)%mod,mod-2);
	}
	for(int d=1;d<=n;d++){
		ll cnt=0;
		for(int i=1;i<=k+2;i++){
			ll tot=1;
			for(int j=1;j<=k+2;j++){
				if(i==j) continue;
				(tot*=(a[d]+1-j)%mod * inv[i-j+k+2] %mod)%=mod;
			}
			(cnt+=sum[i]*tot%mod)%=mod;	
		}
		(ans*=cnt)%=mod;
	}printf("%lld\n",(ans+mod)%mod);
}
int main(){
	scanf("%lld%lld",&n,&k);
	for(int i=1;i<=n;i++){
		scanf("%lld",&a[i]);
		maxn=max(maxn,a[i]);
	}
	if(k>12) solve1();
	else solve2();
	return 0;
}

---

-BZOJ3453. tyvj 1858 XLkxc

有如下函数：

$$\begin{array}{rl}& f(x)=\sum _{i=1}^x{i}^k\\ & g(x)=\sum _{i=1}^{x}f(x)\\ & h(x)=\sum _{i=0}^{x}g(a+id)\end{array}$$

现在给定 $1\le k,a,n,d\le 123456789$，求 $h(n)mod1234567891$。

我们先把要求的这个东西直接拆开！

$$\sum _{u=0}^n\sum _{i=1}^{a+ud}\sum _{j=1}^i{j}^k$$

cy言：「积分多一次，求导少一次」

我们发现 $j^k$ 是个 $k$ 次的，而求和符号 $\sum$ 可以看做特殊的积分符号 $\int$，

那么这个式子相当于在 $k$ 次的式子上套了三个积分，那就多三次，所以这是个关于 $n$ 的 $k+3$ 次多项式。

于是我们只要知道这个式子在 $n\in [1,k+4]$ 时的取值，就可以拉格朗日插值插出来这个式子。

我们考虑怎么求值。

第一层，枚举 $n\in [1,k+4]$，复杂度 $O(k)$；

第二层，枚举 $u\in [0,n]$，复杂度 $O(k)$；（严格来说外面两层一起算是 $O(n\log n)$ 的）

第三层，枚举 $i\in [1,a+ud]$，发现太大了，没法枚举。但如果我们记这么一个函数：

$$t(n)=\sum _{i=1}^n\sum _{j=1}^i{j}^k$$

，那么这一层我们要求的就是 $t(a+ud)$。

$t$ 函数中 $n$ 的取值是 ${10}^{16}$，很大，于是我们可以考虑再套一层拉格朗日插值，把 $t$ 函数插出来。

根据「积分多一次，求导少一次」，我们发现，$t$ 函数是个 $k+2$ 次的多项式，需要求出 $t(n),n\in [1,k+3]$ 这 $k+3$ 个点值，才能插出来。

我们把此时的预处理放到最后说，我们假定现在可以插出来这个 $t$ 了。

于是第三层循环的复杂度就是一个拉格朗日插值的复杂度，我给它做到了 $O(k\log k)$。

所以这么一段求原式在 $n\in [1,k+4]$ 处的取值的处理，总复杂度是 $O((k\log k{)}^2)$，非常小，可以接受。（不确定有没有分析错，先放着吧，等人来 hack）

---

我们现在考虑一下 $t(n),n\in [1,k+3]$ 的取值咋求。

预处理 $i^n,n\in [1,k+3]$，然后直接枚举就行了，$O(k\log k)$。

---

现在我们已经求出了 $n\in [1,k+4]$，原式的取值。

那么我们直接跑拉格朗日插值硬插就行了，嫌麻烦直接写 $O(k^2)$ 的就行。

这样我们这个题就做完了，代码实现的细节蛮多，但好爽！

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1010,mod=1234567891,del=300;
ll T=1;
ll k,a,n,d;
ll aud[N],pd[N],nd[N];
ll inv[N];
ll qpow(ll a,ll n){
	ll res=1ll;
	while(n){
		if(n&1ll) res=(res*a)%mod;
		a=(a*a)%mod;n>>=1ll;
	}return res;
}
void Init_inv(){
	for(ll i=-230;i<=230;i++){
		inv[i+del]=qpow((i+mod)%mod,mod-2);
	}
}
void Init_aud(){//m=O(100)
	for(ll i=1;i<=k+100;i++){
		pd[i]=qpow(i,k)%mod;
	}
	for(ll n=1;n<=k+4;n++){
		aud[n]=0;
		for(ll i=1;i<=n;i++){
			for(ll j=1;j<=i;j++){
				(aud[n]+=pd[j]%mod)%=mod;
			}
		}
	}
}
void Init_n(){
	for(ll n=1;n<=k+4;n++){
		nd[n]=0;
		for(ll u=0;u<=n;u++){
			ll now=(a+u*d%mod)%mod;
			if(now<=k+3){
				(nd[n]+=aud[now])%=mod;
				continue;
			}
			ll M=1,Z=1,tot=0;
			for(ll j=1;j<=k+3;j++) (M*=(now-j)%mod)%=mod;
			for(ll j=-1;j>=-k-2;j--) (Z*=inv[j+del])%=mod;
			for(ll i=1;i<=k+3;i++){
				ll zi=M*qpow(((now-i)%mod+mod)%mod,mod-2) %mod;
				ll tmp=(-k-4+i)%mod+mod;
				if(i>1) (Z*=tmp%mod * inv[i-1+del]%mod)%=mod;
				(tot+=aud[i]*zi%mod*Z%mod)%=mod;
			}
			(nd[n]+=tot)%=mod;
		}
	}
}
int main(){
	Init_inv();
	scanf("%lld",&T);
	for(int Case=1;Case<=T;Case++){
		scanf("%lld%lld%lld%lld",&k,&a,&n,&d);
		Init_aud();
		Init_n();
		if(n<=k+4){
			printf("%lld\n",(nd[n]+mod)%mod);
			continue;
		}
		ll ans=0;
		for(ll i=1;i<=k+4;i++){
			ll tmp=1ll;
			for(ll j=1;j<=k+4;j++){
				if(i==j) continue;
				(tmp*=(n-j)%mod*inv[i-j+del]%mod)%=mod;
			}
			(ans+=tmp%mod*nd[i]%mod)%=mod;
		}printf("%lld\n",(ans+mod)%mod);
	}
	return 0;
}

---

- P6271 [湖北省队互测2014] 一个人的数论

给定：

$$1\le Q\le 100,1\le w\le 1000,N=\prod _{i=1}^w{p}_i^{{\alpha }_i}$$

求：

$$\sum _{i=1}^N[gcd(i,N)=1]i^Q$$

直接推式子：

$$\sum _{i=1}^N\sum _{d\mid gcd(i,N)}\mu (d)\cdot i^Q$$

枚举 $d$：

$$\sum _{d\mid N}\mu (d)\sum _{i=1}^{\frac{N}{d}}(id{)}^Q$$

也就是：

$$\sum _{d\mid N}\mu (d)\cdot d^Q\sum _{i=1}^{\frac{N}{d}}{i}^Q$$

给出了 $N$ 的质因数分解，且 $d$ 是 $N$ 的因数，所以直接枚举 $d$ 的质因数分解：

$$\sum _{i_1=0}^{{\alpha }_1}\sum _{i_2=0}^{{\alpha }_2}\dots \sum _{i_w=0}^{{\alpha }_w}\mu (\prod _{k=1}^w{p}_k^{i_k})\cdot (\prod _{k=1}^w{p}_k^{i_k}{)}^Q\cdot \sum _{j=1}^{{\displaystyle (\prod _{k=1}^w{p}_k^{{\alpha }_k-i_k})}}{j}^Q$$

注意到只要存在 $i_k>1$，$\mu$ 就是 $0$，所以只需要枚举每个 $i_k=[0,1]$：

$$\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^1\mu (\prod _{k=1}^w{p}_k^{i_k})\cdot (\prod _{k=1}^w{p}_k^{i_k}{)}^Q\cdot \sum _{j=1}^{{\displaystyle (\prod _{k=1}^w{p}_k^{{\alpha }_k-i_k})}}{j}^Q$$

借助 $\mathfrak{zzafanti}$ 之力，我们发现后面那一大坨是一个 $Q+1$ 次的关于那坨上界的多项式，我们写成多项式形式：

$$\sum _{j=0}^{Q+1}{a}_j(\prod _{k=1}^w{p}_k^{{\alpha }_k-i_k}{)}^j$$

于是乎，我们考虑枚举 $j$：

$$\sum _{j=0}^{Q+1}\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^1\mu (\prod _{k=1}^w{p}_k^{i_k})\cdot (\prod _{k=1}^w{p}_k^{i_k}{)}^Q\cdot a_j(\prod _{k=1}^w{p}_k^{{\alpha }_k-i_k}{)}^j$$

这几个函数都是积性函数，所以直接把里面的求积提出来：

$$\sum _{j=0}^{Q+1}\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^1{a}_j\prod _{k=1}^w\mu (p_k^{i_k})\cdot p_k^{Q{i}_k}\cdot p_k^{j{\alpha }_k-j{i}_k}$$

这里头有一个 $a_j\cdot \prod _{k=1}^w{p}_k^{j{\alpha }_k}$ 是与 $i$ 无关的，于是提到最前头：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^1\prod _{k=1}^w\mu (p_k^{i_k})\cdot p_k^{Q{i}_k}\cdot p_k^{-j{i}_k}$$

里面的东西稍微合并一下：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^1\prod _{k=1}^w\mu (p_k^{i_k})\cdot p_k^{(Q-j)i_k}$$

我们记 $f(k)=\mu (p_k^{i_k})\cdot p_k^{(Q-j)i_k}$，则原式：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\sum _{i_1=0}^1\sum _{i_2=0}^1\dots \sum _{i_w=0}^{1}f(1)\cdot f(2)\dots f(w)$$

发现 $f(1)$ 只与 $i_1$ 有关，$f(2)$ 只与 $i_2$ 有关……于是这些都可以提出去，原式变成：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\sum _{i_1=0}^{1}f(1)\sum _{i_2=0}^{1}f(2)\dots \sum _{i_w=0}^{1}f(w)$$

也就是：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\prod _{k=1}^w\sum _{i=0}^{1}f(k)$$

我们观察 $i_k$ 分别取 $0$ 和 $1$ 时，$f(k)$ 为多少：

1. $i_k=0$

$$f(k)=\mu (p_k^{i_k})\cdot p_k^{(Q-j)i_k}=\mu (p_k^0)\cdot p_k^0=1$$

2. $i_k=1$

$$f(k)=\mu (p_k^{i_k})\cdot p_k^{(Q-j)i_k}=\mu (p_k^1)\cdot p_k^{Q-j}=-p_k^{Q-j}$$

此时，原式就变成了：

$$\sum _{j=0}^{Q+1}(a_j\cdot \prod _{v=1}^w{p}_v^{j{\alpha }_v})\prod _{k=1}^w(1-p_k^{Q-j})$$

用高斯消元或者拉格朗日插值把 $a$ 插出来，复杂度用高斯消元的话是 $O(Q^3)$；

直接暴算这个式子，预处理 $\prod _{v=1}^w{p}_v^{{\alpha }_v}$是 $O(w\log max(\alpha ))$，预处理出 $p_k^j,k\in [1,w],j\in [1,Q]$ 是 $O(w\log Q)$，而整个式子在预处理的情况下能跑到 $O(wQ)$

总复杂度就是 $O(Q^3+wQ+w\log max(\alpha )+w\log Q)$，能过。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1010,mod=1e9+7;
ll Q,w,p[N],alp[N],P[N][N];
ll x[N],y[N],a[N],f[N],g[N],h[N];
ll qpow(ll a,ll n){
	if(n==-1) n=mod-2;
	ll res=1;
	while(n){
		if(n&1) res=(res*a)%mod;
		a=(a*a)%mod;n>>=1;
	}return res;
}
void lagrange(ll n){//时刻记住：给出了 Q+2 个点值！ 
	for(int i=0;i<n;i++) x[i]=i+1;
	for(int i=0;i<n;i++) y[i]=(y[i-1]+qpow(i+1,Q))%mod;//求点值
//	for(int i=0;i<n;i++){
//		printf("i:%d x:%lld y:%lld\n",i,x[i],y[i]); 
//	}
	for(int i=0;i<n;i++){
		a[i]=1;
		for(int j=0;j<n;j++){
			if(i==j) continue;
			(a[i] *= (x[i]-x[j])%mod )%=mod;
		}a[i] = qpow(a[i],-1) * y[i]%mod;
	}//求a
	g[0]=1;//最初，是个0次的多项式 "1" 
	for(int i=0;i<n;i++){//给它乘上一个 (m-x[i]) 
		for(int j=i+1;j>=1;j--){
			g[j] = ( g[j-1] - g[j] * x[i]%mod ) %mod;
		}g[0]= g[0] * (-x[i])%mod;
	}//求g
	for(int i=0;i<n;i++){
		ll inv=qpow(-x[i],-1);
		if(!inv){//x[i]=0
			for(int j=0;j<n;j++) h[j]=g[j+1];
		}else{
			h[0]=g[0] * inv%mod;
			for(int j=1;j<n;j++){
				h[j] = (g[j] - h[j-1])%mod *inv%mod;
			}
		}
		for(int j=0;j<n;j++){
			(f[j] += a[i] * h[j]%mod) %=mod; 
		}
	}
	for(int i=0;i<n;i++) f[i]=(f[i]+mod)%mod;
}
ll calc(ll n){
	ll res=0;
	for(int i=Q+1;i>=0;i--) res= (res * n%mod + f[i])%mod;
	return (res+mod)%mod;
}
int main(){
	scanf("%lld%lld",&Q,&w);
	for(int i=1;i<=w;i++){
		scanf("%lld%lld",&p[i],&alp[i]);
	}
	lagrange(Q+2);
	ll A=1,NOWA=1,ans=0;
	for(int k=1;k<=w;k++){
		for(int i=-1;i<=Q;i++){
			P[k][i +1]=qpow(p[k],i);
		}
	}
	for(int v=1;v<=w;v++) (A*=qpow(p[v],alp[v]))%=mod;
	for(int j=0;j<=Q+1;j++){
		if(j) (NOWA*=A)%=mod;
		ll tmp=1;
		for(int k=1;k<=w;k++){
			(tmp*=(1-P[k][Q-j +1])%mod)%=mod;
		}
		(ans+=f[j] * NOWA%mod *tmp%mod)%=mod;
	}printf("%lld\n",(ans+mod)%mod);
	return 0;
}

这里附一张我在没求助 $\mathfrak{zzafanti}$ 之前盯着式子发呆画出的《一个「人」的数论》：

---

四.同余题目

- BZOJ2219. 数论之神

太困难了，于是看了题解。

给定正整数 $A,B,K$，求方程：

$$x^A\equiv B(\mathrm{mod}2K+1)$$

在 $[0,2K)$ 以内的解的个数。

首先我们将 $2K+1$ 分解质因数为 $\prod _{i=1}^k{p}_i^{c_i}$，我们索求的答案即为如下 $k$ 个方程：

$$x^A\equiv B(\mathrm{mod}{p}_i^{c_i})$$

的乘积。

这个东西能用中国剩余定理得出（待会再写）。

于是我们着手于解决问题：

$$x^A\equiv B(\mathrm{mod}{p}^{\alpha })$$

的解的个数（我们令 $B<p^{\alpha }$）。

分三种情况讨论：

1.$B=0$

即为求方程：

$$x^A\equiv 0(\mathrm{mod}{p}^{\alpha })$$

的解的数量。

这个方程的含义就是：$x^A$ 为 $p^{\alpha }$ 的倍数。我们先找到 $x$ 为 $p$ 的幂的最小整数解。

设 $x=p^k$，则有：

$$p^{kA}\equiv 0(\mathrm{mod}{p}^{\alpha })$$

即：

$$kA\ge \alpha$$

此时 $k$ 的最小整数解为 $\lfloor \frac{\alpha -1}{A}\rfloor +1$。

于是我们便有了 $x$ 的最小整数解：

$$x_0=p^{\lfloor \frac{\alpha -1}{A}\rfloor +1}$$

我们欲求解的个数，即为 $x<p^{\alpha }$ 时 $x_0$ 的倍数的个数，显然为 ${\displaystyle \frac{p^{\alpha }}{x_0}}$ 个。

于是这种情况下，解的个数即为：

$$p^{\alpha -(\lfloor \frac{\alpha -1}{A}\rfloor +1)}$$

2.$gcd(B,p^{\alpha })>1$

我们令 $B=b\cdot p^{cnt}$，其中 $b$ 与 $p$ 互质，则有：

$$x^A\equiv b\cdot p^{cnt}(\mathrm{mod}{p}^{\alpha })$$

我们将同余式化为等式，即：

$$x^A-b\cdot p^{cnt}=k\cdot p^{\alpha }$$

再令等式左右两边同时除以 $p^{cnt}$：

$$\frac{x^A}{p^{cnt}}-b=k\cdot p^{\alpha -cnt}$$

我们发现等式右边为 $p$ 的倍数，而 $b$ 不是 $p$ 的倍数，则 ${\displaystyle \frac{x^A}{p^{cnt}}}$ 也一定不为 $p$ 的倍数。

故当 $cnt$ 不为 $A$ 的倍数，即 $cntmodA\ne 0$ 时，${\displaystyle \frac{x^A}{p^{cnt}}}$ 一定是 $p$ 的倍数，此时方程无解。

我们再把方程改写为同余式的形式：

$$\frac{x^A}{p^{cnt}}\equiv b(\mathrm{mod}{p}^{\alpha -cnt})$$

由于 $cnt$ 为 $A$ 的倍数，于是可以写作：

$${\left(\frac{x}{p^{\frac{cnt}{A}}}\right)}^A\equiv b(\mathrm{mod}{p}^{\alpha -cnt})$$

此时 $gcd(b,p^{\alpha -cnt})=1$，我们可以通过第三种情况求出此方程在 $[0,p^{\alpha -cnt})$ 范围中的解的数量。

但这种情况下，由于 $x$ 的取值范围是 $[0,p^{\alpha })$，则 ${\displaystyle \frac{x}{p^{\frac{cnt}{A}}}}$ 的取值范围就是 $[0,p^{\alpha -\frac{cnt}{A}})$，显然这个的范围要比我们刚才求出来的解的范围要大。

由于在模意义下，上面方程的解每 $p^{\alpha -cnt}$ 个一循环，我们就能知道，最终要求的答案 $ans$，就是 $[0,p^{\alpha -\frac{cnt}{A}})$ 内，长度为 $p^{\alpha -cnt}$ 的块的数量，再乘以每个块内解的个数。

而这个块数很显然是 ${\displaystyle \frac{p^{\alpha -\frac{cnt}{A}}}{p^{\alpha -cnt}}}$，也就是 $p^{cnt-\frac{cnt}{A}}$。

所以这种情况下的答案，就是上述方程的解的数量，乘以 $p^{cnt-\frac{cnt}{A}}$。

3.$gcd(B,p^{\alpha })=1$

我们找到模 $p^{\alpha }$ 意义下的某一原根 $g$，由于 $B$ 与 $p^{\alpha }$ 互质，故一定存在一个 $y$，使：

$$g^y\equiv B(\mathrm{mod}{p}^{\alpha })$$

我们用 BSGS 求解出 $y$。

设 $x=g^t$，则有：

$$g^{At}\equiv g^y(\mathrm{mod}{p}^{\alpha })$$

通过指标的性质（其实可以感性理解为通过欧拉定理 $a^{\phi (p)}\equiv 1(\mathrm{mod}p)$，左右两边的指数关于 $\phi (p^{\alpha })$ 同余），有：

$$At\equiv y(\mathrm{mod}\phi (p^{\alpha }))$$

这个线性同余方程解的数量，即为原方程解的数量。

当 $ymodgcd(A,\phi (p^{\alpha }))\ne 0$ 时，无解；否则解的个数为 $gcd(A,\phi (p^{\alpha }))$。

这样，我们这个超级困难，超级牛逼，超级复杂的逆天题，终于做完了……

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1e6+10;
const long long INF=0x3f3f3f3f3f3f3f3f;
ll T;
ll A,B,p;
ll pp[N],c[N],t[N],num;
ll qpow(ll a,ll n,ll mod){
	ll res=1%mod;
	while(n){
		if(n&1) res=(res*a)%mod;
		a=(a*a)%mod;n>>=1;
	}return res;
}
ll gcd(ll __m, ll __n){
    while (__n != 0){
		ll __t = __m % __n;
		__m = __n;
		__n = __t;
	}return __m;
}
void divide(ll x){
	num=0;ll tmp=x;
	for(int i=2;i*i<=x;i++){
		if(x%i) continue;
		pp[++num]=i;c[num]=0;
		while(x%i==0) c[num]++,x/=i;
	}
	if(x!=1) pp[++num]=x,c[num]=1;
	for(int i=1;i<=num;i++) t[i]=tmp/pp[i];
	return;
}
ll getgen(ll phi,ll n){
	for(int i=1;i<=n;i++){
		if(gcd(i,n)!=1) continue;
		if(qpow(i,phi,n)!=1) continue;
		bool is=1;
		for(int j=1;j<=num;j++){
			if(qpow(i,t[j],n)==1) {
				is=0;break;
			}
		}
		if(is) return i;
	}return -114514;
}
ll BSGS(ll a,ll b,ll p){
	if(a%p==0){
		if(b%p==1&&a!=0) return 0;
		else if(b%p==0) return 1;
		else return -1;
	}
	map <ll,int> mp; mp.clear();
	b%=p;ll t=ceil(sqrt((double)p)),now=b;
	for(int j=0;j<t;j++){
		if(j) (now*=a)%=p;
		mp[now]=j;
	}
	now=1;a=qpow(a,t,p);
	for(int i=0;i<=t;i++){
		if(i) (now*=a)%=p;
		if(mp.find(now)!=mp.end()&&i*t-mp[now] >=0) return i*t-mp[now];  
	}return -1;
}
ll calc(ll p,ll alp,ll B){
	ll now=qpow(p,alp,INF);
	ll b=B%now;
	if(b==0){
		ll t=(alp-1)/A + 1;
		return qpow(p,alp-t,INF);
	}
	ll cnt=0;
	while(b%p==0) cnt++,b/=p;
//	printf("p:%lld alp:%lld now:%lld b:%lld cnt:%lld\n",p,alp,now,b,cnt);
	if(cnt==0){
		ll phi=now - qpow(p,alp-1,INF);
		divide(phi);
		ll ming=getgen(phi,now);
		ll y=BSGS(ming,b,now);
		ll ty;
		ll tmp=gcd(A,phi);
		if(y%tmp) return 0;
		return tmp;
	}
	if(cnt%A) return 0;
	return calc(p,alp-cnt,b) * qpow(p,cnt-cnt/A,INF);
}
int main(){
	scanf("%lld",&T);
	for(int Case=1;Case<=T;Case++){
		scanf("%lld%lld%lld",&A,&B,&p);
		p=2*p+1;ll res=1;
		for(int i=2;i*i<=p;i++){
			if(p%i) continue;
			ll now=1,alp=0;
			while(p%i==0) p/=i,now*=i,alp++;
//			ll tmp=calc(i,alp,now);printf("p:%d alp:%lld now:%lld res:%lld\n",i,alp,now,tmp);
			res *= calc(i,alp,B);
		}
		if(p!=1){
			res *= calc(p,1,B);
		}
		printf("%lld\n",res);
	}
	return 0;
}

---

- P4139 上帝与集合的正确用法

要求这么一个东西：

$$2^{2^{2^{\dots }}}modp$$

的值。

我们有扩展欧拉定理：

$$a^b\equiv a^{bmod\phi (p)+\phi (p)}(\mathrm{mod}p)$$

故：

$$2^{2^{2^{\dots }}}\equiv 2^{\left(2^{2^{\dots }}\right)mod\phi (p)+\phi (p)}(\mathrm{mod}p)$$

然后递归求解 $2^{2^{2^{\dots }}}mod\phi (p)$ 即可，复杂度 $O(\log p)$。

当然可以加个记忆化，可能跑的更快。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=1e7+10;
int prime[N],tot,v[N],phi[N];
ll T=1,p;
void shai(int n){
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(v[i]==0){
			v[i]=i;prime[++tot]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=tot;j++){
			if(prime[j] > v[i] || prime[j] > n/i) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0){
				phi[i*prime[j]]=prime[j] * phi[i];
			}else{
				phi[i*prime[j]]=(prime[j]-1)*phi[i];
			}
		}
	}
}
ll qpow(ll a,ll n,ll mod){
	ll res=1%mod;
	while(n){
		if(n&1) res=(res*a)%mod;
		a=(a*a)%mod,n>>=1;
	}return res;
}

ll calc(ll p){
	if(p==1||p==2) return 0;
	ll tmp=calc(phi[p])+phi[p];
	return qpow(2,tmp,p);
}
int main(){
	shai(N-10);//0.3s
	scanf("%lld",&T);
	for(int Case=1;Case<=T;Case++){
		scanf("%lld",&p);
		printf("%lld\n",calc(p));
	}
	return 0;
}

---

- BZOJ4454. C Language Practice

给定两个长度分别为 $1\le n,m\le 2\cdot {10}^3$ 的数组 $a$ 和 $b$，求：

$$\sum _{i=0}^{n-1}\sum _{j=0}^{m-1}gcd(a_i,b_j)\mathrm{xor}i\mathrm{xor}j$$

满足 $0\le a_i,b_j\le {10}^6$。

式子推不了一点，遂直接枚举。使用黑科技 $O(1)gcd$，卡常，如过。

如过，遂无提交记录。

点击查看代码

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef unsigned int ll;
const int N=2010,M=1020,maxm=1e6+10,INF=0x3f3f3f3f;
ll T=1,t=1010;
ll a[N],b[N],n,m;
ll prime[maxm],tot;
ll v[maxm],ggcd[N][N],split[maxm][3];
void shai(ll n){
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			v[i]=i,prime[++tot]=i;
		}for(int j=1;j<=tot;j++){
			if(prime[j] > v[i] || prime[j] > n/i) break;
			v[i*prime[j]] = prime[j];
		}
	}return;
}
void Init(){
	for(int i=1;i<=t;i++) ggcd[i][0]=ggcd[0][i]=i;
	for(int i=1;i<=t;i++){
		for(int j=1;j<=t;j++){
			ggcd[i][j]=ggcd[j%i][i];
		}
	}
	split[1][0]=split[1][1]=split[1][2]=1;
	for(int i=2;i<=maxm-10;i++){
		ll minn=INF,mini;
		for(int j=0;j<3;j++){
			split[i][j]=split[i/v[i]][j];
			if(split[i][j] < minn) minn=split[i][j],mini=j;
		}
		split[i][mini] *= v[i];
	}
}
int Gcd(int x,int y){
	if(x<=1010 && y<=1010) return ggcd[x][y];
	if(!x || !y) return x+y;
	ll res=1;
	for(int i=0;i<3;i++){
		if(split[x][i]==1) continue;
		ll tmp;
		if(split[x][i] <= 1010) tmp=ggcd[split[x][i]][y%split[x][i]];
		else if(y%split[x][i]==0) tmp=split[x][i];
		else tmp=1;
		res *= tmp;y/=tmp;
	}return res;
}
int main(){
	shai(maxm-10);
	Init();
	cin>>T;
	ll res;
	ll ans[N],cnt=0;
	while(T--){
		scanf("%u%u",&n,&m);
		for(int i=0;i<n;i++) scanf("%u",&a[i]);
		for(int j=0;j<m;j++) scanf("%u",&b[j]);
		res=0;
		for(int i=0;i<n;i++){
			for(int j=0;j<m;j++){
				res += Gcd(a[i],b[j]) ^ i ^ j;
			}
		}
		printf("%u\n",res);
	}
	return 0;
}
//0.36->6.54

- hdu6209 The Intersection

给定两函数 $f(x)=\sqrt{x},g(x)={\displaystyle \frac{k}{x}}$，多次询问给定 $k$，求两函数交点的 $x$ 坐标。

我们联立解方程，有 $x=k^{\frac{2}{3}}$。

用 SB 树去逼近就行了。

注意开long double。

提交记录

---

- P5172 Sum

给定 $n,r$，求：

$$\sum _{i=1}^n(-1{)}^{i\lfloor \sqrt{r}\rfloor }$$

首先，我们用 SB 树找到一个最逼近 $\sqrt{r}$ 的分数，设为 ${\displaystyle \frac{a}{b}}$。

于是原式变为：

$$\sum _{i=1}^n(-1{)}^{\lfloor \frac{ai}{b}\rfloor }$$

这个东西的指数看起来和万欧很像，于是我们考虑使用万欧求这个东西。

每个操作序列上维护两个东西：

$$y,sum$$

我们考虑对 $l,r$ 进行合并，有：

$$y=y_l+y_r$$

$$sum=su{m}_l+\{\begin{array}{ll}su{m}_r& a_{l}mod2=0\\ -su{m}_r& a_{l}mod2=1\end{array}$$

我们跑一遍万欧，输出答案的 $sum$ 即可。

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define double long double
#define DEBUG
using namespace std;
typedef long long ll;
const long long N=1e16;
ll T=1;
struct node{
	ll y,sum;
	node() : y(0), sum(0) {};
	node operator *(node b){
		node a=*this,c;
		c.y=a.y + b.y;
		if(a.y &1) c.sum = a.sum - b.sum;
		else c.sum = a.sum + b.sum;
		return c;
	}
}U,R,res;
node qpow(node a,ll n){
	node res;
	while(n){
		if(n&1) res = res * a;
		a=a*a,n>>=1;
	}return res;
}
ll div(ll a,ll b,ll c,ll d){
	return ((double)1.0*a*b+c)/d;
}
node calc(ll p,ll q,ll r,ll n,node U,node R){
	if(!n) return node();
	if(r >= q) return qpow(U,r/q) * calc(p,q,r%q,n,U,R);
	if(p >= q) return calc(p%q,q,r,n,U,qpow(U,p/q) * R);
	ll m = div(p,n,r,q);
	if(!m) return qpow(R,n);
	return qpow(R,(q-r-1)/p) * U * calc(q,p,(q-r-1)%p,m-1,R,U) * qpow(R,n-div(q,m,-r-1,p));
}
pair <ll,ll> sbt(ll x){
	ll a=0,b=1,c=1,d=0;
	ll ra,rb;
	ll tmp=sqrt(x);
	if(tmp*tmp==x) return MP(tmp,1);
	while(1){
		ll e=a+c,f=b+d;
		if(f>N) break;
		if(e*e < f*f*x) a=e,b=f;
		else c=e,d=f;
		ra=e,rb=f;
	}
	return MP(ra,rb);
}
ll n,r;
int main(){
	scanf("%lld",&T);
	U.y=1;R.sum=1;
	for(int Case=1;Case<=T;Case++){
		cin>>n>>r;
		pair<ll,ll> qr=sbt(r);
//		cout<<qr.first<<"/"<<qr.second<<endl;
		node res=calc(qr.first,qr.second,0,n,U,R);
		printf("%lld\n",res.sum);
	}
	return 0;
}

---

- PE372 光锥

参考了 beginendzrq 的题解，_ANIG_ 的提示，Qcfff 写的代码。我抄抄抄抄抄抄

给定 $m,n$，求：

$$\sum _{x=m+1}^n\sum _{y=m+1}^n(\lfloor \frac{y^2}{x^2}\rfloor mod2)$$

保证 ${\displaystyle \frac{n}{m}}\le 500$。

我们枚举 $t=\lfloor {\displaystyle \frac{y^2}{x^2}}\rfloor$，有：

$$\sum _{t=1}^{\lfloor \frac{n^2}{m^2}\rfloor }(-1{)}^{t+1}\sum _{x=m+1}^n\sum _{y=m+1}^n[\sqrt{t}x\le y]$$

看一下它的几何含义：

我们要求红绿线围成的正方形里，黑色阴影中的整点数量。

$t$ 其实是在枚举斜率，每条直线都是一个以 $\sqrt{t}$ 为斜率的直线。

我们可以考虑求出每条直线与正方形上边，左边围成的三角形中有多少个点。

以 $y=\sqrt{2}x$ 为例：

我们分别用万欧求出 $\mathrm{△}OAB$ 中整点个数，记为 $u$，再求出 $\mathrm{△}OCD$ 内的整点个数，记为 $v$。

$u-v$ 即为蓝色梯形 $ABCD$ 中整点个数。

我们再用大长方形 $BECD$ 中的整点个数，减去 $u-v$，就能得到 $\mathrm{△}AEC$ 中的整点个数。

当然我们需要考虑线段 $AC$ 上的整点个数。

当 $t$ 不为平方数，上述做法完全可以。

当 $t$ 为平方数，我们万欧时会计算 $AC$ 上的整点，而用长方形一减就减没了。于是我们还需要加上直线上的整点数。

具体看代码：

for(int k=1;k<=lim;k++){
	ll now=0;
	pair <ll,ll> tmp=sbt(k);
	ll l=m,r=(double)n/sqrt((double)k);
	now += calc(tmp.first,tmp.second,0,r,U,R).sumy;
	now -= calc(tmp.first,tmp.second,0,l,U,R).sumy;
	now = (r - l ) * (n) - now;
	ll qk=sqrt(k);
	if(qk*qk==k){
		now += (r-l);
	}
	if(k%2==1) res += now;
	else res -= now;
}printf("%lld\n",res);

r 是点 $C$ 的横坐标，即解方程 $\{\begin{array}{l}y=\sqrt{t}x\\ y=n\end{array}$ 得 $x={\displaystyle \frac{n}{\sqrt{t}}}$。

应该很直观了。我们万欧时套个 SB 树板子，把根号变成分数就能算了。

另外，特别鸣鸣鸣鸣鸣鸣鸣鸣鸣鸣鸣鸣谢 Qcfff 想出的做法和实现的代码，我都是抄他的，他真的我哭死呜呜呜

提交记录

点击查看代码

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const long long N=1e15;
ll m,n;
struct node{
	ll x,y,sumy;
	node() : x(0) , y(0), sumy(0) {};
	node operator *(node b){
		node a=*this,c;
		c.x=a.x + b.x;
		c.y=a.y + b.y;
		c.sumy = a.sumy + b.sumy + a.y * b.x;
		return c;
	}
}U,R,res;
node qpow(node a,ll n){
	node res;
	while(n){
		if(n&1) res = res * a;
		a=a*a,n>>=1;
	}return res;
}
ll div(ll a,ll b,ll c,ll d){
	return ((long double)1.0*a*b+c)/d;
}
node calc(ll p,ll q,ll r,ll n,node U,node R){
	if(!n) return node();
	if(r >= q) return qpow(U,r/q) * calc(p,q,r%q,n,U,R);
	if(p >= q) return calc(p%q,q,r,n,U,qpow(U,p/q) * R);
	ll m = div(p,n,r,q);
	if(!m) return qpow(R,n);
	return qpow(R,(q-r-1)/p) * U * calc(q,p,(q-r-1)%p,m-1,R,U) * qpow(R,n-div(q,m,-r-1,p));
}
pair <ll,ll> sbt(ll x){
	ll a=0,b=1,c=1,d=0;
	ll ra,rb;
	ll tmp=sqrt(x);
	if(tmp*tmp==x) return MP(tmp,1);
	while(1){
		ll e=a+c,f=b+d;
		if(f>N) break;
		if(e*e < f*f*x) a=e,b=f;
		else c=e,d=f;
		ra=e,rb=f;
	}
	return MP(ra,rb);
}
int main(){
	scanf("%lld%lld",&m,&n);
	ll lim=n*n/(m+1)/(m+1),res=0;
	cout<<lim<<endl;
	U.y=1,R.x=1;
	for(int k=1;k<=lim;k++){
		ll now=0;
		pair <ll,ll> tmp=sbt(k);
		ll l=m,r=(double)n/sqrt((double)k);
		now += calc(tmp.first,tmp.second,0,r,U,R).sumy;
		now -= calc(tmp.first,tmp.second,0,l,U,R).sumy;
		now = (r - l ) * (n) - now;
		ll qk=sqrt(k);
		if(qk*qk==k){
			now += (r-l);
		}
		if(k%2==1) res += now;
		else res -= now;
	}printf("%lld\n",res);
	return 0;
}