题意

题意：求∑i=1n∑j=1m[i=j](nmodi)⋅(mmodj)。

数据范围：1⩽n,m⩽109。

分析

这是一道简单的推式子题，但是实现比较恶心。

首先不妨设n⩽m（如果n>m交换一下就好了）

然后可以用容斥将题目拆成两个部分：

i=1∑nj=1∑m(nmodi)⋅(mmodj)−i=1∑n(nmodi)⋅(mmodi)

将两个求和拆开：

i=1∑n(nmodi)⋅j=1∑m(mmodj)−i=1∑n(nmodi)⋅(mmodi)

因为取模很难搞，所以可以用一个性质amodb=a−⌊ba⌋⋅b（取模的定义），将上式转换为：

i=1∑n(n−⌊in⌋⋅i)⋅j=1∑m(m−⌊jm⌋⋅j)−i=1∑n(n−⌊in⌋⋅i)⋅(m−⌊im⌋⋅i)

将括号拆开，可以得到：

(n2−i=1∑ni⋅⌊in⌋)⋅(m2−i=1∑mi⋅⌊im⌋)−i=1∑n(nm−mi⋅⌊in⌋−ni⋅⌊im⌋+i2⋅⌊in⌋⋅⌊im⌋)

此时我们的复杂度已经是O(n)了，然而这个复杂度仍然不足以通过本题。

要做这道题需要一个简单的技巧——整除分块。

我们很容易发现很多⌊in⌋的值是一样的，且所有的⌊in⌋为一个不下降子序列，呈块状分布

通过简单的计算可以得到，对于每个起点为i的块，它的值为⌊in⌋，终点为⌊⌊in⌋n⌋，然后我们就可以用O(n

)的算法计算∑i=1ni⋅⌊in⌋了：

主要步骤，对于每一个块[l,r]，它乘号前面前缀和差分得到，即(1+2+⋯+r)−(1+2+⋯+(l−1))，乘号后面的就是⌊ln⌋。

代码：

inline long long sum1(long long x){

return x*(x+1)%mod*inv2%mod;

}

l=1,sum=n*n%mod;

while(l<=n){

r=n/(n/l);

sum=(sum-(sum1(r)-sum1(l-1)+mod)%mod*(n/l)%mod+mod)%mod;

l=r+1;

}

同样，∑i=1mi⋅⌊im⌋的值也可以用相同的方法求得。

考虑求∑i=1n(nm−mi⋅⌊in⌋−ni⋅⌊im⌋+i2⋅⌊in⌋⋅⌊im⌋)，发现用整除分块完成它需要使区间[l,r]中所有的i都满足⌊in⌋与⌊im⌋相同。

可以用类似的方法，将上面代码中的r=n/(n/l);改为r=min(n/(n/l),m/(m/l));，然后就可以直接求值了！

且慢，虽然前三项nm,mi⋅⌊in⌋,ni⋅⌊im⌋都很容易求得，但是i2⋅⌊in⌋⋅⌊im⌋不是很容易求，因为(12+22+⋯+i2)不好处理。

这里有一个简单的结论：12+22+⋯+n2=6n(n+1)(2n+1)，会在文后证明，现在先给出这一部分的代码：

inline long long sum2(long long x){

return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;

}

l=1,tmp3=0;

while(l<=n){

long long a,b,c;

r=min(n/(n/l),m/(m/l));

a=(r-l+1)*n%mod*m%mod;

b=(sum1(r)-sum1(l-1)+mod)%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;

c=(sum2(r)-sum2(l-1)+mod)%mod*(n/l)%mod*(m/l)%mod;

tmp3=(tmp3+a-b+c+mod)%mod;

l=r+1;

}

代码

记得开long long取模很恶心，如果错了多半是这种原因因为三个整数相乘会爆long long，因此除法用逆元实现（逆元提前求得）

#include<stdio.h>

const long long mod=19940417,inv2=9970209,inv6=3323403;

long long i,j,k,m,n,l,r,ans,tmp1,tmp2,tmp3;

inline long long min(long long a,long long b){

return a<b? a:b;

}

inline void swp(long long &a,long long &b){

a+=b,b=a-b,a-=b;

}

inline long long sum1(long long x){

return x*(x+1)%mod*inv2%mod;

}

