HDU 6833 莫比乌斯反演 + 数论分块

给定下列式子:
∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)f(gcd⁡(a1,a2...,ax))⋅gcd⁡(a1,a2...,ax)\sum_{a_1=1}^n\sum_{a_2=1}^n...\sum_{a_x=1}^n(\prod_{j=1}^xa_j^k)f(\gcd(a_1,a_2...,a_x))⋅ \gcd(a_1,a_2...,a_x)a1=1∑na2=1∑n...ax=1∑n(j=1∏xajk)f(gcd(a1,a2...,ax))⋅gcd(a1,a2...,ax)
       其中f(x)=[∃d>1:d2∣x]f(x)=[\exist d>1:d^2\mid x]f(x)=[∃d>1:d2∣x]，给定 t,k,x(1≤t≤104,1≤k≤109,1≤x≤109)t, k, x\ (1\le t \le 10^4,1\le k\le 10^9,1\le x\le 10^9)t,k,x (1≤t≤104,1≤k≤109,1≤x≤109)，下面 ttt 组查询，每次输入一个 n(1≤n≤2×105)n\ (1\le n\le 2\times 10^5)n (1≤n≤2×105)，输出对应公式结果。

       首先 f(x)=∣μ(x)∣f(x)=|\mu(x)|f(x)=∣μ(x)∣，令公式为 F(x)F(x)F(x)，化简如下（要实现单次询问至多 n\sqrt{n}n 的复杂度，所以当化简后的公式时间复杂度过高时，要想想能否通过更换枚举顺序来使得能预处理的内容更多）：

       令 t=gcd⁡(a1,a2...,ax)t=\gcd(a_1,a_2...,a_x)t=gcd(a1,a2...,ax)，S(n)=∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)S(n)=\sum_{a_1=1}^n\sum_{a_2=1}^n...\sum_{a_x=1}^n(\prod_{j=1}^xa_j^k)S(n)=∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)，根据乘法分配率，S(n)=(∑i=1nik)xS(n)=(\sum_{i=1}^ni^k)^xS(n)=(∑i=1nik)x。我们枚举 ttt，因为 gcd⁡\gcdgcd 为 ttt，所以gcd⁡(a1t,a2t,...,axt)\gcd(\frac{a_1}{t},\frac{a_2}{t},...,\frac{a_x}{t})gcd(ta1,ta2,...,tax) 一定等于1，则有
