zhugezy的完美数学教室(v2.0)

你找到了这个菜鸡的板子。。。
2019/09/20:马上要ccpc了，更新一下，过十几个小时应该还会加一点东西（？）

zhugezy的完美数学教室

一些结论

\(2C_n^2 + n = n^2,C_n^m=C_{n - 1}^m+C_{n-1}^{m-1}\)

\(\sum_{i = 1}^niC_n^i=n2^{n-1},\sum_{i=1}^ni^2C_n^i=n(n+1)2^{n-2}\)

\(C_n^{x+1}+2\sum_{i=0}^xC_n^i=\sum_{i=0}^{x+1}C_{n+1}^i\)

\(\sum_{i=1}^n(-1)^{i-1}\frac{C_n^i}{i}=\sum_{i=1}^n\frac{1}{i}\)

\(\sum_{i=0}^n(C_n^i)^2=C_{2n}^n\)

1~n与n互质的数的和：\(\frac{n*\phi(n)+[n==1]}{2}\)

\(\sum_{i=1}^ngcd(i,n)=\sum_{d|n}d\phi(n/d)\)

1_{p模p的所有逆元值对应1}p中所有的数，例如p=7时，1~6分别对应1,4,5,2,3,6

$1 \leq i \leq n, \lfloor\frac{n}{i}\rfloor $最多只有\(2\sqrt{n}\)种取值

\(\lfloor {\frac{n}{\lfloor{\frac{n}{i}}\rfloor}}\rfloor\)是\(\lfloor\frac{n}{i}\rfloor\)取值相同的一段区间的右端点

\(\lfloor\frac{n}{ab}\rfloor=\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{\lfloor\frac{n}{b}\rfloor}{a}\rfloor\)(对\(\lceil\rceil\)也成立)​

n的所有因子之和(对n的质因数分解式)\(\sigma(n) = \frac{p^{e_1 + 1}_1 - 1}{p_1-1} *\frac{p^{e_2 + 1}_2 - 1}{p_2-1}*...*\frac{p^{e_k + 1}_k - 1}{p_k-1}\)

n个数的全排列的逆序对数和:\(n!\frac{n(n-1)}{4}\)

\(Cayley\)公式：\(n\)个带编号的结点组成\(s\)棵树的森林，其中指定\(s\)个不同的结点属于不同的树，这样的组成方案一共有\(F(n,s)=sn^{n-s-1}\)种。

特别地，\(n\)个结点带编号的完全图一共有\(n^{n-2}\)个生成树，即\(n\)个结点带编号的无根树有\(n^{n-2}\)个。

威尔逊定理：\(p\)为质数\(\ \ iff \ \ p|(p-1)!+1\)

\(gcd(a^n−1,a^m−1)=a^{gcd(n,m)}−1\)

\(a>b,gcd(a,b)=1\)时有\(gcd(a^m-b^m,a^n-b^n)=a^{gcd(m,n)}-b^{gcd(m,n)}\)

设\(G(n)=gcd(C_n^1,C_n^2,...,C_n^{n-1}),\)则对任意质数\(p\)和自然数\(i\),有\(G(p^i)=p\),否则\(G(n)=1\).

\((n+1)lcm(C_n^0,C_n^1,...,C_n^n)=lcm(1,2,...,n+1)\)

对质数\(p,(x_1+x_2+...+x_n)^p\equiv(x_1^p+x_2^p+...+x_n^p)\ (mod\ p)\).

斐波那契数列

\(F_i=\frac{1}{\sqrt5}[(\frac{1+\sqrt5}{2})^i-(\frac{1-\sqrt5}{2})^i]\)

\(gcd(F_i,F_j) = F_{gcd(i,j)}\)

\(F_n^2-F_{n-1}F_{n+1}=(-1)^{n-1}\)

\(F_1+F_3+...+F_{2n-1}=F_{2n}-F_2+F_1\)

\(F_2+F_4+...+F_{2n}=F_{2n+1}-F_1\)

\(F_1^2+F_2^2+...+F_n^2=F_nF_{n+1}\)

\(F_{2n-2m-2}(F_{2n}+F_{2n+2})=F_{2m+2}+F_{4n-2m}(n > m \geq-1,n\geq1)\)

\(\frac{F_{2n}}{F_n}=F_{n-1}+F_{n+1}\)

中国剩余定理

\[x\equiv{a_i}\pmod{m_i}(1\leq i \leq n)(方程组S) \]

且\(m_1, m_2, ...,m_n\)两两互质,则对任意的\(a_1, a_2, ...,a_n\)同余方程组有解，

设\(M = m_1*m_2*...*m_n\) 是整数\(m_1, m_2, ...,m_n\)的乘积，并设 \(M_i\)是除了\(m_i\)以外的n-1个整数的乘积。

设 \(t_i\)为\(M_i\) 模 \(m_i\) 意义下的逆元)
​ \(M_i * t_i\equiv{1}\pmod{m_i}(1\leq i \leq n)\)
​ 方程组 (S) 的通解形式为
​ $x = kM + ∑(a_it_i*M_i),i∈{1,2,3,…,n},k∈Z $
​ 在模 M 的意义下，方程组(S)只有一个解如下图所示

\[x = (\sum_{i = 1}^{n}a_i * t_i * M_i) \% M \]

//求M%A=a,M%B=b,...中的M，其中A,B,C...互质
int CRT(int a[],int m[],int n){  
    int M = 1;  
    int ans = 0;  
    for(int i=1; i<=n; i++)  
        M *= m[i];  
    for(int i=1; i<=n; i++){  
        int x, y;  
        int Mi = M / m[i];  
        ex_gcd(Mi, m[i], x, y);  
        ans = (ans + Mi * x * a[i]) % M;  
    }  
    if(ans < 0) ans += M;  
    return ans;  
}  