inline long long sum2(long long x){

return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;

}

int main(){

scanf("%lld%lld",&n,&m);

if(n>m)

swp(n,m);

l=1,tmp1=n*n%mod;

while(l<=n){

r=n/(n/l);

tmp1=(tmp1-(sum1(r)-sum1(l-1)+mod)%mod*(n/l)%mod+mod)%mod;

l=r+1;

}

l=1,tmp2=m*m%mod;

while(l<=m){

r=m/(m/l);

tmp2=(tmp2-(sum1(r)-sum1(l-1)+mod)%mod*(m/l)%mod+mod)%mod;

l=r+1;

}

l=1,tmp3=0;

while(l<=n){

long long a,b,c;

r=min(n/(n/l),m/(m/l));

a=(r-l+1)*n%mod*m%mod;

b=(sum1(r)-sum1(l-1)+mod)%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;

c=(sum2(r)-sum2(l-1)+mod)%mod*(n/l)%mod*(m/l)%mod;

tmp3=(tmp3+a-b+c+mod)%mod;

l=r+1;

}

ans=(tmp1*tmp2%mod-tmp3+mod)%mod;

printf("%lld\n",ans);

return 0;

}

简单的证明

证明：

12+22+⋯+n2=6n(n+1)(2n+1)

（这个式子可以用简单构造法与数学归纳法证明，由于大部分题解都是用的数学归纳法，且用数学归纳法读者自证不难，因此这里使用简单构造法）

由于有这样一个式子：i2=i2−i+i(i−1)⋅i+i，因此我们可以将这个拆开：

12+22+⋯+n2=i=1∑ni2=i=1∑n(i−1)⋅i+i=1∑ni

然后乘上31：

i=1∑ni2=31i=1∑n(i−1)⋅i⋅3+i=1∑ni

构造一下：

i=1∑ni2=31i=1∑n(i−1)⋅i⋅((i+1)−(i−2))=31i=1∑n(−(i−2)⋅(i−1)⋅i+(i−1)⋅i⋅(i+1))+i=1∑ni

把它们全部展开：

i=1∑ni2=31(−(−1)⋅0⋅1+0⋅1⋅2−0⋅1⋅2+1⋅2⋅3−1⋅2⋅3+⋯(n−1)⋅n⋅(n+1))+2n⋅(n+1)

发现括号里的很多项都可以抵消：

i=1∑ni2=3(n−1)⋅n⋅(n+1)+2n⋅(n+1)=6n⋅(n+1)⋅(2(n−1)+3)=6n⋅(n+1)⋅(2n+1)

证毕。

题目大意

按照格雷码的生成方式，输出 n 位格雷码的第 k 号二进制串。（从 0 开始编号）

思路

填完第 k 个 n 位格雷码的第 1 位后，把它转化成某一个 n−1 位格雷码继续填这个 n−1 位格雷码的第 1 位，以此类推。

如何确定转化成第几个 n−1 位格雷码：

当 k<2n−1，因为这个 n 位格雷码的前 2n−1 个二进制串的后 n−1 位是由 n−1 位格雷码正序排列而成的，所以 k 保持不变。

当 k≥2n−1，因为这个 n 位格雷码的后 2n−1 个二进制串的后 n−1 位是由 n−1 位格雷码逆序排列而成的，所以 k 要先减去 2n−1，变成 n−1 位格雷码，再翻转才能保证是逆序的，即 k=r-k+l

不妨设 \( n \leq m \)。 那原式：

$= \sum_{i = 1} ^ n (n \bmod i) \times \sum_{j = 1} ^ m (m \bmod j) - \sum_{i = 1} ^ n (n \bmod i) \times (m \bmod i)$

$= \sum_{i = 1} ^ n \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \times \sum_{j = 1} ^ m \left( m - \left\lfloor \frac{m}{j} \right\rfloor \times j \right) - \sum_{i = 1} ^ n \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( m - \left\lfloor \frac{m}{i} \right\rfloor \times i \right)$

