数论杂谈

· · 算法·理论

数论杂谈

欧拉函数

::::info[定义]{open} 欧拉函数，即 \varphi(n)，表示的是小于等于 n 和 n 互质的数的个数． :::: 此时易得 \varphi(1) = 1 ，当 n 为质数时，有 \varphi(n) = n-1.

性质

单个欧拉函数

如果只求单个欧拉函数，完全只需要依照质因数分解来求就可以了就可以了. ::::info[引理] 设 n = \prod_{i=1}^{s}p_i^{k_i}，其中 p_i 是质数，有 \varphi(n) = n \times \prod_{i = 1}^s{\dfrac{p_i - 1}{p_i}}． :::: ::::success[代码]

#include <cmath>

int euler_phi(int n) {
  int ans = n;
  for (int i = 2; i * i <= n; i++)
    if (n % i == 0) {
      ans = ans / i * (i - 1);
      while (n % i == 0) n /= i;
    }
  if (n > 1) ans = ans / n * (n - 1);
  return ans;
}

::::

线性求多个欧拉函数

但是会发现如果直接求多个欧拉函数都用单个欧拉函数求法，时间复杂度就是 O(n\sqrt{n})，但是这太劣了，这是思考线性做法.

显而易见，你会发现线性筛有一个非常好的性质，每个数都会被它最小的质因数筛掉.

for(int i=2;i<=m;i++){
    if(!b[i]) prime[++p]=i;
    for(int j=1;j<=p&&i*prime[j]<=m;j++){
        b[i*prime[j]]=1;
        if(i%prime[j]==0) break;
    }
}

那么分别考虑 prime[j] \mid i 和 prime[j] \nmid i 两种情况.

当 prime[j] \nmid i 时:
很明显 \varphi(prime[j]\times i) =\varphi(prime[j]) \times \varphi(i)

当 prime[j] \mid i 时:

\begin{aligned} \varphi(i \times prime[j]) &= i \times prime[j] \times \prod_{i=1}^{s}\frac{p_i-1}{p_i} \\ &= prime[j] \times \varphi(i) \end{aligned}

那么我们就知道如何用线性筛来求多个欧拉函数了.

::::success[代码]

phi[1]=1;
for(int i=2;i<=n;i++){
    if(!b[i]){
        phi[i]=i-1;
        sum++;
        p[sum]=i;
    }
    for(int j=1;j<=sum;j++){
        if(p[j]*1ll*i>n) break;
        b[i*p[j]]=1;
        if(i%p[j]==0){//当整除时
            phi[i*p[j]]=p[j]*phi[i];
            continue;
        }
        phi[i*p[j]]=(p[j]-1)*phi[i];//当不整除时
    }
}

::::

P2398 GCD SUM

首先我们看到这道题时，应该会想到一个朴素的想法，就是把求 \gcd(i,j) = d 转化为每个 d 有多少个 \gcd(i,j).

\sum_{i=1}^n \sum_{j=1}^n \gcd(i, j) \sum_{d=1}^n \sum_{1 \le i,j \le n} [\gcd(i, j)=d] \sum_{d=1}^n \sum_{1 \le i,j \le \frac{n}{d} } [\gcd(i, j)=1]

由于 i,j 的地位是平等的，不妨设 i < j.

\sum_{d=1}^n \sum_{1 \le j \le \frac{n}{d} }\sum_{1 \le i < j } [\gcd(i, j)=1]

你会发现一个很神奇的事\sum_{1 \le j \le \frac{n}{d} }\sum_{1 \le i < j } [\gcd(i, j)=1] 不就是 \varphi(\frac{n}{d}) 吗（注意 \color{red}j=1），那么直接就可以用线性时间复杂度求出 \sum_{i=1}^n \sum_{j=1}^n \gcd(i, j)了（当然也可以用数论分块来优化成 O(\sqrt{n}). ::::success[代码]

#include<bits/stdc++.h>
using namespace std;
int n,p[20005],phi[100005],sum;
long long q[100005];
long long ans;
bool b[100005];
int main(){
    cin>>n;
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!b[i]){
            phi[i]=i-1;
            sum++;
            p[sum]=i;
        }
        for(int j=1;j<=sum;j++){
            if(p[j]*1ll*i>n) break;
            b[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=p[j]*phi[i];
                continue;
            }
            phi[i*p[j]]=(p[j]-1)*phi[i];
        }
    }
    q[2]=phi[2];
    for(int i=3;i<=n;i++){
        q[i]=q[i-1]+phi[i];//把j=1的情况排除
    }
    for(int i=1;i<=n;i++){
        ans+=q[n/i]*i;
    }
    ans*=2;
    ans+=(1+n)*1ll*n/2;//加上i=j的情况
    cout<<ans;
    return 0;
}