#include<bits/stdc++.h>/*模不互质*/
using namespace std;
int a[10],m[10],tong[103];
int exgcd(int a,int &x,int b,int &y)
{
    if(b==0){x=1;y=0;return a;}
    int x2,y2;
    int gcd=exgcd(b,x2,a%b,y2);
    x=y2;
    y=x2-a/b*y2;
    return gcd;
}
int main()
{
    freopen("crt.in","r",stdin);
    freopen("crt.out","w",stdout);
    int T;
    scanf("%d",&T);
    while(T--)
    {
        int n;
        scanf("%d",&n);
        scanf("%d%d",&a[1],&m[1]);
        //用lcm求当前到达过的方程中，所有模数的lcm,sum是最后的结果,但是每一次都会用
        //解出的x去更新它,fail判断是否有解 
        int lcm=m[1],sum=a[1],fail=0;
        for(int i=2;i<=n;i++)
        {
            int x,y;
            scanf("%d%d",&a[i],&m[i]);
            a[i]=((a[i]-sum)%m[i]+m[i])%m[i];//sum就是要求的x,不断更新 
            int d=exgcd(lcm,x,m[i],y);//lcm*x+m[i]*y+sum=a[i]才能保证x%m[i]=a[i] 
            if(a[i]%d==0) x=x*(a[i]/d)%m[i];
            else fail=1;
            sum+=x*lcm;//注意是乘lcm哦！不然可能会出现前面的方程冲突 
            lcm=lcm/d*m[i];//到现在这一个方程了,所有模数的最小公倍数 
            sum=(sum%lcm+lcm)%lcm;
        }
        if(fail) printf("No\n");
        else printf("%d\n",sum); 
    }
}

lucas定理 求 \(C_n^k\% p\)(p为素数，p规模小于1e5)

表达式：C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p。（可以递归）

递归方程：(C(n%p, m%p)*Lucas(n/p, m/p))%p。（递归出口为m==0，return 1）

扩展lucas

输入三个整数\(n,k,p\),输出\(C_n^k\%p\).\(1\leq k \leq n \leq 1e18, 2 \leq p \leq 1e6\),\(p\)不一定是质数。

#include<bits/stdc++.h>
using namespace std;
const int maxn = 100000 + 10;
typedef long long LL;

LL Pow(LL n, LL m, LL mod) {
    LL ans = 1;
    while(m > 0) {
        if(m & 1) ans = (LL)ans * n % mod;
        n = (LL)n * n % mod; m >>= 1;
    }
    return ans;
}
LL Pow(LL n,LL m) {
    LL ans = 1;
    while(m > 0) {
        if(m & 1) ans = ans * n;
        n = n * n; m >>= 1;
    }
    return ans;
}
LL x, y;
LL exgcd(LL a, LL b) {
    if(a == 0) {
        x = 0, y = 1;
        return b;
    }LL r = exgcd(b%a, a);
    LL t = x; x = y - (b/a)*x; y = t;
    return r;
}
LL rev(LL a, LL b) { exgcd(a, b); return ((x % b) + b) % b; }
LL Calc(LL n, LL p, LL t) {
    if(n == 0) return 1;

    LL s = Pow(p, t), k = n / s, tmp = 1;
    for(LL i=1; i<=s; i++) if(i % p) tmp = (LL)tmp * i % s;

    LL ans = Pow(tmp, k, s);
    for(LL i=s*k + 1; i<=n; i++) if(i % p) ans = (LL)ans * i % s;

    return (LL)ans * Calc(n / p, p, t) % s;
}
LL C(LL n, LL m, LL p, LL t) {
    LL s = Pow(p, t), q = 0;
    for(LL i=n; i; i/=p) q += i / p;
    for(LL i=m; i; i/=p) q -= i / p;
    for(LL i=n-m; i; i/=p) q -= i / p;

    LL ans = Pow(p, q);
    LL a = Calc(n, p, t), b = Calc(m, p, t), c = Calc(n-m, p, t);
    return (LL)(ans * a * rev(b, s) * rev(c, s)) % s;
}
LL China(LL A[], LL M[], LL cnt) {
    LL ans = 0, m, n = 1;
    for(LL i=1; i<=cnt; i++) n *= M[i];
    for(LL i=1; i<=cnt; i++) {
        m = n / M[i];
        exgcd(M[i], m);
        ans = (ans + (LL)y * m * A[i]) % n;
    }
    return (ans + n) % n;
}
LL A[maxn], M[maxn], cnt;
LL Lucas(LL n, LL m, LL mod) {
    for(LL i=2; i*i <= mod; i++) if(mod % i == 0) {
        LL t = 0;
        while(mod % i == 0) t++, mod /= i;
        M[++cnt] = Pow(i, t);
        A[cnt] = C(n, m, i, t);
    }if(mod > 1) {
        M[++cnt] = mod;
        A[cnt] = C(n, m, mod, 1);
    }
    return China(A, M, cnt);
}
LL n, k, p;
int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
    freopen("output.out", "w", stdout);
#endif
    cin >> n >> k >> p;
    cout << Lucas(n, k, p) << endl;
    return 0;
}

一些常见的组合数学数

卡特兰数\(C_n\)

1, 2, 5, 14, 42,132, 429, 1430, 4862, 16796,58786, 208012, 742900, 2674440, 9694845, 35357670, 129644790, 477638700, 1767263190, 6564120420, 24466267020, 91482563640, 343059613650, 1289904147324, 4861946401452

通项\(C_n = \frac{1}{n + 1} C^{n}_{2n}\)

性质\(C_n = C^{n}_{2n} - C^{n - 1}_{2n}\) \(C_0 = 1, C_{n+1} = \frac{2(2n+1)}{n+2}C_n\)

\(C_{n+1} = \sum_{i = 0}^nC_iC_{n - i}\) 若\(C_n\)是奇数，那么\(n=2^k-1\)

斯特林数

第一类Stirling数

\(s(n,m)\)表示将 n 个不同元素构成m个圆排列的数目。

无符号第一类斯特林数\(S_u(n,m)\) 带符号第一类斯特林数\(S_s(n, m)\)

\(x^{n\downarrow}=\prod_{i=0}^{n-1}(x-i)=\sum_{k=0}^ns_s(n,k)x^k\)

\(x^{n\uparrow}=\prod_{i=0}^{n-1}(x+i)=\sum_{k=0}^ns_u(n,k)x^k\)

有\(s_s(n,m)=(-1)^{n+m}s_u(n,m)\).

有\(s_u(n+1,m)=s_u(n,m-1)+ns_u(n,m)\)

\(s_s(n+1,m)=s_s(n,m-1)-ns_s(n,m)\)

\(s_u(0,0)=1,s_u(n,0)=0(n \geq 1), s_u(n,n)=1,s_u(n,1)=(n-1)!\)

