生成树总结

生成树

一个无向图的生成树指的是从图中选若干边构成边集,全部点构成点集,使得这个边集加上点集恰好是一棵树。

生成树计数

一个无向无权图(允许重边不允许自环)的邻接矩阵为g,显然这是一个对称矩阵,g[u][v]代表边(u,v)的重数,即若存在一条边(u,v)则g[u][v]的值为1,若存在k条,则g[u][v]的值为k。
一个无向无权图(允许重边不允许自环)的度数矩阵为deg,显然这是一个对角矩阵,deg[u][u]代表点u的度数。
一个无向无权图(允许重边不允许自环)的基尔霍夫矩阵(拉普拉斯矩阵)为hoff,是度数矩阵减去邻接矩阵。
矩阵树定理说一个无向图的生成树的个数刚好等于基尔霍夫矩阵的行列式的任意n-1阶主子式(代数余子式)的行列式的绝对值。
生成树计数复杂度$O(V^3+E)=O(V^3)$
黑暗前的幻想乡
我们利用矩阵树定理就能轻松解决

黑暗前的幻想乡代码

最小生成树

有一个无向带权图,每条边有权$x_i$,需要求出一个生成树T,并最小化$\begin{aligned}\sum_{i\in T}x_i\end{aligned}$
kruskal算法:贪心从小到大枚举边合并森林即可。这里就不证明此算法了。
最小生成树复杂度$O(V+ElgE)=O(ElgE)$
最小生成树

最小生成树代码

最小生成树计数

由于最小生成树各自边权构成的多重集合是一样的,并且易证不同的边权对最小生成树的影响是独立的,所以我们可以通过将边按照边权分类,分别求出每一种边权各自对联通块的贡献,然后利用计数的乘法原理合并即可。我们需要先求出任意一个最小生成树,当我们对某一种边权进行讨论的时候,我们需要将这个生成树中这一边权的边全删掉,然后对剩余联通块进行缩点并重新编号,将待选集合中的边映射到联通块间的边,并去掉自环。这样以后待选集合中的边的边权就相等了。这时我们可以借助矩阵树定理来求解。
最小生成树计数复杂度$O(ElgE+V^3)=O(V^3)$
最小生成树计数

最小生成树计数代码

严格次小生成树

严格次小生成树和最小生成树的边权多重集合中只有一个边权不一样,这样我们就有了一个简单的算法,先求出任意一个最小生成树,然后枚举没有被选中为构建最小生成树的边,假设这个边为$(u,v,w_1)$,我们在最小生成树上求出点$u$到点$v$这条路径上的最大边权$w_2$和严格次大边权$w_3$,若$w_1=w_2$则我们用此边换掉次大边,若$w_1>w_2$则我们用此边换掉最大边,这样我们最终能得到很多非最小生成树,从中选出最小的那个,他就是次小生成树,这个过程需要维护树上的路径信息,有算法都可以树剖、树上倍增、lct等等,我们这里使用树上倍增的办法来解决。
严格次小生成树时间复杂度$O(ElgE+ElnV)=O(ElgE)$
严格次小生成树

严格次小生成树代码

最小乘积生成树

有一个无向带权图(权为正数),每条边有权$x_i$和权$y_i$,需要求出一个生成树T,记$\begin{aligned}X=\sum_{i\in T}x_i,Y=\sum_{i\in T}y_i\end{aligned}$,要求最小化乘积$XY$
我们假设已经求出了所有生成树,他们的权为$(X_i,Y_i)$我们把这个二元组看做二维平面上的点,则最小乘积生成树一定在凸包上。进一步分析,整个凸包都在第一象限,那么我们可以锁定出两个点了,他们一定在凸包上。分别求出最小的$X_i$对映点$A$,和最小的$Y_i$对映点$B$,那么答案就在$AB$左下方,我们求出一个点$C$,若能让叉积$AC*AB$最大且为正数,则$C$一定也在凸包上。我们递归处理$AC$和$CB$即可。凸包上的点期望为lg级别。
最小乘积生成树复杂度$O(ElgElg(V!))=O(ElgElgV)$
最小乘积生成树

最小乘积生成树代码

最小差值生成树

有一个无向带权图,每条边有权$x_i$,需要求出一个生成树T,让T的最大边权和最小边权的差值尽可能小。
我们对边集排序后,等价于求最短的一段区间,这个区间内部的边能够生成一棵树,这种连通性维护的问题,直接使用lct就可以了,
最小差值生成树时间复杂度$O(ElgE+Elg(E+V))=O(ElgE)$
最小差值生成树

最小差值生成树代码

k度限制最小生成树

在最小生成树的要求下,多一个条件: 有一个定点的度数不能超过k。
k度限制最小生成树与k-1度限制最小生成树最多有一条边的区别。
时间复杂度$O(ElgE+kV)$
k度限制最小生成树

k度限制最小生成树代码

最小直径生成树

给无向连通图,求一个直径最小的生成树。
以图的绝对中心为根的最短路树,是一个最小直径生成树。先用floyd求多源最短路,然后对每一条边,假设绝对中心在此边上,求出该绝对中心的偏心率,可以考虑从大到小枚举最短路闭包来实现,汇总得到绝对中心,最终以绝对中心为根,求最短路树。
时间复杂度$O(n^3)$

最小直径生成树代码

最小比值生成树

有一个无向带权图(权为正数),每条边有权$x_i$和权$y_i$,需要求出一个生成树T,记$\begin{aligned}X=\sum_{i\in T}x_i,Y=\sum_{i\in T}y_i\end{aligned}$,要求最小化比值$\frac{X}{Y}$.
我们设$r=\frac{X}{Y}$则有$rY-X=0$我们定义函数$f(r)=rY-X$,为当以$ry_i-x_i$作为权值的时候的最大生成树的值,这里显然f(r)关于r单调增,当我们找到一个r使得f(r)等于0的时候,r就是我们分数规划要的答案。
时间复杂度$O(lgn)*O(MST)$

最小比值生成树代码

hdu6607

name

Easy Math Problem

descirption

One day, Touma Kazusa encountered a easy math problem. Given n and k, she need to calculate the following sum modulo $1e9+7$.
$$∑_{i=1}^n∑^n_{j=1}gcd(i,j)^klcm(i,j)[gcd(i,j)∈prime]%(1e9+7) $$
However, as a poor student, Kazusa obviously did not, so Touma Kazusa went to ask Kitahara Haruki. But Kitahara Haruki is too busy, in order to prove that he is a skilled man, so he threw this problem to you. Can you answer this easy math problem quickly?

