夏晓光(loveJY)的博客

仍要向下扎根,向上结果

首页 归档 标签 关于 友链

赫斯提亚的反演篇

2020-12-29 26 min read # ダンジョンに出会いを求めるのは間違っているだろうか # 莫比乌斯反演 # 杜教筛 # 数论

赫斯提亚!!/se

首先:

∑i=1ni=n∗(n+1)/2\sum_{i=1}^n i=n*(n+1)/2

∑i=1ni2=n∗(n+1)∗(2∗n+1)/6\sum_{i=1}^n i^2=n*(n+1)*(2*n+1)/6

∑i=1ni3=(n∗(n+1)/2)2\sum_{i=1}^n i^3=(n*(n+1)/2)^2

杜教筛

积性函数定义

f1=1f_1=1

fxy=fx∗fy,(x,y)=1f_{xy}=f_x*f_y,(x , y)=1

Si=∑i=1nfiS_i=\sum_{i=1}^n f_i

对于

∑i=1n∑d∣ifd∗gid\sum_{i=1}^n\sum_{d|i}f_d*g_{\frac{i}{d}}

∑d=1ngd∑i=1ndfi\sum_{d=1}^ng_{d}\sum_{i=1}^{\frac{n}{d}}f_i

∑d=1ngdS(⌊nd⌋)=∑(f∗g)\sum_{d=1}^n g_{d} S(\lfloor\frac{n}{d}\rfloor)=\sum(f*g)

∑(f∗g)−∑d=2ngdS(⌊nd⌋)=g1Sn\sum(f*g)-\sum_{d=2}^n g_{d} S(\lfloor\frac{n}{d}\rfloor)=g_1S_n

∑(f∗g)−∑d=2ngdS(⌊nd⌋)=Sn\sum(f*g)-\sum_{d=2}^n g_{d} S(\lfloor\frac{n}{d}\rfloor)=S_n

能够快速求出f*g的前缀和和g的前缀和即可快速求出f的前缀和

预处理前n2/3n^{2/3}可以实现O(n2/3)O(n^{2/3})复杂度

∑i=1nϕi,∑i=1nμi\sum_{i=1}^n \phi i,\sum_{i=1}^n \mu i

都有

gn=1g_n=1

可以加速求和gcd和lcm

∑i=1nϕi∗i\sum_{i=1}^n \phi i * i

gn=ng_n=n

可以做多快啊

卷起来之后那个常数就是n了,然后剩下的是phi可以直接等于n

