简单的数学题

题目连接

https://www.luogu.org/problemnew/show/P3768

题目描述

输入一个正整数n,n≤1010n,n\le 10^{10}n,n≤1010和p,p≤1.1×109p,p \le 1.1 \times 10^9p,p≤1.1×109.且ppp为质数.

计算∑i=1n∑j=1nijgcd(i,j)\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)∑i=1n∑j=1nijgcd(i,j)对ppp取模.

题解

方法一.暴力推公式

考虑把枚举gcd(i,j)gcd(i,j)gcd(i,j)
原式=∑d=1nd3∑i=1n/d∑j=1n/dij×[gcd(i,j)=1]\sum_{d=1}^nd^3\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij\times [gcd(i,j)=1]∑d=1nd3∑i=1n/d∑j=1n/dij×[gcd(i,j)=1]
莫比乌斯反演处理后边的式子
设f(x)=∑i=1n/d∑j=1n/dij×[gcd(i,j)=x]f(x)=\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij\times [gcd(i,j)=x]f(x)=∑i=1n/d∑j=1n/dij×[gcd(i,j)=x]
g(x)=∑x∣df(d)=x2∑i=1n/dx∑j=1n/dxijg(x) = \sum_{x|d}f(d)=x^2\sum_{i=1}^{n/dx}\sum_{j=1}^{n/dx}ijg(x)=∑x∣df(d)=x2∑i=1n/dx∑j=1n/dxij
故f(x)=∑x∣pμ(p/x)g(p)=∑x∣pμ(p/x)p2∑i=1n/dp∑j=1n/dpijf(x)=\sum_{x|p}\mu(p/x)g(p)=\sum_{x|p}\mu(p/x)p^2\sum_{i=1}^{n/dp}\sum_{j=1}^{n/dp}ijf(x)=∑x∣pμ(p/x)g(p)=∑x∣pμ(p/x)p2∑i=1n/dp∑j=1n/dpij
因此f(1)=∑p=1nμ(p)p2∑i=1n/dp∑j=1n/dpijf(1)=\sum_{p=1}^n\mu(p)p^2\sum_{i=1}^{n/dp}\sum_{j=1}^{n/dp}ijf(1)=∑p=1nμ(p)p2∑i=1n/dp∑j=1n/dpij
合并两部分
原式=∑d=1nd3∑p=1nμ(p)p2∑i=1n/dp∑j=1n/dpij\sum_{d=1}^nd^3\sum_{p=1}^n\mu(p)p^2\sum_{i=1}^{n/dp}\sum_{j=1}^{n/dp}ij∑d=1nd3∑p=1nμ(p)p2∑i=1n/dp∑j=1n/dpij
转而枚举a=dpa=dpa=dp.我们得到:
原式=∑a=1n∑d∣aa2dμ(a/d)∑i=1n/ai∑j=1n/aj=\sum_{a=1}^n\sum_{d|a}a^2d\mu(a/d)\sum_{i=1}^{n/a}i\sum_{j=1}^{n/a}j=∑a=1n∑d∣aa2dμ(a/d)∑i=1n/ai∑j=1n/aj
=∑a=1n(⌊na⌋(⌊na⌋+1)2)2a2∑d∣adμ(a/d)=\sum_{a=1}^n (\frac{\lfloor \frac{n}{a} \rfloor(\lfloor \frac{n}{a} \rfloor+1)}{2})^2a^2\sum_{d|a}d\mu(a/d)=∑a=1n(2⌊an⌋(⌊an⌋+1))2a2∑d∣adμ(a/d)
我们发现Id∗μ=ϕId*\mu=\phiId∗μ=ϕ是常见的狄利克雷卷积.
因此原式=∑a=1n(⌊na⌋(⌊na⌋+1)2)2a2ϕ(a)=\sum_{a=1}^n (\frac{\lfloor \frac{n}{a} \rfloor(\lfloor \frac{n}{a} \rfloor+1)}{2})^2a^2\phi(a)=∑a=1n(2⌊an⌋(⌊an⌋+1))2a2ϕ(a)
化简到这一步的时候,前面部分可以分块计算,后面的部分a2ϕ(a)a^2\phi(a)a2ϕ(a)要快速计算前缀和.
于是我们想到了杜教筛:
记f(x)=x2ϕ(x)f(x)=x^2\phi(x)f(x)=x2ϕ(x),找到积性函数g(x)=x2g(x)=x^2g(x)=x2与它做卷积使得g(x)g(x)g(x)的前缀和可以快速求出,而(f∗g)(x)=∑d∣xxd2×d2ϕ(d)=x2∑d∣xϕ(x)=x3(f*g)(x)=\sum_{d|x}\frac{x}{d}^2\times d^2\phi(d)=x^2\sum_{d|x}\phi(x)=x^3(f∗g)(x)=∑d∣xdx2×d2ϕ(d)=x2∑d∣xϕ(x)=x3的前缀和也可以快速求出.
记S(x)=∑i=1xf(x)S(x)=\sum_{i=1}^xf(x)S(x)=∑i=1xf(x)
根据杜教筛公式g(1)S(n)=∑x=1n(f∗g)(x)−∑x=2ng(x)S(nx)g(1)S(n)=\sum_{x=1}^n(f*g)(x) - \sum_{x=2}^ng(x)S(\frac{n}{x})g(1)S(n)=∑x=1n(f∗g)(x)−∑x=2ng(x)S(xn).
S(n)=∑x=1nx3−∑x=2nx2S(nx)S(n)=\sum_{x=1}^nx^3-\sum_{x=2}^nx^2S(\frac{n}{x})S(n)=∑x=1nx3−∑x=2nx2S(xn)
写到这就完了.