input

There are multiple test cases.$(T=5)$ The first line of the input contains an integer$T$, indicating the number of test cases. For each test case:
There are only two positive integers n and k which are separated by spaces.
$1≤n≤1e10$
$1≤k≤100$

output

An integer representing your answer.

sample input

1
10 2

sample output

2829

toturial

$$
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j)^{k-1} gcd is prime
\=&\sum_{d\in prime} \sum_{i=1}^n\sum_{j=1}^nijd^{k-1}[gcd(i,j)=d]
\=&\sum_{d\in prime} \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^{k+1}[gcd(i,j)=1]
\=&\sum_{d\in prime}d^{k+1} \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[gcd(i,j)=1]
\=&\sum_{d\in prime}d^{k+1} \sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}j^2\phi(j)
\end{aligned}
$$
我们可以对n分块了,前面可以min25筛
$\begin{aligned}f(j)=j^2\phi(j)\end{aligned}$
$\begin{aligned}g(j)=j^2\end{aligned}$
$\begin{aligned}f\ast g(j)=\sum_{i|j}i^2\phi(i)(\frac{j}{i})^2=j^2\sum_{i|j}\phi(i)=j^2(\phi\ast 1)(j)=j^3\end{aligned}$
于是后面可以杜教筛

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

// 模意义
const ll mod=1e9+7;
ll qpow(ll a,ll b){
assert(a<mod);
ll res=1;
while(b){
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}return res;
}
const ll inv2=qpow(2,mod-2),inv3=qpow(3,mod-2);
inline ll reduce(ll x){return x<0?x+mod:x;}
inline ll A(ll a,ll b){assert(a<mod&&b<mod); return reduce(a+b-mod);}
inline ll M(ll a,ll b){assert(a<2*mod&&b<2*mod); return a*b%mod;}
inline ll M(ll a,ll b,ll c){return M(M(a,b),c);}

//线性筛
// 3e7 int = 120mb
const ll maxn=2.5e7;
bitset<maxn>vis;
int siiphi[maxn];
ll p[1565927+100];
void f_ini(){
siiphi[1]=1;
for (ll i=2;i<maxn;i++){
if(!vis[i]) p[++p[0]]=i,siiphi[i]=i-1;
for (ll j=1;i*p[j]<maxn;j++){
vis[i*p[j]]=true;
if(i%p[j])siiphi[i*p[j]]=siiphi[i]*(p[j]-1);//由积性函数性质推
else{siiphi[i*p[j]]=siiphi[i]*p[j];break;}
}
}
for(ll i=1;i<maxn;i++) siiphi[i]=A(siiphi[i-1],M(i,i,siiphi[i]));
}

// 分块
const ll sqr=3e5;
ll id1[sqr],id2[sqr],w[sqr],idn,idm;// w[x] 第几大的分块值是多少
inline ll& id(ll x){return x<sqr?id1[x]:id2[idn/x];}//返回x是第几大的整除分块值
void ini(ll n){
idn=n;idm=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
id(n/l)=++idm;
w[idm]=n/l;
}
}

namespace min25shai{
ll g[sqr],sp[sqr];
ll getsum(ll x,ll n){// O(n) n次多项式有n+1项 y[0]...y[n] -> y[x]
static ll prepre[1000],suf[1000],r[1000]={1,1},y[1000],*pre=prepre+1;
if(y[999]!=++n) {//这里非常重要
y[999]=n;
for(ll i=1;i<=n;i++) y[i]=A(y[i-1],qpow(i,n-1));
for(ll i=2;i<=n;i++) r[i]=M(mod-mod/i,r[mod%i]);
for(ll i=2;i<=n;i++) r[i]=M(r[i],r[i-1]);
}
pre[-1]=suf[n+1]=1;
for(ll i=0;i<=n;++i) pre[i]=M(pre[i-1],x%mod-i+mod);//这个地方爆掉了
for(ll i=n;i>=0;i--) suf[i]=M(suf[i+1],i-x%mod+mod);//这个地方爆掉了
ll b=0;
for(ll i=0;i<=n;++i) {
ll up=M(pre[i-1],suf[i+1]);
ll down=M(r[i],r[n-i]);
b=A(b,M(y[i],up,down));
}
return b;
}
void min25(ll*g,ll n,ll k,ll(*f)(ll,ll),ll(*s)(ll,ll)){
for(ll i=1;i<=idm;++i) g[i]=A(s(w[i],k),mod-1);
for(ll j=1;p[j]*p[j]<=n;j++){
ll t=f(p[j],k);
sp[j]=A(sp[j-1],t);
for(ll i=1;w[i]>=p[j]*p[j];++i) g[i]=A(g[i],M(sp[j-1]-g[id(w[i]/p[j])]+mod,t));
// w[i]从大到小 当i等于m的时候 w[i]>=p[j]*p[j]恒不成立
}
}
}

namespace dujiaoshai{
// g(1)S(n)=(1≤i≤n)h(i)+(2≤d≤n)g(d)S(n/d)
// f(n)=n*n*phi(n)
// g(n)=n*n
// h(n)=n*n*n
ll s[sqr];// 前缀和
inline ll s1(ll n){return M(n,n+1,inv2);}
inline ll s2(ll n){return M(s1(n),2*n+1,inv3);}
inline ll s3(ll n){return M(s1(n),s1(n));}
void ini(){for(ll i=1;i<=idm;i++)s[i]=0;}
ll dujiao(ll n){
if(n<maxn) return siiphi[n];
if(s[id(n)]!=0) return s[id(n)];
s[id(n)]=s3(n%mod);
for(ll l=2,r;l<=n;l=r+1){
r=n/(n/l);
s[id(n)]-=(s2(r%mod)-s2((l-1)%mod))*dujiao(n/l)%mod;
}
return s[id(n)]=(s[id(n)]%mod+mod)%mod;
}
}

ll solve(ll n,ll k){
ini(n);
dujiaoshai::ini();
#define F(M) [](ll n,ll k){return ll(M);}
min25shai::min25(min25shai::g,n,k+1,F(qpow(n%mod,k)),F(min25shai::getsum(n,k)));
#undef F
ll res=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
ll t1=dujiaoshai::dujiao(n/l);
ll t2=min25shai::g[id(r)];
if(l!=1) t2+=mod-min25shai::g[id(l-1)];
res+=M(t1,t2);
}
return res%mod;
}

inline ll read(){ll x;cin>>x;return x;}
int main() {
f_ini();
for(ll t=read();t>=1;t--){
ll n=read(),k=read();
cout<<solve(n,k)<<endl;
}
}