\(s_u(n,n-1)=C_n^2,s_u(n,2)=(n-1)!*\sum_{i=1}^{n-1}\frac{1}{i},s_u(n,n-2)=2C_n^3+3C_n^4\)

\(\sum_{k=0}^ns_s(n,k)=0^n\),注意\(0^0=1\).

第二类Stirling数

\(n\)个不同元素拆分成\(m\)个无区别的集合的方案数

\(S(n,m)=\frac{1}{m!}\sum_{k=0}^m(-1)^kC_m^k(m-k)^n\)

\(S(n+1,m)=S(n,m-1)+mS(n,m)\)

两类数的转化\(\sum_{k=0}^nS(n,k)s(k,m)=\sum_{k=0}^ns(n,k)S(k,m)\)

贝尔数

集合划分的数目，如\(B_3=5,\)表示\(\{a,b,c\}\)的五个划分为\(\{\{a\},\{b\}, \{c\}\},\{\{a\}, \{b, c\}\},\{\{b\}, \{a, c\}\},\{\{c\}, \{a, b\}\},\{\{a,b,c\}\}\)

\(B_0=B_1=1\)

$1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 190899322, 1382958545, 10480142147, 82864869804, 682076806159, 5832742205057, $$51724158235372, 474869816156751, 4506715738447323, 44152005855084346, 445958869294805289, 4638590332229999353, 49631246523618756274......$

\(B_{n+1}=\sum\limits_{k=0}^nC_n^kB_k\)

\(B_n=\frac{1}{e}\sum\limits_{k=0}^{\infin}\frac{k^n}{k!}\)

\(Touchard\)同余：若\(p\)为质数，则\(B_{p+n}\equiv B_n + B_{n+1}(mod\ p)\).

\(B_n=\sum\limits_{k=1}^nS(n,k)\),\(S\)是第二类Stirling数.

常用数论函数

欧拉函数

\[\phi(n) = n\prod^{k}_{i = 1}(1 - \frac{1}{p_i})(p_i为n的质因子，相同的质因子只计算一次) \]

\(\phi(p) = p - 1, \phi(p^k) = (p - 1)*p^{k - 1}\)

\(时，gcd(m,n) = 1时， \phi(mn) = \phi(m)\phi(n)\)

\(\sum_{d|n}\phi(d) = n\)

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

欧拉定理(n为素数即为费马小定理)

\(则有gcd(a,n) = 1,则有a^{\phi(n)}\equiv{1}\mod{n}\)

推广:欧拉降幂

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

LL Mod(LL b, LL mod){return b < mod?b:b%mod+mod;}

莫比乌斯函数

n质因数分解后，\(如果有幂次大于的为质因子个数\mu(1) = 1;\mu(n) = 0如果有幂次大于1的;\mu(n) = (-1)^k,k为质因子个数\)

n>2时，n的所有约数对应的莫比乌斯函数值之和为0 即\(\sum_{d|n}\mu(d) = 0\)

\(\sum_{d|n}\frac{\mu(d)}{d} = \frac{\phi(n)}{n}\)

等价于\(\phi(n)=\sum_{d|n}\frac{n}{d}\mu(d)=\sum_{d|n}d\mu(\frac{n}{d})\)

\([gcd(a,b) == 1] = \sum_{d|gcd(a,b)}\mu(d)=\sum_{d|a,d|b}\mu(d)\) 容斥

例题：求1-n中有多少个数和x互质

\(\sum_{d|x}\mu(d)*\lfloor\frac{n}{d}\rfloor\)容斥系数即为莫比乌斯函数

int mobius(int x)
{
    int res=1;
    int i,j;
    for(i=2;i*i<=x;i++)
        if(x%i==0)
        {
            x/=i;
            if(x%i==0) mu[i]=0;
            else mu[i]*=-1;
            while(x%i==0) x/=i;
        }
    return res;
}

exgcd

#include <iostream>
#define LL long long
using namespace std;

LL exgcd(LL a, LL b, LL &p, LL &q)
{
    if (a == 0 && b == 0) return -1;
    if (b == 0)
    {
        p = 1; q = 0; return a;
    }
    LL d = exgcd(b, a % b, q, p);
    q -= a / b * p;
    return d;
}

逆元

1.费马小定理----mod为素数

qp(a, mod - 2);

2.exgcd

给定模数m，求a的逆相当于求解ax=1(mod m),这个方程可以转化为ax-my=1 , 然后套用求二元一次方程的方法，用扩展欧几里得算法求得一组x0,y0和gcd
检查gcd是否为1 ,gcd不为1则说明逆元不存在
若为1，则调整x0到0~m-1的范围中即可

#define LL long long
void extgcd(LL a,LL b,LL &d,LL &x,LL &y)
{
    if(!b){ d=a; x=1; y=0;}
    else{ extgcd(b,a%b,d,y,x); y-=x*(a/b); }
}
LL inv(LL a,LL n)
{
    LL d,x,y;
    extgcd(a,n,d,x,y);
    return d==1?(x+n)%n:-1;
}

3.b|a

\[\frac{a}{b}\%mod = \frac{a\%(mod*b)}{b} \]

4.欧拉定理(m可以不为质数，但要a和m互质)

\[inv(a) = a^{\phi(mod) - 1}\%mod \]

5.逆元打表

\[inv[i] = (M - \lfloor\frac{M}{i}\rfloor) * inv[M\%i]\%M \]

#define LL long long
const int N = 1e5 + 5;
int inv[N];
void getinv(int n, int p) 
{
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) 
    {
        inv[i] = (LL) (p - p / i) * inv[p%i] % p;
    }
}

线性筛(素数筛，欧拉函数，莫比乌斯函数，约数个数，最小因子)