$= \sum_{i = 1} ^ n \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \times \sum_{j = 1} ^ m \left( m - \left\lfloor \frac{m}{j} \right\rfloor \times j \right) - \sum_{i = 1} ^ n \left( nm - n \times i \times \left\lfloor \frac{m}{i} \right\rfloor - m \times i \times \left\lfloor \frac{n}{i} \right\rfloor + i ^ 2 \times \left\lfloor \frac{n}{i} \right\rfloor \times \left\lfloor \frac{m}{i} \right\rfloor \right)$

显然对这 3 坨数论分块就可以了

 https://www.luogu.me/article/v919l2og

最小割必刷题！

首先我们先把棋盘看作二分图黑白染色，比如我们可以把横纵坐标加起来为奇数的看作黑，偶数看作白。

考虑到一个黑色节点只会对它四周的节点产生影响，不难想到我们直接把把源点往所有黑点连边权为这个点在棋盘上的数的边，白点同理，只是往汇点连边，然后相邻的方格之间连无穷大的边，这样这些边就不会被算在割边中。

那我们对于舍弃一个点，就是把这个点往源点或汇点的边割掉，dinic跑最小割即可。（ISAP不会写QAQ)

一道好题。

这是一个树形结构，我们可以先拍成 $dfn$ 序，然后对于每个星球，本质上就是在其子树内区间再 删去一些小子树 所构成的一些区间，因为最多 $n$ 个操作，所以区间最多 $\mathcal{O}(n)$ 个，线段树分治，把这些区间插入，然后我们考虑如何求答案。

$y,z$ 是没用的，最终答案即为 $(x_0 - x)^2 + c$，拆开得 $- 2x_0x + {x_0}^2 + x^2 + c$。

我们考虑斜率优化，在看这个式子 $s = - 2x_0x + {x_0}^2 + x^2 + c$，移项得 $x^2 + c = 2x_0x - {x_0}^2 + s$，即点 $(x,x^2 + c)$ 斜率为 $2x_0$，维护下凸包即可，若二分则复杂度是 $\mathcal{O}(n\log^2{n})$，只需将 $x$ 与斜率 $x_0$ 都从小到大排序，我们可以利用单调栈达成 $\mathcal{O}(n\log{n})$ 的复杂度。

2156. [BZOJ 4407] 于神之怒加强版 - 题解

题目链接

题意

给定 $n, m, k$，求 $\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k$ 模 $10^9+7$ 的结果。$T$ 组数据。（$n,m,k\sim5\times10^6$，$T\sim2\times10^3$）

推导过程

令 $n\le m$

$$\begin{aligned}&\quad\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k\\&=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}d^k\cdot[\gcd(i,j)=1]\\&=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}d^k\sum_{p\mid i,j}\mu(p)\\&=\sum_{d=1}^nd^k\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}\sum_{p\mid i,j}\mu(p)\\&=\sum_{d=1}^nd^k\sum_{p=1}^{\lfloor\frac n d\rfloor}\mu(p)\left\lfloor\frac n{dp}\right\rfloor\left\lfloor\frac m{dp}\right\rfloor\\&\text{令}D=dp\\&=\sum_{D=1}^n\left\lfloor\frac n D\right\rfloor\left\lfloor\frac m D\right\rfloor\sum_{d\mid D}d^k\mu\left(\frac D d\right)\end{aligned}$$

前面的 $\displaystyle\sum_{D=1}^n\left\lfloor\frac n D\right\rfloor\left\lfloor\frac m D\right\rfloor$ 使用整除分块，仅需求后面 $\displaystyle\sum_{d\mid D}d^k\mu\left(\frac D d\right)$；设 $f(D)=\displaystyle\sum_{d\mid D}d^k\mu\left(\frac D d\right)$，显然 $f(D)$ 是积性函数，设 $D=\prod p_i^{c_i}$，则有：

$$\begin{aligned}f(D)&=f\left(\prod p_i^{c_i}\right)\\&=\prod f\left(p_i^{c_i}\right)\\&=\prod\left[1^k\mu(p_i^{c_i})+p_i^k\mu(p_i^{c_i-1})+p_i^{2k}\mu(p_i^{c_i-2})+\cdots+p_i^{(c_i-1)k}\mu(p_i)+p_i^{c_ik}\mu(1)\right]\\&=\prod\left[p_i^{c_ik}-p_i^{(c_i-1)k}\right]\\&=\prod\left[p_i^{(c_i-1)k}(p_i^k-1)\right]\end{aligned}$$

