『你从未见过的船新数论』

Preface

\(\mathbf{References}:\)\(\text{zqy's blogs : cnblogs},\)云烟万象但过眼\(.\)

前排提示：本文数学公式较多，加载\(\LaTeX\)需要一定时间，可能会导致浏览器暂时卡顿，请耐心等待数学公式正常显示.

快速幂和快速乘

原理：二进制分解和取模的定义

\[a\times b\bmod p=a\times b-p\times \left\lfloor\frac{a\times b}{p} \right\rfloor \]

inline int fastpow(int a,int b) { int c = 1; for (; b; Mul(a,a), b>>=1) if (1&b) Mul(c,a); return c; }
inline ll fastmul(ll a,ll b) { return ( a * b - (ll)( (long double) a * b / p ) * p + p ) % p; }

整除和约数

整除

如果\(a\)是\(b\)的约数，则称\(a\)整除\(b\)，记作\(a|b\)。

整除的性质：

\(1.\) \(a|b,b|c\Rightarrow a|c\)

\(2.\) \(a|b,c|d\Rightarrow ac|bd\)

\(3.\) \(a|b,a|c\Rightarrow a|(bn+cm)\)

\(4.\) \(ma|mb\Rightarrow a|b\)

约数

\[d|a,d|b\Rightarrow d|(a-b) \]

所以有\(\gcd(a,b)=\gcd(a-b,b)=\gcd(a\bmod b,b)\)，每次取模\(a\)减小一半，时间复杂度\(\mathcal O(\log n)\)。

\[\mathrm{lcm}(a,b)=\frac{a\times b}{\gcd(a,b)} \]

inline int gcd(int a,int b) { return b == 0 ? a : gcd( b, a % b ); }
inline int lcm(int a,int b) { return a / gcd(a,b) * b; }

素数和筛法

素数定理

\[\pi(n)=\sum_{i=1}^n[i\in \mathrm{Prime}]\sim \frac{n}{\ln n} \]

Eratosthenes筛法

\[\sum_{i=1}^n\frac{[i\in \mathrm{Prime}]}{i}\sim \ln\ln n \]

线性筛

用每个数字的最小质因子来筛掉它，时间复杂度\(\mathcal O(n)\)。

inline void EulerSieve(int Lim) {
    for (int i = 2; i <= Lim; i++) {
        if ( !flag[i] ) Prime[++cnt] = i;
        for (int j = 1; j <= cnt && i * Prime[j] <= Lim; j++) {
            flag[ i * Prime[j] ] = true;
            if ( i % Prime[j] == 0 ) break;
        }
    }
}

算术基本定理和推论

每个正整数\(n\)都有如下的唯一分解：

\[n=\prod_{i=1}^k p_i^{a_i},\forall i\in[1,k]\ p_i\in\mathrm{Prime} \]

有如下推论成立：

\[\tau (n)=\prod_{i=1}^k (a_i+1),\sigma(n)=\prod_{i=1}^k\left(\sum_{j=0}^{a_i}p_i^j\right) \\ a=\prod_{i=1}^kp_i^{a_i},b=\prod_{i=1}^kp_i^{b_i}\Longrightarrow \gcd(a,b)=\prod_{i=1}^k p_i^{\min(a_i,b_i)},\mathrm{lcm}(a,b)=\prod_{i=1}^kp_i^{\max(a_i,b_i)} \]

同余

同余的定义和性质

若\(m|(a-b)\)，则称\(a,b\)模\(m\)同余，记作\(a\equiv b\pmod m\)。

同余的基本性质：自反性，对称性，传递性，同加性，同乘性，同幂性。

同余的其他性质：

\(1.\) \(a\equiv b\pmod n,a\equiv b\pmod n\Rightarrow a\equiv b\pmod {\mathrm{lcm}(n,m)}\)

\(2.\) \(n|m,a\equiv b\pmod n\Rightarrow a\equiv b\pmod {n}\)

\(3.\) \(\gcd(c,m)=d,ac\equiv bc\pmod m\Rightarrow a\equiv b\pmod{\frac{m}{d}}\)

Fermat定理

若\(p\)是质数，则对于任意整数\(a\)，\(a^p\equiv a\pmod p\)。

Wilson定理

若\(p\)是素数，则\((p-1)!\equiv -1\pmod p\)。

同余类和剩余系

对于任意\(a\in[0,p-1]\)，集合\(\{a+kp|k\in\mathbb{Z}\}\)中所有数模\(p\)都同余，余数为\(a\)，称该集合为模\(p\)意义下的一个同余类，简记为\(\overline a\)。

模\(p\)意义下的同余类一共有\(p\)个，分别为\(\overline 0,\overline 1,\overline 2,\cdots,\overline {p-1}\)，他们构成模\(p\)意义下的完全剩余系，简称完系。

\(1\sim p\)中与\(p\)互质的数有\(\varphi(p)\)个，这些数字对应的同余类组成了\(p\)的化简剩余系，也称既约剩余系或缩系。

既约剩余系的重要性质是关于模\(p\)乘法封闭，也就是\(a,b\)与\(p\)互质，显然\(a\times b, a\times b\bmod p\)都与\(p\)互质，所以\(a\times b\bmod p\)也在既约剩余系里面。

Euler定理

\[a^b\equiv \begin{cases} a^{b\bmod \varphi(p)},& \gcd(a,p)=1 \\ a^b,& \gcd(a,p)\not = 1, b<\varphi(p) \\ a^{b\bmod \varphi(p)+\varphi(p)},&\gcd(a,p)\not= 1, b\geq \varphi(p) \end{cases} \pmod p \]

可以用来缩小指数规模。

对于\(\gcd(a,p)=1\)的情况，可以用剩余系的知识来证明：

设\(p\)的既约剩余系为\(\{\overline{a_1},\overline{a_2},\cdots,\overline{a_{\varphi(p)}}\}\)，对于\(a_i,a_j\)，若\(a\times a_i\equiv a\times a_j\pmod p\)，则\(a(a_i-a_j)\equiv 0\pmod p\)，由于\(\gcd(a,p)=1\)，所以\(a_i\equiv a_j\pmod p\)。故对于\(a_i\not =a_j\)，都有\(a\times a_i\)和\(a\times a_j\)表示不同的同余类。

又因为既约剩余系乘法封闭，所以\(\overline{aa_i},\overline{aa_j}\)都在既约剩余系中，因此\(\{\overline{aa_1},\overline{aa_2},\cdots,\overline{aa_{\varphi(p)}}\}\)也表示\(p\)的既约剩余系，故：

\[\prod_{i=1}^{\varphi(p)}aa_i\equiv a^{\varphi(p)}\prod_{i=1}^{\varphi(p)}a_i\equiv\prod_{i=1}^{\varphi(p)}a_i\pmod p \]

所以\(a^{\varphi(p)}\equiv 1\pmod p\)，证毕。

同余方程

裴蜀定理

方程\(ax+by=c\)有整数解，当且仅当\(c|\gcd(a,b)\)，此时有无数多组整数解。

证明：

显然要证方程\(ax+by=\gcd(a,b)\)有解，那么\(x'=x\frac{c}{\gcd(a,b)},y'=y\frac{c}{\gcd(a,b)}\)就是原方程的整数解。

不妨设\(ax_1+by_1=\gcd(a,b),bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)\)，所以：

\[ax_1+by_1=bx_2+(a\bmod b)y_2=bx_2+\left(a-b\times \left\lfloor\frac{a}{b}\right\rfloor\right)y_2 \\ ax_1+by_1=ay_2+b\left(x_2-\left\lfloor\frac{a}{b}\right\rfloor y_2\right) \]

显然只要方程\(bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)\)有解，就可以构造一组原方程的解。