min25筛

min25筛是什么

min25筛是一种筛法,他能以亚线性的时间复杂度筛出一类函数的前缀和

定义一部分符号

$M(x) x\gt1$代表$x$的最小质因子

我们再设$P_j$为第$j$小的质数, $P_1=2,P_2=3,P_3=5…$

先看最简单的第一类函数

$$
\begin{aligned}
f(x)=\left{\begin{matrix}
x^k&x\in primes\
0&x \notin primes
\end{matrix}\right.
\end{aligned}
$$
对于这个函数我们可以利用min25筛来达到$O(\frac{n^\frac{3}{4}}{lg(n)})$的时间复杂度,我们没有办法直接求这个函数的前缀和,但是我们可以另外设一个相对好求的函数$h(x)=x^k$,通过h来求f,因为$\begin{aligned}\sum_{i=2}^nh(i)[i\in primes]=\sum_{i=2}^nf(i)[i\in primes]\end{aligned}$

设$\begin{aligned}
g(n,j)=\sum_{i=2}^nh(i)[i \in primes或M(i)\gt P_j]
\end{aligned}$

即 i要么是质数,要么i的最小质因子大于$P_j$。对g函数的理解,我们甚至可以回忆埃式筛,每一轮都会选择一个最小的质数,然后筛掉他的所有倍数,最终唯有所有的质数不会被筛掉,我们的这个函数就是那些没有被筛掉的数的函数值的和。
$$
\begin{aligned}
g(n,j)=\left{\begin{matrix}
g(n,j-1)-x&M(n)\le P_j\
g(n,j-1)& M(n)\gt P_j
\end{matrix}\right.
\end{aligned}
$$
x处是什么呢?第j-1次的结果和第j次的结果有什么不同呢?第j次埃式筛筛掉了$P_j$的倍数,他们的最小质因子都是$P_j$,所以
$$
\begin{aligned}
x&=\sum_{i=2P_j}^nh(i)[M(i)=P_j]
\&=\sum_{i=2}^{\frac{n}{P_j}}h(iP_j)[M(iP_j)=P_j]
\&=h(P_j)\sum_{i=2}^{\frac{n}{P_j}}h(i)[M(i)\ge P_j]
\&=h(P_j)\sum_{i=2}^{\frac{n}{P_j}}h(i)[M(i)\gt P_{j-1}]
\&=h(P_j)(\sum_{i=2}^{\frac{n}{P_j}}h(i)[M(i)\gt P_{j-1}或i \in primes]-\sum_{i=1}^{j-1}h(P_i))
\&=h(P_j)(g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}h(P_i))
\end{aligned}
$$

最后就成了这个
$$
\begin{aligned}
g(n,j)=\left{\begin{matrix}
g(n,j-1)-h(P_j)(g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}h(P_i))&M(n)\le P_j\
g(n,j-1)& M(n)\gt P_j
\end{matrix}\right.
\end{aligned}
$$
到这里就已经可以记忆化递归解决了,但是递归比较慢,我们考虑把它变成非递归,我们观察这个式子。

我们发现我们可以先枚举j因为$g(n,j)$是由$g(?,j-1)$推导过来的,然后从大到小枚举n,来更新数组,又因为n的前一项可能与$\frac{n}{P_j}$有关,所以我们可以把他们都用map映射一下,再进一步分析,根据整除分块的传递性,$\frac{\frac{a}{b}}{c}=\frac{a}{bc}$我们可以得出所有$g(x,y)$中x构成的集合,恰好是集合${x|x=\frac{n}{t},t\in [1,n]}$,最后预处理一下$\sum^{j-1}_{i=1}h(P_i)$即可,对于整除分块的映射,我们又可以加速到O(1)此处不做过多分析。

最后我们得到了这个$O(\frac{n^{\frac{3}{4}}}{lg(n)})$算法

再看复杂一些的第二类函数

第二类函数是抽象的积性函数$f$。

如果我们能够通过一些方法求出$\sum_{i=1}^{n}f(P_i)$和$f(P_i^k)$,那么我们就能够很简单得推出f的前缀和。

对于$\sum_{i=1}^{n}f(P_i)$, 我们这样来求,比如说f(x)在x是一个质数的时候能表示为某个简单多项式,那么我们就可以通过将多项式函数看做某些幂函数的线形组合,先求出幂函数各自的质数前缀和,然后相加就可以得到f的质数前缀和。

而对于另外一个$f(P_i^k)$则需要通过函数的定义来求了。

现在假设我们已经预处理出了$\sum_{i=1}^xf(P_i)(x \in n的数论分块即x=\frac{n}{?})其实就是g(x,\infty)$。

我们设$\begin{aligned}S(n,j)=\sum_{i=2}^nf(i)[M(i)\ge P_j]\end{aligned}$注意和$g(n,j)$对比。
$$
\begin{aligned}
&S(n,j)
\=&\sum_{i=j}^{P_i\le n}f(P_i)+f(P_i)S(\frac{n}{P_j},j+1)+f(P_i^2)S(\frac{n}{P_i^2},j+1)+f(P_i^3)S(\frac{n}{P_i^3},j+1)+…
\end{aligned}
$$
这里已经可以了,第一项可以通过两个前缀和相减得到,后边的递归。这就是min25筛的灵魂所在。

我们现在好好来分析一下这个递归式子。我们发现第一项是最好求的,就是第一类函数,但是后边的几项必须要求积性函数。这也是min25筛只能对积性函数起作用的地方。

min25筛能处理更多的函数吗?

我们暂定这些函数为f,显然我们必须要能够求出g和s,这就是min25筛,对于g,这里不对此作过多分析,没有这个必要,我们假定都是一类与幂函数线形组合有关的函数,抑或是某几项与幂函数有关,反正只要能够找到完全积性函数h在质数自变量和f函数存在相等关系即可。s的话,第一项简单差分,后边的看似要求f是积性函数,其实不然,我们仔细分析,其实他要求的是这样的要求: 假定y是x的最小质因子,$z=y^k且z|x且k最大$,我们只要能够通过$f(z)和f(\frac{x}{z})$这两项能够推出f(x)即可,这里并没有强制要求$f(x)=f(z)*f(\frac{x}{z})即f(x)=f(M(x))$。举个例子,若$f(x)=f(z)=f(y)=y$,我们也是可以求的。

贴一个求$f(a^b)=a \bigotimes b$和$f(x)=M(x)$的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const ll sqr=2e5+10;// 2sqrt(n)
ll p[sqr],np[sqr]={1,1};
void prime_ini(){// 素数不可能筛到longlong范围
for(int i=2; i<sqr; i++){
if(!np[i])p[++p[0]]=i;//把质数收录
for(int j=1; 1ll*i*p[j]<sqr; j++){
np[i*p[j]]=1;//对每一个数字枚举他的最小因子
if(i%p[j]==0)break;//在往后的话就不是最小因子了
}
}
}

const ll mod=1e9+7;
ll w[sqr],g[sqr][2],sp[sqr][2],id1[sqr],id2[sqr],mpn;
inline ll& mp(ll x){return x<sqr?id1[x]:id2[mpn/x];}
void min25(ll n){// 计算质数位置之和的时候 只用到了f(x,1) 和 oddsum(x)
mpn=n;
ll m=0;
for(ll l=1,r;l<=n;l=r+1){// i从小到大 n/i从到小
r=n/(n/l);
mp(n/l)=++m;
w[m]=n/l;
g[m][0]=(w[m]-1)%mod;// f(x)=1, s(x)=x
g[m][1]=(__int128(w[m])*(w[m]+1)/2-1)%mod; // f(x)=x, s(x)=x(x+1)/2 这里的int128非常关键,因为n是1e10级别的
}//assert(w[m]==1);
for(ll j=1;p[j]*p[j]<=n;j++){
sp[j][0]=sp[j-1][0]+1;// f(x)=1
sp[j][1]=(sp[j-1][1]+p[j])%mod;// f(x)=x
for(ll i=1;w[i]>=p[j]*p[j];++i){// w[i]从大到小 当i等于m的时候 w[i]>=p[j]*p[j]恒不成立
g[i][0]-=(g[mp(w[i]/p[j])][0]-sp[j-1][0])*1;// f(x)=1
g[i][1]-=(g[mp(w[i]/p[j])][1]-sp[j-1][1])*p[j];// f(x)=x
g[i][0]=(g[i][0]%mod+mod)%mod;
g[i][1]=(g[i][1]%mod+mod)%mod;
}
}
}

// f(pow(a,b))=a^b f为积性函数
inline ll f(ll a,ll b){return a^b;} // 当且仅当a是一个素数
ll s(ll n,ll j){// sum of f(x) x<=n minfac(x)>=p[j]
ll res=(g[mp(n)][1]-g[mp(n)][0])-(sp[j-1][1]-sp[j-1][0])+2*mod;// 减掉p[j]前面的质数 : [p[j],n]上的质数的函数的和
if(n>=2&&j==1) res+=2;
for(ll k=j;p[k]*p[k]<=n;k++){// 枚举的最小质因子
for(ll x=p[k],e=1;x*p[k]<=n;x*=p[k],e++){//枚举该因子出现次数
res+=s(n/x,k+1)*f(p[k],e)%mod+f(p[k],e+1);// 每次增加2mod res不可能超过 long long
}
}
return res%mod;
}

// f(x)=minfac(x) f不为积性函数 但我们用积性函数来做他
typedef pair<ll,ll> pll;
pll s2(ll n,ll j){//
ll res1=g[mp(n)][0]-sp[j-1][0]+2*mod;
ll res2=g[mp(n)][1]-sp[j-1][1]+2*mod;// 减掉p[j]前面的质数 : [p[j],n]上的质数的函数的和
for(ll k=j;p[k]*p[k]<=n;k++){// 枚举的最小质因子
for(ll x=p[k],e=1;x*p[k]<=n;x*=p[k],e++){//枚举该因子出现次数
pll tmp=s2(n/x,k+1);
res1+=tmp.first*1%mod+1;
res2+=tmp.first*p[k]%mod+p[k];// 每次增加2mod res不可能超过 long long
}
}
return pll(res1%mod,res2%mod);
}

int main() {
prime_ini();
ll n;
while(cin>>n){
min25(n);
if(n==1) cout<<1<<endl;
else cout<<(s(n,1)+1)%mod<<endl;
}
}

常用数论函数变换

中间结论

$$
f_n = \sum_{i=0}^n (-1)^i {n \choose i} g_i
\Leftrightarrow
g_n = \sum_{i=0}^n (-1)^i {n \choose i} f_i
$$

$$
f_n = \sum_{i=0}^n {n \choose i} g_i
\Leftrightarrow
g_n = \sum_{i=0}^n (-1)^{n-i} {n \choose i} f_i
$$

杜教筛

$$
g(1)\sum_{i=1}^nf(i)=\sum_{i=1}^{n}(f*g)(i)-\sum_{d=2}^{n}g(d) \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)
$$