code:

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MAXN = 4e6 + 7;
map<ll, ll> mp1, mp2;
int mu[MAXN], isp[MAXN], pri[MAXN];
ll phi[MAXN];
int N, T, tot;
inline void init() {
	N = 4e6;
	mu[1] = 1;
	phi[1] = 1;
	for(int i = 2; i <= N; ++i) {
		if(!isp[i]) {
			pri[++tot] = i;
			phi[i] = i - 1;
			mu[i] = -1;
		}
		for(int j = 1; j <= tot && i * pri[j] <= N; ++j) {
			isp[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
			phi[i * pri[j]] = phi[i] * (pri[j] - 1);
		}
	}
	for(int i = 1; i <= N; ++i) {
		mu[i] += mu[i - 1];
		phi[i] += phi[i - 1];
	}
}

inline ll solve1(ll x) {
	if(x <= N)return mu[x];
	if(mp1.find(x) != mp1.end())return mp1[x];
	ll ans = 1;
	for(ll l = 2, r; l <= x; l = r + 1) {
		r = x / (x / l);
		ans -= solve1(x / l) * (r - l + 1);
	}
	mp1[x] = ans;
	return ans;
}

inline ll solve2(ll x) {
	if(x <= N)return phi[x];
	if(mp2.find(x) != mp2.end())return mp2[x];
	ll ans = x * (x + 1) / 2;
	for(ll l = 2, r; l <= x; l = r + 1) {
		r = x / (x / l);
		ans -= solve2(x / l) * (r - l + 1);
	}
	mp2[x] = ans;
	return ans;
}

ll n;
int main() {
	scanf("%d", &T);
	init();
	while(T-- > 0) {
		scanf("%lld", &n);
		printf("%lld %lld\n", solve2(n), solve1(n));
	}
	return 0;
}

P3172 [CQOI2015]选数

问题显然可以转化为从n个数里面选数,一定要包括至少一次最小的数

那么这个问题可以容斥,我们可以钦点最小的数出现一次,然后减去他出现两次,加上他出现三次....

又因为每个集合是本质相同的,所以可以用组合数加速qwq

怎么求[L,H]中有多少个数是k的倍数啊

有⌊Hk⌋\lfloor\frac{H}{k}\rfloor个上界

⌊Lk⌋\lfloor\frac{L}{k}\rfloor

如果L不是k的倍数的话,否则是-1qwq

然后么得了,因为次数可能达到O(N)级别

考虑变成gcd容斥,就是fif_i表示最大公约数为i的倍数的方案数gig_i表示最大公约数为i的方案数

不难发现fi=gi+g2i+g3i....f_i=g_i+g_{2i}+g_{3i}....

不难发现gi=fi−g2i−g3i−...g_i=f_i-g_{2i}-g_{3i}-...

不难发现fn=gn=1f_n=g_n=1

所以说我们只需要求出一个数的最后gcd出现他的倍数的方案数

应该是这个数的倍数个数的n次幂

然后这个值fi=xnf_i=x^n难道就是了吗?

当然不是,因为我们是枚举了一个约数,这个约数只可能是H-L量级的,任何一个x大于H−LH-L,都不可能有两个以上的数使得他们的gcd等于x了(y和y+x一定超过了这个H,L的跨度吧)

然鹅如果我们所有数都是一样的话,这意味着我们的gcd一定会大于H-L这个范围,也就是说我们转换之后,能够保证操作序列随便给出最后的gcd都是H-L范围内的,但是如果所有数都一样就是例外,此时一定不会合法了,他们的gcd>H−Lgcd>H-L!!,所以就是要钦定

所以fi=xn−xf_i=x^n-x

于是我们就做完了这道题/yun

注意特判k是否在L,R内,此时我们要算上一次全部相同的qwq

code:

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 3e5 + 7;
const int P = 1e9 + 7;

ll f[MAXN];
ll n, k, l, h, m, ans;
inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void solve() {
	for(int i = m; i >= 1; --i) {
		ll t = h / i - l / i - (l % i ? 1 : 0) + 1;
		if(t <= 0)continue;
		f[i] = (ksm(t, n) - t + P) % P;
		// printf("%d %lld\n", t, f[i]);
		for(int j = i + i; j <= m; j += i) {
			f[i] += (P - f[j]) % P;
			f[i] %= P;
		}
	}
	printf("%lld\n", ans + f[1]);
	return;
}

int main() {
	scanf("%lld%lld%lld%lld", &n, &k, &l, &h);
	ans = (k >= l && k <= h);
	h /= k;
	l = l / k + (l % k ? 1 : 0);
	m = h - l;
	// printf("%lld %lld\n", l, h);
	// printf("%lld %lld\n", m, ans);
	solve();
	return 0;
}

杜教筛做法

(∑ai=LR)n[gcd(a1,a2,.....,an)==K](\sum_{a_i=L}^R)^n[gcd(a_1,a_2,.....,a_n)==K]

这个等价于直接枚举倍数,然鹅枚举倍数要小心转换下界哦

l′=L/k+[L%k],r′=R/kl'=L/k+[L\%k],r'=R/k

(∑ai=l′r′)n[gcd(a1,a2,.....,an)==1](\sum_{a_i=l'}^{r'})^n [gcd(a_1,a_2,.....,a_n)==1]

然后我们可以考虑反演

∑d=1r′(∑ai=l′r′)n[d∣ai]\sum_{d=1}^{r'}(\sum_{a_i=l'}^{r'})^n[d|a_i]

l′′=l′/d+[l′%d],r′′=r′/dl''=l'/d+[l' \% d],r''=r'/d

∑d=1r′μ(d)(∑ai=l′′r′′)n\sum_{d=1}^{r'}\mu (d)(\sum_{a_i=l''}^{r''})^n

∑d=1r′μ(d)(r′′−l′′+1)n\sum_{d=1}^{r'}\mu (d)(r''-l''+1)^n

杜教筛前面的部分,后面整除分块即可

这个整除分块有毒/吐血

code:


#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MAXN = 3e6 + 7;
const int P = 1e9 + 7;
map<ll, ll> mp1;
int mu[MAXN], isp[MAXN], pri[MAXN];
int N, tot;

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void init() {
	N = 3e6;
	mu[1] = 1;
	for(int i = 2; i <= N; ++i) {
		if(!isp[i]) {
			pri[++tot] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= tot && i * pri[j] <= N; ++j) {
			isp[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= N; ++i) {
		mu[i] += mu[i - 1];
	}
}

inline ll solve1(ll x) {
	if(x <= N)return mu[x];
	if(mp1.find(x) != mp1.end())return mp1[x];
	ll ans = 1;
	for(ll l = 2, r; l <= x; l = r + 1) {
		r = x / (x / l);
		ans -= solve1(x / l) * (r - l + 1);
	}
	mp1[x] = ans;
	return ans;
}


ll n, k, L, H, ans;
int main() {
	scanf("%lld%lld%lld%lld", &n, &k, &L, &H);
	init();
	L = L / k + (L % k ? 1 : 0) - 1;
	H = H / k;
	// printf("%lld %lld\n", L, H);
	ll ans = 0;
	for(int l = 1, r; l <= H; l = r + 1) {
		int tmpl = L / l + 1;
		int tmpr = H / l;
		// printf("%d %d %d %d\n", l, r, tmpr, tmpl);
		if(l <= L)r = min(H / tmpr, L / (L / l));
		else r = H / tmpr;
		// printf("%d %d %d %d\n", l, r, tmpr, tmpl);
		ans = ans + (solve1(r) - solve1(l - 1)) * ksm(tmpr - tmpl + 1, n) % P;
		ans %= P;
	}
	ans = (ans + P) % P;
	printf("%lld\n", ans);
	return 0;
}

51nod1220 约数之和

如果求约数个数和?

有一个可爱的sdoi2015反演原题

于是我们套用结论

∑i=1n∑j=1nd(ij)\sum_{i=1}^n\sum_{j=1}^nd(ij)

∑i=1n∑j=1n∑x∣i∑y∣i[gcd(x,y)==1]\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|i}[gcd(x,y)==1]

∑x=1n∑y=1n[gcd(x,y)==1]∑i=1nx∑j=1ny\sum_{x=1}^n\sum_{y=1}^n[gcd(x,y)==1]\sum_{i=1}^{\frac{n}{x}}\sum_{j=1}^{\frac{n}{y}}

∑d=1nμ(d)∑x=1nd∑y=1ndnxdnyd\sum_{d=1}^n\mu(d)\sum_{x=1}^{\frac{n}{d}}\sum_{y=1}^\frac{n}{d}\frac{n}{xd}\frac{n}{yd}

前面的我会杜教筛

后面的?其实直接做是O(n3/4)O(n^{3/4})的

但是如果我们预处理前O(n^{2/3),那么再做还是O(n2/3)O(n^{2/3}),因为本身就是一个整除分块的形式qwq

求约数和???

我的天哪这个也可以构造

di∗j=∑x∣i∑y∣jx∗j/y∗[gcd(x,y)==1]d_{i*j}=\sum_{x|i}\sum_{y|j}x*j/y*[gcd(x,y)==1]

原理其实理解了之前那个构造这个也能想懂

我们不妨只看一个约数的,多个约数互相独立也就一样

之前我们对于pa+bp^{a+b}计算方式是:单独x一个人计算a次,单独y一个人计算b次,然后两个加起来就计算了a+b次,两个都没有计算了一次

而这个相当于给每个计算方案分配一个权值:单独x一个人计算的是pi(i<a)+bp^{i(i<a)+b},单独y一个人计算的是pj(j<b)p^{j(j<b)},两个加起来就几乎都算了,两个都没有就是0次

于是剩下的就是

∑d=1nμ(d)d∑x=1ndnxd∗x∑y=1nyS(nyd)\sum_{d=1}^n\mu(d)d\sum_{x=1}^{\frac{n}{d}}\frac{n}{xd}*x\sum_{y=1}^{\frac{n}{y}}S(\frac{n}{yd})

其中S(x)=x∗(x+1)/2S(x)=x*(x+1)/2

你仔细品一品后面这两个东西都可以整除分块/se

于是直接上,就有O(n3/4)O(n^{3/4})优劣复杂度可以过2500ms

但为人不能这样,我们再一看,后面两个式子是一样的

∑d=1nμ(d)d(∑x=1ndS(nxd))2\sum_{d=1}^n\mu(d)d(\sum_{x=1}^{\frac{n}{d}}S(\frac{n}{xd}))^2

而!!再一看后面那个东西本质上就是约数个数和的前缀和!!

所以说我们可以对于前O(n2/3)O(n^{2/3})预处理出来

然后后面的直接整除分块做,复杂度还是O(n2/3)O(n^{2/3})的复杂度证明和杜教筛是一样的qwq

下面代码是O(n3/4)O(n^{3/4})的,因为懒惰

code:


#include<bits/stdc++.h>
const int MAXN = 2e6 + 7;
using namespace std;
#define ll long long
const int P = 1e9 + 7;
const int inv2 = (P + 1) / 2;
int N, isp[MAXN], pri[MAXN], tot;
ll mu[MAXN];
inline void add(int &x, ll y) {
	x += y;
	if(x >= P)x -= P;
}
inline void add(ll &x, ll y) {
	x += y;
	if(x >= P)x -= P;
}
inline void init() {
	N = 1e6;
	mu[1] = 1;
	for(int i = 2; i <= N; ++i) {
		if(!isp[i]) {
			pri[++tot] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= tot && i * pri[j] <= N; ++j) {
			isp[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= N; ++i) {
		mu[i] = (mu[i] * i % P + mu[i - 1] + P) % P;
		// printf("%lld ?\n", mu[i]);
	}
	// for(int m = 1; m <= N; ++m) {
	// 	for(int l = 1, r; l <= m; l = r + 1) {
	// 		r = m / (m / l);
	// 		add(f[m], 1ll * (r - l + 1) * (m / l) % P);
	// 	}
	// 	// if(m % 10000 == 0)printf("%d ?%lld\n", m, f[m]);
	// }
	return;
}

map<ll, ll>  mp;

inline ll solve1(ll x) {
	if(x <= N)return mu[x];
	if(mp.find(x) != mp.end())return mp[x];
	ll ans = 1;
	for(ll l = 2, r; l <= x; l = r + 1) {//222!!!
		r = x / (x / l);
		ll Len = (r * (r + 1) % P * inv2 % P - l * (l - 1) % P * inv2 % P + P) % P;
		add(ans, (P - solve1(x / l) * Len % P) % P);
	}
	mp[x] = ans;
	return ans;
}

inline ll solve2(ll n) {
	ll ans = 0;
	for(ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ll Len = (r * (r + 1) % P * inv2 % P - l * (l - 1) % P * inv2 % P + P) % P;
		// printf("%d %d %lld\n", l, r, Len);
		add(ans, Len * (n / l) % P);
	}
	// printf("%lld ?\n", ans);
	return ans;
}

inline ll solve3(ll n) {
	ll ans = 0;
	for(ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		add(ans, 1ll * (n / l) * ((n / l) + 1) % P * inv2 % P * (r - l + 1) % P);
	}
	// printf("%lld ?!\n", ans);
	return ans;
}


int n;
int main() {
	scanf("%d", &n);
	init();
	// puts("qwq");
	ll ans = 0;
	// printf("%lld\n", solve1(11));
	for(int l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		// printf("%d %d!!\n", l, r);
		ans += (solve1(r) - solve1(l - 1)) * solve2(n / l) % P * solve3(n / l) % P;
		ans %= P;
	}
	ans = (ans + P) % P;
	printf("%lld\n", ans);
	return 0;
}

P3700 [CQOI2017]小Q的表格

其实还是很妙的

b×f(a,a+b)=(a+b)×f(a,b)

盯....

f(a,b)(a∗b)=f(a,a+b)a∗(a+b)\frac{f(a,b)}{(a*b)} = \frac{f(a,a+b)}{a*(a+b)}

相当于欧几里得算法,辗转相减

f(a,b)a∗b=f(a,a%b)a∗(a%b)\frac{f(a,b)}{a*b}=\frac{f(a,a\%b)}{a*(a\%b)}

f(a,b)a∗b=f(gcd(a,b),gcd(a,b))gcd(a,b)2\frac{f(a,b)}{a*b}=\frac{f(gcd(a,b),gcd(a,b))}{gcd(a,b)^2}

f(a,b)=f(gcd(a,b),gcd(a,b))∗(a∗b)gcd(a,b)2\frac{f(a,b)}=\frac{f(gcd(a,b),gcd(a,b))*(a*b)}{gcd(a,b)^2}

坏事了

∑i=1n∑j=1nfi,j\sum_{i=1}^n\sum_{j=1}^nf_{i,j}

直接带入那个gcd,d=gcd

∑i=1n∑j=1ni∗j∗fd,d/d2\sum_{i=1}^n\sum_{j=1}^ni*j*f_{d,d}/d^2

∑dnfd,d∑i=1nd∑j=1ndi∗j[(i,j)==1]\sum_d^n f_{d,d} \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}} i*j[(i,j)==1]

∑dnfd,d∑t=1n/dμ(t)t2∑i=1ndt∑j=1ndti∗j\sum_d^n f_{d,d} \sum_{t=1}^{n/d} \mu(t)t^2 \sum_{i=1}^{\frac{n}{dt}}\sum_{j=1}^{\frac{n}{dt}} i*j

∑dnfd,d∑t=1n/dμ(t)S(n/dt)\sum_d^n f_{d,d} \sum_{t=1}^{n/d} \mu(t) S(n/dt)

S是约数的和

∑dnfd,dF(n/d)\sum_d^n f_{d,d} F(n/d)

F是F(N)=∑t=1Nμ(t)S(N/t)F(N)=\sum_{t=1}^{N} \mu(t) S(N/t)

gg了?带修!

修改的是每个f(d,d)f(d,d)qaq

然鹅这个查询一看就有mnm\sqrt n组

一看这个m就只有1e4

所以根号平衡一下就做 完 啦 !!

P4464 [国家集训队]JZPKIL

黑题之间的经典套路(雾)

∑i=1ngcd(i,n)xlcm(i,n)ymod  (1e9+7)\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y \mod (1e9+7)

你随意反演一下,注意没有两个西格玛意味着d直接是n的约数

∑d∣ndx∑i=1nd(iy)ny[gcd(i,n/d)==1]\sum_{d|n}d^x\sum_{i=1}^{\frac{n}{d}}(i^y)n^y[gcd(i,n/d)==1]

这里注意一下我们不能随便替换n/d和d,因为这意味着我们式子中其他地方涉及到的都要改变

∑d∣ndx∑t∣ndμ(t)∑i=1ndt(it)yny\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)\sum_{i=1}^{\frac{n}{dt}}(it)^yn^y

∑d∣ndx∑t∣ndμ(t)ty∗ny∑i=1ndtiy\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y*n^y\sum_{i=1}^{\frac{n}{dt}}i^y

注意到后面是自然幂数和,所以我们直接自闭,套用一个人的数论里面的结论有

∑d∣ndx∑t∣ndμ(t)ty∗ny∑i=0y+1(ndt)ifi\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y*n^y\sum_{i=0}^{y+1}(\frac{n}{dt})^if_i

提前

ny∗∑i=0y+1fi∑d∣ndx∑t∣ndμ(t)ty(ndt)in^y*\sum_{i=0}^{y+1}f_i\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y(\frac{n}{dt})^i

然 后 开 始 神 乎 其 技

我们注意到,后面那一坨东西:(忽略常数x,y,i)

∑d∣ndx∑t∣ndμ(t)ty(ndt)i\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y(\frac{n}{dt})^i

这个东西就是直接是n的一个函数啊!设为Z

所以说就是说这个东西可以观察得到(别问我为啥可以)

fn=nx,gx=μ(x)∗xy,hxi=nif_n=n^x,g_x=\mu(x)*x^y,h^i_x=n^i

Z=f∗g∗hZ=f*g*h

!!!也就是说他是积性函数

对于积性函数,单点求值最好的解决方法就是分成所有质因数然后把质因数卷起来(或者我只会这个)

打表法

Z(pk)Z(p^k)

  1. k=0

Z1=1Z1=1

  1. k=1

Zp=pi+px−pyZp=p^i+p^x-p^y

因为此时相当于随意分配一下

  1. k>=2

你观察到我们有多余两个指数分给中间那一项就可以变成0了

所以我们中间那个要么分0要么分1

对应了:

f∗hi(pk)−f∗hi(pk−1)∗pyf*h^i(p^k)-f*h^i(p^{k-1})*p^y

考虑卷积是啥

f∗h=∑l+j=k(plx∗pji)f*h=\sum_{l+j=k}(p^{l^x}*p^{j^i})

真不错,所以求值就是把所有这样的计算一下然后加起来即可

综上,我们还差一步求质因数,这一步可以使用MR暴力质因数分解,然后暴力算次幂qwq

自然幂数和多项式?

正当您信心满满的想去写高消的时候,发现y<=2000y<=2000

说明我们要伯努利数或者斯特林数了,插值就要实现多项式除法还毒瘤qaq

待填坑

Powered by Gridea RSS