此时归纳推理，问题转为求方程\(\gcd(a,b)x=\gcd(a,b)\)，显然\(x=1,y\in \mathbb{Z}\)就是一组解，那么原命题得证。

拓展欧几里得算法

把上述迭代过程实现一下，时间复杂度\(\mathcal O(\log n)\)。

inline int Exeuclid(int a,int b,int &x,int &y) {
    if ( b == 0 ) return x = 1, y = 0, a;
    int p = Exeuclid(b,a%b,y,x); return y -= a / b * x, p;
}

设\(\gcd(a,b)=p\)，原方程的一组特解为\((x_0,y_0)\)，则原方程的通解为：

\[x=x_0+\frac{b}{p}\times k,y=y_0-\frac{a}{p}\times k, k\in \mathbb{Z} \]

乘法逆元

求方程\(ax\equiv 1\pmod p\)的解。

方法\(1\)：当\(p\in \mathrm{Prime}\)时，\(x\equiv a^{-1}\equiv a^{p-2}\pmod p\)，快速幂即可。

方法\(2\)：原方程有解等价于\(p|(ax-1)\)，不妨设\(ax-1=-py\)，那么\(ax+py=1\)，显然原方程有解等价于\(\gcd(a,p)=1\)，那么可以用扩展欧几里得算法求解。

方法\(3\)：不妨设\(p=ka+t(t< a)\)，那么\(ka+t\equiv 0\pmod p\)，等式两边同时乘\(a^{-1}t^{-1}\)，那么\(a^{-1}\equiv -kt^{-1}\equiv -\left\lfloor\frac{p}{a}\right\rfloor(p\bmod a)^{-1}\pmod p\)，递归下去，时间复杂度\(\mathcal O(\log n)\)。如果改成递推，可以\(\mathcal O(n)\)求\(1\sim n\)所有数字的逆元。

线性同余方程组

求解方程组

\[\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\\ \ \ \cdots \\ x\equiv a_n\pmod {m_n} \end{cases} \]

当所有\(m_i\)互质的时候，可以令\(M=\prod m_i,M_i=\frac{M}{m_i},t_i=M_i^{-1}\pmod {m_i}\)，则\(x=\sum a_iM_it_i\)是原方程组的一个解，这个结论叫做中国剩余定理。