#define MAXN 1000000
int prime[MAXN + 8];
bool notprime[MAXN + 8];
//int phi[MAXN + 8];
///int mu[MAXN + 8];
////int facnum[MAXN + 8], d[MAXN + 8];
int ind = 0;
void getprime()
{
    //phi[1] = 1;
    ///mu[1] = 1;
    ////facnum[1] = 1;
    for (int i = 2; i <= MAXN; ++i)
    {
        if(!notprime[i])
        {
            prime[++ind] = i;
            //phi[i] = i - 1;
            ///mu[i] = -1;
            ////facnum[i] = 2, d[i] = 1;
            /////minfactor[i] = i;
        }
        for (int j = 1; j <= ind && i * prime[j] <= MAXN; ++j)
        {
            notprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                //phi[i * prime[j]] = phi[i] * prime[j];
                ///mu[i * prime[j]] = 0;
                ////facnum[i * prime[j]] = facnum[i] / (d[i] + 1) * (d[i] + 2), d[i * prime[j]] = d[i] + 1;
                break;
            }
            //phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            ///mu[i * prime[j]] = -mu[i];
            ////facnum[i * prime[j]] = facnum[i] * 2, d[i * prime[j]] = 1;
            /////minfactor[i*prime[j]] = prime[j];
        }
    }
}

例子：前n个数的约数个数筛

facnum[i]表示i的约数个数

通过素数筛得到前n个数的约数个数非常巧妙，首先根据约数个数定理：

对于一个大于1正整数n可以分解质因数：\(n = p_1^{d_1} * p_2^{d_2} *... * p_k^{d_k}\)，其中pi为素数

则n的正约数的个数就是

：facnum[n] = (1 + d1) * (1 + d2) * ... * (1 + dk)
我们需要一个辅助数组d[i]，表示i的最小质因子的次幂，(最小的原因是素数筛里每次都是用最小的质因子来筛合数的)，还是三种情况：
1.当i为素数时，facnum[i] = 2；d[i] = 1，很好理解
2.当i % p[j] != 0时，gcd(i, p[j]) =1，由积性函数的性质可得facnum[i * p[j]] = facnum[i] * facnum[p[j]] = facnum[i] * 2
​ d[i * p[j]] = 1(无平方因子)
3.当i % p[j] == 0时，出现平方因子，最小质因子的次幂加1，因此有

facnum[i * p[j]] = facnum[i] / (d[i] + 1) * (d[i] + 2)
​ d[i * p[j]] = d[i] + 1

康托展开

康托展开

康托展开是一个全排列到一个自然数的双射，常用于构建hash表时的空间压缩。设有n个数（1，2，3，4,…,n），可以有组成不同(n!种)的排列组合，康托展开表示的就是是当前排列组合在n个不同元素的全排列中的名次。

\[X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! \]

在（1，2，3，4，5）5个数的排列组合中，计算 34152的康托展开值。

首位是3，则小于3的数有两个，为1和2，a[5]=2，则首位小于3的所有排列组合为 a[0]*(5-1)!

第二位是4，则小于4的数有两个，为1和2，注意这里3并不能算，因为3已经在第一位，所以其实计算的是在第二位之后小于4的个数。因此a[4]=2

第三位是1，则在其之后小于1的数有0个，所以a[3]=0

第四位是5，则在其之后小于5的数有1个，为2，所以a[2]=1

最后一位就不用计算啦，因为在它之后已经没有数了，所以a[1]固定为0

根据公式：
X = 2 * 4! + 2 * 3! + 0 * 2! + 1 * 1! + 0 * 0! = 2 * 24 + 2 * 6 + 1 = 61
所以比 34152 小的组合有61个，即34152是排第62。

static const int FAC[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};   // 阶乘
int cantor(int *a, int n)
{
    int x = 0;
    for (int i = 0; i < n; ++i) {
        int smaller = 0;  // 在当前位之后小于其的个数
        for (int j = i + 1; j < n; ++j) {
            if (a[j] < a[i])
                smaller++;
        }
        x += FAC[n - i - 1] * smaller; // 康托展开累加
    }
    return x;  // 康托展开值
} 

逆康托展开

一开始已经提过了，康托展开是一个全排列到一个自然数的双射，因此是可逆的。即对于上述例子，在（1，2，3，4，5）给出61可以算出起排列组合为 34152。由上述的计算过程可以容易的逆推回来，具体过程如下：

用 61 / 4! = 2余13，说明a[5]=2,说明比首位小的数有2个，所以首位为3。

用 13 / 3! = 2余1，说明a[4]=2，说明在第二位之后小于第二位的数有2个，所以第二位为4。

用 1 / 2! = 0余1，说明a[3]=0，说明在第三位之后没有小于第三位的数，所以第三位为1。

用 1 / 1! = 1余0，说明a[2]=1，说明在第二位之后小于第四位的数有1个，所以第四位为5。

最后一位自然就是剩下的数2啦。

通过以上分析，所求排列组合为 34152。

static const int FAC[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};   // 阶乘
void decantor(int x, int n)
{
    vector<int> v;  // 存放当前可选数
    vector<int> a;  // 所求排列组合
    for(int i=1;i<=n;i++)
        v.push_back(i);
    for(int i=m;i>=1;i--)
    {
        int r = x % FAC[i-1];
        int t = x / FAC[i-1];
        x = r;
        sort(v.begin(),v.end());// 从小到大排序 
        a.push_back(v[t]);      // 剩余数里第t+1个数为当前位
        v.erase(v.begin()+t);   // 移除选做当前位的数
    }
}

莫比乌斯反演+数论分块

如果存在函数F(x)和f(x),满足 \(F(n) = \sum_{d|n}f(d)\),那么

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

或如果\(F(n) = \sum_{n|d}f(d)\),那么

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

例题:求对于区间[a，b]内的整数x和[c, d]内的y，满足gcd(x,y)=k的数对的个数。

\(的对数F(t):gcd(x, y) \% t == 0 的(x,y)对数\) \(的对数f(t):gcd(x, y) == t 的(x, y)对数\)

考虑 $x \in [1,b] $ and $ y \in [1,d]$ : \(F(t) = \lfloor\frac{b}{t}\rfloor\lfloor\frac{d}{t}\rfloor\)

$Ans = f(k), gcd(x,y) = k \leftrightarrow gcd(x/k, y/k) = 1 \leftrightarrow $

\(gcd(x,y) = 1(x\in[1,b/k]\) and \(y\in[1,d/k])\)

因此直接把b，d都除上k,\(f(1) = \sum^{min(b/k,d/k)}_{i = 1}\mu(i)F(i)\)就是结果

当t接近最大值时，对一段连续的t它们的\(F(t)\)值都相同，因此可以打出莫比乌斯前缀和然后分块求解。

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