线性筛即可。

详细解释

原式如下：

$$\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k$$

可令 $n\le m$，若 $n>m$ 直接交换即可。

首先枚举 $\gcd(i,j)$ 可能等于的数，即 $1$ 到 $n$（因为 $n\le m$）。也就是枚举 $d$，当 $\gcd(i,j)=d$ 时才会产生贡献，即：

$$\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^md^k[\gcd(i,j)=d]$$

因为 $i$ 是从 $1$ 到 $n$ 的数字，$d$ 也是从 $1$ 到 $n$ 的数字，所以可以用 $d$ 乘上一个 $1$ 到 $n/d$ 的数来表示 $1$ 到 $n$ 的数，$m$ 同理：

$$\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}d^k[\gcd(i\times d,j\times d)=d]$$

既然 $i\times d$ 与 $j\times d$ 已经有了公共的因数 $d$，那么要使这俩最大公约数为 $d$ 的话，一定有 $\gcd(i,j)=1$：

$$\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}d^k[\gcd(i,j)=1]$$

因为 $\mu*1=\varepsilon$，即 $\sum_{d\mid n}\mu(d)=[n=1]$，替换掉上式的艾弗森括号：

$$\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}d^k\sum_{p\mid i,j}\mu(d)$$

$d$ 的取值只与第一层求和有关，与第二、三层求和无关，故可以把 $d^k$ 提到前面：

$$\sum_{d=1}^nd^k\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}\sum_{p\mid i,j}\mu(p)$$

与其枚举 $i,j$，不如只枚举 $\mu(p)$ 会出现多少次，如果固定了 $p$，那么在 $i$ 在遍历 $1$ 到 $\lfloor\frac n d\rfloor$ 时，会出现 $\lfloor\frac n{dp}\rfloor$ 次 $p$，或者说 $1$ 到 $\lfloor\frac n d\rfloor$ 中有 $\lfloor\frac n{dp}\rfloor$ 个数是 $p$ 的倍数（$j$ 同理）。因为 $p$ 是 $i,j$ 的因数且 $i$ 的枚举范围小于 $j$ 的枚举范围（因为 $n\le m$），所以 $p$ 的枚举范围和 $i$ 相同：

$$\sum_{d=1}^nd^k\sum_{p=1}^{\lfloor\frac n d\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac n{dp}\rfloor}\sum_{j=1}^{\lfloor\frac m{dp}\rfloor}1$$

后面两个求和已经能直接算出来了：

$$\sum_{d=1}^nd^k\sum_{p=1}^{\lfloor\frac n d\rfloor}\mu(p)\left\lfloor\frac n{dp}\right\rfloor\left\lfloor\frac m{dp}\right\rfloor$$

令 $D=dp$，转为枚举 $D$ 之后 $\lfloor\frac n{dp}\rfloor\lfloor\frac m{dp}\rfloor$ 就可以直接提出来；$p=D/d$，所以 $\mu(p)$ 被替换成 $\mu(D/d)$；第二层循环与 $d$ 无关，所以 $d^k$ 可以提到第二层循环；此时第二层循环里面就变成了 $d^k\mu(D/d)$，因为 $d=D/p$，所以 $d$ 的枚举范围就是 $d\mid D$：

$$\sum_{D=1}^n\left\lfloor\frac n D\right\rfloor\left\lfloor\frac m D\right\rfloor\sum_{d\mid D}d^k\mu\left(\frac D d\right)$$

观察前面的 $\sum_{D=1}^n\lfloor\frac n D\rfloor\lfloor\frac m D\rfloor$ 可以想到用整除分块做，但是想要用整除分块，需要知道后面 $\sum_{d\mid D}d^k\mu(\frac D d)$ 的前缀和（整体的，指在 $D$ 取 $1,2,3\cdots$ 时整个求和结果的前缀和），那就单独来研究它。

令

$$f(D)=\sum_{d\mid D}d^k\mu\left(\frac D d\right)$$

$d^k$ 其实就是恒等函数 $\text{id}_k(d)$，$\mu(\frac D d)$ 就是欧拉函数 $\mu(\frac D d)$，那 $f(D)$ 不就是 $\text{id}_k$ 和 $\mu$ 的狄利克雷卷积嘛！