当\(m_i\)不满足两两互质的时候，可以归纳求解，不妨设前\(i-1\)个方程的解为\(x'\)，\(\mathrm{lcm}_{j=1}^{i-1}\{m_j\}=m\)，显然\(x'+km(k\in \mathbb{Z})\)都是前\(i-1\)个方程的解，只要让\(x'+km\equiv a_i\pmod {m_i}\)成立就求出了前\(i\)个方程的解，这时候解一个单个的同余方程即可，如果无解，那么原方程组就无解。

这样解方程很容易有整型溢出的风险，要注意数据范围，即使\(\mathrm{lcm}\)在\(\mathrm{long\ long}\)范围内，也用用快速乘。

inline ll ExCRT(void) {
    ll res = a[1] % m[1], l = m[1], p, x, y, v;
    for (int i = 2; i <= n; i++) {
        v = ( ( a[i] - res ) % m[i] + m[i] ) % m[i];
        if ( v % gcd(l,m[i]) ) return -1; p = Exeuclid(l,m[i],x,y);
        x = ( x % (m[i]/p) + (m[i]/p) ) % (m[i]/p), x = fastmul( x, v/p, m[i]/p );
        res = res + fastmul( x, l, lcm(l,m[i]) ), res %= ( l = lcm(l,m[i]) );
    }
    return res;
}

离散对数

求方程\(a^{x}\equiv b\pmod p\)的解，\(\gcd(a,p)=1\)。

根据欧拉定理，\(a^{\varphi(p)}\equiv 1\pmod p\)，所以\(x\)一定在区间\([1,\varphi(p)]\)里面，直接枚举的话时间复杂度是\(\mathcal O(\varphi(p))\)的。

考虑分块，令\(x=kT-r(r\leq T)\)，那么\(a^{kT-r}\equiv b\pmod p\)，所以\(a^{kT}\equiv a^{r}b\pmod p\)。此时可以\(\mathcal O(T)\)的时间枚举\(a^{r}b\)的取值，存在一个\(\mathrm{hash}\)表里面，那么可以枚举\(k\in\left[1,\left\lceil\frac{\varphi(p)}{T}\right\rceil+1\right]\)，然后查表即可，取\(T=\sqrt{\varphi(p)}\)，时间复杂度\(\mathcal O(\sqrt {\varphi(p)})\)。

inline ll BSGS(ll a,ll b,ll p) {
    ll T = sqrt(p), v1 = 1, v2; map<ll,ll> h;
    for (int i = 0; i < T; v2 = v1 = Mul(v1,a,p), i++) h[Mul(v1,b,p)] = i;
    for (int i = 1; i <= T + 2; v2 = Mul(v2,v1,p), i++)
        if ( h.count(v2) ) return i * T - h[v2];
    return -1;
}

Lucas定理

如果\(p\)是质数，那么

\[{n\choose m}\equiv {n\bmod p\choose m\bmod p}\times {\left\lfloor\frac{n}{p}\right\rfloor\choose\left\lfloor\frac{m}{p}\right\rfloor}\pmod p \]

可以在模数较小时的求组合数。

inline int C(int n,int m) { return n < m ? 0 : mul( fac[n], mul( finv[m], finv[n-m] ) ); }
inline int Lucas(int n,int m) { return n < p && m < p ? C(n,m) : mul( C(n%p,m%p), Lucas(n/p,m/p) ); }

素数判定和分解质因数

MillerRabin测试

首先我们可以用费马小定理来逆向测试，也就是随机一些\(a\)，判定\(a^{x-1}\bmod x\)是否等于\(1\)，如果不是\(1\)，显然\(x\)是合数。

但是这样做仍然有些强伪素数会通过测试，也就是说一个数\(x\)现在满足对于随机的\(a\)，都有\(a^{x-1}\equiv 1\pmod x\)，现在我们要判定\(x\)是不是素数。

我们知道当\(p\)是素数的时候\(x^2\equiv 1\pmod p\)的解是\(x=1/x=p-1\)，所以我们就把\(a^{x-1}\)开方，如果不是\(1/p-1\)就可以判定它不是素数，如果是开方后是\(1\)的话可以继续判定。

这样做在\(\mathrm{long\ long}\)范围内不会误判。

inline bool MillerRabin(int x) {
    if ( x == 2 || x == 3 || x == 5 || x == 7 ) return true;
    if ( x == 1 || x % 2 == 0 ) return false;
    int a = x - 1, b = 0; while ( not( a & 1 ) ) a>>=1, b++;
    for (int i = 1, j; i <= 8; i++) {
        int base = rand() % ( x - 2 ) + 2, v = fastpow(base,a,x);
        if ( v == 1 || v == x - 1 ) continue;
        for (j = 0; j < b; j++) if ( ( v = 1LL * v * v % x ) == x - 1 ) break;
        if ( j == b ) return false;
    }
    return true;
}

Pollard-Rho分解质因数

我们可以在\([1,n]\)间随机\(x_1,x_2,\cdots,x_k\)这\(k\)个数字，然后试图判定\(\gcd(|x_i-x_j|,n)>1\)，如果存在，那么我们就找到了一个\(n\)的非平凡质因子，根据生日悖论，其成功率很可观。

随机数可以用\(x_1=c,x_i=(x_{i-1}^2+c)\bmod n\)的方式生成，分布比较均匀，但是可能会出现\(\rho\)型环，需要特判。

一个很有效的判环策略是：令一个变量\(x'\)以每次迭代两下的速度同时计算，如果出现\(x'=x\)，那么就说明出现环，此时只要更改常数\(c\)的值重新计算即可。

这样做的期望复杂度是\(\mathcal O(n^{\frac{1}{4}}\mathrm{polylog}(n))\)的，要配合\(\mathrm{Miller\ Rabin}\)测试，基本不会误判。处理\(10^{18}\)范围内的数据时，要注意配合快速乘。

inline ll PollardRho(ll x) {
    if ( MillerRabin(x) ) return x;
    for (ll c = 3, x1, x2, i, j; x1 = x2 = i = 1, j = 2; ++c)
        for (ll d; x1 = (Mul(x1,x1,x)+c) % x, d = gcd(abs(x1-x2),x);) {
            if ( 1 < d && d < x ) return max( PollardRho(d), PollardRho(x/d) );
            if ( x1 == x2 ) break; if ( ++i == j ) j <<= 1, x2 = x1;
        }
}

阶和原根

阶的定义和性质

若\(\gcd(a,p)=1\)，称方程\(a^x\equiv 1\pmod p\)的最小正整数解\(x\)为\(a\)模\(p\)的阶，记作\(\delta_p(a)\)。

性质\(1\)：对于任意的\(n\)满足\(a^n\equiv 1\pmod p\)，都有\(\delta_p(a)|n\)，证明如下：

若不成立，则一定有\(n=k\times \delta_p(a)+r(r<\delta_p(a))\)，那么：\(a^{n}\equiv a^{k\times \delta_p(a)}\times a^r\equiv 1\pmod p\)，所以\(a^r\equiv 1\pmod p\)，这与\(\delta_p(a)\)的最小性矛盾，故假设不成立。

根据欧拉定理可得推论：\(\delta_p(a)|\varphi(p)\)。

性质\(2\)：设\(x=\delta_p(a)\)，则\(a^0,a^1,\cdots,a^{x-1}\)两两不同余，证明如下：

若不成立，不妨设\(0\leq i<j<x\)满足\(a^i\equiv a^j\pmod p\)，那么\(a^{j-i}\equiv 1\pmod p\)，显然\(0<j-i<x\)，与\(\delta_p(a)\)的最小性矛盾，故假设不成立。

性质\(3\)：设\(x=\delta_p(a)\)，则\(\delta_p(a^t)=\frac{x}{\gcd(x,t)},t\in \mathbb{N}^+\)，证明如下：

不妨设\(\delta_p(a^t)=y\)，根据阶的定义\(a^{ty}\equiv 1\pmod p\)，根据阶的性质\(1\)，有\(x|ty\)，可以导出\(\frac{x}{\gcd(x,t)}|\frac{t}{\gcd(x,t)}\times y\)，于是\(\frac{x}{\gcd(x,t)}|y\)，根据阶的最小性，显然\(y=\frac{x}{\gcd(x,t)}\)。

原根的定义和性质

若\(\gcd(g,p)=1,\delta_p(g)=\varphi(p)\)，则称\(g\)为模\(p\)意义下的一个原根。

根据阶的性质和既约剩余系的乘法封闭性，我们还可以知道\(g\)为模\(p\)意义下的一个原根，当且仅当\(\{\overline {g},\overline{g^2},\cdots,\overline{g^{\varphi(p)}}\}\)构成一个模\(p\)意义下的既约剩余系。

最小原根的求法：我们可以直接在\(2\sim p-1\)范围内枚举\(g\)，只要验证\(\forall i\in[2,\varphi(p)-1]\)，\(g^{i}\equiv 1\pmod p\)恒不成立即可。而根据推论，如果存在\(g^i\equiv 1\pmod p\)，一定有\(i|\varphi(p)\)，所以只要枚举\(\varphi(p)\)的约数验证即可。再进一步，我们可以对于每个\(\varphi(p)\)的质因数\(p_i\)用\(\frac{\varphi(p)}{p_i}\)来验证，因为他们涵盖了\(\varphi(p)\)所有真因子的倍数。

可以证明，一个正整数\(p\)的最小原根是\(O(n^{\frac{1}{4}})\)级别的，那么结合\(\text{Pollard-}\rho\)算法，可以实现\(O(p^{\frac{1}{4}}\log^2 p)\)求一个最小原根。

inline void PollardRho(ll x) {
    if ( MillerRabin(x) ) return fac.push_back(x), void();
    for (ll c = 3, x1, x2, i, j; x1 = x2 = i = 1, j = 2; ++c)
        for (ll d; x1 = (Mul(x1,x1,x)+c) % x, d = gcd(abs(x1-x2),x);) {
            if ( 1 < d && d < x ) return PollardRho(d), PollardRho(x/d);
            if ( x1 == x2 ) break; if ( ++i == j ) j <<= 1, x2 = x1;
        }
}
inline ll Phi(ll p) {
    if ( MillerRabin(p) ) return p-1; ll val = p;
    PollardRho(p), sort( fac.begin(), fac.end() );
    fac.erase( unique(fac.begin(),fac.end()), fac.end() );
    for (int i = 0; i < fac.size(); i++)
        val /= fac[i], val *= (fac[i]-1);
    return fac.clear(), val;
}
inline ll PrimitiveRoot(ll p,ll val = 0) {
    if ( p == 2 ) return 1; PollardRho( val = Phi(p) );
    sort( fac.begin(), fac.end() );
    fac.erase( unique(fac.begin(),fac.end()), fac.end() );
    for (ll i = 2, flag; flag = 1, i < p; i++) {
        for (int j = 0; flag && j < fac.size(); j++)
            if ( Pow(i,val/fac[j],p) == 1 ) flag = 0;
        if ( flag ) return fac.clear(), i;
    }
}

原根存在判定定理：在模\(m\)意义下存在原根，当且仅当\(m=2,4,p^k,2p^k\)，其中\(p\)为质数，\(k\)为正整数。

求出一个数的最小原根\(g\)以后，我们可以进一步求出这个数的所有原根：集合\(G=\{g^s|s\in[1,\varphi(p)],\gcd(s,\varphi(p))=1\}\)给出了模\(p\)意义下的所有原根，一共有\(\varphi(\varphi(p))\)个，证明如下：

根据原根的性质，集合\(G_0=\{g^s|s\in[1,\varphi(p)]\}\)生成了模\(p\)意义下的既约剩余系，已经考察了所有可能的原根，所以只要说明\(g^s\)是模\(p\)意义下的原根，一定满足\(\gcd(s,\varphi(p))=1\)即可。根据阶的性质：\(\delta_p(g^s)=\frac{\varphi(p)}{\gcd(\varphi(p),s)}\)，如果\(g^s\)是原根，显然\(\delta_p(g^s)=\varphi(p)\)，那么就有\(\gcd(\varphi(p),s)=1\)，进而可知原根一共有\(\varphi(\varphi(p))\)个。

原根的应用

假设\(g\)是模\(p\)意义下的原根，那么称方程\(g^x\equiv a\pmod p\)的最正小整数解为离散对数，记作\(\text{ind}_g(a)\)。

离散对数具有和连续对数一样的性质：

\[\text{ind}_g(ab)\equiv \text{ind}_g(a)+\text{ind}_g(b)\pmod{\varphi(p)} \\ \text{ind}_g(a^t)\equiv t\times \text{ind}_g(a)\pmod{\varphi(p)} \]

离散对数可以在\(O(\sqrt p)\)的时间内求解\(k\)次剩余，即方程\(x^{k}\equiv a\pmod p\)的解：

\[x^{k}\equiv a\pmod p \\ k\times \text{ind}_g (x)\equiv \text{ind}_g (a)\pmod {\varphi(p)} \\ \left(g^k\right)^{\text{ind}_g(x)}\equiv g^{\text{ind}_g(a)}\equiv a\pmod p \]

处理出原根\(g\)之后，\(\text{ind}_g (x)\)的值就可以用\(\text{BSGS}\)算法求出，然后再快速幂还原出\(x\)即可。

不妨设\(x_0\equiv g^c\pmod p\)是原问题的一个解，那么\(x^k\equiv g^{kc+t\varphi(p)}\equiv a\pmod p\)，显然解\(x\equiv g^{c+\frac{t\varphi(p)}{k}}\pmod p\)，参数\(t\)显然要满足\(\frac{k}{\gcd(k,\varphi(p))}|t\)，不妨设\(t=i\times \frac{k}{\gcd(k,\varphi(p))},i\in\mathbb{N}\)，那么原问题的所有解为：\(x\equiv g^{c+\frac{\varphi(p)}{\gcd(\varphi(p),k)}\times i}\pmod p,i\in \mathbb{N}\)。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll; vector<ll> fac;
inline ll Mul(ll a,ll b,ll p) { return ( a * b - (ll)( (long double) a * b / p ) * p + p ) % p; }
inline ll Pow(ll a,ll b,ll p) { ll c = 1; for (; b; a = Mul(a,a,p), b >>= 1) if (1&b) c = Mul(c,a,p); return c; }
inline ll gcd(ll a,ll b) { return b == 0 ? a : gcd( b, a % b ); }
inline ll lcm(ll a,ll b) { return a / gcd(a,b) * b; }
inline bool MillerRabin(ll x) {
    if ( x == 2 || x == 3 || x == 5 || x == 7 ) return true;
    if ( x == 1 || x % 2 == 0 ) return false;
    ll a = x - 1, b = 0; while ( not( a & 1 ) ) a>>=1, ++b;
    for (int i = 1, j; i <= 8; i++) {
        ll Base = rand() % ( x - 2 ) + 2, v = Pow(Base,a,x);
        if ( v == 1 || v == x - 1 ) continue;
        for (j = 0; j < b; j++) if ( ( v = Mul(v,v,x) ) == x - 1 ) break;
        if ( j == b ) return false;
    }
    return true;
}
inline void PollardRho(ll x) {
    if ( MillerRabin(x) ) return fac.push_back(x), void();
    for (ll c = 3, x1, x2, i, j; x1 = x2 = i = 1, j = 2; ++c)
        for (ll d; x1 = (Mul(x1,x1,x)+c) % x, d = gcd(abs(x1-x2),x);) {
            if ( 1 < d && d < x ) return PollardRho(d), PollardRho(x/d);
            if ( x1 == x2 ) break; if ( ++i == j ) j <<= 1, x2 = x1;
        }
}
inline ll Phi(ll x) {
    if ( MillerRabin(x) ) return x - 1; ll val = x;
    PollardRho(x), sort( fac.begin(), fac.end() );
    fac.erase( unique(fac.begin(),fac.end()), fac.end() );
    for (int i = 0; i < fac.size(); i++) val /= fac[i], val *= (fac[i]-1);
    return fac.clear(), val;
}
inline ll PrimitiveRoot(ll x,ll val = 0) {
    if ( x == 2 ) return 1; PollardRho( val = Phi(x) );
    sort( fac.begin(), fac.end() );
    fac.erase( unique(fac.begin(),fac.end()), fac.end() );
    for (ll i = 2, flag; flag = 1, i < x; i++) {
        if ( gcd(i,x) != 1 ) continue;
        for (int j = 0; flag && j < fac.size(); j++)
            if ( Pow(i,val/fac[j],x) == 1 ) flag = 0;
        if ( flag ) return fac.clear(), i;
    }

}
inline ll BSGS(ll a,ll b,ll p) {
    ll T = sqrt(p), v1 = 1, v2; map<ll,ll> h;
    for (int i = 0; i < T; v2 = v1 = Mul(v1,a,p), i++) h[Mul(v1,b,p)] = i;
    for (int i = 1; i <= T + 1; v2 = Mul(v2,v1,p), i++)
        if ( h.count(v2) ) return i * T - h[v2];
    return -1;
}
int main(void)
{
    ll k,a,p,g,x0,phi,t; vector<ll> x; int T = 0;
    while ( ~scanf( "%lld%lld%lld", &k, &p, &a ) ) {
        printf( "case%d:\n", ++T );
        g = PrimitiveRoot(p), x0 = BSGS(Pow(g,k,p),a,p), phi = Phi(p);
        if ( x0 == -1 ) { puts("-1"); continue; }
        t = phi / gcd(k,phi), x0 %= t, x.clear();
        for (; x0 < phi; x0 += t) x.push_back(x0);
        for (int i = 0; i < x.size(); i++) x[i] = Pow(g,x[i],p);
        sort( x.begin(), x.end() );
        for (int i = 0; i < x.size(); i++) printf( "%lld\n", x[i] );
    }
    return 0;
}
// HDU3930

基础数论例题

GCD Table

假设符合要求的位置起点为\((x,y)\)，那么显然\(\forall i\in[1,k],a_i|x\)，也就是\(\mathrm{lcm}(a_1,a_2,\cdots ,a_k)|x\)，显然\(x\)可以取\(k\times \mathrm{lcm}(a_1,\cdots,a_k)\)，但是\(k>1\)还要求\(\gcd(k,\frac{y+i-1}{a_i})=1\)，显然\(k>1\)有解，\(k=1\)也有解。

注意到我们还要求\(\forall i\in[1,k],a_i|(y+i-1)\)，所以\(\forall i\in[1,k],y\equiv 1-i\pmod{a_i}\)，用\(\mathrm{ExCRT}\)解这个同余方程组，注意要让\(y\)是正整数。

即使\((x,y),\cdots,(x,y+k-1)\)都落在给定的表格内，也要再次检验是否符合题意，因为保证下标都是对应\(a_i\)的倍数，有可能其\(\gcd\)也是\(a_i\)的倍数，不是\(a_i\)。

仪仗队

\[2+\sum_{i=1}^{n-1}\sum_{j=1}^{n-1}e(\gcd(i,j)) \\ =3+2\sum_{i=1}^{n-1}\sum_{j=1}^{i-1}e(\gcd(i,j))=3+2\sum_{i=1}^{n-1}\varphi(i) \]

Longge的问题

\[\sum_{i=1}^n\gcd(i,n)=\sum_{d|n}d\sum_{i=1}^{\frac{n}{d}}e\left(\gcd\left(i,\frac{n}{d}\right)\right) \\ =\sum_{d|n}d\varphi\left(\frac{n}{d}\right) \]

沙拉公主的困惑

考虑到\(\gcd(x,y)=\gcd(x+y,y)=\gcd(x+ky,y)\)，所以对于\(i\in[1,m!]\)，当一个\(\gcd(i,m!)\)取到\(1\)时，总共会有\(\frac{n!}{m!}\)个\(\gcd(i+k(m!),m!)\)取到\(1\)，所以：

\[\sum_{i=1}^{n!}e(\gcd(i,m!))=\frac{n!}{m!}\sum_{i=1}^{m!}e(\gcd(i,m!))=\frac{n!}{m!}\varphi(m!) \]

外星人

首先要发现只有\(\varphi(2)=1\)，然后考虑唯一分解式取\(\varphi\)之后会发生什么：

\[\varphi\left(\prod p_i^{a_i}\right)=\prod(p_i-1)p_i^{a_i-1} \]

我们发现，每次取\(\varphi\)后原数中的因子\(2\)会因为\(p_i=2\)的部分恰好减少一个，但是会被其他\(p_i\)是奇质数的部分增加若干个。当这个数字\(n\)没有任何因子\(2\)的时候，\(n\)就是\(1\)了，所以取\(\varphi\)函数的次数，就恰好是\(n\)所有质因子的幂会产生的\(2\)的个数之和。

设\(f(x)\)表示\(x\)会产生多少个因子\(2\)，显然\(f(p)=f(p-1)\)，\(f(a\times p)=f(a)\times f(p)(p\in\mathrm{Prime})\)，这一过程在线性筛里面处理即可，那么根据输入就可以直接计算答案了。

要注意，当\(n\)一开始为奇数的时候，第一次取\(\varphi\)没有因子\(2\)可以消去，所以答案要\(+1\)。

上帝与集合的正确用法

\[2^{2^{2^{\cdots}}}\equiv 2^{\left(2^{2^{2^{\cdots}}}\right)\bmod \varphi(p)+\varphi(p)}\pmod p \]

直接递归，当\(\varphi(p)=1\)时，返回\(0\)即可。

数论函数

符号及约定

\(1.\) \([P]\)是指正则表达式，当\(P\)为\(\mathrm{true}\)时，\([P]=1\)，当\(P\)为\(\mathrm{false}\)时，\([P]=0\)。

\(2.\) 约数个数函数定义为\(\tau(n)=\sum_{d|n}1\)。

\(3.\) 约数和函数定义为\(\sigma(n)=\sum_{d|n}d\)。

\(4.\) 元函数定义为\(e(n)=[n=1]\)。

\(5.\) 恒等函数定义为\(I(n)=1\)。

\(6.\) 单位函数定义为\(\epsilon_k(n)=n^k\)。

\(7.\) 欧拉函数定义为\(\varphi(n)=\sum_{i=1}^n[\gcd(i,n)=1]\)。

\(8.\) 设\(n=\prod p_i^{c_i}\)，则莫比乌斯函数定义为\(\mu(n)=\begin{cases}1&n=1\\(-1)^k& \forall c_i=1\\0&\exist c_i>1\end{cases}\)。

\(9.\) 对于数论函数\(f\)，若满足\(\gcd(a,b)=1\)时，有\(f(ab)=f(a)f(b)\)，则称函数\(f\)为积性函数。

\(10.\) 对于两个函数\(f,g\)，定义他们的\(\mathrm{dirichlet}\)卷积为：\((f\times g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})\)，函数\(f,g\)不必要是积性函数。

