@
如有错误，忘大佬不吝赐教，及时指出。
积性函数好难啊，求大佬教我积性函数，或者推荐点博客。。
如果感觉有点卡，请点这里。
莫比乌斯函数
参考：peng-ym
(sum_{d|n} mu(d) = [n = 1]) (（拆成二项式定理就很易证明了)
(sum_{d|n}phi(d) = n) (（wiki上的证明挺易理解,就是化简{frac 1n,frac 2n...frac nn})
(phi(n) = sum_{d|n} mu(d)frac{n}{d}) ((容斥过程))
证明：
莫比乌斯反演
(式子C: sum_{d|gcd(i,j)} mu(d) = [gcd(i,j) = 1])
(式子B: F(n)=sum_{d|n}f(d) => f(n)=sum_{d|n}u(d)F(⌊frac{n}{d}⌋))
(式子A: F(n)=sum_{n|d}f(d) => f(n)=sum_{n|d}u(frac{d}{n})F(d))
从莫比乌斯到欧拉：here
【CJOJ2512】gcd之和 -反演
求：(ans = sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j);;;nleq mleq1e7)
转换一下：
法1：
法2：(无敌易理解)
参考：这是一个神奇的Blog
AC_CODE： 和下面类似，就不写了。
P2257 YY的GCD -反演
题目：传送门
分析：
把(prime)提到前面去：
莫比乌斯反演一下：
把(k)提到前面来：
这样写会超时，因为多组数据，还得优化，就是交换求和顺序，令(T=kd：)
(O(nlog(n)))预处理这个式子(sum_{d|T}mu(frac{T}{d}))：枚举每个素数，给它的(x)倍数加上(mu(x))，然后前缀和。
这个式子也可以(O(n))预处理：
AC_CODE：
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e7+ 7;
const int mod = 998244353;
LL n, m;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN];
LL pre_mu[MXN];
void init_prime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
//phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
if(i % pp[j] == 0) {
//phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
break;
}
}
}
}
void init_ans() {
for(int i = 0; i < pcnt; ++i) {
for(int j = 1; j*pp[i] < MXN; ++j) {
pre_mu[j*pp[i]] += mu[j];
}
}
for(int i = 2; i < MXN; ++i) pre_mu[i] += pre_mu[i-1];
}
int main() {
init_prime();
init_ans();
int tim; scanf("%d", &tim);
while(tim --) {
scanf("%lld%lld", &n, &m);
if(n > m) swap(n, m);
LL ans = 0, tmp;
for(LL L = 2, R; L <= n; L = R + 1) {
R = min(n/(n/L),m/(m/L));
tmp = pre_mu[R]-pre_mu[L-1];
ans = (ans + tmp*(n/L)*(m/L));
}
printf("%lld
", ans);
}
return 0;
}
小清新数论 -杜教筛
题目：wannafly winter camp day3 F
这个式子与(gcd)之和那题非常相似，如果数据范围是(1e7)的话，这两题代码就基本一样了。可这题数据范围(1e11)，就用杜教筛咯。
杜教筛学习博客：here
方法1：利用欧拉函数前缀和
化简一下：$$ans = sum_{d=1}^{n}mu(d) sum_{i=1}{n}sum_{j=1}{m}[gcd(i,j)=d]=sum_{d=1}{n}mu(d)sum_{i=1}{frac{n}{d}}sum_{j=1}^{frac{n}{d}}[gcd(i,j)=1]$$
反演一下：$$ans=sum_{d=1}{n}mu(d)sum_{i=1}{frac{n}{d}}sum_{j=1}^{frac{n}{d}}sum_{k|gcd(i,j)}mu(k)$$
把(k)提出来：$$ans=sum_{d=1}{n}mu(d)sum_{k=1}{frac{n}{d}}mu(k)(frac n{kd})^2$$
又到了熟悉的步骤：令(T=kd)：
接着问题变成求下面这个式子的前缀和：$$sum_{i=1}^n (mu*mu)(i)$$
AC_CODE1：
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e7+ 7;
const LL mod = 998244353;
const LL inf = 1e18;
LL n, m, inv2;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN];
LL pre_mu[MXN], pre_phi[MXN];
//unordered_map<LL, LL> mp1, mp2;
void init_prime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
if(i % pp[j] == 0) {
phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) {
pre_mu[i] = pre_mu[i-1] + mu[i];
pre_phi[i] = (pre_phi[i-1] + phi[i])%mod;
}
}
struct Hash_map{
static const int mask=0x7fffff;
LL p[mask+1],q[mask+1];
void clear(){
memset(q,0,sizeof(q));
}
LL& operator [](LL k){
LL i;
for(i=k&mask;q[i]&&p[i]!=k;i=(i+1)&mask);
p[i]=k;
return q[i];
}
}mp1, mp2;
LL solve_mu(LL n) {
if(n < MXN) return pre_mu[n];
if(mp1[n] == inf) return 0;
if(mp1[n]) return mp1[n];
LL ans = 1;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans -= (R-L+1LL)%mod*solve_mu(n/L);
}
mp1[n] = ans;
if(mp1[n] == 0) mp1[n] = inf;
return ans;
}
LL solve_phi(LL n) {
if(n < MXN) return pre_phi[n];
if(mp2[n]) return mp2[n];
LL ans = n%mod*(n+1)%mod*inv2%mod;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1LL)%mod*solve_phi(n/L)%mod)%mod;
}
ans = (ans + mod) % mod;
mp2[n] = ans;
return ans;
}
int main() {
inv2 = 499122177;
init_prime();
scanf("%lld", &n);
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_mu(R)-solve_mu(L-1))%mod*(solve_phi(n/L)*2LL%mod-1LL)%mod+mod)%mod;
}
printf("%lld
", (ans+mod)%mod);
return 0;
}
AC_CODE2：
#include<bits/stdc++.h>
#define fi first
#define se second
#define iis std::ios::sync_with_stdio(false);cin.tie(0)
#define pb push_back
#define o2(x) (x)*(x)
using namespace std;
typedef long long LL;
typedef pair<int, int> pii;
typedef pair<LL, LL> pll;
const int INF = 0x3f3f3f3f;
const LL MOD = 998244353;
const int MXN = 5e6 + 6;
LL n, q, m;
int noprime[MXN], pp[MXN/2], pcnt;
int mu[MXN], phi[MXN];
int pre_mu[MXN];//莫比乌斯函数的前缀和
LL mumu[MXN];//莫比乌斯函数卷积莫比乌斯函数的前缀和
unordered_map<LL,LL>mp1,mp2;
void init_rime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i]=-1;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) pre_mu[i] = pre_mu[i-1] + mu[i], mumu[i] = mu[i];
for(int i = 2; i < MXN; ++i) {//预处理莫比乌斯卷积莫比乌斯
for(int j = i; j < MXN; j += i) {
mumu[j] += mu[i]*mu[j/i];
if(mumu[j] >= MOD) mumu[j] %= MOD;
}
}
for(int i = 2; i < MXN; ++i) mumu[i] = (mumu[i]+mumu[i-1]+MOD)%MOD;
}
LL solve_u(LL n) {//求莫比乌斯前缀和
if(n < MXN) return pre_mu[n];
if(mp1.count(n)) return mp1[n];
LL ans = 1;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1)%MOD*solve_u(n/L)%MOD+MOD)%MOD;
}
mp1[n] = ans;
return ans;
}
LL solve_uu(LL n) {//求莫比乌斯卷积莫比乌斯前缀和
if(n < MXN) return mumu[n];
if(mp2.count(n)) return mp2[n];
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_u(R)-solve_u(L-1)+MOD)%MOD*solve_u(n/L)%MOD);
if(ans >= MOD) ans %= MOD;
}
mp2[n] = (ans+MOD)%MOD;
return ans;
}
int main(int argc, char const *argv[]) {
init_rime();
scanf("%lld", &n);
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (n/L)%MOD*(n/L)%MOD*(solve_uu(R)-solve_uu(L-1)+MOD)%MOD)%MOD;
}
printf("%lld
", (ans+MOD)%MOD);
return 0;
}
/*
void init_rime() {//因为莫比乌斯卷积莫比乌斯是积性函数，所以可以On预处理
noprime[0] = noprime[1] = 1;
mu[1] = 1;mumu[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i]=-1, mumu[i]=-2;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
mumu[i*pp[j]] = mumu[i]*mumu[pp[j]];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
if((i/pp[j])%pp[j]) mumu[i*pp[j]] = mumu[i/pp[j]];
else mumu[i*pp[j]] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) pre_mu[i] = pre_mu[i-1] + mu[i];
for(int i = 2; i < MXN; ++i) mumu[i] = (mumu[i]+mumu[i-1]+MOD)%MOD;
}
*/
(求sum_{i=1}^nsum_{j=1}^{m}lcm(i,j))
2154 化简一下：
经过一系列化简可得：
这个式子是(O(n))，但是通过改变枚举方法可以(O(sqrt n))求解。
令(T=kd,h(n)=sum_{i=1}^ni)可得：
然后，你发现(f(n)=sum_{d|n}d imes mu(d))是积性函数。为什么呢？
证明1:
(inv(x)=x^{-1})是积性函数。
(f(n)=n imes sum_{d|n}frac 1d mu(frac nd)=n imes inv*mu(n))
证明2:
(s(x)=x imes mu(x))是积性函数
(f(x)=s*I(x)=sum_{d|n}d imes mu(d))
很自然的可以得到：(f(n))是一个积性函数。
当((a,b)=1)时，(a)的因子(j)和(b)的因子(k)当然也是互质的，然后(mu)是积性函数。
所以$$j imes mu(j) imes k imes mu(k)=jk imes mu(jk)$$
也就是说(f(a))的某一项乘上(f(b))的某一项一定是(f(ab))的某一项且是完备的。
当(gcd(a,b)=1)时：$$f(ab)=f(a) imes f(b)$$
当(gcd(a,b)!=1)且(b)是(a)的最小素因子时：$$f(ab)=f(a)$$
为什么呢？
我们现在有: (a)是(b)的倍数, (f(ab)=sum_{d|(ab)}d imes mu(d))：
如果(ab)的因数(d)素因数分解后(b)的系数小于等于(1)，那么这一项一定在(f(a))中也会出现；
如果(ab)的因数(d)素因数分解后(b)的系数大于(1)，那么(mu(d)=0)。故上式成立。
证明完(f(n)=sum_{d|n}d imes mu(d))是积性函数和上述式子，(f(n)和f(n) imes n)都可以线性筛的同时预处理出来：
void init_rime() {
noprime[0] = noprime[1] = 1;
F[1] = sum_F[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, F[i] = 1 - i;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
F[i*pp[j]] = F[pp[j]]*F[i];
if(i % pp[j] == 0) {
F[i*pp[j]] = F[i];
break;
}
}
sum_F[i] = sum_F[i-1] + F[i]*i;
}
}
所以(O(sqrt n))即可求解：$$ans=sum_{T=1}^nh(frac nT)h(frac mT) imes f(T) imes T，O(sqrt n)$$
(求sum_{i=1}^nsum_{j=1}^{m}d(i imes j))
3994 有一个公式：$$d(i imes j)=sum_{x|i}sum_{y|j}[gcd(x,y)=1]$$
因为(frac ni)表示的是(1,..,n)有几个数是(i)的倍数。当然也可以(nsqrt n)预处理。
(O(n))筛法：
bool noprime[MXN];
int pp[MXN/2], pcnt, mu[MXN], pre_mu[MXN];
int a1[MXN];
LL g[MXN];
void init_prime(int MXN) {
noprime[0] = noprime[1] = mu[1] = pre_mu[1] = 1;
g[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i] = -1, g[i] = 2, a1[i] = 1;
for(int j = 0; j < pcnt && i * pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1; mu[i*pp[j]] = -mu[i];
g[i*pp[j]] = g[i] * 2;
a1[i*pp[j]] = 1;
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
a1[i*pp[j]] = a1[i] + 1;
g[i*pp[j]] = g[i]*(a1[i*pp[j]]+1)/a1[i*pp[j]];
break;
}
}
pre_mu[i] = pre_mu[i-1] + mu[i];
}
for(int i = 2; i < MXN; ++i) g[i] += g[i-1];
}
预处理(2^{63})范围内(d(n))的前缀和：洛谷SP26073
小D的demo
牛客练习赛40：$$x=prod_{i=1}qp_i{k_i}，g(x)=sum_{i=1}^qk_i，g(1)=1$$
所以线性筛的时候，可以顺便把(g(n))给筛出来。
待求解式子：$$prod_{i=1}nprod_{j=1}m g(gcd(i,j))$$
提取(gcd)：$$prod_{d=1}^n g(d){sum_{i=1}nsum_{j=1}^m[gcd(i,j)=d]}$$
反演一下指数那个式子：$$sum_{i=1}nsum_{j=1}m[gcd(i,j)=d]=sum_{i=1}^{frac nd}sum_{j=1}^{frac md}[gcd(i,j)=1]=sum_{i=1}^{frac nd}sum_{j=1}^{frac md}sum_{k|gcd(i,j)}mu(k)$$
得到这个式子：$$prod_{d=1}^n g(d){sum_{k=1}{frac nd}mu(k)frac{n}{kd}frac{m}{kd}}$$
上面这个式子可以(O(n))求值了，但是还不够，因为这题有多组数据。
想到一个常用的优化方式：
先把指数提下来：$$prod_{d=1}^n prod_{k=1}^{frac nd} g(d)^{mu(k)frac{n}{kd}frac{m}{kd}}$$
令(T=kd)，换元求积：$$prod_{T=1}nprod_{d|T}g(d){mu(frac Td)frac nT frac mT}$$
再变换一下：$$prod_{T=1}n(prod_{d|T}g(d){mu(frac Td)})^{frac nT frac mT}$$
发现这个式子(prod_{d|T}g(d)^{mu(frac Td)})可以(O(nlog(n)))预处理出来，(O(1))求值。
所以最后总的复杂度是(O(nlog(n)+Tsqrt n log(MOD)))
代码：小D的demo
2440: [中山市选2011]完全平方数
(T=50)，求第(k(1e9))个非完全平方数倍数的数字。
显然二分是一个不错的选择，当你列出求(n)以内完全平方数倍数的个数的式子后，你会发现这东西不是可以用莫比乌斯函数求吗？都不需要反演的。复杂度就是枚举因数的复杂度：(O(sqrt n))。
HDU 6134
(nleq 1e6)，求这个式子的值：$$f(n)=sum_{d=1}nmu(d)sum_{x=1}{frac nd}sum_{y=1}^x ⌈frac xy⌉$$
牛批网友说：$$g(n)=sum_{i=1}^n⌈frac ni⌉，phi(n)sum_{i=1}^n g(i)$$
BZOJ 3529 [Sdoi2014]数表
(T(2e4) n,m(1e5)，A(1e9))，求：$$sum_{i=1}nsum_{j=1}msigma(gcd(i,j))$$
若不考虑(A)，可以很容易得到：$$sum_{d=1}nsigma(d)sum_{k=1}{frac nd}mu(k)frac n{kd}frac n{kd}$$
再优化一下：$$sum_{T=1}^nfrac nTfrac mTsum_{d|T}sigma(d) imes mu(frac Td)$$
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e5 + 20;
const LL mod = 1LL<<31;
/*
T(2e4) n,m(1e5)
不考虑A的话，式子很好化简，可以sqrt(n)算答案
这题不离线下来，处理巨奇怪，也不清楚正确性
离线下来的话，主要复杂度是O(n*sqrt(n)*log(n))吧？
先把莫比乌斯函数和因数和函数给线性筛出来，离线询问按A排序
不考虑A的话，我们正常做法就是枚举gcd，用的是g(i)算答案
考虑了A，我们还要按g(i)从小到大的顺序来枚举gcd，算出卷积的贡献
因为要求卷积的前缀和，所以要用一个数据结构维护，树状数组就够了，所以复杂度多了log。
取模也可以不用取模，因为模数是1LL<<31，利用int的自然溢出就行了
*/
bool noprime[MXN];
int pp[MXN/2], pcnt;
int mu[MXN], pre_mu[MXN];
int num1[MXN];//p1^a1
int g[MXN];//因数和函数
int p[MXN];
struct lp {
int n, m, A, id;
}cw[MXN];
int ANS[MXN];
int bit[MXN];
int lowbit(int x) {return x&(-x);}
void add_bit(int x, int val, int N) {for(;x <= N; x += lowbit(x)) bit[x] = (bit[x]+val);}
int query_bit(int x) {int ans = 0; for(; x; x -= lowbit(x)) ans = (ans+bit[x]); return ans;}
void init_prime(int MXN) {
noprime[0] = noprime[1] = mu[1] = pre_mu[1] = 1;
g[1] = 1; p[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i] = -1, g[i] = 1+i, num1[i] = i;
for(int j = 0; j < pcnt && i * pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
g[i*pp[j]] = g[i] * (1+pp[j]);
num1[i*pp[j]] = pp[j];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
num1[i*pp[j]] = num1[i] * pp[j];
g[i*pp[j]] = (g[i]/g[num1[i]])*(g[num1[i]]*pp[j]+1);
break;
}
}
pre_mu[i] = pre_mu[i-1] + mu[i]; p[i] = i;
}
}
bool cmp(const lp &a, const lp &b) {
return a.A < b.A;
}
bool cmp1(const int &a, const int &b) {
return g[a] < g[b];
}
int solve(int n, int m) {
int ans = 0, lst = 0, now;//除法分块的一点点优化，还是可以快一点哒
for(int L = 1, R; L <= n; L = R + 1) {
R = min(n/(n/L), m/(m/L));
now = query_bit(R);
ans = (ans+(n/L)*(m/L)*(now-lst));
lst = now;
}
return (ans);
}
int main(){
int tim; scanf("%d", &tim);
int ret = 0;
for(int i = 1; i <= tim; ++i) {
scanf("%d%d%d", &cw[i].n, &cw[i].m, &cw[i].A);
if(cw[i].n > cw[i].m) swap(cw[i].n, cw[i].m);
cw[i].id = i;
ret = max(ret, cw[i].n);
}
init_prime(ret+10);
sort(cw + 1, cw + 1 + tim, cmp);
int j = 1;
sort(p + 1, p + 1 + ret, cmp1);
for(int i = 1; i <= tim; ++i) {
for(; j <= ret && g[p[j]] <= cw[i].A; ++j) {
for(int k = p[j]; k <= ret; k += p[j]) {
if(mu[k/p[j]]) add_bit(k, g[p[j]]*mu[k/p[j]], ret);
}
}
ANS[cw[i].id] = solve(cw[i].n, cw[i].m);
}
for(int i = 1; i <= tim; ++i) {
if(ANS[i] < 0) ANS[i] += 2147483647, ++ ANS[i];
printf("%d
", ANS[i]);
}
return 0;
}
2019西安邀请赛B. Product
显然是一个欧拉降幂套上一个莫比乌斯反演，看这数据范围还得杜教筛。
由：$$prod_{i=1}{n}prod_{j=1}{n}prod_{k=1}{n}m{gcd(i,j)[k|gcd(i,j)]}$$
可化简为：$$m{sum_{i=1}nsum_{j=1}^ngcd(i,j) imes D(gcd(i,j))}$$
意思转化为求：$$sum_{d=1}^n d imes D(d)sum_{i=1}{n/d}sum_{j=1}{n/d}[(i,j)==1]$$
参考小清新数论那题：$$sum_{d=1}^n d imes D(d) imes (2 imes P(frac nd)-1)$$
又因为：$$sum_{i=1}^n i imes D(i)=sum_{i=1}^n i imes sum_{d|i}1=sum_{d=1}^n d imes sum_{i=1}^frac nd 1=sum d imes frac{frac nd(frac nd+1)}{2}$$
LL n, m;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN], a1[MXN];
LL pre_phi[MXN], pre_g[MXN], g[MXN];
void init_prime() {//线性筛预处理
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1; g[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1, g[i] = 2, a1[i] = 1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
g[i*pp[j]] = g[i] * 2;
a1[i*pp[j]] = 1;
if(i % pp[j] == 0) {
phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
a1[i*pp[j]] = a1[i] + 1;
g[i*pp[j]] = g[i]*(a1[i*pp[j]]+1)/a1[i*pp[j]];
break;
}
}
}
}
unordered_map<LL, LL> mp1, mp2;
LL solve_g(LL n) {
if(n < MXN) return pre_g[n];
if(mp1.find(n) != mp1.end()) return mp1[n];
//if(mp1[n]) return mp1[n];
//if(mp1.count(n)) return mp1[n];
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + ((n/L)*(n/L+1)/2)%mod*((R-L+1)*(L+R)/2)%mod)%mod;
}
ans = (ans + mod) % mod;
mp1[n] = ans;
return ans;
}
LL solve_phi(LL n) {//杜教筛欧拉函数前缀和
if(n < MXN) return pre_phi[n];
if(mp2.find(n) != mp2.end()) return mp2[n];
//if(mp2[n]) return mp2[n];
//if(mp2.count(n)) return mp2[n];
LL ans = n*(n+1)/2;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1)*solve_phi(n/L)%mod)%mod;
}
ans = (ans + mod) % mod;
mp2[n] = ans;
return ans;
}
int main() {
init_prime();
scanf("%lld%lld%lld", &n, &m, &mud);
mod = mud - 1;
for(int i = 1; i < MXN; ++i) {
pre_phi[i] = (pre_phi[i-1] + phi[i])%mod;
pre_g[i] = (pre_g[i-1] + i * g[i] % mod)%mod;
}
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_g(R)-solve_g(L-1))%mod*(solve_phi(n/L)*2LL%mod-1LL)%mod)%mod;
}
ans = (ans+mod)%mod;
printf("%lld
", ksm(m, ans, mud));
return 0;
}
// 当x≥ϕ(p)时，有a^x ≡ a^(x mod ϕ(p)+ϕ(p)) mod p
// gcd(a,mod)=1 a^(x) = a^(x%ϕ(mod)) % mod
南昌邀请赛和湘潭邀请赛
计蒜之道2019 1-D. 商汤AI园区的n个路口（困难）
插曲：二项式反演
插曲：集合反演
(max(S)=sum_{Tin S}(-1)^{|T|-1} imes min(T))
(min(S)=sum_{Tin S}(-1)^{|T|-1} imes max(T))