数论变换

$$
\begin{aligned}
&e(n)=[n=1]
\&id(n)=n
\&I(n)=1
\&d(n)=\sum_{x|n}1 (就是因子个数)
\&σ(n)=\sum_{x|n}x (就是因子和)
\&\mu(n) 莫比乌斯函数
\&\phi(n) 欧拉函数
\&I\ast I=d
\&I\ast id=σ
\&I\ast \phi=id
\&I\ast \mu=e
\end{aligned}
$$

自变量互质的前缀和函数分析

$$
\begin{aligned}
&\sum_{i=1}^n{f(i)[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i*d)}\
&\sum_{i=1}^n{[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)\lfloor\frac{n}{d}\rfloor}\
&\sum_{i=1}^n{[gcd(i,n)=1]}=\phi(n)\
&\sum_{i=1}^n{i[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)d\frac{\lfloor\frac{n}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)}{2}}\
&\sum_{i=1}^n{i[gcd(i,n)=1]}=\frac{n}{2}(\phi(n)+e(n))\
&\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=1]=\sum_{j=1}^nj^2\phi(j)\
\end{aligned}
$$

牛客挑战赛33D

###name
种花家的零食

###descirption
在很久以前,有一颗蓝星,蓝星上有一个种花家。
种花家有1到n共n包零食,同时种花家的兔子又有1到n共n个朋友(比如毛熊,鹰酱,脚盆鸡等)。
昨天,兔子的n个朋友都到他家来玩了。他的n个朋友瓜分了他的n包零食,每个人都恰好吃了一包零食,没有两个人吃了同一包零食。
兔子发现,第i个朋友吃第j包零食能获得的愉悦值是$i\mod j$。
今天,兔子想回忆起每个朋友吃的是哪包零食,他想不起来了,但是他却记得了所有人的愉悦值之和s。于是,兔子找上了你,请你构造出一种可能的方案。
由于兔子记忆力不好,他有可能记错了s,所以可能会存在无解的情况。