欧拉筛法

前置知识里已经提到过，以下给出本文可能会使用的线性筛代码：

int flag[N],Prime[N],cnt;
inline void EularSieve(void)
{
    for (int i = 2; i <= Lim; i++) {
        if ( !flag[i] ) Prime[++cnt] = i;
        for (int j = 1; j <= cnt && i * Prime[j] <= Lim; j++) {
            flag[ i * Prime[j] ] = true;
            if ( i % Prime[j] == 0 ) break;
        }
    }
}

这是最朴素的欧拉筛法，不过当我们要筛一些比较复杂的积性函数的时候，线性筛方程的推导可能会非常复杂。

因此我们考虑引入一种更好的线性筛法，目的是能够方便的筛出积性函数的值。

首先，我们注意到积性函数的性质：当\(\gcd(a,b)=1\)时，有\(f(ab)=f(a)f(b)\)。那么对于一个任意的正整数\(n\)，我们可以用算术基本定理进行分解：

\[n=\prod_{i=1}^kp_i^{c_i} \]

此时，任意的\(i,j\in[1,k]\)都满足\(\gcd(p_i^{c_i},p_j^{c_j})=1\)。那么我们就可以得到：

\[f(n)=f\left(\prod_{i=1}^kp_i^{c_i}\right)=\prod_{i=1}^kf(p_i^{c_i}) \]