方法二.ϕ\phiϕ卷积

根据公式:

gcd(i,j)=∑d∣gcd(i,j)ϕ(d)=∑d∣i,d∣jϕ(d)gcd(i,j) = \sum_{d|gcd(i,j)}\phi(d)=\sum_{d|i,d|j}\phi(d)gcd(i,j)=∑d∣gcd(i,j)ϕ(d)=∑d∣i,d∣jϕ(d)

因此

∑i=1n∑j=1nijgcd(i,j)=∑i=1n∑j=1nij∑d∣i,d∣jϕ(d)=∑d=1nϕ(d)d2∑i=1n/d∑j=1n/dij=∑d=1nϕ(d)d2(⌊nd⌋(⌊nd⌋+1)2)2\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)=\sum_{i=1}^n\sum_{j=1}^nij\sum_{d|i,d|j}\phi(d)=\sum_{d=1}^n\phi(d)d^2\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij=\sum_{d=1}^n\phi(d)d^2(\frac{\lfloor \frac{n}{d} \rfloor(\lfloor \frac{n}{d} \rfloor+1)}{2})^2∑i=1n∑j=1nijgcd(i,j)=∑i=1n∑j=1nij∑d∣i,d∣jϕ(d)=∑d=1nϕ(d)d2∑i=1n/d∑j=1n/dij=∑d=1nϕ(d)d2(2⌊dn⌋(⌊dn⌋+1))2

同样也得到了我们的式子,是不是推导方法简单多了?

遇见gcd(i,j)gcd(i,j)gcd(i,j)的时候,多想想是不是可以用ϕ\phiϕ卷积来做?

可能会减少很多不必要的推导过程.

代码

#include <iostream>
#include <algorithm>
#include <cstring>
#include <unordered_map>
#define pr(x) std::cout << #x << ':' << x << std::endl
#define rep(i,a,b) for(LL i = a;i <= b;++i)
const int N = 1e7;
typedef long long LL;
LL n,p;
LL phi[N+10];
int prime[N+10],zhi[N+10],low[N+10],pcnt;
LL mod_pow(LL x,LL n) {LL res = 1;while(n) {if(n&1) res = res * x % p;x = x * x % p;n >>= 1;}return res;
}
LL inv6,inv2,inv4;
void sieve() {pcnt = 0;low[1] = phi[1] = zhi[1] = 1;rep(i,2,N) {if(!zhi[i]) {phi[i] = i-1;prime[pcnt++] = i;low[i] = i;}for(LL j = 0;j < pcnt && prime[j]*i <= N;++j) {zhi[i*prime[j]] = 1;if(i % prime[j] == 0) {low[i*prime[j]] = low[i] * prime[j];if(i == low[i]) {phi[i*prime[j]] = phi[i]*prime[j];}else {phi[i*prime[j]] = phi[i/low[i]] * phi[low[i]*prime[j]];}break;}else{low[i*prime[j]] = prime[j];phi[i*prime[j]] = phi[i] * phi[prime[j]];}}}
}
LL sum3(LL n) {n %= p;LL res = n*(n+1)%p*(2*n+1)%p*inv6%p;return res;
}
LL sum4(LL n) {n %= p;LL x = n*(n+1)%p;return x*x%p*inv4%p;
}
std::unordered_map<LL,LL> vis,rec;
LL F(LL n) {if(n <= N) return phi[n];if(vis[n]) return rec[n];LL res = sum4(n);for(LL x = 2,last;x <= n;x = last) {last = n/(n/x)+1;res = (res - ((sum3(last-1)-sum3(x-1)+p)%p*F(n/x)%p) + p) % p;}vis[n] = 1;return rec[n] = res;
}signed main() {sieve();std::ios::sync_with_stdio(false);std::cin >> p >> n;inv6 = mod_pow(6,p-2);inv4 = mod_pow(4,p-2);inv2 = mod_pow(2,p-2);rep(i,1,N) phi[i] = ((phi[i]*i%p*i%p) + phi[i-1]) % p;LL last,ans = 0;for(LL x = 1;x <= n;x = last+1) {last = n/(n/x);LL y = n/x%p;LL z = (1+y)*(y)/2%p;ans = (ans + (z*z%p*((F(last)-F(x-1)+p)%p))) % p;}std::cout << ans << std::endl;return 0;
}