###input
一行两个整数$n(1\leq n\leq 10^6)$和$s(1\leq s\leq10^{18})$

###output
如果不存在满足条件的方案,输出一行-1。
否则输出n行,每行一个整数,第i行的整数表示第i个朋友吃的是哪包零食。

###sample input
5 7

###sample output
1
4
3
5
2

###sample input
5 100

###sample output
-1

###toturial
分析出上界为$\frac{n(n-1)}{2}$后,分类讨论,用数学归纳法证明特例即可

###code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include<bits/stdc++.h>
#define rep(i,j,k) for(ll i=(j);i<=(k);++i)

using namespace std;
typedef long long ll;

const int maxn=1e6+10;
ll ans[maxn];
vector<ll>vec;
void solve(ll n,ll s){
if(n==1) {
vec.push_back(1);
return;
}
// choose n-1
s-=n-1;
n--;
if(s!=n*(n-1)/2-1&&s!=2&&s!=0&&s>0){
vec.push_back(n);
solve(n,s);
return;
}
n++;
s+=n-1;


solve(n-1,s);
}

int main() {
ll n,s;
cin>>n>>s;
if(s>n*(n-1)/2) ans[0]=-1;
else if(n==1){
if(s==0) ans[1]=1;
else ans[0]=-1;
}
else if(n==2){
if(s==0) ans[1]=1,ans[2]=2;
else if(s==1) ans[1]=2,ans[2]=1;
else ans[0]=-1;
}
else if(s==0) rep(i,1,n) ans[i]=i;
else if(s==2) {
ans[1]=3;;
ans[2]=1;
ans[3]=2;
rep(i,4,n) ans[i]=i;
}
else if(s==n*(n-1)/2-1){
if(n%2==0){
ans[1]=1;
rep(i,2,n-1) ans[i]=i+1;
ans[n]=2;
}
else{
ans[1]=3;
ans[2]=1;
rep(i,3,n-1) ans[i]=i+1;
ans[n]=2;
}
}
else {
vec.push_back(n);
solve(n,s);
rep(i,1,n) ans[i]=i;
int sz=vec.size();
// rep(i,0,sz-1) cout<<vec[i]<<endl;
// cout<<endl;
rep(i,0,sz-1)ans[vec[i]]=vec[(i-1+sz)%sz];
}
if(ans[0]==-1) printf("-1\n");
else {
ll ss=0;
rep(i,1,n) ss+=i%ans[i],printf("%lld\n",ans[i]);
assert(ss==s);
// cout<<s<<endl;
}
}

自变量互质的前缀和函数分析

有这样一类问题,他们的形式常常是这个样子

$$
\begin{aligned}
\sum_{i=1}^n{f(i)[gcd(i,j)=1]}
\end{aligned}
$$

我们来对他进行变形

$$
\begin{aligned}
&\sum_{i=1}^n{f(i)[gcd(i,j)=1]}\
=&\sum_{i=1}^n{f(i)e(gcd(i,j))}\
=&\sum_{i=1}^n{f(i)(\mu1)(gcd(i,j)}\
=&\sum_{i=1}^n{f(i)\sum_{d|gcd(i,j)}\mu(d)}\
=&\sum_{i=1}^n{f(i)\sum_{d|i,d|j}\mu(d)}\
=&\sum_{d|j}{\mu(d)\sum_{d|i,1<=i<=n}f(i)}\
=&\sum_{d|j}{\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i
d)}\
\end{aligned}
$$

如果$f(i)=1$ 则

$$
\begin{aligned}
\sum_{i=1}^n{[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)\lfloor\frac{n}{d}\rfloor}\
\end{aligned}
$$

更加特殊的 如果$j=n$ 则

$$
\begin{aligned}
\sum_{i=1}^n{[gcd(i,n)=1]}=\sum_{d|j}{\mu(d)\frac{n}{d}}=(\mu*id)(n)=\phi(n)\
\end{aligned}
$$

如果$f(i)=i$ 则

$$
\begin{aligned}
&\sum_{i=1}^n{i[gcd(i,j)=1]}\
=&\sum_{d|j}{\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i*d}\
=&\sum_{d|j}{\mu(d)d\frac{\lfloor\frac{n}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)}{2}}\
\end{aligned}
$$

更加特殊的 如果$j=n$ 则

$$
\begin{aligned}
&\sum_{i=1}^n{i[gcd(i,n)=1]}\
=&\sum_{d|n}{\mu(d)d\frac{\frac{n}{d}(\frac{n}{d}+1)}{2}}\
=&\frac{n}{2}\sum_{d|n}{\mu(d)(\frac{n}{d}+1)}\
=&\frac{n}{2}(\sum_{d|n}{\mu(d)\frac{n}{d}}+\sum_{d|n}{\mu(d)})\
=&\frac{n}{2}(\phi(n)+e(n))\
\end{aligned}
$$

总结

$
\begin{aligned}
&\sum_{i=1}^n{f(i)[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i*d)}\
&\sum_{i=1}^n{[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)\lfloor\frac{n}{d}\rfloor}\
&\sum_{i=1}^n{[gcd(i,n)=1]}=\phi(n)\
&\sum_{i=1}^n{i[gcd(i,j)=1]}=\sum_{d|j}{\mu(d)d\frac{\lfloor\frac{n}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)}{2}}\
&\sum_{i=1}^n{i[gcd(i,n)=1]}=\frac{n}{2}(\phi(n)+e(n))\
\end{aligned}
$

hdu6703

###name
array

###descirption
You are given an array $a_1,a_2,…,a_n(∀i∈[1,n],1≤a_i≤n)$. Initially, each element of the array is unique.

Moreover, there are m instructions.

Each instruction is in one of the following two formats:

  1. (1,pos),indicating to change the value of $a_{pos}$ to $a_{pos}+10,000,000$;
  2. (2,r,k),indicating to ask the minimum value which is not equal to any $a_i$ ( 1≤i≤r ) and not less than k.

Please print all results of the instructions in format 2.

###input
The first line of the input contains an integer T(1≤T≤10), denoting the number of test cases.

In each test case, there are two integers n(1≤n≤100,000),m(1≤m≤100,000) in the first line, denoting the size of array a and the number of instructions.

In the second line, there are n distinct integers $a_1,a_2,…,a_n (∀i∈[1,n],1≤a_i≤n)$,denoting the array.
For the following m lines, each line is of format $(1,t_1) or (2,t_2,t_3)$.
The parameters of each instruction are generated by such way :

For instructions in format 1 , we defined $pos=t_1⊕LastAns$ . (It is promised that 1≤pos≤n)