F(x)=∑t=1ntxk∑a1=1⌊nt⌋∑a2=1⌊nt⌋...∑ax=1⌊nt⌋(∏j=1xajk)f(t)⋅t[gcd⁡(a1,a2,...,ax)=1]=∑t=1ntxkS(⌊nt⌋)f(t)⋅t[gcd⁡(a1,a2,...,ax)=1]=∑t=1ntxkS(⌊nt⌋)f(t)⋅t∑d∣tμ(d)=∑t=1ntxk+1f(t)∑d=1⌊nt⌋μ(d)dxkS(⌊ntd⌋)\begin{aligned}F(x)&=\sum_{t=1}^nt^{xk}\sum_{a_1=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\sum_{a_2=1}^{\left\lfloor\frac{n}{t}\right\rfloor}...\sum_{a_x=1}^{\left\lfloor\frac{n}{t}\right\rfloor}(\prod_{j=1}^xa_j^k)f(t)⋅ t[\gcd(a_1,a_2,...,a_x)=1] \\ &= \sum_{t=1}^nt^{xk}S(\left\lfloor\frac{n}{t}\right\rfloor)f(t)⋅ t[\gcd(a_1,a_2,...,a_x)=1] \\ &=\sum_{t=1}^nt^{xk}S(\left\lfloor\frac{n}{t}\right\rfloor)f(t)⋅ t\sum_{d|t}\mu(d)\\&= \sum_{t=1}^nt^{xk+1}f(t)\sum_{d=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\mu(d)d^{xk}S(\left\lfloor\frac{n}{td}\right\rfloor)\end{aligned}F(x)=t=1∑ntxka1=1∑⌊tn⌋a2=1∑⌊tn⌋...ax=1∑⌊tn⌋(j=1∏xajk)f(t)⋅t[gcd(a1,a2,...,ax)=1]=t=1∑ntxkS(⌊tn⌋)f(t)⋅t[gcd(a1,a2,...,ax)=1]=t=1∑ntxkS(⌊tn⌋)f(t)⋅td∣t∑μ(d)=t=1∑ntxk+1f(t)d=1∑⌊tn⌋μ(d)dxkS(⌊tdn⌋)
       化简到这一步更换枚举对象，令 T=tdT=tdT=td，则有（这里枚举 ttt 或 ddd 都可以，若剩下除法，则计算对应的逆元）
F(x)=∑T=1nS(⌊nT⌋)∑t∣Tf(t)μ(Tt)txk+1(Tt)xk=∑T=1nS(⌊nT⌋)Txk∑t∣Tf(t)μ(Tt)t=∑T=1nS(⌊nT⌋)Txk∑t∣T∣μ(t)∣μ(Tt)t\begin{aligned} F(x)&=\sum_{T=1}^nS(\left\lfloor\frac{n}{T}\right\rfloor)\sum_{t|T}f(t)\mu(\frac{T}{t})t^{xk+1}(\frac{T}{t})^{xk}\\ &=\sum_{T=1}^nS(\left\lfloor\frac{n}{T}\right\rfloor)T^{xk}\sum_{t|T}f(t)\mu(\frac{T}{t})t\\&= \sum_{T=1}^nS(\left\lfloor\frac{n}{T}\right\rfloor)T^{xk}\sum_{t|T}|\mu(t)|\mu(\frac{T}{t})t\end{aligned}F(x)=T=1∑nS(⌊Tn⌋)t∣T∑f(t)μ(tT)txk+1(tT)xk=T=1∑nS(⌊Tn⌋)Txkt∣T∑f(t)μ(tT)t=T=1∑nS(⌊Tn⌋)Txkt∣T∑∣μ(t)∣μ(tT)t
       令 B(T)=Txk∑t∣T∣μ(t)∣μ(Tt)tB(T)=T^{xk}\sum_{t|T}|\mu(t)|\mu(\frac{T}{t})tB(T)=Txk∑t∣T∣μ(t)∣μ(tT)t，利用倍数枚举（调和级数复杂度）计算 μ(t)\mu(t)μ(t) ，μ(Tt)\mu(\frac{T}{t})μ(tT) 和 ttt 对 B(T)B(T)B(T) 的贡献即可。

       则单次查询的答案即为 ∑T=1nS(⌊nT⌋)B(T)\sum_{T=1}^nS(\left\lfloor\frac{n}{T}\right\rfloor)B(T)∑T=1nS(⌊Tn⌋)B(T)。采用数论分块，单次查询时间复杂度为 n\sqrt{n}n （实测 405ms405ms405ms）。

       然而，大佬还有更猛的做法（%%%）

       对于 ij≤n<(i+1)jij\le n< (i+1)jij≤n<(i+1)j，S(i)B(j)S(i)B(j)S(i)B(j) 都会对其中的 nnn 产生贡献。我们考虑构造一个差分数组，在 ijijij 的位置加，在 (i+1)j(i+1)j(i+1)j 的位置减 S(i)B(j)S(i)B(j)S(i)B(j)，然后求个前缀和就得到了答案。于是我们实现了 O(1)O(1)O(1) 查询（实测 249ms249ms249ms）。

#include<bits/stdc++.h>
#define endl '\n'
using namespace std;
typedef long long LL;
const int maxn = 2e5 + 5;
const LL mod = 1e9 + 7;
int cnt;
LL global_x, global_k;
int S[maxn], B[maxn], mu[maxn], vis[maxn], prime[maxn], ans[maxn];inline LL powmod(LL a, LL b){LL ans = 1;while(b){if(b&1) ans = ans*a%mod;a = a*a%mod;b >>= 1;}return ans;
}inline void init(){mu[1] = 1;for(int i = 2; i < maxn; i++){if(!vis[i]){prime[vis[i] = ++cnt] = i;mu[i] = -1;}for(int j = 1, v; j <= vis[i] && (v = i*prime[j]) < maxn; j++){vis[v] = j;mu[v] = -mu[i];if(j == vis[i]) mu[v] = 0;}}for(int i = 1; i < maxn; i++)for(int j = 1, tol = i; tol < maxn; j++, tol += i)B[tol] = (B[tol] + i*abs(mu[i])*mu[j])%mod;for(int i = 1; i < maxn; i++){S[i] = (S[i - 1] + powmod(i, global_k))%mod;B[i] = B[i]*powmod(i, global_x*global_k)%mod;}for(int i = 1; i < maxn; i++)S[i] = powmod(S[i], global_x);for(int i = 1; i < maxn; i++)for(int j = 1, tol = i; tol < maxn; j++, tol += i){ans[tol] = (ans[tol] + (LL)S[i]*B[j])%mod;if(tol + j < maxn) ans[tol + j] = (ans[tol + j] - (LL)S[i]*B[j]%mod + mod)%mod;}for(int i = 1; i < maxn; i++) ans[i] = (ans[i] + ans[i - 1])%mod;
}int main(){cin.tie(0);cout.tie(0);ios::sync_with_stdio(false);int t, n;cin >> t >> global_k >> global_x;init();while(t--){cin >> n;cout << ans[n] << endl;}
}

HDU 6833 莫比乌斯反演 + 数论分块相关推荐

BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+数论分块)
BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+数论分块) 手动博客搬家: 本文发表于20171216 13:34:20, 原地址https://blog.csd ...
P2522 HAOI2011 Problem b [莫比乌斯反演,数论分块]
P2522 HAOI2011 题意对于给出的n个询问,每次求有多少个数对(x,y)(x,y)(x,y),满足a≤x≤ba≤x≤ba≤x≤b,c≤y≤dc≤y≤dc≤y≤d,且gcd(x,y)=kgc ...
小A的数学题(莫比乌斯反演数论)
小A的数学题(莫比乌斯反演&数论) 1.容斥原始化简为 ∑ d = 1 n d 2 ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = d ] \sum\limi ...
luogu P3455 [POI2007]ZAP-Queries （莫比乌斯反演 + 整除分块）
整理的算法模板合集: ACM模板题目传送门本题中数据为5e4,我们只需要筛一次5e4就行了. 双倍经验的P4450 双亲数中数据达到了1e6,我们直接筛1e6的莫比乌斯函数有点不可取,因为只有一组 ...
欧拉函数+狄利克雷卷积+莫比乌斯函数+莫比乌斯反演+整除分块+杜教筛
Powered by:NEFU AB-IN 文章目录欧拉函数狄利克雷卷积莫比乌斯函数莫比乌斯反演 P3455 [POI2007]ZAP-Queries 整除分块 P2522 [HAOI2011 ...
HDU 6134 莫比乌斯反演
题目链接题意: 已知 n < = 1 e 6 n<=1e6 n<=1e6,求: f ( n ) = ∑ i = 1 n ∑ j = 1 i ⌈ i j ⌉ [ g c d ( i ...
HDU - Mophues(莫比乌斯反演)
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=4746 Time Limit: 2000/1000 MS (Java/Others) Memory Li ...
bzoj 3994 约数个数和 —— 反演+数论分块
题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3994 推导过程和这里一样:https://www.cnblogs.com/MashiroSk ...
数论 —— 莫比乌斯反演
[反演] 假设我们手头有个数列 F,通过某种变换 H,可以得到函数 G.,即: 但现在只有函数 G,需要求 F,那么我们就需要寻找一种变换 ,使得 G 在经过这种变换后能够获得 F,这个过程即为反演, ...

HDU 6833 莫比乌斯反演 + 数论分块

HDU 6833 莫比乌斯反演 + 数论分块相关推荐

最新文章

热门文章