P3768 简单的数学题 [狄利克雷卷积,杜教筛,莫比乌斯反演]相关推荐

【学习小记】狄利克雷卷积+杜教筛
Preface 这东西分明就是玄学暴力用来求简单的数论函数的前缀和,像φ,μφ,\mu这类的东西当然,约数和,约数个数之类的也是可以的 Text 数论函数是指定义域是整数,陪域是复数的函数 Dir ...
huntian oy (数论卷积杜教筛)
huntian oy (数论&卷积&杜教筛) 思路: f(n,a,b)=∑i=1n∑j=1igcd(ia−ja,ib−jb)[gcd(i,j)=1](mod109+7)f(n,a,b) ...
[学习笔记] 初次见面，请多关照 (公式推导+题集)——杜教筛
筛积性函数的前缀和常见积性函数公式推导狄利克雷卷积杜教筛实现常见卷积习题集 Sum 神犇和蒟蒻简单的数学题常见积性函数 μ\muμ φφφ d(n)d(n)d(n):nnn的约数个数 ...
P3768 简单的数学题（杜教筛）
P3768 简单的数学题推式子 ∑i=1n∑j=1mijgcd(i,j)=∑d=1nd∑i=1n∑j=1mij(gcd(i,j)=d)=∑d=1nd3∑i=1nd∑j=1ndij∑k∣gcd(i,j ...
（每日一题）P3768 简单的数学题（确信）（莫反 + 欧拉反演 + 杜教筛）
整理的算法模板合集: ACM模板点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划每日一题(莫反 / 多项式 / 母函数 / 群论) 2021.4.20 莫反和上午的那道题比较类似的 ...
#6229. 这是一道简单的数学题（反演 + 杜教筛）
#6229. 这是一道简单的数学题推式子 ∑i=1n∑j=1ilcm(i,j)gcd(i,j)=(∑i=1n∑j=1nlcm(i,j)gcd(i,j)+n)∗inv2所以重点求∑i=1n∑j=1nl ...
matlab狄利克雷函数,数论入门1——莫比乌斯函数,欧拉函数,狄利克雷卷积,线性筛,莫比乌斯反演,杜教筛...
数论入门1 一个菜鸡对数论的一点点理解... 莫比乌斯函数定义函数$\mu(n)$为: 当n有平方因子时,$\mu(n)=0$. 当n没有平方因子时,$\mu(n)=(-1)^{\omega(n)} ...
欧拉函数+狄利克雷卷积+莫比乌斯函数+莫比乌斯反演+整除分块+杜教筛
Powered by:NEFU AB-IN 文章目录欧拉函数狄利克雷卷积莫比乌斯函数莫比乌斯反演 P3455 [POI2007]ZAP-Queries 整除分块 P2522 [HAOI2011 ...
【洛谷3768】简单的数学题【莫比乌斯反演】【杜教筛】【小学奥数】
传送门题意:给定p,Np,Np,N,求 ∑i=1N∑j=1Nijgcd(i,j)modp\sum_{i=1}^{N}\sum_{j=1}^{N}ijgcd(i,j)\text{ }mod \text ...

P3768 简单的数学题 [狄利克雷卷积,杜教筛,莫比乌斯反演]