For instructions in format 2 , we defined $r=t_2⊕LastAns,k=t_3⊕LastAns$. (It is promised that 1≤r≤n,1≤k≤n )

(Note that ⊕ means the bitwise XOR operator. )

Before the first instruction of each test case, LastAns is equal to 0 .After each instruction in format 2, LastAns will be changed to the result of that instruction.

(∑n≤510,000,∑m≤510,000 )

###output
For each instruction in format 2, output the answer in one line.

###sample input
3
5 9
4 3 1 2 5
2 1 1
2 2 2
2 6 7
2 1 3
2 6 3
2 0 4
1 5
2 3 7
2 4 3
10 6
1 2 4 6 3 5 9 10 7 8
2 7 2
1 2
2 0 5
2 11 10
1 3
2 3 2
10 10
9 7 5 3 4 10 6 2 1 8
1 10
2 8 9
1 12
2 15 15
1 12
2 1 3
1 9
1 12
2 2 2
1 9

###sample output
1
5
2
2
5
6
1
6
7
3
11
10
11
4
8
11

###hint
note:
After the generation procedure ,the instructions of the first test case are :
2 1 1, in format 2 and r=1 , k=1
2 3 3, in format 2 and r=3 , k=3
2 3 2, in format 2 and r=3 , k=2
2 3 1, in format 2 and r=3 , k=1
2 4 1, in format 2 and r=4 , k=1
2 5 1, in format 2 and r=5 , k=1
1 3 , in format 1 and pos=3
2 5 1, in format 2 and r=5 , k=1
2 5 2, in format 2 and r=5 , k=2

the instructions of the second test case are :
2 7 2, in format 2 and r=7 , k=2
1 5 , in format 1 and pos=5
2 7 2, in format 2 and r=7 , k=2
2 8 9, in format 2 and r=8 , k=9
1 8 , in format 1 and pos=8
2 8 9, in format 2 and r=8 , k=9

the instructions of the third test case are :
1 10 , in format 1 and pos=10
2 8 9 , in format 2 and r=8 , k=9
1 7 , in format 1 and pos=7
2 4 4 , in format 2 and r=4 , k=4
1 8 , in format 1 and pos=8
2 5 7 , in format 2 and r=5 , k=7
1 1 , in format 1 and pos=1
1 4 , in format 1 and pos=4
2 10 10, in format 2 and r=10 , k=10
1 2 , in format 1 and pos=2

###toturial1
先不考虑修改,若只有查询,我们发现每次都是前缀的查询,这里显然是可以使用主席树用log的复杂度完成的,然后我们考虑修改,我们发现修改等价于删除数字,那么这样一来,又因为每个数都是独一无二的,删除只会让答案变得更小,且恰好变成删掉的数字,我们可以尝试用一个集合记录所有删掉的数字,然后用lower_bound来查询,和主席树得到的答案取得最小值,就是真正的答案。证明过程很简单,分类证明即可。

###code1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// 主席树+set
#include<bits/stdc++.h>
using namespace std;
#define rep(i,j,k) for(int i=j;i<=int(k);++i)

inline int read(){int x;scanf("%d",&x);return x;}

const int maxn = 1e5+5;
int ls[maxn*20*1],rs[maxn*20*1],siz[maxn*20*1],tot,rt[maxn];//update用了几次,就要乘以多少
void update(int pre,int&u,int l,int r,int pos,int val){//把u按照pre复制,然后更新pos
u=++tot;
ls[u]=ls[pre];rs[u]=rs[pre];
siz[u]=siz[pre]+val;
if(l==r)return ;
int mid=(l+r)>>1;
if(pos<=mid) update(ls[pre],ls[u],l,mid,pos,val);
else update(rs[pre],rs[u],mid+1,r,pos,val);
}

int query(int u,int l,int r,int ql,int qr){
int mid=(l+r)>>1,res=1e9;
if(ql<=l&&r<=qr){
if(l==r)return siz[u]==0?l:1e9;
if(siz[ls[u]]!=mid-l+1) return query(ls[u],l,mid,ql,qr);
else return query(rs[u],mid+1,r,ql,qr);
}
if(ql<=mid)res=min(res,query(ls[u],l,mid,ql,qr));
if(res!=1e9)return res;
if(qr>=mid+1)res=min(res,query(rs[u],mid+1,r,ql,qr));
return res;
}

int a[maxn];
int main(){
int T=read();
rep(times,1,T){
tot=0;
set<int>se;
se.insert(1e9);
int n=read(),m=read();
rep(i,1,n) update(rt[i-1],rt[i],1,n+1,a[i]=read(),1);
int lastans=0;
rep(i,1,m){
if(read()==1) se.insert(a[read()^lastans]);
else{
int r=read()^lastans,k=read()^lastans;
printf("%d\n",lastans=min(*se.lower_bound(k),query(rt[r],1,n+1,k,n+1)));
}
}
}
}

###toturial2
逆向思维,反转键值,题目让我们在键区间[1,r]上找到最小的不小于k的值,我们反转后变成了在值区间[k,n+1]上找到值最小的键,其键不小于k,修改操作就成了把值所对的键修改为无穷大,这个问题用普通最值线段树很轻松就能解决

###code2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
// 逆向思维 键值颠倒
#include<bits/stdc++.h>
using namespace std;
#define rep(i,j,k) for(int i=j;i<=int(k);++i)

inline int read(){int x;scanf("%d",&x);return x;}

#define ml ((l+r)>>1)
#define mr (ml+1)
const int maxn = 1e5+20;
int ls[maxn*2],rs[maxn*2],mx[maxn*2],a[maxn],pos[maxn],tot;//update用了几次,就要乘以多少
void build(int&u,int l,int r){
u=++tot;
if(l==r) {mx[u]=pos[l];return;}
build(ls[u],l,ml);
build(rs[u],mr,r);
mx[u]=max(mx[ls[u]],mx[rs[u]]);
}

void update(int&u,int l,int r,int q,int d){
if(l==r) {mx[u]=d;return;}
if(q<=ml) update(ls[u],l,ml,q,d);
else update(rs[u],mr,r,q,d);
mx[u]=max(mx[ls[u]],mx[rs[u]]);
}

int query(int u,int l,int r,int ql,int qr,int x){// >x
int ans=1e9;
if(ql<=l&&r<=qr){
if(mx[u]<=x) return 1e9;
if(l==r) return l;
ans=query(ls[u],l,ml,ql,qr,x);
if(ans!=1e9) return ans;
return query(rs[u],mr,r,ql,qr,x);
}
if(ml>=ql) ans=min(ans,query(ls[u],l,ml,ql,qr,x));
if(ans!=1e9) return ans;
if(mr<=qr) ans=min(ans,query(rs[u],mr,r,ql,qr,x));
return ans;
}

