0x31 质数

定义

若一个数，只有一和它本身两个因子，那么这个数就是一个质数

在自然数集中，小于$n$的质数约有$ln(n)$个

试除法

试除法是常用判断素数的方法

inline bool is_prime( int x )
{
    if( x < 2 ) return 0;
    for( register int i = 2 ; i * i <= x ; i ++ )
    {
        if( x % i == 0 ) return 0;
    }
    return 1;
}

素数的筛选

Eratosthenes 筛法（埃拉托色尼筛法）

每次只用素数去筛复杂度$O(nlog_{log_{n}})$

const int N = 1005;
int n,prime[N];

inline void primes()
{
    for(register int i = 2;i <= n;i ++)
    {
        if(prime[i]) continue;
        for(register int j = 2; i * j <= n; j++) prime[i*j] = 1;
    }
    return ;
}

线性筛（欧拉筛法）

每次只用一个数用小于当前这个数最小质因子的质数去筛其他数，即保证每个数都被自己的最小质因子筛掉

const int N = 10005;
int n,prime[N],cnt;
bool vis[N];

inline void primes()
{
    for(register int i = 2;i <= n;i ++)
    {
        if(!vis[i]) prime[++cnt] = i;
        for(register int j = 1;j <= cnt && i * prime[j] <= n; j ++)
        {
            vis[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

质因数分解

算数基本定理

任何一个大于1的数都可以被分解成有限个质数乘积的形式

试除法

类似埃式筛，我们直接枚举影子然后把当前因子全部除尽即可

分解成 $p_{1}\times _{2}\times p_{3}\times \cdots \ p_{n}$这种形式

const int N = 1005;
int p[N];

inline int factorize(int x)
{
    register int cnt = 0;
    for(register int i = 2;i * i <= x;i ++)
    {
        while(x % i == 0)
        {
            p[++cnt] = i;
            x /= i; 
        }
    }
    if(x > 1) p[++cnt] = x;
    return cnt;
}

分解成 $p_{1}^{k_{1}} \times p_{2}^{k_{2}} \times p_{3}^{k_{3}} \times \cdots \ p_{n}^{k_{n}}$

const int N = 1005;
int p[N],power[N];


inline int factorize(int x)
{
    register int cnt = 0;
    for(register int i = 2;i *i <= x;i ++)
    {
        if(x%i) continue;
        p[++cnt] = i;
        while(x%i == 0) x /= i,power[cnt] ++;
    }
    if(x == 1) goto end;
    p[++cnt] = x;
    power[cnt] = 1;
    end : return cnt
}

AcWing 196. 质数距离

这道题数据的范围非常的大，我们没有办法在一秒内求出所有质数

但是我们知道一个合数$x$在一定是一个小与$\sqrt{x}$的质数的倍数

扫描二维码关注公众号，回复： 7861304 查看本文章

所以我们可以求出$(1\cdots \sqrt{U})$的所有质数，然后对于每个区间做埃式筛

然后暴力遍历一遍区间即可

#include <bits/stdc++.h>
#define LL long long

using namespace std;

const int N = 1000010;

int prime[N], cnt , p[N] , tot;
bitset< N > v;

inline void primes()
{
    int n = 50000;
    for( register int i = 2 ; i <= n ; i ++ )
    {
        if( !v[i] ) prime[ ++ cnt ] = i;
        for( register int j = 1 ; j <= cnt && i * prime[j] <= n ; j ++ )
        {
            v[ i * prime[j] ] = 1;
            if( i % prime[j] == 0 ) break;
        }
    }
}


int main()
{
    long long l, r;
    primes();
    while (cin >> l >> r)
    {
        v.reset();
        for (int i = 1; i <= cnt; i ++ )
        {
            int t = prime[i];
            // 把[l, r]中所有t的倍数筛掉
            for (long long j = max((l + t - 1) / t * t, 2ll * t); j <= r; j += t)
                v[j - l] = true;
        }

        tot = 0;
        for (int i = 0; i <= r - l; i ++ )
            if (!v[i] && i + l > 1)
                p[tot ++ ] = i + l;

        if (tot < 2) puts("There are no adjacent primes.");
        else
        {
            int minp = 0, maxp = 0;
            for (int i = 0; i + 1 < tot; i ++ )
            {
                int d = p[i + 1] - p[i];
                if (d < p[minp + 1] - p[minp]) minp = i;
                if (d > p[maxp + 1] - p[maxp]) maxp = i;
            }

            printf("%d,%d are closest, %d,%d are most distant.\n", p[minp], p[minp + 1], p[maxp], p[maxp + 1]);
        }
    }
    return 0;
}

AcWing 197. 阶乘分解

$N!$中质因子$p$的个数等于$1\cdots N$中每个数的质因子$p$的个数之和

包含一个质因子p的数显然都是$p$的倍数，所以有$\left \lfloor \frac{N}{p}\right \rfloor$个质因子

同理包含第二个（不是包含两个）质因子$p$的数显然都是$p^2$的倍数，所以有所以有$\left \lfloor \frac{N}{p^2}\right \rfloor$个质因子，注意不是$2\times\left \lfloor \frac{N}{p^2}\right \rfloor$个，因为第一个已经统计过了

所以$N!$中一共包含$\sum_{k=1}^{p^k\le N}\lfloor\frac{N}{p^k}\rfloor$个质因子

所以我们先求出所有的质因子，再用以上方法求出所有的质因子数即可

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


const int N = 1000010;
int n , prime[N] , tot;
bool v[N];

inline void primes()
{
    for( register int i = 2 ; i <= n ; i ++ )
    {
        if( !v[i] ) prime[ ++ tot ] = i;
        for( register int j = 1 ; j <= tot && i * prime[j] <= n ; j ++ )
        {
            v[ i * prime[j] ] = 1;
            if( i % prime[j] == 0) break;
        }
    }
    return ;
}


int main()
{
    cin >> n;
    primes();

    for (int i = 1; i <= tot; i ++ )
    {
        register int p = prime[i] , cnt = 0;
        for ( p ; p <= n ; p *= prime[i] )
        {
            cnt += n / p;
            if( p > n / prime[i] ) break;
            //等价于 if( p * prime[i] > n ) break; 防止溢出
        }
        cout << prime[i] << ' ' << cnt << endl;
    }

    return 0;
}

0x32 约数

定义

若整数$n$除以整数$x$的余数为$0$，即$d$能整除$n$，则称$d$是$n$的约数，$n$是$d$的倍数，记为$d|n$

算术基本定理推论

由算数基本定理得正整数N可以写作$N=p_1^{C_1}\times p_2^{C_2} \times p_3^{C_3} \cdots \times p_m^{C_m}$

N的正约数个数为($\Pi$是连乘积的符号，类似$\sum$)

\[ (c_1+1)\times (c_2+1)\times \cdots (c_m+1)=\Pi_{i=1}^{m}(ci+1) \]

$N$的所有正约数和为
\[ (1+p_1+p_1^2+\cdots +p_1^{c_1})\times\cdots\times(1+p_m+p_m^2+\cdots +p_m^{c_m})=\prod_{i=1}^{m}(\sum_{j=0}^{c_i}(p_i)^j) \]

求$N$的正约数的集合

对于任意的$d|n$则$(\frac{n}{d})|n$

所以只要扫描$1\cdots\sqrt n$就能找到n的所有正约数，复杂度$O(\sqrt n)$

int factor[N] , tot = 0;
for( int i = 1 ; i * i <= n ; i ++ )
{
    if( n % i ) continue;
    factor[ ++ tot] = i;
    if( i != n / i ) factor[ ++ tot ] = n / i;
}

AcWing 198. 反素数

引理1

$1\cdots N$中最大的反素数就是约数个数最多的数中最小的一个

证明：

设$m$是$1\cdots N$中约数个数最多的数中最小的一个。根据$m$的定义，$m$满足

$\forall x < m ,g(x)\le g(m) $
$\forall x>m,g(x) \le m$

第一条性质说明$m$是反素数，第二条性质说明大于$m$的都不是反素数

引理2

$1\cdots N$中任何数的不同质因子都不超过$10$个，任何质因子的指数总和不超过$30$

证明：

最小的$11$个质因子乘积$2\times3\times5\times7\times11\times13\times17\times19\times23\times29\times31>2\times10^9$

最小的质数的$31$次方$2^{31}>2\times10^9$

引理3

$\forall x \in [1,N]$，x为反素数的必要条件是：

x分解质因数后可写作$2^{c_1}\times3^{c_3}\times5^{c_3}\times7^{c_4}\times11^{c_5}\times13^{c_6}\times17^{c_7}\times19^{c_8}\times23^{c_9}\times29^{c_{10}}$

并且$c_{1}\geq c_{2}\geq c_{3}\geq c_{4}\geq c_{5}\geq c_{6}\geq c_{7}\geq c_{8}\geq c_{9}\geq c_{10}\geq0 $

证明：

反正法，由引理$2$，$x$的质因数分解式中存在一项$p^k(p>29)$，则必定有一个不超过29的质因子${p}’$不能整除$x$。根据算数基本定理的推论，$\frac{x}{p^k}\times{p}'^k$的约数个数与$x$的约数个数相等，但前者更小，所以，这与反质数的定义矛盾。故$x$只包含$29$以内的质因子

同理，若$x$的质因子不是连续若干最小的，或者质数不单调递减，我们可以通过上述交换质因子的方法，的到一个比$x$更小的、但约数个数相同的整数。因此假设不成立，原命题成立

综上所述，我们可以用$DFS$确定前十个质数的指数，并满足指数单调递减，总乘积不超过N，同时记录约数个数

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const int N = 2e9+5 , t = 10 , p[t] = { 2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 21 , 29 };
int n , sum = 0 , minx;


inline void dfs( int u , int last , int pi , int s )
{
    if( s > sum || s == sum && pi < minx ) sum = s , minx = pi;
    for( register int i = 1 ; i <= last ; i ++ )
    {
        if( (LL)pi * p[u] > n ) break;
        pi *= p[u];
        dfs( u + 1 , i , pi , s * ( i + 1 ) );
    }
    return ;
}


int main()
{
    cin >> n;
    dfs( 0 , 30 , 1 , 1 ); 
    cout << minx << endl;
}

AcWing 199. 余数之和

首先注意到$k \mod i = k- \left \lfloor \frac{k}{x} \right \rfloor$,故可以转化为$n\times k - \sum_{i=1}^{n}\left \lfloor \frac{k}{i} \right \rfloor\times i$

对于任意的$x\in[1,k]$，设$g(x)= \left \lfloor k/ \left \lfloor \frac{k}{x}\right \rfloor \right \rfloor$

因为
\[ \left\lfloor\frac{k}{x} \right\rfloor \le \frac{k}{x}\Rightarrow g(x)\ge \left\lfloor \frac{k}{k /x}\right\rfloor = x \]
所以可得
\[ \left\lfloor\frac{k}{g(x)}\right\rfloor\le\left\lfloor\frac{k}{x}\right\rfloor \]

又因为
\[ g(x)\le\frac{k}{\lfloor k/x \rfloor}\Rightarrow\left\lfloor \frac{k}{g(x)}\right\rfloor\ge\left\lfloor \frac{k}{k/\lfloor\frac{k}{x}\rfloor}\right\rfloor = \left\lfloor \frac{k}{x}\right\rfloor \]
所以可得
\[ \left\lfloor \frac{k}{g(x)}\right\rfloor=\left\lfloor \frac{k}{x}\right\rfloor \]
所以可知$\forall i \in[x,g(x)]$，$\lfloor \frac{k}{i}\rfloor$恒为定值

所以在区间[x,g(x)]中$k\mod i$的值是个等差数列，利用高斯公式$O(1)$求和

对于$\forall x\in[1,n]$ ，$\left\lfloor \frac{k}{x}\right\rfloor$只有$2\sqrt{n}$中结果，相当于把$[1,n]$分成了$2\sqrt{n}$段，每段$\left\lfloor \frac{k}{x}\right\rfloor$都相同

每一段的余数和可以直接用公式求得，所以最终的复杂度就是$O(\sqrt{n})$

int main()
{
    long long n , k , ans;
    cin >> n >> k;
    ans = n * k;
    for( register int x = 1 , gx ; x <= n ; x = gx + 1 )
    {
        gx = ( k / x ) ? min( k / ( k / x ) , n ) : n;
        ans -= (k/x) * ( x + gx ) * ( gx - x + 1 ) / 2;
    }
    cout << ans << endl;
    return 0;
}

最大公约数

定义

若自然数$d$同时是$a$和$b$的约数，则称$d$是$a$和$b$的公约数

在所有的公约数中最大的一个就是最大公约数，记作$gcd(a,b)$

若自然数$m$同时是$a$和$b$的倍数，则成$m$是$a$和$b$的公约数

在所有的公倍数中最小的一个就是最小公倍数，记住$lcm(a,b)$

同理，我们也可以定义三个数以及更多数的最大公约数或最小公倍数

定理

\[ \forall a,b \in N,gcd(a,b)\times lcm(a,b)=a\times b \]

证明：

设$d=gcd(a,b),a_0=\frac{a}{d},b_0=\frac{b}{d}$

根据定义得$gcd(a_0,b_0)=1,lcm(a_0,b_0)=a_o\times b_0$

所以$lcm(a,b)=lcm(a_0\times d,b_0\times d)=lcm(a_0,b_0)\times d=a_0\times b_0 \times d = \frac{a \times b }{d} $

更相减损术

\[ \forall a,b \in N , a > b,有gcd(a,b)=gcd(b,a-b)=gcd(a,a-b) \\ \forall a,b \in N , 有gcd(2a,2b) = 2gcd(a,b) \]

证明

根据最大公约数的定义，后者显然成立，我们主要证明前者

对于任意的公约数$d$，因为$d|a,d|b$,所以$d|(a-b)$因此d也是$b,a-b$的公约数

反之亦然成立，故$a,b$的公约数集合与$b,a-b$的公约数集合相同

欧几里得算法

\[ \forall a , b \in N , b \ne 0, gcd(a,b) = gcd(b , a \mod b) \]

原理其实和更相减损术是一样的，做多次减法直到减不下实际上就等价于做一次取摸

inline int gcd( int a , int b ) { return b ? gcd( b , a % b ) : a; }

欧几里得算法是常用的求最小公倍数的算法，但是如果遇到需要取模和高精度的时候，考虑到高精度实现取模运算比较复杂，可以考虑用更相减损术来代替欧几里得，并且时间复杂度是相同的

Luogu P1072 Hankson 的趣味题

为了方便叙述我们把$a_0,a_1,b_0,b_1$重新定义为$a,b,c,d$

首先可以最小公倍数的定义可知$x|d$

所以我们可以直接枚举$d$的因子即可

对于最小公倍数的判断可以进行一些变形
\[ \because lcm(a,b)\times gcd(a,b)=a\times b\\ \therefore lcm(x,c)=d\Leftrightarrow d\times gcd(x,c)= x \times c \]

int main()
{
    T = read();
    while( T -- )
    {
        a = read() , b = read() , c = read() , d = read() , ans = 0;
        for( register int i = 1 , j ; i * i <= d ; i ++ )
        {
            if( d % i ) continue;
            if( gcd( i , a ) == b && d * gcd( i , c ) == i * c ) ans ++;
            j = d / i;
            if( gcd( j , a ) == b && d * gcd( j , c ) == j * c ) ans ++;
        }
        printf( "%d\n" , ans );
    }
    
    return 0;
}

互质与欧拉函数

定义
\[ \forall a,b\in N 若gcd(a,b)=1,则称a,b互质 \]
对于三个数或更多的数，$gcd(a,b,c)=1 \Leftrightarrow gcd(a,b) = gcd(a,c)=gcd(b,c)=1$

欧拉函数

$1\cdots N$中与N互质的数的个数，被称为欧拉函数，记作$\phi(N)$

由算数基本定理得$N= p_{1}^{k_{1}} \times p_{2}^{k_{2}} \times p_{3}^{k_{3}} \times \cdots \ p_{m}^{k_{m}}$

所以$\phi(N) = N \times \frac{p_1-1}{p_1}\times \frac{p_2-1}{p_2}\times \cdots \times \frac{p_m-1}{p_m} = N \times \Pi_{p|N}(1-\frac{1}{p})$

证明

假设$N$有两个质因子$p,q$

在$1\cdots N$中$p$的倍数有$\frac{N}{p}$个，同理$q$的倍数也有$\frac{N}{q}$个

我们把$\frac{N}{p}$和$\frac{N}{q}$排除掉，显然$p\times q$ 的倍数被排除了两次所以要加上$\frac{N}{p\times q}$
\[ \phi(N) = N - \frac{N}{p}-\frac{N}{q}+\frac{N}{p\times q} = N\times(1-\frac{1}{p}-\frac{1}{q}+\frac{1}{p\times q}) = N\times(1-\frac{1}{p})\times(1-\frac{1}{q}) \]
根据数学归纳法可以把上诉结论扩展至任意个质因子

所以我们可以在做质因数的同时，求出欧拉函数

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

性质1~2

$\forall n > 1 , 1\cdots n中与互质的数的和为n\times \phi(n) / 2 $
若a,b互质，则$\phi(ab)=\phi(a)\phi(b)$

积性函数

如果当a,b互质时，满足$f(ab)=f(a)f(b)$的函数$f$称为积性函数

性质3~6

若f是积性函数，且在算数基本定理中$n=\Pi_{i=1}^{m} p_i^{c_i}$，则$f(n) = \Pi_{i=1}^{m} f(p_i^{c_i})$
设$p$为质数，若$p|n$且$p^2|n$ , 则$\phi(n) = \phi( \frac{n}{p})\times p$
设$p$为质数，若$p|n$且$p^2 \not|n$,则$\phi(n) = \phi(\frac{p}{n})\times (p-1)$
$\sum_{d|n}\phi(d)=n$

性质4~6，只适用于欧拉函数，并非所有的积性函数都适用

$O(n)$递推求欧拉函数

const int N = 1e3+5;
int prime[N] , phi[N] , tot ;
bool v[N];

for( register int i = 2 ; i <= n ; i ++ )
{
    if( !v[i] ) prime[ ++ tot ] = i , phi[i] = i - 1;
    for( register int j = 1 ; j <= tot && i * prime[j] <= n ; j ++ ) 
    {
        v[ i * prime[j] ] = 1;
        phi[ i * prime[j] ] = phi[i] * ( i % prime[j] ? prime[j] - 1 : prime[j] );
        if( i % prime[j] == 0 ) break; 
    }   
}

根据性质4、性质5得

补充证明

若$i|p_j$ ，则$i=k\times p_j$

所以$n=i\times p_j = k\times p_j^2$

故$p_j|n$,且$p_j^2|n$

[SDOI2008]仪仗队

JZYZOJ P1377

AcWing 201. 可见的点和本题基本相同，不单独讲

首先我们将原图从$(1,1)$到$(n,n)$重新标号$(0,0)$到$(n-1，n-1)$

分析题目可知，当$n=1$时，显然答案是$0$，特判断即可

我们考虑 $n \ne1$的情况，首先$(0,1),(1,0),(1,1)$这三个点是一定能看到的

考虑剩余的点

对于任意的点$(x,y)，2\le x,y\le n-1$，若不被其它的点挡住必定满足$gcd(x,y) = 1$

证明

若点$(x,y)$ 不满足$gcd(x,y)=1$，令$d = gcd(x,y)$

则$(x,y)$与$(0,0)$连线的斜率$k=\frac{y}{x} = \frac{y/d}{x/d}$

所以$(x,y)$会被$(x/d,y/d)$挡住

然后我们把这个图按照对角线分成两个等腰直角三角形

若点$(x,y)$在一个三角形中，并满足$gcd(x,y)=1$

则必有一个点$(y,x)$在另一个三角形中,并满足$gcd(y,x) = 1$

所以只统计其中一个三角形即可

现在我们在$y<x$的三角形中考虑

对于任意一个$x$,满足条件的$y$的个数就是$\phi(x)$

所以答案就是
\[ 3+2\times\sum_{i=2}^{n-1}\phi(i) \]
所以们$O(n)$的地推欧拉函数的同时求和即可

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

const int N = 40005;
 
int n , sum , prime[N] , phi[N] , tot ;
bool v[N];

int main()
{
    cin >> n;
    if( n == 1 ) puts("0") , exit(0);
    n --;
    for( register int i = 2 ; i<= n ; i ++ )
    {
        if( !v[i] ) prime[ ++ tot ] = i , phi[i] = i - 1;
        sum += phi[i]; 
        for( register int  j = 1 ; j <= tot && i * prime[j] <= n ; j ++ )
        {
            v[ i * prime[j] ] = 1;
            phi[ i * prime[j] ] = phi[i] * ( i % prime[j] ? prime[j] - 1 : prime[j] );
            if( i % prime[j] == 0 ) break;
        }
    } 
    cout << sum * 2 + 3 << endl;
}

0x33同余

定义

若正整数$a$和$b$除以$m$的余数相等，则称$a$，$b$模$m$同余，记作$a \equiv b(mod \ m)$

性质

自反性：$a\equiv a (mod \ m)$
对称性：若$a\equiv b(mod \ m)$，则$b \equiv a (mod\ m)$
传递性：若$a \equiv b(mod\ m) , b\equiv c(mod\ m)$，则$a\equiv c(mod\ m)$
同加性：若$a\equiv b(mod\ m)$，则$a+c\equiv b+c(mod\ m)$
同乘性：若$a\equiv b(mod\ m)$，则$a\times c \equiv b\times c(mod\ m)$，若$a\equiv b(mod\ m) , c\equiv d (mod\ m)$，则$a\times c \equiv b\times d (mod\ m)$
同幂性：若$a\equiv b(mod\ m)$，则$a^c\equiv b^c(mod\ m)$
$a\times b\%\ k=(a\%k)\times(b\%k)\%k$
若$a\%p=x,a\%q=x,gcd(p,q)=1$，则$a\%(p\times q) = x$
不满足同除，若$a\equiv b(mod\ m)$不满足$a\div c\equiv b\div c(mod\ m)$

费马小定理

若$p$是质数，则对于任意正整数$a$都有 $a^p\equiv a(mod\ p)$或$a^{p-1}\equiv 1(mod\ p)$

欧拉定理

若正整数$a,n$互质，则$a^{\phi(n)}\equiv 1(mod\ n)$其中$\phi(n)$是欧拉函数

费马小定理是欧拉定理的一种特殊情况

欧拉定理推论

若正整数$a,n$互质，则对于任意的正整数$b$，有$a^b\equiv a^{b\ mod \ \phi(n)}(mod \ n)$

在一些计数的问题中，常常要求对结果取模，但面对乘方时，无法直接取模，可以先把底数对$p$取模，指数对$\phi(p)$取模，再计算乘方

特别的，当不保证$a,n$互质时，若$b>\phi(n)$，有$a^b\equiv a^{b\ mod \ \phi(n)+\phi(n)}(mod \ n)$

这也就意味着即使不保证互质，我们也可把乘方运算缩小到容易计算的规模

AcWing 202. 最幸运的数字

$x$个8连在一起组成的整数可写作$8(10^x-1)/9$，题目就是要求我们求最小的$x$满足$L|8(10^x-1)/9$

设$d = gcd(L,8)$
\[ L|\frac{8(10^x-1)}{9}\Leftrightarrow 9L|8(10^x-1)\Leftrightarrow \frac{9L}{d}|10^x-1\Leftrightarrow 10^x\equiv1(mod\ \frac{9L}{d}) \]
引理：

若正整数$a,n$互质，则满足$a^x\equiv 1(mod n)$的最小正整数$x_0$是$\phi(n)$的约数

证明

反证法：假设满足$a^x\equiv 1(mod n)$的最小正整数$x_0$不是$\phi(n)$的约数

设$\phi(n)=qx_0+r,(0<r<x_0)$，因为$a^{x_0}\equiv1(mod\ n)$，所以$a^{qx_0}\equiv 1(mod\ n)$

又因欧拉定理得$a^{\phi(n)} \equiv 1(mod\ n)$，所以可得$a^{r}\equiv 1(mod \ n)$

假设不成立，所以假设错误

证毕

根据以上定理，我们只需求出$\phi(\frac{9L}{d})$,枚举他所有的约数，用快速幂逐一检查条件即可

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const LL INF = 1e18;
LL l, mod , ret , p ;

inline int read()
{
    register int x = 0;
    register char ch = getchar();
    while( ch < '0' || ch > '9' ) ch = getchar();
    while( ch >= '0' && ch <= '9' ) x = ( x << 3 ) + ( x << 1 ) + ch - '0' , ch = getchar();
    return x;
}

inline int gcd( int a , int b ) { return b ? gcd( b , a % b ) : a ; }

inline LL power(LL x , LL  p )
{
    register LL ans = 1;
    while( p )
    {
        if( p & 1 ) ans = ( ans * x ) % mod;
        x = ( x * x ) % mod;
        p >>= 1; 
    }
    return ans;
}

inline int phi( int x )
{
    register int ans = x;
    for( register int i = 2 ; i * i <= x ; i ++ )
    {
        if( x % i ) continue;
        ans = ans / i * ( i - 1 );
        while( x % i == 0 ) x /= i;
    }
    if( x > 1 ) ans = ans / x * ( x - 1 );
    return ans;
}

inline void work()
{
    mod = l / gcd( 8 , l ) * 9 , p = phi( mod ) ,ret = INF;
    for( register LL i = 1 ; i * i <= p ; i ++)
    {
        if( p % i ) continue;
        if( power( 10 , i  ) ==  1 ) ret = min( ret , i );
        if( power( 10 , p / i ) == 1 )  ret = min( ret , p / i ); 
    }
    return ;
}

int main()
{
    for( register int i = 1 ; i ; i ++  )  
    {
        l = read();
        if( l == 0 ) break;
        work() , printf( "Case %d: %d\n" , i , (ret == INF ? 0 : ret ) );
    } 
    return 0;
}

裴蜀定理

对于任意的$a$,$b$存在无穷多组整数对$(x,y)$满足不定方程$ax+by=d$,其中$d=gcd(a,b)$

$ax+by=d \ (d=gcd(a,b))$这个方程叫丢番图方程

扩展欧几里得算法

在求$gcd(a,b)$的同时，可以求出关于$x,y$的不定方程$ax+by=d$的一组解

考虑递归计算：假设已经算出了$(b\%a,a)$的一组解$(x_{0} , y_{0})$满足$(b\%a)x_0+ay_0=d$

可以得到$(b-a\lfloor \frac{b}{a} \rfloor)x_0+ay_0 = d$

整理得到$a(y_0 - \lfloor \frac{b}{a} \rfloor x_0) + bx_0 = d$

为了方便表示$swap(x_0,y_0)$得$a(x_0 - \lfloor \frac{b}{a} \rfloor y_0) + by_0 = d$

就可以得到一组解$((x_0 - \lfloor \frac{b}{a} \rfloor y_0),y_0)$

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

扩展欧几里得一元线性同余方程

解关于$x$的一元线性同余方程$ax\equiv b(mod\ m)$

已知扩展欧几里得可以解$ax_0+by_0=d,(d = gcd(a,m))$

将原方程变形得$ax + my = b$

通过扩展欧几里得可解出$ax^\prime + my^\prime = d$

当$b\%d==0$时$ax + my = b$有解,反之无解

方程$ax^\prime + my^\prime = d$两边同乘$\frac{b}{d}$得$a\frac{b}{d}x^\prime+m\frac{b}{d}y^\prime = b$

我们就把原方程解出来$x = \frac{b}{d}x^\prime$

乘法逆元

在之前我们提到过，除法不能直接取模，但是可以用逆元

若$x$满足$ax\equiv 1(mod \ n )$，则称$x$是$a$在模$n$意义下的逆元记作$a^{-1}$

当然满足条件的逆元不止一个，通常情况下我们只用最小正整数的逆元，并且逆元也不在任何条件下都有逆元

求逆元

下面介绍三种最常用的求逆元的方法

费马小定理

费马小定理 $a^{p-1} \equiv 1(mod\ p)$,其中$p$为素数

根据费马小定理变形得$aa^{p-2}\equiv 1(mod\ p)$

则$a^{p-2}\% p$就是逆元，直接快速幂求得

inline int invers(int x,int p)
{
    return ksm(x,mod-2,p) % p;
}

扩展欧几里得

根据逆元的定义$ax \equiv 1(mod\ m)$得$ax+my=1$

解这个丢番图方程可得$x\%m$就是逆元，同时可得逆元存在的条件$gcd(a,m) = 1$

inline int invers(int a,int mod)
{
    register int x,y;
    exgcd(a,mod,x,y);
    return (x % mod + mod) % mod;
}

线性递推

$a^{-1} =\displaystyle -\left \lfloor \frac{m}{a}\right\rfloor\times(m\%a)^{-1}$

const int N = 1005;
int inv[N] = {};

inline void invers(int n,int mod)
{
    inv[1] = 1;
    for(register int i = 2;i <= n;i ++)
            inv = (long long)(mod - mod / i) * inv[mod % i] % mod;
    return ;
}

线性同余方程

给定整数$a,b,m$，求一个整数$x$满足$a\times x\equiv b(mod \ m)$，或者无解。因为未知数的指数为一，所以成为一次同余方程或线性同余方程

解关于$x$的一元线性同余方程$ax\equiv b(mod\ m)$

已知扩展欧几里得可以解$ax_0+by_0=d,(d = gcd(a,m))$

将原方程变形得$ax + my = b$

通过扩展欧几里得可解出$ax^\prime + my^\prime = d$

当$b\%d==0$时$ax + my = b$有解,反之无解

方程$ax^\prime + my^\prime = d$两边同乘$\frac{b}{d}$得$a\frac{b}{d}x^\prime+m\frac{b}{d}y^\prime = b$

我们就把原方程解出来$x = \frac{b}{d}x^\prime$

中国剩余定理CRT

又称中国余数定理或孙子定理，最早可见于中国南北朝时期（公元5世纪）的数学著作《孙子算经》卷下第二十六题，叫做“物不知数”问题，原文如下：

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

简单来说就求一元线性同余方程组
$\left\{\begin{matrix} x \equiv a_1(mod\ m_1)\\x \equiv a_2(mod\ m_2)\\\vdots\\x \equiv a_n(mod\ m_n) \end{matrix}\right.$

对于这样的一个方程组如果任意的两个$m$都互素，则一定有解
所以中国剩余定理就是解但$m$互素的情况

我们尝试来构造这样一组解

设$M = m_1 \times m_2 \times \cdots \times m_n,M_i = \frac{M}{m_i}$

$M_i$只有在模$m_i$时不为零，模其他$m$都是零

由于我们希望在模$m_i$时得到$a_i$所以乘上$a_it_i$

其中$t_i$是$M_i$在模$m_i$时的逆元

所以最小正整数解就是$x = \sum_{i=1}^{n}a_it_iM_i$，通解就是$x+kM$

如果题目说明了的$m$都是素数则可以用费马小定理

inline int crt()
{
    register int ans = 0,M = 1;
    for(register int i = 1;i <= n;i ++) M *= m[i];
    
    for(register int i = 1;i <= n;i ++)
    {
        register int Mi = M/m[i],invers = ksm(Mi,m[i]-2,m[i]);
        ans = (ans + a[i]*invers*Mi) % M;
    }
    return ans;
}

如果题目没有说明我们还是用扩展欧几里得来求逆元

inline int crt()
{
    register int M = 1, ans = 0;
    for(register int i = 1;i <= n; i ++) M *= m[i];
    for(register int i = 1;i <= n; i ++)
    {
        register LL Mi = M/m[i],x,y;
        exgcd(Mi,m[i],x,y);
        ans = ((ans + a[i] * Mi * x) % M + M) % M;
    }
    return ans;
}

证明中国剩余定理一定有解

我们举这个同余方程为例

$\left\{\begin{matrix} x \equiv a_1(mod\ m_1)\\x \equiv a_2(mod\ m_2)\\ \end{matrix}\right.$

假设我们分别解出

$\left\{ \begin{matrix} x_1 = a_1+k_1 \times m_1 \\ x_2 = a_2 + k_2 \times m_2 \end{matrix}\right.$

因为 $x_1 = x_2$

所以$a_1+k_1\times m_1 = a_2 + k_2 \times m_2$

整理得$k_2\times m_2 - k_1\times m_1 = a_1-a_2$

这就是$ax+by=d$的形式

所以当$gcd(m_1,m_2)|a_1-a_2$ 时一定有解

又因为$m_1$与$m_2$互质

所以$gcd(m_1,m_2) = 1$

显然中国剩余定理一定有解

扩展中国剩余定理EXCRT

扩展中国剩余定理就是解一元线性同余方程组

$\left\{\begin{matrix} x \equiv a_1(mod\ m_1)\\x \equiv a_2(mod\ m_2)\\\vdots\\x \equiv a_n(mod\ m_n) \end{matrix}\right.$

但是这里的$m$不要求互素，这是与中国剩余定理不同的地方

假设们我已经求出前$k-1$个等式的最小整数解$x$

令$M = LCM(m_1,m_2,\cdots,m_{k-1})$

则前$k-1$个等式的通解为$x+tM\ (t\in Z)$

我们要找$x+t^\prime M \equiv a_k (mod\ m_k)$

新的解$x^\prime = x + t^\prime M$

所以解$n$个方程就可以求出最后的解

怎么解$x+t^\prime M \equiv a_k (mod\ m_k)$呢?

我们可以化成$t^\prime M \equiv a_k - x(mod\ m_k)$

我们可以用扩展欧几里得解出$t^\prime$,这样我们就构造出呢前$k$个方程的一组特解

inline int excrt()
{
    register int x = a[1],M = m[1],t,y,c,d,mod; 
    for(register int i = 2;i <= n;i ++)
    {
        c = ((a[i] - x) % m[i] + m[i]) % m[i]; d = exgcd(M,m[i],t,y);
        if(c%d) return -1;
        mod = m[i] / d;
        t = ((t * (c/d)) % mod + mod) % mod;
        x += t * M;
        M *= mod;
        x = (x + M) % M;
    }
    return x;
}

对于第一个方程$x=a_1 , M= m_1$
$c = a_i - x$
$x \% \frac{M\times m_i}{d} = (x\%\frac{m_i}{d})\%M$ 所以$mod = \frac{m_i}{d}$

Luogu P1082 同余方程

非常非常落的模板题，记得用long long

#define LL long long

int main()
{
    LL a , b ;
    cin >> a >> b;
    LL x , y;
    exgcd( a , b , x , y );
    cout << ( x % b + b ) % b << endl;
    return 0;
}

Luogu P4777 扩展中国剩余定理

本题有弱化版AcWing 204. 表达整数的奇怪方式

非常裸的$excrt$，记住用快速乘

0x34 矩阵乘法

矩阵

数学上，一个$n\times m$的矩阵是一个有$m$行$n$列的元素排列成的矩形阵列。矩阵里的元素可以是数字、符号或数学式。一下是一个由$6$个元素构成的$2$行$3$列的矩阵
\[ \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \]
大小相同（行列数都相同）的矩阵可以相加减，具体做法是每个对应位置上的元素相加减

特殊的矩阵

方阵

如果一个行数和列数都相同的矩阵，我们称作方阵，下面是一个$3\times3$的方阵
\[ \begin{bmatrix}1 & 2 & 3\\4&5&6\\7&8&9\end{bmatrix} \]

单位矩阵

对于$n\times n$的方阵，如果主对角线上的元素为$1$，其余元素全为$0$我们称为$n$阶单位矩阵，一般用$In$或$En$表示，通常情况下简称$I$或$N$，下面是一个$4$阶单位矩阵
\[ \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix} \]
主对角线：从左上到右下的对角线

反对角线：从右上到左下的对角线

矩阵乘法

详细的矩阵乘法的运算过程可以参考这个视频
\[ \begin{bmatrix} 3 & 1 & -1\\ 2 & 0 & 3 \end{bmatrix} \times \begin{bmatrix} 1&6\\ 3&-5\\ -2&4\\ \end{bmatrix} = \begin{bmatrix} 3\times 1+1\times3+(-1)\times(-2)&2\times1+0\times3+3\times(-2)\\ 3\times6+1\times(-5)+(-1)\times4&2\times6+0\times(-5)+3\times4 \end{bmatrix} =\begin{bmatrix} 8&9\\-4&24\end{bmatrix} \]
定义

设$A$是$n\times m$的矩阵$B$是$m\times p$ 的矩阵，则$C=A\times B$是$n\times p$的矩，并且
\[ \forall i\in [1,n], \forall j\in[1,p]：c_{i,j}=\sum_{k=1}^{m}A_{i,k}\times B_{k,j} \]
也就是说，参与矩阵乘法的矩阵第一个矩阵的列数必须等于第二个矩阵的行数，否则不能运算

矩阵乘法满足结合律$(A\times b)\times C = A \times (B \times C)$，满足分配律$(A+B)\times C = A\times C + B \times C $，但不满足交换律$A\times B \ne B \times A$

方阵乘幂

方阵乘幂是指，$A$是一个方阵，将$A$连乘$n$次，即$C = A^n$

若不是方阵则不能进行方阵乘幂

由于矩阵乘法满足结合律，所以可以用快速幂的实现来实现矩阵快速幂

矩阵乘法的应用

矩阵乘法奥妙在于：

很容易将有用的状态储存在一个矩阵中
通过状态矩阵和状态转移矩阵相乘快速的得到一次$DP$的值（这个$DP$的状态转移方程必须是一次的）
求矩阵乘法的结果是要做很多次乘法，这样的效率非常慢甚至比不上直接递推进行$DP$的状态转移，但由于矩阵乘法满足结合律，可以先算后面的转移矩阵，即用快速幂，迅速的处理好后面的转移矩阵，再用初始矩阵乘上后面的转移矩阵得到结果，复杂度是$O( log(n) ) $的

Loj 10219.矩阵 A×B

矩阵乘法的模板题

const int N = 105;
int n , m , q , a[N][N] , b[N][N] , c[N][N];

int main()
{
    n = read() , m = read();
    for( register int i = 1 ; i <= n ; i ++ )
    {
        for( register int j = 1 ; j <= m ; j ++ ) a[i][j] = read();
    }
    q = read();
    for( register int i = 1 ; i <= m ; i ++ )
    {
        for( register int j = 1 ; j <= q ; j ++ ) b[i][j] = read();
    }
    //代码的核心矩阵乘法
    for( register int i = 1 ; i <= n ; i ++ )
    {
        for( register int j = 1 ; j <= q ; j ++ )
        {
            for( register int k = 1 ; k <= m ; k ++ ) c[i][j] += a[i][k] * b[k][j];
        }
    }
    for( register int i = 1 ; i <= n ; i ++ )
    {
        for( register int j = 1 ; j <= q ; j ++ ) printf( "%d " , c[i][j] );
        puts("");
    }
    return 0;
}

Loj 10220.Fibonacci 第 n 项

本题看似简单，但巨大的数据范围开数组肯定开不下，并且线性求解效率也很低

我们知道
\[ f[i]=1\times f[i-1] + 1\times f[i-2]\ (1)\\ f[i-1] = 1\times f[i - 1]+0\times f[i-2]\ (2) \]
观察上面两个式子，我们可以把他改写成矩阵乘法的形式
\[ \begin{bmatrix}f[i]\\f[i-1]\end{bmatrix}=\begin{bmatrix}1&1\\1&0\end{bmatrix}\times\begin{bmatrix}f[i-1]\\f[i-2]\end{bmatrix}=\begin{bmatrix}1&1\\1&0\end{bmatrix}\times \begin{bmatrix}1&1\\1&0\end{bmatrix}\times\begin{bmatrix}f[i-2]\\f[i-3]\end{bmatrix} \]
根据数学归纳法们就能得到
\[ \begin{bmatrix}f[n]\\f[n-1]\end{bmatrix}=\begin{bmatrix}1&1\\1&0\end{bmatrix}^{n-2}\times\begin{bmatrix}f[2]\\f[1]\end{bmatrix} \]
之前我们提到过矩阵可以快速幂，所以我们就可以在$O(log(n))$的复杂度内求出$f[n]$

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const int N = 5;
LL n , m ;

struct matrix
{
    int x , y ;
    LL v[N][N];
    inline friend matrix operator * ( matrix a , matrix b )
    {
        register matrix cur;
        cur.x = a.x , cur.y = b.y;
        for( register int i = 1 ; i <= a.x ; i ++ )
        {
            for( register int j = 1 ; j <= b.y ; j ++ )
            {
                cur.v[i][j] = 0;
                for( register int k = 1 ; k <= a.y ; k ++ )
                {
                    cur.v[i][j] += a.v[i][k] * b.v[k][j];
                    cur.v[i][j] %= m; 
                }
            }
        }
        return cur;
    }
}f , res , cus ;


inline void matrix_I( matrix & k , int t )
{
    k.x = k.y = t;
    for( register int i = 1 ; i <= t ; i ++ )
    {
        for( register int j = 1 ; j <= t ; j ++ )
        {
            if( i == j ) k.v[i][j] = 1;
            else k.v[i][j] = 0;
        }
    }
    return ;
}

inline matrix matrix_power( matrix x , int k )
{
    register matrix t; matrix_I( t , 2 );
    while( k )
    {
        if( k & 1 ) t = t * x;
        x = x * x;
        k >>= 1;
    }
    return t ;
}

inline LL fib()
{
    if( n <= 2 ) return 1;
    f.x = f.y = 2;
    f.v[1][1] = f.v[1][2] = f.v[2][1] = 1 , f.v[2][2] = 0;
    res.x = 2 , res.y = 1;
    res.v[1][1] = res.v[2][1] = 1;
    f = matrix_power( f , n - 2 );
    res = f * res;
    return res.v[1][1] ;
}


int main()
{
    cin >> n >> m;
    cout << fib() << endl;
    return 0;
}

Loj 10221.Fibonacci 前 n 项和

本题在上一题的基础上增加了求和，我们依旧要构造状态转移矩阵
\[ s[n] = 1 \times s[n-1] + 1 \times f[n] + 0 \times f[n-1]\\ f[n+1] = 0 \times s[n-1] + 1 \times f[n] + 1 \times f[n-1]\\ f[n] = 0 \times s[n-1] + 1 \times f[n] + 0 \times f[n-1] \]
所以我么可以根据上面的三个式子推出状态转移矩阵
\[ \begin{bmatrix} s[n] \\ f[n+1] \\ f[n] \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} s[n-1] \\ f[n] \\ f[n-1] \end{bmatrix} \]
然后就依旧是矩阵快速幂

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const int N  = 5;
int n , m ;

struct matrix
{
    LL v[N][N];
    int x , y;
    inline friend matrix operator * ( matrix a , matrix b )
    {
        register matrix cur;
        cur.x = a.x , cur.y = b.y;
        for( register int i = 1 ; i <= cur.x ; i ++ )
        {
            for( register int j = 1 ; j <= cur.y ; j ++ )
            {
                cur.v[i][j] = 0;
                for( register int k = 1 ; k <= a.y ; k ++ ) cur.v[i][j] = ( cur.v[i][j] + a.v[i][k] * b.v[k][j] ) % m;
            }
        }
        return cur;
    }
} res , f;


inline void matrix_I( matrix & k , int t )
{
    k.x = k.y = t;
    for( register int i = 1 ; i <= t ; i ++ )
    {
        for( register int j = 1 ; j <= t ; j ++ ) k.v[i][j] = ( i == j ) ? 1 : 0;
    }
    return ;
}

inline matrix matrix_power( matrix x , int k )
{
    register matrix t ; matrix_I( t , 3 );
    while( k )
    {
        if( k & 1 ) t = t * x;
        x = x * x;
        k >>= 1;
    }
    return t;
}

inline LL fib()
{
    if( n == 1 ) return 1;
    if( n == 2 ) return 2;
    f.x = f.y = res.x = 3 , res.y = 1; 
    f.v[1][1] = f.v[1][2] = f.v[2][2] = f.v[2][3] = 1 , f.v[3][2] = res.v[1][1] = res.v[2][1] = res.v[3][1] = 1;
    f.v[1][3] = f.v[2][1] = f.v[3][1] = f.v[3][3] = 0 ;
    f = matrix_power( f , n - 1 );
    res = f * res;
    return res.v[1][1];
}


int main()
{
    cin >> n >> m;
    cout << fib() << endl;
}

Lougu P1939 矩阵加速

矩阵快速幂的模板题，推到过程很简单
\[ a[i] = 1\times a[i-1] + 0\times a[i-2]+1\times a[i-3]\\ a[i-1] = 1\times a[i-1] + 0\times a[i-2]+0\times a[i-3]\\ a[i-2] = 0\times a[i-1] + 1\times a[i-2]+0\times a[i-3]\\ \Downarrow \\ \begin{bmatrix}f[i]\\f[i-1]\\f[i-2]\end{bmatrix}=\begin{bmatrix}1&0&1\\1&0&0\\0&1&0\end{bmatrix}\times \begin{bmatrix}f[i-1]\\f[i-2]\\f[i-3]\end{bmatrix} \]
直接暴力矩阵快速幂是可以过的，但有没有优化呢？

我们考虑离线来做这道题

首先我们把所有的数据读入，然后离散化，并排序

对与$ans_i$我们可以计算$d=n_i-n_i$得到一个差值

然后只需将$ans_{i-1}$乘上状态转移矩阵$f$的$d$次方，就是$ans_i$

可以通过这种操作来减少多次重复的运算

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const int N = 105 , M = 5 , mod = 1e9 + 7;
int n , key[N] , val[N] , ans[N];

struct matrix
{
    int x , y;
    LL v[M][M];
    
    matrix() { memset( v , 0 , sizeof(v) ); }
    
    inline friend matrix operator *  (matrix a , matrix b )
    {
        register matrix cur;
        cur.x = a.x , cur.y = b.y;
        for( register int i = 1 ; i <= cur.x ; i ++ )
        {
            for( register int j = 1 ; j <= cur.y ; j ++ )
            {
                cur.v[i][j] = 0;
                for( register int k = 1 ; k <= a.y ; k ++ ) cur.v[i][j] = ( cur.v[i][j] + a.v[i][k] % mod * b.v[k][j] % mod ) % mod;
            }
        }
        return cur;
    }
} f , res;


inline LL read()
{
    register LL x = 0;
    register char ch = getchar();
    while( ch < '0' || ch > '9' ) ch = getchar();
    while( ch >= '0' && ch <= '9' )
    {
        x = ( x << 3 ) + ( x << 1 ) + ch - '0';
        ch = getchar();
    }
    return x;
}

inline void matrix_I( matrix &t , int k)//构造单位矩阵
{
    t.x = t.y = k;
    for( register int i = 1 ; i <= k ; i ++ ) t.v[i][i] = 1;
}

inline matrix matrix_power( matrix x , int k )//矩阵快速幂
{
    if( k == 1 ) return x;
    register matrix t; matrix_I( t , 3 );
    while( k )
    {
        if( k & 1 ) t = t * x ;
        x = x * x;
        k >>= 1;
    }
    return t;
}

inline void init()
{
    f.x = f.y = 3;
    memset( f.v , 0 , sizeof( f.v ) );
    f.v[1][1] = f.v[1][3] = f.v[2][1] = f.v[3][2] = 1 ;
}


int main()
{
    n = read();
    for( register int i = 1 ; i <= n ; i ++ ) key[i] = val[i] = read();
    sort( val + 1 , val + 1 + n ) ;
    for( register int i = 1 ; i <= n ; i ++ ) key[i] = lower_bound( val + 1 , val + 1 + n , key[i] ) - val;
    res.x = 3 , res.y = res.v[1][1] = res.v[2][1] = res.v[3][1] = 1;
    val[0] = 3;
    for( register int i = 1 ; i <= n ; i ++ )
    {
        if( val[i] <= 3 ) { ans[i] = 1 , val[i] = 3 ; continue; }
        register int d = val[i] - val[ i - 1 ];//计算差值
        if( d == 0 ) //特殊情况
        {
            ans[i] = ans[ i - 1 ] ;
            continue;
        }
        init();//初始化 矩阵f 
        f = matrix_power( f , d );
        res = f * res;
        ans[i] = res.v[1][1] % mod;
    }
    for( register int i = 1 ; i <= n ; i ++ ) printf( "%d\n" , ans[ key[i] ] );
    return 0;
}

0x35组合数学

加法原理

若完成一件事有$n$类，第$i$类有$a_i$中方法，那么完成这件事共有$a_1+a_2+\cdots+a_n$种方法

乘法原理

若完成一件事$n$个步骤，第$i$个步骤有$a_i$种方法，那么完成这件事共有$a_1\times a_2\times \cdots\times a_n$种方法

排列数

假设有这样一个问题：现在有甲乙丙丁四个小朋友，老师想从中挑出两个小朋友排成一列，有几种排法。
很容易枚举出来，有12种排法：

甲乙，甲丙，甲丁，乙甲，乙丙，乙丁，丙甲，丙乙，丙丁，丁甲，丁乙，丁丙

但问题假如变成从$66666$个小朋友中挑出$2333$个小朋友有几种排法，显然枚举是不现实的。我们需要一个通用的计算方法。

排列就是指这类问题。它的定义是：从$m$个不同的元素中取出$n$ ($0\leq n\leq m$)个元素排成一列的方法数。记作$A_{m}^{m}$。
我们可以一步一步的考虑这个问题，首先从$m$个小朋友中抽取一个放在第一个有$m$种方法。之后从剩下的$m-1$个中挑一个放在第二个即有$m-1$中方法。以此类推，根据乘法原理我们可以得到
\[ A_{m}^{n} = m \times (m-1) \times (m-2) \times \cdots \times (m-n+1) \]
通过化简我们可以得到
\[ A_{m}^{n} = \frac{m!}{(m-n)!} \]
PS：阶乘$n! = 1 \times 2 \times 3 \times \cdots \times n$。特别的我们规定$0! = 1$。

组合数

如果老师只是先从4个小朋友中挑出2人组成一组，不用考虑排列顺序呢么有几种呢？很显然通过枚举得知有6种

甲乙，甲丙，甲丁，乙丙，乙丁，丙丁

组合数就是指这类问题。它的定义是：从$m$个不同的元素中取出$n$ ($0\leq n\leq m$)个元素组成一组的方法数。记作$C_{m}^{m}$。注意与排列不同的是不考虑顺序。
呢么我们就要重新推到，类比刚才的结论我们得知排列只是将组合的每种结果进行了不同的排序。反而言之对于组合来说，排序是重复的。重复的情况是每个组合的排序。呢么对于每个组合重复了多少种情况呢。相当于从$n$个小朋友种抽出$n$个小朋友排一对，也就是
\[ A_{n}^{n} = \frac{n!}{(n-n)!} = \frac{n!}{0!} = n!\]
\[ 排序是对组合的每种情况进行排序，根据乘法原理 \]
A_{m}^{n} = C_{m}^{n} \times A_{n}^{n}
\[ 所以我们得到 \]
C_{m}^{n} = \frac{m!}{(m-n)! \times n!}
$$

性质：

$C_n^m=C_n^{n-m}$
$C_n^m=C_{n-1}^{m}+C_{n-1}^{m-1}$
$C_n^0+C_n^1+\cdots +C_n^n=2^N$
$C_n^0+C_n^2+C_n^4+\cdots = C_n^1+C_n^3+C_n^5\cdots=2^{n-1}$

组合数的计算

加法递推$O(n^2)$

利用性质2，边界条件$C_n^0=C_n^n=1$，其实就是杨辉三角

for( register int i = 0 ; i <= n ; i ++ )
{
    c[i][0] = 1;
    for( register int j = 1 ; j <= i ; j ++) c[i][j] = c[i-1][j] + c[i-1][j-1];
}

乘法地推$O(n)$

$C_n^m=\frac{n-m+1}{n}\times C_n^{m-1}$，边界条件$C_n^0=1$。必须先乘后除，不然可能除不开

可以利用性质$1$，减少部分运算

c[0] = 1;
for( register int i = 1 ; i * 2 <= n ; i ++ ) c[i] = c[n-i] = ( n - i + 1 ) * c[ i - 1 ] / i;

生成全排列

一些题目可能会给定一串数字，要求找到其中的某一个顺序，暴力的做法就是生存这个序列的全排列，然后暴力判断

这里提供一个$STL$来生成全排列next_permutation

next_permutation是生成当前序列的下一个序列，直到整个序列单调递减

int a[3];

do
{
    for( register int i = 0 ; i < 3 ; i ++ ) cout << a[i] << ' ';
    cout << endl;
}while( next_permutation( a , a + 3 ) );

当然next_permutation也可以只针对某一段生成排列，修改参数即可，同时也支持自定义排序方式，类似sort

注意next_permutation操作的序列必须是单调递增，不然只能生成当前序列的下一个序列，无法生成所有的序列

枚举子集

如果用一个数的二进制来表示一个集合，那么就可以用下面的方法生成当前集合的所有子集，并能做到不重不漏

for( register int i = s ; i ; i = ( i - 1 ) & s );

二项式定理

\[ (a+b)^n=\sum_{k=0}^{n}C_n^ka^{n-k}b^{k} \]

NOIP2011计算系数

二项式定理的模板题直接计算即可

#include <bits/stdc++.h>
#define LL long long
using namespace std;


const int N = 1e6+5;
const LL mod = 10007;
LL a,b,k,n,m,result = 1,t[N];


inline LL gcd(LL a,LL b) {return !b?a:gcd(b,a%b);}

inline LL ksm(LL a,LL b)
{
    register LL ans = 1;
    while(b)
    {
        if(b & 1) ans = (ans * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return ans % mod;
}


int main()
{
    cin >> b >> a >> k >> n >> m;
        
    result = (result * ksm(b,n)) % mod;
    result = (result * ksm(a,m)) % mod;
    
    n = min(n,m);
    
    for(register int i = 1,j = k - n + 1;i <= n;i ++,j ++) t[i] = j;
    
    for(register int i = 2;i <= n;i ++)
    {
        register LL x = i,y;
        for(register int j = 1;j <= n;j ++)
        {
            if(x > t[j]) y = gcd(x,t[j]);
            else y = gcd(t[j],x);
            if(y == 1) continue;
            x /= y; t[j] /= y;
            if(x == 1) break;
        }
    }
    
    for(register int i = 1;i <= n;i ++) result = (result * t[i]) % mod;
    
    cout << result << endl;
}

Catalan数列

给定$n$个$0$和$n$个$1$,他们按照某种顺序排列成长度为$2n$的序列，满足任意前缀中$0$的个数都不少于$1$的个数的序列的数量为
\[ Cat_n=\frac{1}{n+1}\times C_{2n}^{n} \]
推论：

$n$个左括号和$n$个右括号组成的合法序列为$Cat_n$
$1,2,3,\cdots n$经过一个栈形成的合法出栈顺序为$Cat_n$
$n$个节点构成不同的二叉树的数量为$Cat_n$
在平面直角

Lucas定理

若$p$是质数，则对于任意的整数$1\le m \le n$，有：
\[ C_n^m=C_{n\mod p}^{m\mod p}\times C_{n/p}^{m/p}(\mod P) \]

Loj 10228. 组合

卢卡斯定理的模板题，组合数根据定义直接求

Luogu P3807卢卡斯定理小细节不同

#include <bits/stdc++.h>
#define LL long long 
using namespace std;

LL p;

inline int read()
{
    register int x = 0;
    register char ch = getchar();
    while( ch < '0' || ch > '9' ) ch = getchar();
    while( ch >= '0' && ch <= '9' )
    {
        x = ( x << 3 ) + ( x << 1 ) + ch - '0';
        ch = getchar();
    }
    return x;
}

inline LL power( LL x , LL y )
{
    register LL res = 1;
    while( y )
    {
        if( y & 1 ) res = ( res * x ) % p;
        x = ( x * x ) % p;
        y >>= 1;
    } 
    return res;
} 

inline LL inv ( LL x ) { return power( x , p - 2 ) % p ; }

inline LL C ( LL n , LL m )
{
    if( m > n )  return 0;
    LL a = 1 , b = 1;
    for( register LL i = n - m + 1 ; i <= n ; i ++ ) a = a * i % p;
    for( register LL i = 2 ;i <= m ; i ++ ) b = b * i % p;
    return a * inv( b );
}

inline LL lucas( LL n , LL m )
{
    if( !m ) return 1;
    return C( n % p , m % p ) * lucas( n / p , m / p ) % p;
}

int main()
{
    for( register int T = read() , n , m ; T >= 1 ; T -- ) n = read() , m = read() , p = read() , printf( "%lld\n" , lucas( n , m ) ); 
    return 0;
}

0x3e补充证明

part 1.

已知$t = x+y,a^x\equiv1(mod\ n),a^t \equiv 1(mod \ n)$

求证$a^y\equiv 1(mod\ n)$

证明
\[ 由已知式子得a^x=k_1 \times n+1\\ 设a^y=k_2\times n + p\\ 则a^t=a^x\times a^y=k_1k_2n^2+(tk_1+k_2)n+t\\ 又因为a^t\% n=1,所以t=1\\ 所以a^y\equiv1(mod\ n) \]

证明素数是无数个

假设数是$n$个，每个素数是$p_i$

令$P = \Pi_{i=1}^{n} pi + 1$

因为任何一个数都可以分解成多个质数相乘

所以$P$除以任何一个质数都余$1$

显然$P$就也是一个质数

与假设矛盾，所以假设错误

所以质数是无限个

$P = \Pi_{i=1}^{n} pi + 1$不是质数的通项公式

举例$ n = 2 $

则$P = 2 \times 3 + 1 = 7$

显然第三个质数是$5$

所以这个式子只能用于证明素数是无限个，并不是质数的通项公式

目前人类还没有发现质数的通项公式，只有近似公式

求斐波那契数列通项公式

已知$F_{n+2} = F_{n+1} + F_n$ 其中$F_1 = F_2 = 1$ ,求$F_n$通项公式

首先令首项为$1$公比为$x$的等差数列$\{ x^n \}$ 满足$F_n$的递推式

代入递推式得$x^{n+2}=x^{n+1} + x^{n}$

整理得$(x^2-x-1)x^n = 0$

因为$x^n\neq 0 $

所以$x^2-x-1=0$

解得$x_1=\frac{1+\sqrt{5}}{2}$或$x_2=\frac{1-\sqrt{5}}{2}$

所以数列${x_1^n}$和${x_2^n}$都满足$F_n$的递推公式，但${x_1^n}$和${x_2^n}$都不是斐波那契数列的通项公式

构造$F_n = ax_1^n+bx_2^n$

证明

$F_{n+2} = ax_1^{n+2} + bx_2^{n+2} =a(x_1^nx_1^2) +b(x_2^nx_2^2)$

因为$x_1^2=x_1+1,x_2^2=x_1+1$

所以$F_{n+2}=ax_1^n(x_1+1)+bx_2^n(x_2+1) = ax_1^{n+1}+bx_2^{n+1}+ax_1^n+bx_2^n =F_{n+1}+F_n$

所以$F_n = ax_1^n+bx_2^n$

已知$F_1=1,F_2=1$，代入得

$\left \{\begin{matrix} ax_1+bx_2=1 \\a^2x_1^2+b^2x_2^2=1 \end{matrix} \right .$

解得 $a=\frac{1}{\sqrt{5}},b=-\frac{1}{\sqrt{5}}$

所以
\[ F_n=\frac{\sqrt{5}}{5}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n] \]

0x30 数学

0x31 质数

定义

试除法

素数的筛选

Eratosthenes 筛法 （埃拉托色尼筛法）

线性筛 （欧拉筛法）

质因数分解

算数基本定理

试除法

0x32 约数

定义

算术基本定理推论

求\(N\)的正约数的集合

最大公约数

定义

定理

更相减损术

欧几里得算法

互质与欧拉函数

\(O(n)\)递推求欧拉函数

0x33同余

定义

性质

费马小定理

欧拉定理

欧拉定理推论

裴蜀定理

扩展欧几里得算法

扩展欧几里得一元线性同余方程

乘法逆元

求逆元

费马小定理

扩展欧几里得

线性递推

线性同余方程

中国剩余定理CRT

证明中国剩余定理一定有解

扩展中国剩余定理EXCRT

0x34 矩阵乘法

矩阵

特殊的矩阵

方阵

单位矩阵

矩阵乘法

方阵乘幂

矩阵乘法的应用

0x35组合数学

加法原理

乘法原理

排列数

组合数

组合数的计算

加法递推\(O(n^2)\)

乘法地推\(O(n)\)

生成全排列

枚举子集

二项式定理

Catalan数列

Lucas定理

0x3e补充证明

part 1.

证明素数是无数个

求斐波那契数列通项公式

猜你喜欢

Eratosthenes 筛法（埃拉托色尼筛法）

线性筛（欧拉筛法）