快速冪和快速乘

原理：二進位制分解和取模的定義

\[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\boldsymbol{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 \boldsymbol{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 \boldsymbol{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 \boldsymbol{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 \boldsymbol{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\boldsymbol{N}\)，那麼原問題的所有解為：\(x\equiv g^{c+\frac{\varphi(p)}{\gcd(\varphi(p),k)}\times i}\pmod p,i\in \boldsymbol{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)\)回答一組詢問了.