int main(){
int T=read();
rep(times,1,T){
tot=0;
int n=read(),m=read(),rt;
rep(i,1,n) a[i]=read(),pos[a[i]]=i;
a[n+1]=n+1,pos[n+1]=n+1;
build(rt,1,n+1);
int lastans=0;
rep(i,1,m){
if(read()==1) {
int val=a[read()^lastans];
update(rt,1,n+1,val,n+1);
// pos[val]=n+1;
}
else{
int r=read()^lastans,k=read()^lastans;
printf("%d\n",lastans=query(rt,1,n+1,k,n+1,r));
// rep(i,k,n+1) if(pos[i]>r) {cout<<" "<<i<<endl;lastans=i;break;}
}
}
}
}

hdu6588

###name

###descirption

###input

###output

###sample input

###sample output

###toturial
先来看一个简单的变形
$$
\begin{aligned}
&\sum_{i=1}^{n}gcd(x,i)\
=&\sum_{d|x}\sum_{i=1}^{n}[gcd(x,i)=d]\
=&\sum_{d|x}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(\frac{x}{d},i)=1]\
=&\sum_{d|x}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y|\frac{x}{d},y|i}\mu(y)\
=&\sum_{y|x}\sum_{d|\frac{x}{y}}\sum_{y|i,i\leq\lfloor\frac{n}{d}\rfloor}\mu(y)\
=&\sum_{y|x}\sum_{d|\frac{x}{y}}\frac{\lfloor\frac{n}{d}\rfloor}{y}\mu(y)\
=&\sum{}\
\end{aligned}
$$

题目让我们求的东西是这个
$$
\begin{aligned}
\sum_{i=1}^{n}gcd(\lfloor\sqrt[3]{i}\rfloor,i)
\end{aligned}
$$

###code

1
2
#include<bits/stdc++.h>
using namespace std;

hdu6586

###name
String

###descirption
Tom has a string containing only lowercase letters. He wants to choose a subsequence of the string whose length is k and lexicographical order is the smallest. It’s simple and he solved it with ease.
But Jerry, who likes to play with Tom, tells him that if he is able to find a lexicographically smallest subsequence satisfying following 26 constraints, he will not cause Tom trouble any more.
The constraints are: the number of occurrences of the ith letter from a to z (indexed from 1 to 26) must in $[L_i,R_i]$.
Tom gets dizzy, so he asks you for help.

###input
The input contains multiple test cases. Process until the end of file.
Each test case starts with a single line containing a string $S(|S|≤10^5)$and an integer k(1≤k≤|S|).
Then 26 lines follow, each line two numbers$ L_i,R_i(0≤L_i≤R_i≤|S|)$.
It’s guaranteed that S consists of only lowercase letters, and $∑|S|≤3×10^5$.

###output
Output the answer string.
If it doesn’t exist, output −1.

###sample input
aaabbb 3
0 3
2 3
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
###sample output
abb

###toturial
遇到这种题一般要想到一位一位去构造,贪心的选择小的字母,从而构造出最小字典序,而这一步我们需要的是验证此字母是否合法。因为是在选子序列,所以我们只需要统计后缀是否满足要求即可,后缀中的字母都满足数量足够即可

###code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include<bits/stdc++.h>
using namespace std;

#define rep(i,j,k) for(int i=(j);i<=(k);++i)
#define per(i,j,k) for(int i=(j);i>=(k);--i)

const int maxn=1e5+5;
char s[maxn];
int k;
int l[maxn],r[maxn];

int suf[maxn][26],ct[maxn][26];
char ans[maxn];

int main(){
while(~scanf("%s%d",s,&k)){
rep(i,0,25) scanf("%d%d",l+i,r+i);
int len=strlen(s);

rep(j,0,25) suf[len][j]=-1,ct[len][j]=0;
per(i,len-1,0){
rep(j,0,25) suf[i][j]=suf[i+1][j],ct[i][j]=ct[i+1][j];
suf[i][s[i]-'a']=i;
ct[i][s[i]-'a']++;
}

int cur=0;// no choose
rep(i,0,k-1){
rep(j,0,25){
if(suf[cur][j]!=-1) {
l[j]--,r[j]--;
int nex=suf[cur][j]+1;
int ok=1,nd=0;
for(int t=0;t<26;t++){
if((nex==len?0:ct[nex][t])<l[t]||r[t]<0) ok=0;
nd+=max(0,l[t]);
}
if(ok==1&&i+1+nd<=k){
ans[i]=j+'a';
cur=nex;
break;
}
l[j]++,r[j]++;
}
if(j==25){
printf("-1\n");
goto failed;
}
}
}
for(int i=0;i<k;i++) printf("%c",ans[i]);
printf("\n");

failed:;
}
}

hdu6584

###name
Meteor

###descirption
hough time passes, there is always someone we will never forget.
“The probability of being hit by a meteor is one in a billion, but it is much more miraculous, to meet you in my life.” said Tom to Jerry with affection.
“One in a billion? I may not agree with you.” answered Jerry proudly, “Let’s do the math.”

Thinking of the days they have happily spent together, Tom can’t help bursting into tears. Though Jerry has been gone for a long time, Tom still misses him every day. He remembers it was a sunny afternoon when Jerry and him were lying in the yard, working on the probability of a man being hit by a meteor.
Unlike Jerry, he was always slow. Jerry got the answer soon, but Tom was stuck as usual. In the end, Tom lost patience and asked Jerry to tell him the answer.
“I can’t be so straightforward,” snickered Jerry, “the only thing I will tell you is that the answer is $\frac{p}{q}$, where p,q≤n,gcd(p,q)=1.”
“Is it $\frac{1}{n}$?”
“Is it $\frac{1}{n-1}$?”

If answered “No” , he would try the minimum larger number that satisfies the requirement.
Tom only remembered n given by Jerry, and k, the times that he tried, but forgot what matters the most: Jerry’s answer. Now, he wants you to help him find it.

###input
The first line contains an integer $T(T≤10^2)$, the number of test cases.
The next T lines, each line contains two number$s, n,k(2≤n≤10^6)$, indicating a query.
The answer is guaranteed to be in (0,1].

###output
T lines, each line contains a fraction in term of p/q ,where gcd(p,q)=1.