那么我们就有一个想法：利用定义求出积性函数\(f\)在所有素数幂出的取值，然后直接递推出所有函数值。

具体地说，可以这样递推：

f[1] = 1;
for (int i = 2; i <= Lim; i++)
    if ( i == Pow[i] ) f[i] = Calc(p[i],e[i]);
    else f[i] = f[i/Pow[i]] * f[Pow[i]];

其中\(p_i\)代表数\(i\)的最小素因子，\(e_i\)代表\(i\)的分解式中\(p_i\)这个素因子的指数，\(\mathrm{Calc}\)就是求素数幂处取值的函数，而\(\mathrm{Pow}_i=p_i^{e_i}\)，这样是不是就符合我们的要求了呢？

那么现在我们的问题就是如何求出\(p,e,\mathrm{Pow}\)这三个数组，幸运的是，他们都可以在线性筛的过程中求。

根据欧拉筛法每次用一个数的最小素因子筛去这个合数，就可以更新这三个数组的值了。

代码如下：

int p[N],e[N],Pow[N],Prime[N],cnt;
inline void EularSieve(void)
{
    Pow[1] = p[1] = f[1] = 1 , e[1] = 0;
    for (int i = 2; i <= Lim; i++) {
        if ( p[i] == 0 ) p[i] = Pow[i] = Prime[++cnt] = i , e[i] = 1;
        for (int j = 1; j <= cnt && i * Prime[j] <= Lim; j++) {
            int Next = i * Prime[j]; p[Next] = Prime[j];
            if ( p[i] == Prime[j] ) {
                e[Next] = e[i] + 1;
                Pow[Next] = Pow[i] * Prime[j];
                break;
            }
            else e[Next] = 1 , Pow[Next] = Prime[j];
        }
    }
    for (int i = 2; i <= Lim; i++)
        if ( i == Pow[i] ) f[i] = Calc(p[i],e[i]);
        else f[i] = f[i/Pow[i]] * f[Pow[i]];
}