bool check[MAXN + 8];
int prime[MAXN + 8];
int mu[MAXN + 8];
int sum[MAXN + 8];
void Mobius()
{
    memset(check, false, sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for (int i = 2; i <= MAXN; ++i)
    {
        if (!check[i])
        {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot; ++j)
        {
            if (i * prime[j] > MAXN)
                break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= MAXN; ++i)
    {
        sum[i] = sum[i - 1] + mu[i];
    }
}

LL calc(int a, int b)
{
    LL ret = 0;
    int last = 0;
    for (int i = 1; i <= min(a, b); i = last + 1)
    {
        last = min(a / (a / i), b / (b / i));
        ret += 1LL * (sum[last] - sum[i - 1]) * (a / i) * (b / i);
    }
    return ret;
}
int main()
{
    Mobius();
    int T, a, b, c, d, k;
    cin >> T;
    for (int t = 1; t <= T; ++t)
    {
        cin >> a >> b >> c >> d >> k;
        a--, c--;
        cout << calc(b / k, d / k) - calc(a / k, d / k) - calc(b / k, c / k) + calc(a / k, c / k);//容斥
        if (t != T) cout << endl;
    }
    return 0;
}

卷积

\(id(n)=n,\epsilon(n)=[n=1],I(n)=1,d(n)=\sum_{d|n}1,\sigma(n)=\sum_{d|n}d,\lambda(n)=(-1)^k\)

\(\mu*I=\epsilon,\phi*I=id,\phi=id*\mu,d=I*I,1=\mu*d,\sigma*\mu=id,\sigma=id*I\)

\(\epsilon\)卷任何\(f\)都为\(f\)

常见求和公式及其推导

1.\(\sum_{a=1}^n\sum_{b=1}^mgcd(a,b)=\sum_{d=1}^{min(n,m)}\phi(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\)

不妨设\(n\leq m\),枚举\(d=gcd(a,b)\),反演得到左边=\(\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}d\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor\),变形，用\(t\)替换\(i*d\),得\(\sum_{t=1}^n\lfloor\frac{n}{t}\rfloor\lfloor\frac{m}{t}\rfloor\sum_{d|t}d\mu(\frac{t}{d})\),后面那个式子是典型卷积，所以原公式成立。

2.\(\sum_{a=1}^n\sum_{b=1}^mlcm(a,b)=\frac{1}{4}\sum_{d=1}^{min(n,m)}d\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)(\lfloor\frac{m}{d}\rfloor+1)\sum_{i|d}i\mu(i)\)

反演：设\(F(t)=\sum_{i=1}^n\sum_{j=1}^mij[t|gcd(i,j)],f(t)=\sum_{i=1}^n\sum_{j=1}^mij[t=gcd(i,j)]\),则有\(F(t)=\sum_{i=1}^n\sum_{j=1}^mij[t|i][t|j]=\frac{\lfloor\frac{n}{t}\rfloor\lfloor\frac{m}{t}\rfloor(\lfloor\frac{n}{t}\rfloor+1)(\lfloor\frac{m}{t}\rfloor+1)}{4}\).因此答案=\(\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)i^2\frac{\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor(\lfloor\frac{n}{id}\rfloor+1)(\lfloor\frac{m}{id}\rfloor+1)}{4}\).\(t=id\)来替换\(d\),套路求解。记\(G(d)=\sum_{i|d}i\mu(i)\),显然是积性函数，对质数\(p\)有\(G(p^k)=1-p\),随便推一推就能用线性筛预处理，前面再一分块，就能\(O(\sqrt n)\)求解。

3.\(\sum_{a=1}^n\sum_{b=1}^nlcm(a,b)=\sum_{i=1}^n(-i+2\sum_{j=1}^ilcm(i,j))\).

预处理,O(1)输出。

4.\(\sum_{i=1}^ngcd(i,n)=\sum_{d|n}\frac{n}{d}\phi(d)\).

枚举\(d\),老套路了。

5.\(\sum_{i=1}^n\sigma(i)=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\).

根据\(\sigma\)的定义轻松得出。左边=\(\sum_{i=1}^n\sum_{j|i}j=\sum_{j=1}^nj\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}1\).

一类数论函数的更快求和方法

杜教筛---\(O(n^{2/3})\)计算一类积性函数的前缀和

我们要求\(\sum_{i=1}^nf(i)\),其中\(f\)为积性函数。

构造两个积性函数\(g,h\)满足\(h=f*g\)

\(\sum_{i=1}^nh(i)=\sum_{i=1}^n\sum_{d|i}g(d)f(i/d)=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)=\sum_{d=1}^ng(d)S_f(\lfloor\frac{n}{d}\rfloor)\)

因此\(g(1)S_f(n)=\sum_{i=1}^{n}h(i)-\sum_{i=2}^ng(d)S_f(\lfloor\frac{n}{d}\rfloor)\)

包含\(\phi(i),\mu(i),i\phi(i),i\mu(i)\)的板子

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

const int mod = 1000000007;
const int inv2 = 500000004;
const int inv4 = 250000002;
const int inv6 = 166666668;
LL qp(LL a, LL b)
{
    LL ret = 1;
    while (b)
    {
        if (b & 1)
            ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}

LL inv(LL a)
{
    return qp(a, mod - 2);
}
int prime[MAXN + 8];
bool notprime[MAXN + 8];
LL phi[MAXN + 8], phisum[MAXN + 8];
LL mu[MAXN + 8], musum[MAXN + 8];
LL iphisum[MAXN + 8];
LL imusum[MAXN + 8];
int ind = 0;
void getprime()
{
    phi[1] = 1;
    mu[1] = 1;
    for (int i = 2; i <= MAXN; ++i)
    {
        if(!notprime[i])
        {
            prime[++ind] = i;
            phi[i] = (i - 1) % mod;
            mu[i] = -1;
        }
        for (int j = 1; j <= ind && i * prime[j] <= MAXN; ++j)
        {
            notprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j] % mod;
                mu[i * prime[j]] = 0;
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1) % mod;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= MAXN; ++i)
    {
        phisum[i] = (phisum[i - 1] + phi[i]) % mod;
        musum[i] = (musum[i - 1] + mu[i] + mod) % mod;
        iphisum[i] = (iphisum[i - 1] + i * phi[i] % mod) % mod;
        imusum[i] = (imusum[i - 1] + i * mu[i] % mod + mod) % mod;
    }
}
unordered_map<LL,LL> _phi, _mu, _iphi, _imu;
LL Calcphi(LL n)
{
    if (n <= MAXN)
        return phisum[n];
    auto it = _phi.find(n);
    if (it != _phi.end())
        return it->second;
    LL last;
    LL ans = (n % mod) * (((n % mod) + 1 + mod) % mod) % mod * inv2 % mod;
    for (LL i = 2; i <= n; i = last + 1)
    {
        last = (n / (n / i));
        ans = (ans - (((last % mod) - (i % mod) + 1 + mod) % mod * (Calcphi(n / i) % mod) % mod) + mod) % mod;
    }
    return _phi[n] = ans;
}