::::

Lucas

::::info[定理]{open} 对于素数 p，有

\binom{n}{k}\equiv \binom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}\binom{n\bmod p}{k\bmod p}\pmod p
其中，当 n<k 时，二项式系数 \dbinom{n}{k} 规定为 0．
::::info[证明]{open}
首先我们知道
$$
(1+x)^p \equiv (1+x^p) \pmod p
$$
:::info[证明]{open}
$$
(1+x)^p= \prod_{i=0}^{p} {p \choose i} x^i
$$
那么很明显当 i \in (0,p) 时 \binom{p}{i} = \dfrac{p!}{i!(p-i)!} 中不可以把分母的 p 约掉
于是只有 i=1 和 i=p 的情况保留了
那么
$$
(1+x)^p \equiv (1+x^p) \pmod p
$$
:::
设 m=ap+b,n=cp+d，其中 (0 \le b,d < p,m \le n).
那么
$
\begin{aligned}
(1+x)^n &=(1+x)^{cp}(1+x)^d \
&=(1+x^p)^c(1+x)^d \pmod p
\end{aligned}
$
考虑 x^m 项的系数得
$
\begin{aligned}
{n \choose m}x^m &\equiv {c \choose a} x^{ap}{d \choose b}x^b \pmod p\
{n \choose m}x^m &\equiv {c \choose a}{d \choose b}x^{ap+b} \pmod p\
{n \choose m}x^m &\equiv {c \choose a}{d \choose b}x^m \pmod p\
{n \choose m} &\equiv {c \choose a}{d \choose b} \pmod p \
{n \choose m} &\equiv {\lfloor \frac{n}{p} \rfloor \choose \lfloor \frac{m}{p} \rfloor}{n \bmod p \choose m \bmod p} \pmod p
\end{aligned}
$
得证
::::
这个定理可以在组合数取模数并不是很大的情况下使用.

P3773 [CTSC2017] 吉夫特

一眼发现是lucas板子题.
使 {n \choose m} \bmod 2=1 时:

{n \choose m} \equiv {\lfloor \frac{n}{2} \rfloor \choose \lfloor \frac{m}{2} \rfloor}{n \bmod 2 \choose m \bmod 2} \pmod 2

那么 {n \bmod 2 \choose m \bmod 2} \equiv 1 \pmod 2.

因为只在 {1 \choose 0} 时为 0.

所以可以说只能 m 在二进制下是 n 的子集，才满足条件.

这很明显直接用枚举子集 dp 就可以了. ::::success[代码]

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
int n,a[211987];
int dp[233334],ans;
inline int read() {
    int x=0,f=1;
    char ch=getchar_unlocked();
    while(ch<'0'||ch>'9') {
        if(ch=='-')
            f=-1;
        ch=getchar_unlocked();
    }
    while(ch>='0'&&ch<='9')
        x=x*10+ch-'0',ch=getchar_unlocked();
    return x*f;
}
int main(){
    n=read();
    for(register int i=1;i<=n;++i){
        a[i]=read();
    }
    for(register int i=1;i<=n;++i){
        for(register int j=a[i]&(a[i]-1);j;j=(j-1)&a[i]){//子集枚举
            dp[j]+=dp[a[i]]+1;
            dp[j]%=mod;
        }
        ans+=dp[a[i]];
        ans%=mod;
    }
    printf("%d",ans);
}

::::

扩展欧几里得

这种算法常用于求 ax+by=\gcd(a,b) 的一组可行解.

::::info[证明]{open} 本证明属于构造法

不妨设：

ax_1 + by_1 = \gcd(a,b) bx_2 + (a \bmod b)y_2 = \gcd(b,a\bmod b)

由 欧几里得定理 可知 \gcd(a,b)=\gcd(b,a\bmod b) 那么 ax_1 + by_1 = bx_2 + (a \bmod b)y_2
又因为 a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b)
所以 ax_1 + by_1 = bx_2 + (a-(\lfloor\frac{a}{b}\rfloor\times b))y_2

ax_1 + by_1 = bx_2 + ay_2 -b\lfloor\frac{a}{b}\rfloor y_2

因为 a=a,b=b，所以 x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2

:::: 函数返回的值为 \gcd，在这个过程中计算 x,y 即可．