Dirichlet卷积

Dirichlet卷积的性质

\(1.\) 两个积性函数\(f\)和\(g\)的\(\mathrm{dirichlet}\)卷积仍为积性函数。

证明：

设有两个积性函数\(f\)和\(g\)，则它们的\(\mathrm{dirichlet}\)卷积为：

\[h=f\times g=\sum_{d|n}f(d)g\left(\frac{n}{d}\right) \]

对于函数\(h\)则可以得到：

\[h(x)h(y)=\left (\sum_{d_1|x}f(d_1)g\left(\frac{x}{d_1}\right)\right)\left(\sum_{d_2|y}f(d_2)g\left(\frac{y}{d_2}\right)\right) \\ \ \\=\sum_{d_1|x,d_2|y}f(d_1d_2)g\left(\frac{xy}{d_1d_2}\right)=\sum_{d|xy}f(d)g\left(\frac{xy}{d}\right)=h(xy) \]

故函数\(h\)为积性函数。

\(2.\) \(\mathrm{dirichlet}\)卷积满足交换律。

证明：

设有数论函数\(f\)和\(g\)，则有：

\[f\times g=\sum_{d|n}f(n)g\left(\frac{n}{d}\right)=\sum_{d|n}g(n)f\left(\frac{n}{d}\right)=g\times f \]

\(3.\) \(\mathrm{dirichlet}\)卷积满足结合律。

证明：

设有数论函数\(f\)，\(g\)和\(h\)，则有：