LL Calciphi(LL n)
{
    if (n <= MAXN)
        return iphisum[n];
    auto it = _iphi.find(n);
    if (it != _iphi.end())
        return it->second;
    LL ret = n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
    LL last;
    for (LL i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ret = (ret - ((last + i) % mod * (last - i + 1) % mod * inv2 % mod + mod) % mod * Calciphi(n / i) % mod + mod) % mod;
    }
    return _iphi[n] = ret;
}

LL Calcmu(LL n)
{
    if (n <= MAXN)
        return musum[n];
    auto it = _mu.find(n);
    if (it != _mu.end())
        return it->second;
    LL ans = 1;
    LL last;
    for (LL i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ans = (ans - (last - i + 1) * Calcmu(n / i) % mod + mod) % mod;
    }
    return _mu[n] = ans;
}

LL Calcimu(int n)
{
    if (n <= MAXN)
        return imusum[n];
    auto it = _imu.find(n);
    if (it != _imu.end())
        return it->second;
    int last;
    LL ret = 1;
    for (int i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ret = (ret - ((1LL * last + i) % mod * (last - i + 1) % mod * inv2 % mod * Calcimu(n / i) % mod) + mod) % mod;
    }
    return _imu[n] = ret;
}

洲阁筛（还没学）

简单的写了个求1~n质数个数 -O2 N=1e11要2.5s

#pragma GCC optimize(2)
#include <bits/stdc++.h>
using namespace std;
#define LL long long
int L0[1000005],L1[1000005],sn,m,P[1000005],x,ga[1000005],gb[1000005];
LL n,a[1000005],b[1000005];
int main(){
    n=100000000000;
//    scanf("%lld",&n);
    sn=sqrt(n);
    for (int i=2;i<=sn;++i){
        if (!P[i]) P[++m]=i;
        for (int j=1;j<=m;++j){
            x=P[j]*i;
            if (x>sn) break;
            P[x]=1;
            if (i%P[j]==0) break;
        }
    }
    for (int i=1;i<=sn;++i) a[i]=i,b[i]=n/i;
    for (int i=1,j=1;i<=sn;++i){
        while (j<=m&&P[j]*P[j]<=i) ++j;
        L0[i]=j;
    }
    for (int i=sn,j=L0[sn];i;--i){
        while (j<=m&&P[j]*P[j]<=b[i]) ++j;
        L1[i]=j;
    }
    for (int i=1;i<=m;++i){
        for (int j=1;j<=sn&&i<L1[j];++j){
            LL y=j*P[i]; gb[j]=i;
            if (y<=sn) b[j]-=b[y]+gb[y]+1-i;
            else b[j]-=a[n/y]+ga[n/y]+1-i;
        }
        for (int j=sn;j&&i<L0[j];--j){
            LL y=j/P[i]; ga[j]=i;
            a[j]-=a[y]+ga[y]+1-i;
        }
    }
    printf("%lld\n",b[1]+m-1);
    return 0;
}

min25筛（学了一半）

待咕

类欧几里得算法

可快速计算:

\(f(a,b,c,n)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor\)

\(g(a,b,c,n)=\sum\limits_{i=0}^ni\lfloor\frac{ai+b}{c}\rfloor\)

\(h(a,b,c,n)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2\)

常用推导公式

\(a\leq \lfloor\frac{b}{c}\rfloor⇔ac\leq b\)

\(a< \lceil\frac{b}{c}\rceil⇔ac< b\)

\(a\geq \lceil\frac{b}{c}\rceil⇔ac\geq b\)

\(a> \lfloor\frac{b}{c}\rfloor⇔ac> b\)

\(\lfloor\frac{b}{c}\rfloor=\lfloor\frac{b+c-1}{c}\rfloor=\lceil\frac{b-c+1}{c}\rceil\)

以f的推导过程为例

1. \(a\geq c\)或\(b \geq c\)

\(\lfloor\frac{ai+b}{c}\rfloor=\lfloor\frac{(a\%c)i+b\%c}{c}\rfloor+\lfloor\frac{a}{c}\rfloor i+\lfloor\frac{b}{c}\rfloor\)

因此\(f(a,b,c,n)=f(a\%c,b\%c,c,n)+\frac{n(n+1)}{2}\lfloor\frac{a}{c}\rfloor+(n+1)\lfloor\frac{b}{c}\rfloor\)

2.\(a<c\)且\(b<c\)

\(a=0:f(a,b,c,n)=0\)

否则：

\[f(a,b,c,n)=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor -1}1 \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor] \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}\sum_{i=0}^n[j<\lceil\frac{ai+b-c+1}{c}\rceil] \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}\sum_{i=0}^n[cj<ai+b-c+1] \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}\sum_{i=0}^n[ai>cj-b+c-1] \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}\sum_{i=0}^n[i>\lfloor\frac{cj-b+c-1}{c}\rfloor] \\=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor -1}n-\lfloor\frac{cj-b+c-1}{c}\rfloor \\=n\lfloor\frac{an+b}{c}\rfloor-f(c,-b+c-1,a,\lfloor\frac{an+b}{c}\rfloor-1) \]

\(或g(a,b,c,n)=g(a\%c,b\%c,c,n)+\frac{n(n+1)(2n+1)}{6}\lfloor\frac{a}{c}\rfloor+\frac{n(n+1)}{2}\lfloor\frac{b}{c}\rfloor\ (a\geq c或b\geq c)\)