$$f(D)=\sum_{d\mid D}d^k\mu\left(\frac D d\right)=(\text{id}_k*\mu)(D)$$

狄利克雷卷积具有积性保持性（狄利克雷生成函数 - OI Wiki），所以 $f(D)$ 也是积性函数！那就好办了！如果想求 $F(D)$ 的话（注意！上面那个式子 $f(D)$ 表示一个关于 $D$ 的函数，$D$ 是变量，而现在说的 $D$ 是一个给定的常数！），不妨设 $D=\prod p_i^{c_i}$，则根据积性函数的定义：

$$f(D)=f\left(\prod p_i^{c_i}\right)=\prod f(p_i^{c_i})$$

把 $f(p_i^{c_i})$ 都带入 $f(D)=\sum_{d\mid D}d^k\mu(\frac D d)$（同样，这里的 $D$ 只是一个函数中的自变量，与上面那个式子中的 $D$ 含义不同！）：

$$\prod\left[1^k\mu(p_i^{c_i})+p_i^k\mu(p_i^{c_i-1})+p_i^{2k}\mu(p_i^{c_i-2})+\cdots+p_i^{(c_i-1)k}\mu(p_i)+p_i^{c_ik}\mu(1)\right]$$

因为 $p_i$ 是质数，根据莫比乌斯函数的定义（[莫比乌斯反演 - OI Wiki](https://oiwiki.com/math/number-theory/mobius/#定义)），当 $x$ 含有平方因子时 $\mu(x)=0$，那么 $\mu(p_i^{c_i}),\mu(p_i^{c_i-1}),\mu(p_i^{c_i-2})\cdots\mu(p_i^2)$ 都等于 $0$！带入到上式：

$$\prod\left[p_i^{c_ik}-p_i^{(c_i-1)k}\right]$$

居然只剩这么两项了！镇棒！把 $p_i^{(c_i-1)k}$ 给提出来：

$$\prod\left[p_i^{(c_i-1)k}(p_i^k-1)\right]$$

所以，我们得到了函数 $f(D)$ 的表达式：

$$f(D)=\prod\left[p_i^{(c_i-1)k}(p_i^k-1)\right]$$

哇哇！这样不就可以线性筛了嘛！为什么呢？看看线性筛的代码：

void solve(int n) {
   for (int i = 2; i <= n; ++i) {
       if (!vis[i]) {
           primes.push_back(i);
           // 更新 f[i] 值
       }
       for (int p : primes) {
           if (i * p > n) break;
           vis[i * p] = 1;
           if (i % p == 0) {
               // 更新 f[i] 值
               break;
           } else {
               // 更新 f[i] 值
           }
       }
   }
}

假如说现在 $i$ 是质数，程序运行到第 5 行，根据上面得到的式子：$f(D)=\prod[p_i^{(c_i-1)k}(p_i^k-1)]$，这个时候 $p_i=i$ 且 $c_i=1$，那 $f(i)$ 就等于 $p_i^k-1$ 也就是 $i^k-1$。

if (!vis[i]) {
   primes.push_back(i);
   f[i] = (kasumi(i, k) - 1 + MOD) % MOD;
}

到了下面筛 $i\times p$ 时，假如 $i\bmod p\neq0$，也就是 $p$ 不是 $i$ 的因数，$i$ 乘上 $p$ 之后为 $i$ 多添加了一个全新的因数，而且这个因数还是个质数（根据线性筛的定义，每个数只被它最小的质因数筛去（筛法 - OI Wiki）），那么根据刚才的式子：$f(D)=\prod[p_i^{(c_i-1)k}(p_i^k-1)]$，乘积中多添加了 $p^{(c-1)k}(p^k-1)$ 这一项。因为 $p$ 是新添加的因数，所以 $c=1$，那 $p^{(c-1)k}=1$，相当于在乘积中只多添加了 $p^k-1$ 这一项（又或者说是添加了 $f(p)$ 这一项，因为 $p$ 是质数，则 $f(p)=p^k-1$；另一种解释是 $\gcd(i,p)=1$，所以根据积性函数的定义也能得证。），则 $f(i\times p)=f(i)(p^k-1)$（或 $f(i\times p)=f(i)f(p)$）。

if (i % p == 0) {
   // ...
   break;
} else {
   f[i * p] = f[i] * (kasumi(p, k) - 1 + MOD) % MOD; // 或 f[i * p] = f[i] * f[p] % MOD;
}

那如果 $i\bmod p=0$ 呢？也就是 $p\mid i$，$p$ 在 $i$ 的质因数中已经出现过了，$i$ 乘上 $p$ 之后相当于把 $p$ 这个质因数的指数加了 $1$，那么根据上式 $f(D)=\prod[p_i^{(c_i-1)k}(p_i^k-1)]$，$i$ 乘 $p$ 相当于把 $c_i$ 加 $1$，相当于在乘积中添加了 $p_i^k$ 这一项（$\displaystyle\frac{p_i^{(c_i+1-1)k}(p_i^k-1)}{p_i^{(c_i-1)k}(p_i^k-1)}=p_i^k$），则 $f(i\times p)=f(i)p_i^k$。

if (i % p == 0) {
   f[i * p] = f[i] * kasumi(p, k) % MOD;
   break;
} else {
   // ...
}

终于解决完 $f(D)$ 的问题啦，回到原来的答案式：

$$\sum_{D=1}^n\left\lfloor\frac n D\right\rfloor\left\lfloor\frac m D\right\rfloor f(D)$$

因为线性筛时涉及到快速幂，那么只能用 $O(n\log_2k)$ 的时间复杂度求出 $f(D)$ 即其前缀和，加上整除分块，可以用 $O(\sqrt n)$ 的时间复杂度求出答案式。那么整体的时间复杂度就是 $\Theta(n\log_2k+T\sqrt n)$，$\Omega(n\log_2k)$。

代码

#include <cstdio>
#include <vector>
#include <algorithm>
typedef long long ll;

const int MAXN = 5e6 + 10;
const int MOD = 1e9 + 7;

ll kasumi(ll x, ll y) {
  ll res = 1;
  while (y) {
    if (y & 1) res = res * x % MOD;
    y >>= 1;
    x = x * x % MOD;
  }
  return res;
}

ll k;

bool vis[MAXN];
ll f[MAXN];
::std::vector <ll> primes, mi;
void solve(ll n) {
  f[1] = 1;
  for (ll i = 2; i <= n; ++i) {
    if (!vis[i]) {
      primes.push_back(i);
      f[i] = (kasumi(i, k) - 1 + MOD) % MOD;
    }
    for (ll p : primes) {
      if (i * p > n) break;
      vis[i * p] = 1;
      if (i % p == 0) {
        f[i * p] = f[i] * kasumi(p, k) % MOD;
        break;
      } else {
        f[i * p] = f[i] * f[p] % MOD; /* kasumi(p, k) - 1*/
      }
    }
  }
  for (ll i = 1; i <= n; ++i) {
    f[i] = (f[i] + f[i - 1]) % MOD;
  }
}

int T;
ll n, m;
ll ans;

int main() {
  freopen("bzoj_4407.in", "r", stdin);
  freopen("bzoj_4407.out", "w", stdout);
  scanf("%d %lld", &T, &k);
  solve(MAXN - 10);
  while (T--) {
    scanf("%lld %lld", &n, &m);
    if (n > m) n ^= m ^= n ^= m;
    ans = 0;
    for (ll i = 1, ed; i <= n; i = ed + 1) {
      ed = ::std::min(n / (n / i), m / (m / i));
      ans = (ans + (n / i) * (m / i) % MOD * (f[ed] - f[i - 1] + MOD) % MOD) % MOD;
    }
    printf("%lld\n", ans);
  }
  return 0;
}

Update

[2025/06/20] 修正了一处代码错误。

if (i % p == 0) {
 // ...
 break;
} else {
 f[i * p] = f[i] * (kasumi(p, k) - 1 + MOD) % MOD; // 或 f[i * p] = > f[i] * f[p] % MOD;
 // - f[i * p] = f[i] * kasumi(p, k) % MOD; // 或 f[i * p] = f[i] * f[p] % MOD;
}

忘记减一了（

[2025/07/04] 修改了几处错别字。