\[(f\times g)\times h=f\times g \times h\\=g\times h \times f=(g\times h)\times f\\=f\times (g\times h) \]

\(4.\) \(\mathrm{dirichlet}\)卷积满足分配律。

证明：

设有数论函数\(f\)，\(g\)和\(h\)，则有：

\[(g+h)\times f=\sum_{d|n}(g(d)+h(d))f\left(\frac{n}{d}\right) \\=\sum_{d|n}g(d)f\left(\frac{n}{d}\right)+\sum_{d|n}h(d)f\left(\frac{n}{d}\right) \\=g\times f+h\times f \]

常见的Dirichlet卷积

\(1.\) \(f\times e=f\)

正确性显然，对于任意函数\(f\)成立.

\(2.\) \(\epsilon=\varphi \times I\)

\[\epsilon=e\times \epsilon=\epsilon\times \mu \times I=\varphi \times I \]

\(3.\) \(\tau=I\times I\)

\[I^2(n)=\sum_{d|n}I(d)I\left(\frac{n}{d}\right)=\sum_{d|n}1=\tau(n) \]

\(4.\) \(\sigma=\epsilon\times I\)

\[(\epsilon\times I)(n)=\sum_{d|n}\epsilon(d)I\left(\frac{n}{d}\right)=\sum_{d|n}d=\sigma(n) \]

\(5.\) \(e=\mu\times I\)

\[(\mu \times I)(n)=\sum_{d|n}\mu(d)I\left(\frac{d}{n}\right)=\sum_{d|n}\mu(d) \]

不妨设\(n=\prod_{i=1}^k p_i^{a_i}\)，\(\mathrm{P}=\{1,p_1,p_2,\cdots,p_{k-1}\}\)，因为当\(x\)有平方因子时\(\mu(x)=0\)，所以可以进行如下转化：

\[(\mu\times I)(n) = \sum_{d|n}\mu(d)=\sum_{S\subseteq \mathrm{P}}\left (\mu\left(\prod_{p\in S}p\right) + \mu\left(p_k\prod_{p\in S}p\right)\right ) \]

对于质数\(p\)和正整数\(a\)满足\(p\not | a\)，显然有\(\mu(p)+\mu(ap)=0\)，所以：

\[(\mu \times I)(n)=\sum_{d|n}\mu(d)=[n=1]=e(n) \]

\(6.\) \(\varphi=\epsilon\times \mu\)

\[(\epsilon \times \mu)(n)=\sum_{d|n}\mu(d)\frac{n}{d}=\varphi(n) \]

可以用欧拉函数的容斥计算来理解.

\(7.\) \(\sigma=\tau \times \varphi\)

\[\sigma =I\times \epsilon=I^2\times \varphi=\tau\times \varphi \]

下取整函数的性质

\(1.\) 对于\(i\in \left [x, \left \lfloor \frac{k}{ \left\lfloor \frac{k}{x}\right \rfloor } \right \rfloor \right ]\)，\(\left\lfloor \frac{k}{i} \right\rfloor\)的值都相等。

证明：

设\(f(x)= \left \lfloor \frac{k}{ \left\lfloor \frac{k}{x} \right\rfloor } \right \rfloor\)，显然有\(f(x)\geq \left \lfloor \frac{k}{ \left( \frac{k}{x} \right) } \right \rfloor=x\)，则可得\(\left \lfloor \frac{k}{f(x)} \right \rfloor\leq \left \lfloor \frac{k}{x} \right \rfloor\)。

从另一方面考虑，则有\(\left \lfloor \frac{k}{f(x)} \right \rfloor\geq\left \lfloor \frac{k}{ \frac{k}{\left \lfloor k/x \right \rfloor } } \right \rfloor=\left\lfloor \frac{k}{x}\right \rfloor\)，则可得：\(\left \lfloor \frac{k}{f(x)} \right \rfloor=\left \lfloor \frac{k}{x} \right \rfloor\)。

\(2.\) \(\left\lfloor\frac{n}{k}\right\rfloor\)最多只有\(2\sqrt n\)种取值。

证明：

当\(k\leq\sqrt n\)时，至多只有\(\sqrt n\)个\(k\)，所以对应的取值只有\(\sqrt n\)种。

当\(k>\sqrt n\)时，有\(1\leq\left\lfloor\frac{n}{k}\right\rfloor\leq\sqrt n\)，所以对应的取值也只有\(\sqrt n\)种。

那么当我们要对某个有关下取整函数的值求和的时候，就可以使用如下的整除分块算法，时间复杂度\(O(\sqrt n)\)。

inline int Calc(void) {
    int res = 0;
    for (int l = 1, r; l <= n; l = r + 1)
        ( r = n/l ? min( n/(n/l) , n ) : n ), res += f(n/l);
    return res;
}

对于二元的情况，其实也是一样的，我们考虑间断点合并，就可以证明其时间复杂度为\(O(\sqrt n + \sqrt m)\)。

inline int Calclim(int n,int k,int l) { return k/l ? min( k/(k/l) , n ) : n; }
inline int Calc(void)
{
    int res = 0;
    for (int l = 1, r; l <= min(n,m); l = r + 1)
        r = min( Calclim( min(n,m), n, l ), Calclim( min(n,m), m, l ) ),
        // 此时，区间[l,r]中所有的 [n/l] 都相等,所有的 [m/l] 也都相等 
        res += f( n/l , m/l );
    return res;
}

\(3.\) 若\(m\)是正整数，则有\(\left\lfloor\frac{\lfloor x \rfloor}{m}\right\rfloor=\left\lfloor\frac{x}{m}\right\rfloor\)

证明：

设\(x=km+r(r<m)\)，则

\[\left\lfloor\frac{\lfloor x \rfloor}{m}\right\rfloor=\left\lfloor k+ \frac{\lfloor r \rfloor}{m}\right\rfloor=k\\ \ \\ \left\lfloor\frac{x}{m}\right\rfloor=\left\lfloor k+ \frac{r}{m}\right\rfloor=k \]

推论 ：

\[\left\lfloor\frac{\left\lfloor k/n \right\rfloor}{m}\right\rfloor=\left\lfloor\frac{k}{nm}\right\rfloor \]

反演原理

一般来说会用到两种反演原理：

\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]=\sum_{i=1}^n\sum_{j=1}^me(\gcd(i,j))=\sum_{i=1}^n\sum_{j=1}^m\sum_{k|\gcd(i,j)}\mu(k) \\ \]