\(否则g(a,b,c,n)=\frac{\lfloor\frac{an+b}{c}\rfloor n(n+1)-f(c,-b+c-1,a,\lfloor\frac{an+b}{c}\rfloor-1)-h(c,-b+c-1,a,\lfloor\frac{an+b}{c}\rfloor-1)}{2}(否则)\)

\(或h(a,b,c,n)=h(a\%c,b\%c,c,n)+\frac{n(n+1)(2n+1)}{6}\lfloor\frac{a}{c}\rfloor^2+(n+1)\lfloor\frac{b}{c}\rfloor^2+2\lfloor\frac{b}{c}\rfloor f(a\%c,b\%c,c,n)\\+2\lfloor\frac{a}{c}\rfloor g(a\%c,b\%c,c,n)+\lfloor\frac{a}{c}\rfloor\lfloor\frac{b}{c}\rfloor n(n+1)(a\geq c或b\geq c)\)

\(否则h(a,b,c,n)=\lfloor\frac{an+b}{c}\rfloor(\lfloor\frac{an+b}{c}\rfloor+1)n-2g(c,-b+c-1,a,\lfloor\frac{an+b}{c}\rfloor-1)\\-2f(c,-b+c-1,a,\lfloor\frac{an+b}{c}\rfloor-1)-f(a,b,c,n)(否则)\)

/*洛谷P5170，输入T组，一开始写了个递归暴力，给我T飞掉了，我人都傻了*/
#include <bits/stdc++.h>
#define LL long long
#define MAXN 2000
#define x first
#define y second
#define pii pair<long long, long long>
using namespace std;
const LL mod = 998244353;
const LL inv2 = 499122177;
const LL inv6 = 166374059;
const LL inv4 = 748683265;
struct Node
{
    LL f, g, h;
    Node(){}
    Node(LL _f, LL _g, LL _h)
    {
        f = _f; g = _g; h = _h;
    }
};

Node solve(LL a, LL b, LL c, LL n)
{
    if (a == 0)
    {
        n %= mod;
        LL bc = (b/c) % mod;
        return Node((n+1)*bc%mod, n*(n+1)%mod*inv2%mod*bc%mod, (n+1)*bc%mod*bc%mod);
    }
    LL t = (a*n+b)/c;
    LL ac=(a/c)%mod, bc=(b/c)%mod;
    if (a >= c || b >= c)
    {
        Node n1 = solve(a%c,b%c,c,n);
        return Node(( n1.f + n*(n+1)%mod*inv2%mod*ac%mod + (n+1)*bc%mod) % mod,
                    ( n1.g + n*(n+1)%mod*(2*n+1)%mod*inv6%mod*ac%mod % mod + n*(n+1)%mod*inv2%mod*bc%mod) % mod,
                    ( n1.h + n*(n+1)%mod*(2*n+1)%mod*inv6%mod*ac%mod*ac%mod + (n+1)*bc%mod*bc%mod + 2*bc%mod*n1.f%mod + 2*ac%mod*n1.g%mod + ac*bc%mod*n%mod*(n+1)%mod) % mod
                    );
    }
    else
    {
        Node n2 = solve(c,c-b-1,a,t-1);
        t %= mod;
        return Node((n * t % mod - n2.f + mod) % mod,
                    inv2 * ((t*n%mod*(n+1)%mod - n2.f - n2.h + 2*mod) % mod) % mod,
                    (t*(t+1)%mod*n%mod - 2*n2.g%mod - 2*n2.f%mod - (n * t % mod - n2.f + mod) % mod + 4*mod) % mod
                    );
    }
}

int main()
{
    LL a, b, c, n;
    int T;
    scanf("%d", &T);
    while (T--)
    {
        scanf("%lld %lld %lld %lld", &n, &a, &b, &c);
        Node ans = solve(a, b, c, n);
        printf("%lld %lld %lld\n", ans.f, ans.h, ans.g);
    }
	return 0;
}

拉格朗日插值法

朴素算法

给定\(k+1\)个点\((x_0,y_0),...,(x_k,y_k)\)表示的\(k\)次多项式\(L(x)\)，给出\(a\)求\(L(a)\).

记\(\ell_j(x)=\prod\limits_{i=0,i\neq j}^k\frac{x-x_i}{x_j-x_i}\),则有\(L(x)=\sum\limits_{j=0}^ky_j\ell_j(x)\). \(O(n^2)\).

/*洛谷P4781 给出n,X,接下来给出n对(x,y),求L(X).*/
#include <bits/stdc++.h>
#define LL long long
#define MAXN 2000
#define x first
#define y second
#define pii pair<long long, long long>
using namespace std;
const int mod = 998244353;
int n;
LL X;
pii arr[MAXN + 8];