long long exgcd(long long a, long long b, long long &x, long long &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

::::info[最小 x 求法] 其实如果想求其中最小的 x ，只需要找到各个 x 之间的关系就行. :::info[推导]{open} 不妨设 ax + by = gcd(a,b) 有两组解，分别为 \{x,y\} 和 \{x_0,y_0\}. 那么 :::align{center}

ax-ax_0+by-by_0=0 a(x-x_0)+b(y-y_0)=0 a(x-x_0)=-b(y-y_0) \frac{a(x-x_0)}{\gcd(a,b)}=-\frac{b(y-y_0)}{gcd(a,b)}

∵\gcd(\frac{a}{\gcd(a,b)},\frac{b}{\gcd(a,b)})=1
∴\frac{b}{\gcd(a,b)} \mid x-x_0
∴\frac{b}{\gcd(a,b)}\times t=x-x_0,t \in \Z
∴对于任意一个 x_i 都有 x_i=x_0+t \times \frac{b}{\gcd(a,b)}

:::

逆元

::::info[定义]{open} 对于非零整数 a , m 如果存在 b 使得 ab \equiv 1 \pmod m，就称 b 是 a 在模 m 意义下的逆元. :::: 这相当于说，b 是线性同余方程 ax \equiv 1 \pmod m 的解．根据 线性同余方程 的性质可知，当且仅当 \gcd(a,m)=1，即 a,m 互素时，逆元 a^{-1}\bmod m 存在，且在模 m 的意义下是唯一的.

快速幂法

::::info[回顾一下费马小定理]{open} 若 p 为素数，a 为正整数，且 a、p 互质。 则有a^{p-1} \equiv 1 \pmod p. ::::

那么可推得 :::align{center}

a\times x \equiv 1 \pmod m a\times x \equiv a^{p-1} \pmod m x \equiv a^{p-2} \pmod m

::: 很明显求 a^{p-2} \bmod m 是容易的，只需要快速幂求就可以. 但是这种代码很明显只能在 p 为质数时候用. ::::success[代码]

long long poww(long long x,int y,int p){//快速幂
    int ans=1;
    while(y){
        if(y&1!=0){
            ans=(ans*x)%p;
        }
        x*=x;
        x%=p;
        y>>=1;
    }
    return ans;
}
long long ny(int a,int p){//逆元
    return poww(a,p-2,p);
}

::::

exgcd法

很容易发现逆元其实是一个特殊的扩展欧几里得. ::::info[推导]{open} :::align{center}

ax \equiv 1 \pmod m ax = 1-my ax+my=1
:::
这种算法不需要考虑 p 的奇偶性.
::::success[代码]
```cpp
long long exgcd(long long a, long long b, long long &x, long long &y) {//扩展欧几里得
if (!b) {
x = 1, y = 0;
return a;
}
long long d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
long long ny(long long a){//逆元
long long d=exgcd(a,(long long)mod,x,y);
return (x%p+p)%p;
}
```
::::

线性求法

当遇到需要连续求 n 个正整数在素数 p 下的逆元时，可以用这个方法. ::::info[推导]{open} 对于 1<i<p 的正整数 i :::align{center}

p =⌊\frac{p}{i}⌋ + (p \bmod i) 0 \equiv ⌊\frac{p}{i}⌋ + (p \bmod i) \pmod p i^{-1} \equiv -⌊\frac{p}{i}⌋ (p \bmod i)^{-1} \pmod p
:::
这样我们只需要从 1^{-1} \bmod p =1 开始，依此把 i^{-1} \bmod p 的问题转化成更小的问题 (p \bmod i)^{-1} .

::::success[代码]

ny[1]=1;//逆元数组
cout<<1<<'\n';
for(int i=2;i<=n;i++){
    ny[i]=(long long)(p-p/i)*ny[p%i]%p;
    cout<<ny[i]<<'\n';
}
::::

### 另一种多个元素的求法

有些会要求快速求出数 $a_1,a_2,\cdots,a_n$ 在模 $m$ 意义下的逆元，保证 $m$ 和 $a$ 里面所有数互素.
::::info[证明]{open}
考虑序列 $\{a_i\}$ 的前缀积：
:::align{center}
$S_0 = 1,~ S_i = a_iS_{i-1},~ i=1,2,\cdots,n$
:::

因 $a$ 里面所有数都与 $m$ 互质，那么 $a$ 的前缀积也与 $m$ 互素。

我们知道求整体乘积的逆元等于所有单个逆元的乘积：
:::align{center}
$S_n^{-1} \times S_n \equiv 1 \pmod m$  
$(a_1\times a_2\times \cdots \times a_n)^{-1}*(a_1\times a_2\times \cdots \times a_n) \equiv 1 \pmod m$
:::
那么我们每次只需要用 $S_i^{-1} \times a_i$ 就可以把 $a_i^{-1}$ 消掉然后得到 $S_{i-1}^{-1}$，从而推出所有的前缀积的逆元.

然后用 $S_{i-1} \times S_i^{-1}$ 就可以把 $S_{i-1}$ 的逆元消掉，只留下 $a_i$ 的逆元.

由此就可以求出所有的逆元.
::::
这样因为只求了一次逆元，时间复杂度为度是 $O(n+\log m)$ .
::::success[代码]
```cpp
q[0]=1;
q[1]=1;
  for(int i=2;i<=p-1;i++){//求前缀和
    q[i]=q[i-1]*i%p;
}
qny[n]=poww(q[n],p-2);
for(int i=n-1;i>=0;i--){//求前缀积逆元
    qny[i]=qny[i+1]*(i+1)%p;
}
for(int i=1;i<=n;i++){//求单个数逆元
    ny[i]=qny[i]*q[i-1];
}

::::

容斥

好(屎)题？

AT_abc465_f Sjeltzer?

P3813 [FJOI2017] 矩阵填数

P1316 Mivik 写书

附录

数论分块

再推式子时，如果遇到了一个形如 \sum_{i=1}^nf(\lfloor \frac{n}{i} \rfloor)时，就可以快速的直接用数论分块做到 O(\sqrt{n}).

想法是直接考虑有多少个不同的 \lfloor \frac{n}{i} \rfloor，很容易发现当 i\le \sqrt{n} 时有 \sqrt{n} 种可能，当 i>\sqrt{n} 时 \lfloor \frac{n}{i} \rfloor < \sqrt{n} 肯定也不超过 \sqrt{n} 种取值.

注意到由于 \lfloor \frac{n}{i} \rfloor 单调递减，那么这些相同的值肯定在连续的下标上，那么只需要求出这个连续下标的 l 和 r 就可以了，并且显然易见，下一个 l 为现在的 r+1.

那么考虑知道 l 之后 r 怎么求，显然 \lfloor \frac{n}{l} \rfloor=\lfloor \frac{n}{r} \rfloor.

我们知道 k \le \frac{n}{r} < k+1， 那么 \frac{n}{k} \ge r > \frac{n}{k+1}，那么 r 最大值为 \lfloor \frac{n}{k} \rfloor ，很明显 k=\lfloor \frac{n}{l} \rfloor.

代码如下：

inline int Calc(int n){
    int Ans=0;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        Ans+=(SF(r)-SF(l-1))*g(n/l);
    }
    return Ans;
}

如果要考虑给更高维的数论分块，相当于把 r 取两个端点的 \min ，可证明 m 维数论分块是 O（m\sqrt{n}）

P2260 [清华集训 2012] 模积和

先钦定 n \le m，把式子化为 \sum_{i=1}^{n}(n \bmod i) \sum_{j=1}^{m} (m \bmod j)-\sum_{i=1}^{n}(n \bmod i)(m \bmod i)

考虑到 n \bmod i 为 n - i{\lfloor\frac{n}{i} \rfloor},那么就可以将式子化为两个数论分块，其中有一个是二维的。

::::success[代码]

#include<bits/stdc++.h>
using namespace std;
const int mod=19940417;
long long sum1,sum2,ans;
long long n,m;
int main(){

    cin>>n>>m;
    if(n>m)
        swap(n,m);
    sum1=n*1ll*n%mod;
    sum1%=mod;
    for(long long l=1,r;l<=n;l=r+1){//一维数论分块
        r=n/(n/l);
        sum1+=-1ll*(r+l)*(r-l+1)/2%mod*(n/r)%mod;
        sum1%=mod;
    }

    sum1=(sum1+mod)%mod;
    sum2=m*1ll*m%mod;
    sum2%=mod;
    for(long long l=1,r;l<=m;l=r+1){//把两个一维数论分块分开求
        r=m/(m/l);
        sum2+=-1ll*(r+l)*(r-l+1)/2%mod*(m/l)%mod;
        sum2%=mod;
    }
    sum2=(sum2+mod)%mod;
    ans=sum1*sum2%mod-n*1ll*m%mod*n%mod;
    ans=(ans+mod)%mod;
    for(long long l=1,r;l<=n;l=r+1){//二位数论分块
        r=min(n/(n/l),(m/(m/l)));
        ans+=(r+l)*1ll*(r-l+1)/2%mod*(n*1ll*(m/r)%mod+m*1ll*(n/r)%mod)%mod;
        ans%=mod;
        ans-=(r*1ll*(r+1)%mod*(2*r+1)%mod*3323403%mod-(l-1)*1ll*l%mod*(2*l-1)%mod*3323403%mod+mod)%mod*(m/l)%mod*(n/l)%mod;
        ans%=mod;
        ans=(ans+mod)%mod;
    }
    ans=(ans+mod)%mod;
    cout<<ans;
    return 0;
}

::::