第一个，本质上是\(\mu\times I=e\)，用于一般的莫比乌斯反演.

\[F(n)=\sum_{d|n}f(d)\Longleftrightarrow f(n)=\sum_{d|n}g(d)\mu\left(\frac{n}{d}\right) \]

第二个，本质上是\(F=f\times I\Leftrightarrow f=F\times \mu\)，又称为因数形式莫比乌斯定理（倍数形式莫比乌斯定理基本上可以用反演原理一代替，这里不提了）.

\[F(n)=\sum_{d|n}f(d)\Longleftrightarrow f(n)=F(n)- \sum_{d|n\and d\not = n}f(d) \]

第三个，没什么好说的，就是直接根据定义变形.

数论函数例题：几种常见的反演

NOI2010 能量采集

\(n,m\leq 10^5\)，求：

\[\sum_{i=1}^n\sum_{j=1}^m\left( 2\gcd(i,j)-1\right) \]

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)=\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d] \\ = \sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor m/d\rfloor}[\gcd(i,j)=1] \\ = \sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor m/d\rfloor}\sum_{k|\gcd(i,j)}\mu(k) \\ = \sum_{d=1}^{\min(n,m)}d\sum_{k=1}^{\min(\lfloor n/d\rfloor,\lfloor m/d\rfloor)}\mu(k) \sum_{k|j}^{\lfloor n/d\rfloor}\sum_{k|j}^{\lfloor m/d\rfloor}1 \\ =\sum_{d=1}^{\min(n,m)}d\sum_{k=1}^{\min(\lfloor n/d\rfloor,\lfloor m/d\rfloor)}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor \\ = \sum_{T=1}^{\min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \sum_{d|T}d\times \mu\left(\frac{T}{d} \right) \\ = \sum_{T=1}^{\min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\varphi(T) \]

显然\(\varphi\)函数可以线性筛，那么时间复杂度就优化到\(O(\sqrt n+\sqrt m)\)回答每一组询问.

看到\(\gcd\)，其实可以用欧拉反演秒杀：

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)=\sum_{i=1}^n\sum_{j=1}^m\sum_{k|\gcd(i,j)}\varphi(k) = \sum_{k=1}^{\min(n,m)}\varphi(k)\left\lfloor\frac{n}{k}\right\rfloor\left\lfloor\frac{m}{k}\right\rfloor \]

Luogu2257 YY的gcd

\(T=10^4\)，\(n,m\leq 10^7\)，求：

\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\ \mathrm{is\ a \ prime\ number}] \]

\[\mathrm{LHS}=\sum_{p\ \mathrm{is\ prime\ number}}\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p] \\ \sum_{p}\sum_{i=1}^{\lfloor n/p\rfloor}\sum_{j=1}^{\lfloor m/p\rfloor}[\gcd(i,j)=1]=\sum_{p}\sum_{i=1}^{\lfloor n/p\rfloor}\sum_{j=1}^{\lfloor m/p\rfloor}\sum_{k|\gcd(i,j)}\mu(k) \\ = \sum_{p} \sum_{k=1}^{\min(\lfloor n/p\rfloor,\lfloor m/p\rfloor)}\mu(k)\sum_{k|i}^{\lfloor n/p\rfloor}\sum_{k|j}^{\lfloor m/p\rfloor}1 \\ = \sum_{p} \sum_{k=1}^{\min(\lfloor n/p\rfloor,\lfloor m/p\rfloor)}\mu(k)\left\lfloor\frac{n}{kp}\right\rfloor\left\lfloor\frac{m}{kp}\right\rfloor \\ =\sum_{T=1}^{\min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{p\mathrm{\ is\ a\ prime\ number},p|T}\mu\left(\frac{T}{p}\right) \]

显然我们可以令

\[f(T)=\sum_{p\mathrm{\ is\ a\ prime\ number},p|T}\mu\left(\frac{T}{p}\right) \]

我们完全可以枚举每一个质数\(p\)然后调和级数地枚举其倍数，累加贡献，时间复杂度是\(O(\frac{n}{\ln n}\times \ln n)\approx O(n)\)，当然也可以分类讨论线性筛，不过没什么必要.

BZOJ2693 jzptab

\(T\leq 10^4\)，\(n,m\leq 10^7\)，求：

\[\sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i,j) \]

\[\mathrm{LHS}=\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{\gcd(i,j)}=\sum_{d=1}^{\min(n,m)}\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{d}[\gcd(i,j)=d] \\ = \sum_{d=1}^{\min(n,m)}\sum_{i=1}^{\lfloor n/d \rfloor}\sum_{j=1}^{\lfloor m/d \rfloor}dij[\gcd(i,j)=1]=\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\lfloor n/d \rfloor}\sum_{j=1}^{\lfloor m/d \rfloor}ij[\gcd(i,j)=1] \\ \]

令\(F(n,m)=\sum_{i=1}^n\sum_{j=1}^mij[\gcd(i,j)=1]\)，则：

\[F(n,m)=\sum_{i=1}^n\sum_{j=1}^mij\sum_{k|\gcd(i,j)}\mu(k)=\sum_{k=1}^{\min(n,m)}\mu(k)\sum_{k|i}^{n}\sum_{k|j}^mij \\ = \sum _{k=1}^{\min(n,m)}\mu(k)k^2\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}ij \]

这里记\(S(n,m)=\sum_{i=1}^n\sum_{j=1}^mij=\frac{nm(n+1)(m+1)}{4}\)，那么：

\[F(n,m)=\sum _{k=1}^{\min(n,m)}\mu(k)k^2S\left(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor\right) \\ \mathrm{LHS}=\sum_{d=1}^{\min(n,m)}d\sum_{k=1}^{\min(\lfloor n/d\rfloor,\lfloor m/d\rfloor)}\mu(k)k^2S\left(\left\lfloor\frac{n}{kd}\right\rfloor,\left\lfloor\frac{m}{kd}\right\rfloor\right) \\ =\sum_{T=1}^{\min(n,m)}S\left(\left\lfloor\frac{n}{T}\right\rfloor,\left\lfloor\frac{m}{T}\right\rfloor\right)\sum_{k|T}T\times k\mu(k) \]

显然后面的函数又是积性函数，可以线性筛，那么就可以\(O(\sqrt n+\sqrt m)\)回答一组询问了.