LL qp(LL a, LL b)
{
    LL ret = 1;
    while (b)
    {
        if (b & 1)
            ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}

LL inv(LL a)
{
    return qp(a, mod - 2);
}

LL calcell(int j, LL X)
{
    LL ret = 1;
    for (int i = 1; i <= n; ++i)
    {
        if (i == j || arr[i].x == arr[j].x)
            continue;
        ret = ret * (X - arr[i].x + mod) % mod * inv((arr[j].x - arr[i].x + mod) % mod) % mod;
    }
    return ret;
}

LL calcF(LL X)
{
    LL ret = 0;
    for (int i = 1; i <= n; ++i)
    {
        ret = (ret + arr[i].y * calcell(i, X) % mod) % mod;
    }
    return ret;
}

int main()
{
    cin >> n >> X;
    for (int i = 1; i <= n; ++i)
    {
        cin >> arr[i].x >> arr[i].y;
    }
    cout << calcF(X) << endl;
	return 0;
}

重心拉格朗日插值法

插入新点时单次插入优化至\(O(n)\):

记\(\ell(x)=\prod\limits_{i=0}^k(x-x_i)\),\(w_j=\frac{1}{\prod_{i=0,i\neq j}^k(x_j-x_i)}\),则\(\ell_j(x)=\ell(x)\frac{w_j}{x-x_j}\).

则\(L(x)=\ell(x)\sum\limits_{j=0}^k\frac{w_j}{x-x_j}y_j\).

对多项式\(g(x)=1\)插值，有\(\forall x,g(x)=\ell(x)\sum\limits_{j=0}^k\frac{w_j}{x-x_j}\),因此\(L(x)\)两边同除以\(g(x)\),有

\(L(x)=\frac{\sum_{j=0}^k\frac{w_j}{x-x_j}y_j}{\sum_{j=0}^k\frac{w_j}{x-x_j}}\).插入新值\(x_{k+1}\)时，只需对\(w_j\)乘上\(x_j-x_{k+1}\)即可更新\(w\),更新所有\(w\)总共需要\(O(n)\).

/*依然是上面那个模板题*/
#include <bits/stdc++.h>
#define LL long long
#define MAXN 2000
#define x first
#define y second
#define pii pair<long long, long long>
using namespace std;
const int mod = 998244353;
int n;
LL X;
pii arr[MAXN + 8];
LL w[MAXN + 8];
LL qp(LL a, LL b)
{
    LL ret = 1;
    while (b)
    {
        if (b & 1)
            ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
}

LL inv(LL a)
{
    return qp(a, mod - 2);
}

void insert(int i)
{
    w[i] = 1;
    for (int j = 1; j < i; ++j)
    {
        if (arr[i].x == arr[j].x)
            continue;
        w[i] = w[i] * (arr[i].x - arr[j].x + mod) % mod;
        w[j] = w[j] * inv(arr[j].x - arr[i].x + mod) % mod;
    }
    w[i] = inv(w[i]);
    return;
}

LL query(int n)
{
    LL fz = 0, fm = 0;
    for (int i = 1; i <= n; ++i)
    {
        fz = (fz + w[i] * inv((X - arr[i].x + mod) % mod) % mod * arr[i].y % mod) % mod;
        fm = (fm + w[i] * inv((X - arr[i].x + mod) % mod) % mod) % mod;
    }
    return fz * inv(fm) % mod;
}

int main()
{
    cin >> n >> X;
    for (int i = 1; i <= n; ++i)
    {
        cin >> arr[i].x >> arr[i].y;
        insert(i);
    }
    cout << query(n) << endl;
	return 0;
}

佩尔方程

\(x^2-Dy^2=1,D\)不是平方数。

设\((x_1,y_1)\)是使\(x_1\)最小的解，则通解为

\(x_k+\sqrt{D}y_k=(x_1+\sqrt{D}y_1)^k\).

\(x^2-Dy^2=1,D\)不是平方数。

设\((x_0,y_0)\)是使\(x_0\)最小的解，则通解为

\(x_k+\sqrt{D}y_k=(x_0+\sqrt{D}y_0)^{2k+1}\).

佩尔方程(\(x^2-Dy^2=1\))的最小正整数解

bool PQA(LL D, LL &p, LL &q) {//来自于PQA算法的一个特例
    const double eps = 1e-6;
    LL d = sqrt(D);
    if (abs(d * d - D) <= eps)
        return false;//这里是判断佩尔方程有没有解
    LL u = 0, v = 1, a = int(sqrt(D)), a0 = a, lastp = 1, lastq = 0;
    p = a, q = 1;
    do {
        u = a * v - u;
        v = (D - u * u) / v;
        a = (a0 + u) / v;
        LL thisp = p, thisq = q;
        p = a * p + lastp;
        q = a * q + lastq;
        lastp = thisp;
        lastq = thisq;
    } while ((v != 1 && a <= a0));//这里一定要用do~while循环
    p = lastp;
    q = lastq;
    //这样求出后的(p,q)是p2 – D * q2 = (-1)k的解，也就是说p2 – D * q2可能等于1也可能等于-1,如果等于1,(p,q)就是解,如果等于-1还要通过(p2 + D * q2,2 * p * q)来求解，如下
    if (p * p - D * q * q == -1) {
        p = lastp * lastp + D * lastq * lastq;
        q = 2 * lastp * lastq;
    }
    return true;
}

正整数的拆分方案数\(O(n\sqrt{n})\)

void init() {
    for(int i = 0; i < 1000; ++i) {
        g[i << 1] = 1LL * i * (i * 3 - 1) / 2;
        g[i << 1 | 1] = 1LL * i * (i * 3 + 1) / 2;
    }
    f[0] = 1, f[1] = 1, f[2] = 2;
    for(int i = 3; i < N; ++i) {
        f[i] = 0;
        int k = 0;
        for(int j = 2; g[j] <= i; ++j) {
            if(k & 2) 
                f[i] = ((f[i] - f[i - g[j]]) % mod + mod) % mod;
            else
                f[i] = ((f[i] + f[i - g[j]]) % mod + mod) % mod;
            k++; k %= 8;
        }
    }
}
/*g是五边形数，f是拆分方案数*/
/*若限制每个数大于等于3，则 ans=f[i]+f[i-3]-f[i-1]-f[i-2]*/

博弈论

常用简单博弈模型

巴什博弈

一堆n个石子，两人轮流取，每次最多取m个，取完最后一颗的人胜。

解:\(n\%(m+1) = 0\)时，先手必败，否则必胜。

斐波那契博弈

一堆n个石子，先手取任意多个，但不能取完，之后每次取的石子数只能在1到对手刚取的石子数的2倍之间（闭区间），取完最后一颗的人胜。

解:\(n\)是斐波那契数时，先手必败，否则必胜。

威佐夫博弈

两堆各若干个石子，每次可以一堆取任意多个，或两堆取同样多个，取完最后一颗的胜。

解:\((0,0),(1,2),(3,5),(4,7),(6,10)....\)先手必败。一般地，对奇异局势(先手必败局势)\((a,b)(a\le b),\)有$a=\lfloor(b-a)\frac{\sqrt 5 + 1}{2}\rfloor $

或者说，第\(k\)组奇异局势\((a_k,b_k)\)有\(a_k=\lfloor\frac{\sqrt 5 + 1}{2}k\rfloor, b_k=a_k +k\).

尼姆博弈

若干堆石子，每堆若干个，每次选一堆取至少一个，取完最后一颗的胜。

解:\(\bigoplus_{i=1}^n a[i] = 0\)时，先手必败，否则必胜。

阶梯博弈

若干堆石子，堆标号依次\(1,2,3,...,n\)每次选一堆(标号为\(i\))的若干个放到标号\(i-1\)的堆(可以放到标号为0的堆，称作地面)，把最后一颗石子放到\(0\)号堆的胜。

解:等价于保留奇数号，忽视偶数号的Nim博弈。