###sample input
5
4 6
5 1
9 9
3 4
7 11

###sample output
1/1
1/5
1/3
1/1
3/5

###toturial
$$
\begin{aligned}
&答案肯定是找到一个分子分母小于n的分数\frac{p}{q}他满足下面的特征\
&(\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=1][\frac{i}{j} \leq \frac{p}{q}]) = k\
&我们对左式化简\
&=\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=1][i \leq \frac{p}{q}j]\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n[gcd(i,j)=1]\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^ne(gcd(i,j))\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n(u*1)(gcd(i,j))\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n\sum_{d|gcd(i,j)}u(d)*1(d)\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n\sum_{d|gcd(i,j)}u(d)\
&=\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n\sum_{d|i,d|j}u(d)\
&=\sum_{d=1}^{n}u(d)\sum_{i=1}^{\lfloor \frac{p}{q}j\rfloor}\sum_{j=1}^n[d|i,d|j]\
&=\sum_{d=1}^{n}u(d)\sum_{xd=1}^{\lfloor \frac{p}{q}(yd)\rfloor}\sum_{yd=1}^n[d|(xd),d|(yd)]\
&=\sum_{d=1}^{n}u(d)\sum_{x=1}^{\lfloor\frac{\lfloor \frac{p}{q}(yd)\rfloor}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{n}{d}\rfloor}1\
&=\sum_{d=1}^{n}u(d)\sum_{y=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{\lfloor \frac{p}{q}(yd)\rfloor}{d}\rfloor\
&=\sum_{d=1}^{n}u(d)\sum_{y=1}^{\lfloor\frac{n}{d}\rfloor}{\lfloor \frac{p}{q}y\rfloor}\
\end{aligned}
$$
这里是可以求出答案的,对d分块,右边的部分采用类欧几里得算法
我们一直往下二分,直到区间足够小,最后用 Stern-Brocot Tree 或 法雷序列找出答案

###code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;

/**** * 超级积性函数线性筛 * ****/
typedef long long ll;
const ll maxn=2e6+10;
ll no_pri[maxn]={0,1,0},pri[maxn],low[maxn];
ll PHI[maxn],DDD[maxn],XDX[maxn],MUU[maxn],SIG[maxn];
void f_ini(){
for(ll i=2;i<maxn;i++){
if(!no_pri[i]) low[i]=pri[++pri[0]]=i;
for(ll j=1;pri[j]*i<maxn;j++){
no_pri[pri[j]*i]=1;
if(i%pri[j]==0) {
low[pri[j]*i]=low[i]*pri[j];
break;
}
else low[pri[j]*i]=pri[j];
}
}

DDD[1]=PHI[1]=MUU[1]=SIG[1]=1;// 改这里
for(ll i=1;i<=pri[0];i++){
for(ll mul=pri[i],ct=1;mul<maxn;mul*=pri[i],ct++){
DDD[mul]=ct+1;// 改这里
SIG[mul]=SIG[mul/pri[i]]+mul;// 改这里
MUU[mul]=ct==1?-1:0;// 改这里
PHI[mul]=mul/pri[i]*(pri[i]-1);// 改这里
}
}

for(ll i=2;i<maxn;i++){
for(ll j=1;pri[j]*i<maxn;j++){
ll x=low[i*pri[j]], y=i*pri[j]/x;
DDD[x*y]=DDD[x]*DDD[y];
MUU[x*y]=MUU[x]*MUU[y];
PHI[x*y]=PHI[x]*PHI[y];
SIG[x*y]=SIG[x]*SIG[y];
if(i%pri[j]==0) break;
}
}

for(ll i=1;i<maxn;i++) {
DDD[i]+=DDD[i-1];
MUU[i]+=MUU[i-1];
PHI[i]+=PHI[i-1];
SIG[i]+=SIG[i-1];
XDX[i]=(DDD[i]-DDD[i-1])*i+XDX[i-1];
}
}


struct frac{
ll x,y;
frac(ll x_=0,ll y_=1){
ll gcd=__gcd(x_,y_);
x=x_/gcd;
y=y_/gcd;
}
frac operator +(const frac&rhs){
ll lcm=y/__gcd(y,rhs.y)*rhs.y;
return frac(x*(lcm/y)+rhs.x*(lcm/rhs.y),lcm);
}
frac operator /(ll k){
ll gcd=__gcd(k,x);
return frac(x/gcd,y*(k/gcd));
}
bool operator <=(const frac&rhs){
ll lcm=y/__gcd(y,rhs.y)*rhs.y;
return x*(lcm/y)<=rhs.x*(lcm/rhs.y);
}
};

// a>=0 b>=0 c>0 n>=0 -> O(lg(a,c))
void calfgh(ll a,ll b,ll c,ll n,ll&f,ll&g,ll&h){
ll A=a/c,B=b/c,s0=n+1,s1=n*(n+1)/2,s2=n*(n+1)*(2*n+1)/6;
f=s1*A+s0*B;
g=s2*A+s1*B;
h=s2*A*A+s0*B*B+2*s1*A*B-2*B*f-2*A*g;// 先减掉
a%=c,b%=c;
ll m=(a*n+b)/c;
if(m!=0) {
ll ff,gg,hh; calfgh(c,c-b-1,a,m-1,ff,gg,hh);
f+=n*m-ff;
g+=(n*m*(n+1)-hh-ff)/2;
h+=n*m*m-2*gg-ff;
}
h+=2*B*f+2*A*g;//再加上
}


ll count(frac k,int n){
ll ret=0;
for(int i=1,ed;i<=n;i=ed+1){
ed=n/(n/i);
ll a[3]; calfgh(k.x,0,k.y,n/i,a[0],a[1],a[2]);
ret+=1ll*(MUU[ed]-MUU[i-1])*a[0];
}
return ret;
}

int main(){
f_ini();
ll t,n,k;
scanf("%lld",&t);
while(t--){
scanf("%lld%lld",&n,&k);
frac l(0,1),r(1,1);// [l,r]
for(int ijk=0;ijk<40;ijk++){
frac mid=(l+r)/2;
ll ct=count(mid,n);//[0,mid]
if(ct>=k)r=mid;
else l=mid;
}
//[l,r]
frac L(0,1),R(1,0);
while(true){
frac mid(L.x+R.x,L.y+R.y);
if(mid.x<=n&&mid.y<=n&&l<=mid&&mid<=r){
printf("%lld/%lld\n",mid.x,mid.y);
break;
}
if(!(l<=mid)){
L=mid;
}
if(!(mid<=r)){
R=mid;
}
}
}
}