bqlp的数学模板笔记 / 开普饭

数学

数学
- 数论
  - 素数
- 素性探测
  - Fermat 素性测试
  - 卡迈克尔数
  - 二次探测定理
  - 实现
- miller-rabin Pollard_rho(kuangbin模板)
  - 短小互推的
    - 欧拉定理
    - 扩展欧拉定理
    - 费马小定理
    - 裴蜀定理
    - 类欧几里得算法
    - 线性同余方程
    - 中国剩余定理
    - 二次剩余
    - N次剩余
    - 卢卡斯定理
    - exLucas定理
    - BSGS
    - 扩展BSGS
    - 原根
    - 杜教筛
  - 乘法逆元
    - 扩展欧几里得：n
    - 线性递推：1到n
    - 线性递推：1!到n!
    - 线性递推：给定随机n个数
  - 常见数论函数
    - 1）u(n), e(n), I(n)
    - 2）d(n) 约数个数
    - 3）Omega(n) 可重复素因子个数
    - 4）omega(n) 不可重复素因子个数
    - 5）sigma_lambda(n) 约数幂和
    - 6）varphi(n) 互质个数
    - 7）mu(n)
    - 8）Lambda(n)
    - 9）lambda(n)
    - 可乘函数（积性函数）
    - 完全（绝对）可乘函数
  - Dirichlet 乘积
    - 数论分块
    - 常见卷积
    - 莫比乌斯反演
- FT
  - 快速傅里叶变换
  - 快速沃尔什变换
- 快速幂
- 快速乘
- 常见数列
  - 斐波那契数
  - 卡特兰数
  - Harmonic(调和级数)
- 常见常数
  - e
  - pi
- 泰勒公式
- 组合数学
  - 排列组合
- 其它数学
  - Cantor
  - 数位dp
  - 计算基础
    - 求区间$[i, j]$的两两组合乘积的和：
    - 求$n^k$的最高$m$位：
    - 求$lcm(i,j)=n,\ i\leq j$的$(i,j)$对数：
    - 求$n!$中$p$的个数：
    - 求$ax bz=c$且$x\geq0,b\geq0$的$(x,y)$对数，$\gcd(a,b)=1$：
    - 求$ax (a 1)y (a 2)z==b$，$x,y,z$为非负整数的$(x,y,z)$对数：
    - 求有限$\pm2^i$表示数：
  - 四边形不等式
- 多项式
  - 拉格朗日插值
- 高精度

数论

素数

素数计数函数：小于或等于$x$的素数的个数，用$\pi(x)$表示。随着$x$的增大，有这样的近似结果：$\pi(x)\sim\frac{x}{ln(x)}$

素数判定：

bool isPrime(a) {
if (a<2) return 0;
for (int i=2; i*i<=a; i  )
if (a%i==0) return 0;
return 1;
}

素数筛：

$O(n\log\log n)$

int m=sqrt(n 0.5);
memset(vis, 0, sizeof(vis));
for(int i=2; i<=m; i  ) if(!vis[i])
for(int j=i*i; j<=n; j =i) vis[j]=1;
for(int i=2; i<=n; i  ){
if(!vis[i]) pr[tot  ]=i;
}

欧拉筛法，优化：取消合数的重复标记。$O(n)$

//phi, MAXN, vis, pri, cnt
void init() {
phi[1]=1;
for (int i=2; i<MAXN; i  ) {
if(!vis[i]) {
phi[i]=i-1;
pri[cnt  ]=i;
}
for(int j=0; j<cnt; j  ) {
if (1ll*i*pri[j]>=MAXN) break;
vis[i*pri[j]]=1;
if (i%pri[j]) {
phi[i*pri[j]]=phi[i]*(pri[j]-1);
} else {
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
}
}
}

求大数小区间$(a,b)$间素数：

int cal(ll a, ll b){
memset(v, 0, sizeof(v));
for(int i=0; i<tot; i  ){
if(pr[i]>b) break;
ll x=a/pr[i] (a%pr[i]!=0);
if(x<=1) x=2;
for(ll j=x*pr[i]; j<=b; j =pr[i])
if(j>=a) v[j-a]=1;
}
int ans=0;
for(ll i=0; i<=b-a; i  ){
if(!v[i]){
ans  ;
}
}
return ans;
}

素性探测

检测$n$是否为素数。

Fermat 素性测试

使用费马小定理，在$[2,n-1]$中不断取随机数$a$，检验$a^{n-1}\equiv1\pmod n$不成立时，一定不为素数。

对所有$a$成立时，可能为素数。

卡迈克尔数

对所有$[2,n-1]$中的$a$都满足$a^{n-1}\equiv1\pmod n$，但$n$为合数。

eg：$561=3\times11\times17$

性质：若$n$为卡迈克尔数，则$m=2^n-1$也是卡迈克尔数，所以有无穷多个。

卡迈克尔数的形式一定不为$p^e$。

二次探测定理

原理：若$p$为奇素数，则$x^2\equiv1\pmod p$的解为$x=1$和$x=p-1$。

实现

分解$n-1=2^t\times u$，不断对$u$进行平方$\dots$

再极端的情况也只会把合数判成素数，不会把素数判成合数。

常用大素数：

18位素数：154590409516822759
19位素数：2305843009213693951 (梅森素数)
19位素数：4384957924686954497

一些伪素数（合数）：

开头中间结束各选了3个，

561
1105
1729
426821473
429553345
434330401
2851896013395343201
3589102252820237401
4954039956700380001

抄完模板后，以上用于验证模板抄错没。

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

ll mul(ll a,ll b,ll mod){
    a%=mod;
    b%=mod;
    ll c=(long double)a*b/mod;
    ll ans=a*b-c*mod;
    return (ans%mod mod)%mod;
}
ll pow_mod(ll x,ll n,ll mod){
    ll res=1;
    while(n){
        if(n&1)
        res=mul(res,x,mod);
        x=mul(x,x,mod);
        n>>=1;
    }
    return (res mod)%mod;
}

bool Miller_Rabbin(ll a,ll n){

    ll s=n-1,r=0;
    while((s&1)==0){
        s>>=1;r  ;
    }

    ll k=pow_mod(a,s,n);

    if(k==1)return true;
    for(int i=0;i<r;i  ,k=k*k%n){
        if(k==n-1)return true;
    }
    return false;
}
bool isprime(ll n){
    if(n==1) return false;
    ll times=7;
    ll prime[100]={2,3,5,7,11,233,331};
    for(int i=0;i<times;i  ){
        if(n==prime[i])return true;
        if(Miller_Rabbin(prime[i],n)==false)return false;
    }
    return true;
}
ll n;
int main(){
int T;
cin>>T;
while(T--){
scanf("%lld", &n);
if(isprime(n)) printf("yes\n");
else printf("no\n");
}
}

miller-rabin Pollard_rho(kuangbin模板)

输入是$<2^{63}$

#include<bits/stdc  .h>
using namespace std;
#define ll long long
#define S 20//随即判定的次数，S越大出错率

ll mult_mod(ll a,ll b,ll c)//a*b%c
{
    a%=c;
    b%=c;
    ll ret=0;
    while(b)
    {
        if(b&1){ret =a;ret%=c;}
        a<<=1;
        if(a>=c)a%=c;
        b>>=1;
    }
    return ret;
}

ll pow_mod(ll x,ll n,ll mod)//x^n%mod
{
    if(n==1)return x%mod;
    x%=mod;
    ll tmp=x;
    ll ret=1;
    while(n)
    {
        if(n&1) ret=mult_mod(ret,tmp,mod);
        tmp=mult_mod(tmp,tmp,mod);
        n>>=1;
    }
    return ret;
}
//一定是合数返回1，不一定返回0
bool check(ll a,ll n,ll x,ll t)
{
    ll ret=pow_mod(a,x,n);
    ll last=ret;
    for(int i=1;i<=t;i  )
    {
        ret=mult_mod(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1) return true;
        last=ret;
    }
    if(ret!=1) return true;
    return false;
}

//是素数或者伪素数为1， 非伪素数的合数为0
bool Miller_Rabin(ll n)
{
    if(n<2)return false;
    if(n==2)return true;
    if((n&1)==0) return false;
    ll x=n-1;
    ll t=0;
    while((x&1)==0){x>>=1;t  ;}
    for(int i=0;i<S;i  )
    {
        ll a=rand()%(n-1) 1;
        if(check(a,n,x,t))
            return false;
    }
    return true;
}

ll factor[100];//存因子，就算是2也不超过63个
//得到的factor是不排序的
int tol;

ll gcd(ll a,ll b)
{
    if(a==0)return 1;
    if(a<0) return gcd(-a,b);
    while(b)
    {
        ll t=a%b;
        a=b;
        b=t;
    }
    return a;
}
//找到一个因子
ll Pollard_rho(ll x,ll c)
{
    ll i=1,k=2;
    ll x0=rand()%x;
    ll y=x0;
    while(1)
    {
        i  ;
        x0=(mult_mod(x0,x0,x) c)%x;
        ll d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k){y=x0;k =k;}
    }
}
//因式分解n
void findfac(ll n)
{
    if(Miller_Rabin(n))
    {
        factor[tol  ]=n;
        return;
    }
    ll p=n;
    while(p>=n)p=Pollard_rho(p,rand()%(n-1) 1);
    //可以分为p,n/p，其中p可以不为素数
    findfac(p);
    findfac(n/p);
}

int main()
{
    //srand(time(NULL));
    ll n;
    while(scanf("%lld",&n)!=EOF)
    {
        tol=0;
        findfac(n);
        for(int i=0;i<tol;i  )printf("%lld ",factor[i]);
        printf("\n");
        if(Miller_Rabin(n))printf("Yes\n");
        else printf("No\n");
    }
    return 0;
}

卡迈克尔数

反素数

短小互推的

欧拉定理

$gcd(a,m)=1$，则$a^{\varphi(m)}\equiv1(\mod m)$

扩展欧拉定理

$a^b\equiv\begin{cases}a^{b{\mod\varphi(p)}},\quad &gcd(a,p)=1 \\ a^b,\quad &gcd(a,p)\neq1,b<\varphi(p) \\ a^{b\mod\varphi(p) \varphi(p)},\quad &gcd(a.p)\neq1,b\geq\varphi(p)\end{cases}\pmod p$

费马小定理

$p$为素数, 且$gcd(a,p)=1$, 则$a^{p-1}\equiv1(modp)$

裴蜀定理

$a,\ b$是不全为零的整数，则存在整数$x,\ y$, 使得$ax by=gcd(a,b)$

证明及求解见扩展欧几里得。

类欧几里得算法

$a,b,c$均为整数，求$f(n,a,b,c)=\sum_{i=0}^n \lfloor\frac{ai b}{c} \rfloor\\ g(n,a,b,c)=\sum_{i=0}^n i\lfloor\frac{ai b}{c} \rfloor\\ h(n,a,b,c)=\sum_{i=0}^n \lfloor\frac{ai b}{c} \rfloor^2$

类似辗转相除法，$O(\log n)$

可单独算$f$

#define M 998244353//由题目给定
#define inv2 499122177
#define inv6 166374059

ll t,n,a,b,c;
struct node {
ll f, g, h;
};
node solve(ll a, ll b, ll c, ll n) {
node ans, tmp;
if(a==0) {
ans.f=(b/c)*(n 1)%M;
ans.g=(b/c)*n%M*(n 1)%M*inv2%M;
ans.h=(b/c)*(b/c)%M*(n 1)%M;
} else if(a>=c || b>=c) {
tmp=solve(a%c, b%c, c, n);
ans.f=(tmp.f n*(n 1)%M*inv2%M*(a/c)%M (n 1)*(b/c)%M)%M;
ans.g=(tmp.g (a/c)*n%M*(n 1)%M*(2*n 1)%M*inv6%M (b/c)*n%M*(n 1)%M*inv2%M)%M;
ans.h=(tmp.h (a/c)*(a/c)%M*n%M*(n 1)%M*(2*n 1)%M*inv6%M
       (n 1)*(b/c)%M*(b/c)%M 2*(a/c)%M*tmp.g%M 2*(b/c)%M*tmp.f%M
       2*(a/c)%M*(b/c)%M*n%M*(n 1)%M*inv2%M)%M;
} else {
long long m=(a*n b)/c;
tmp=solve(c, c-b-1, a, m-1);
ans.f=(n*(m%M)%M-tmp.f)%M;
ans.g=(n*(n 1)%M*(m%M)%M-tmp.f-tmp.h)%M*inv2%M;
ans.h=(n*(m%M)%M*((m 1)%M)%M-2*tmp.g-2*tmp.f-ans.f)%M;
}
return ans;
}
int main() {
    int T;
cin>>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 M)%M,(ans.g M)%M,(ans.h M)%M);
}
return 0;
}

线性同余方程

形如$ax\equiv c\pmod b$的方程。

原方程等价于$ax by=c$，有整数解的充要条件为$\gcd(a,b)\mid c$。

特殊：$\gcd(a,b)=1,\ x_0,y_0$为方程一组解，通解为$x=x_0 bt,\ y=y_0-at$，$t$为任意整数。

int Exgcd(int a, int b, int &x, int &y) {
if (!b){
x=1, y=0;
return a;
}
int d=Exgcd(b, a%b, x, y);
int t=x;
x=y;
y=t-(a/b)*y;
return d;
}
bool liEu(int a, int b, int c, int& x, int& y){
  int d=ex_gcd(a, b, x, y);
  if(c%d!=0) return 0;
  int k=c/d;
  x*=k;
  y*=k;
  return 1;
}
ll x, y;
bool ans=liEu(a, b, c, x, y);

若$gcd(a,b)=1$，则$ax=c\pmod b$在$[0,b-1]$上有唯一解。

若$gcd(a,b)=d$，则$ax=c\pmod b$在$[0,b/d-1]$上有唯一解。

$\therefore$我的求$ax\equiv c\pmod b$：可以但没必要，这样算了两次，慢一些

ll x, y, ste=(a,b);
if(c%ste) return -1;
a/=ste, b/=ste, c/=ste;
exgcd(a, c, x, y);//求c=1
x=(x*c)%b;
x=(x%b b)%b;//求最小非负x

中国剩余定理

中国剩余定理（Chinese Remainder Theorem, CRT ）

求解形如$\begin{cases}x&\equiv &a_1\pmod {n_1}\\x&\equiv &a_2\pmod{n_2}\\x&\equiv &a_3\pmod{n_3}\\&\vdots&\\x&\equiv &a_n\pmod{n_k}\end{cases}$的一元线性同余方程，其中$n_1, n_2, n_3, \dots, n_k$两两互质，求最小非负$x$。

求解：

$n= \prod _{i=1} ^kn_i$

第$i$个方程：$m_i:m_i=\frac n{n_i},\\ m_i^{-1}:m_im_i^{-1}\equiv 1\pmod {n_i},\\ c_i: c_i=m_im_i^{-1}.$

$c_i$不对$n_i$取模。

方程组唯一解：$ans=\sum_{i=1}^ka_ic_i\pmod n$

证明：不会证，不想证，懒得证。

//共r行，b[maxn]存n，a[maxn]存a
//辅助 Exgcd
ll b[11], a[11];//b=n, b>a
int r;
ll ans;
int main()
{
    scanf("%d", &r);
    for(int i=0; i<r; i  )//先读M，再读余数
        scanf("%lld%lld", &b[i], &a[i]);
    ll sb=1;
    for(int i=0;i<r;i  )
        sb*=b[i];
    for(int i=0;i<r;i  )
    {
        ll m=sb/b[i], x, y;
        Exgcd(m, b[i], x, y);
        ans=((ans m*x*a[i])%sb sb)%sb;
    }
    printf("%lld\n",ans);
}

//防爆ll
ll mul(ll a,ll b, ll mod){
bool f=0;
if(a<0) f^=1, a=-a;
if(b<0) f^=1, b=-b;
    a%=mod, b%=mod;
    ll ans=0;
    while(b){
    if(b&1) ans=(ans a)%mod;
    b>>=1;
    a=(a a)%mod;
}
if(f) ans=-ans;
    return ans;
}

ll Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x=1, y=0;
return a;
}
ll d=Exgcd(b, a%b, x, y);
ll t=x;
x=y;
y=t-(a/b)*y;
return d;
}

//b[maxn]->n,a[maxn]->a
ll b[11], a[11];//b=n, b>a
int r;
ll zg()
{
ll ans=0;
    ll sb=1;
    for(int i=0;i<r;i  )
        sb*=b[i];
    for(int i=0;i<r;i  )
    {
        ll m=sb/b[i], x, y;
        Exgcd(m, b[i], x, y);
        ll t1=mul(m, x, sb);
        ll t2=mul(t1, a[i], sb);
        ans=((ans t2)%sb sb)%sb;
    //  ans=((ans m*x*a[i])%sb sb)%sb;
    }
    return ans;
}

应用：一些计数问题给出的模不是质数，但发现该模数没有平方因子，即由两两互质的的质数相乘得到，则可分别对素因子进行计算，再用CRT求解答案。

模数不互质解决方法？？

二次剩余

a为p的二次剩余：$p\nmid a$(a不是p的倍数)，且$a\mod p$为平方数。

a为p的非二次剩余：$p\nmid a$(a不是p的倍数)，且$a\mod p$不为平方数。

对二次剩余求解，即解方程：$x^2\equiv a\pmod p$

当p为奇素数时，使用Cipolla算法求解。（若$p$换为非素数，使用N次剩余）

满足$x^2\equiv n\pmod p$的$n$有$\frac {p-1}2$个(不含0)，非二次剩余也是有$\frac {p-1}2$个。

显然，若$x^2\equiv n\pmod p$，对$x kp,\ k\in {N^ }$有相同$n$。

证：对不同$n$给出一对解$x_1,\ x_2,\ x_1<p, x_2<p$（若$x_i>p$先模$p$即可同$n$），则有$x_1^2\equiv x_2^2\pmod p$，则$p\mid (x_1-x_2)(x_1 x_2)$，因为$\mid x_1-x_2\mid<p$且$p\nmid x_1-x_2$，所以$p\mid x_1 x_2$，满足的$x_i$对共$\frac {p-1}2$种，由于为充要条件，所以$x_1 x_2\neq p$是，对应的$n$必不同。

故只需令$x=1,2,3\dots \frac{p-1}2$则可得除全部$x_1,x_2$与$n$的对应。

勒让德符号

$(\frac np)=\begin{cases} 1,&p\nmid n 且n是p的二次剩余\\ -1,&p\nmid n 且n不是p的二次剩余\\ 0,&p\mid n\end{cases}$

欧拉判别准则

$(\frac np)\equiv n^{\frac{p-1}2}\pmod p$

证明：由费马小定理$n^{p-1}\equiv 1\pmod p$则有$(n^{p-1}2-1)(n^{p-1}2 1)\equiv0\pmod p$，两者中必有一个为$p$的倍数，所以$n^{\frac{p-1}2}\mod p$必为$ 1$或$-1$，然后还没看懂。

Cipolla算法

解方程$x^2\equiv n\pmod p$：

找一个数$a$满足$a^2-n$为非二次剩余，随机数生成$a$，再用欧拉判别判断是否为非二次剩余。

自建一个类似复数域，令$i^2=a^2-n$，就可将所有数表示为$a bi$形式，$a,b$是模意思上的数。

方程的解为$x\equiv n^\frac12\equiv(a i)^{\frac{p 1}2}\pmod p$，$i^2$开方后有正负两个值，带入对应两个解。

特殊的：$n=0$，因为$p$为奇素数，所以只有一个解为0。其他，因为$a$随机选取，随机返回一个解$x_1$，另一个解为$p-x_1$。

特殊的：$p=2$，只有一个解为$n\mod p$。

复杂度为：先测试是否为二次剩余$O(\log(\frac{p-1}2))$。若为二次剩余再随机数生成一个非二次剩余，期望为$O(1.5)$。再建立复数域，计算复数的$\frac{p-1}2$次方，复杂度为$O(\log(\frac{p-1}2))$。故总复杂度，若为非二次剩余为$O(\log(\frac{p-1}2))$；若为二次剩余为$O(\log(\frac{p-1}2) 1.5 \log(\frac{p-1}2))=O(2\log(\frac{p-1}2) 1.5)$。

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

int t;
ll n, p;
ll w;

struct num {  //建立一个复数域
ll x, y;//实部单位为1，虚部单位为i
};

num mul(num a, num b, ll p) {  //复数乘法
num ans = {0, 0};
ans.x = ((a.x * b.x % p   a.y * b.y % p * w % p) % p   p) % p;
ans.y = ((a.x * b.y % p   a.y * b.x % p) % p   p) % p;
return ans;
}

ll binpow_real(ll a, ll b, ll p) {  //实数快速幂
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans % p;
}

ll binpow_imag(num a, ll b, ll p) {  //复数快速幂
num ans = {1, 0};
while (b) {
if (b & 1) ans = mul(ans, a, p);
a = mul(a, a, p);
b >>= 1;
}
return ans.x % p;//最后虚部会等于0
}

ll cipolla(ll n, ll p) {//返回x_1, x_2=p-x_1, 返回大数或小数是随机的，
//因为随机获得a和虚部i^2=a^2-n
n %= p;
if (n == 0) return 0;
if (p == 2) return n;
if (binpow_real(n, (p - 1) / 2, p) == p - 1) return -1;
ll a;
while (1) {  //生成随机数再检验找到满足非二次剩余的a
a = rand() % p;
w = ((a * a % p - n) % p   p) % p;
if (binpow_real(w, (p - 1) / 2, p) == p - 1) break;
}
num x = {a, 1};
return binpow_imag(x, (p   1) / 2, p);
}

int main(){
srand(time(0));//不加也行
int T;
cin>>T;
while(T--){
scanf("%lld%lld", &n, &p);
if(n==0){
printf("0\n");
continue;
}
if(p==2){
printf("%lld\n", n%p);
continue;
}
//求p为奇素数,且n%p!=0,无解或解有两个
ll x1=cipolla(n, p);
if(x1==-1){
printf("Hola!\n");
continue;
}
ll x2=p-x1;
if(x1>x2) swap(x1, x2);
printf("%lld %lld\n", x1, x2);
}

return 0;
}

N次剩余

求解$x^A\equiv B\pmod m$，$m$可以为任意正整数，不限于素数。大体思路：

完全分解$m$，可得到$\prod p_i^{a_i}$，由CRT可求出$x^A\equiv B\pmod{p_i^{a_i}}$的解。

在每个方程中选出一个解，就可对应一个方程组（原方程）的解。

各方程的解数相乘，就是方程组解的个数。

要求所有的解，就要枚举每个方程选取每个解时，计算方程组解。（类似dfs）

下面求解$x^A\equiv B\pmod M$，其中$M=p^a$：

当$B\equiv0\pmod M$，则令$x=p^t$，$t\times A\geq a$即可。

当$\gcd(B,M)=1$，求$M$原根$g$（$p^a$形式的$M$一定有原根），$\dots$

当$\gcd(B,M)\neq1$，则$\gcd(B,M)$必为$p^q$，$0<q<a$的形式，则令$x=p^t$。要求$p^{t\times A}\equiv p^q\pmod{p^a}$，等价$t\times A\equiv q\pmod a$，求$t$，若有解则必$A\mid q$。$\dots$

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

int A, B, mod;
//快速幂
int pow(int x, int y, int mod = 0, int ans = 1) {
    if (mod) {
        for (; y; y >>= 1, x = (LL) x * x % mod)
            if (y & 1) ans = (LL) ans * x % mod;
    } else {
        for (; y; y >>= 1, x = x * x)
            if (y & 1) ans = ans * x;
    }
    return ans;
}

//分解m，求x^A=B mod m，转化为x^A=B mod p_i^a_i，
//p_i^a_i为m的完全分解，从而转换为中国剩余定理可解决的问题
//tot为本质不同p_i个数
//prime为p_i从小到大，expo为p_i对应a_i，pk为对应p_i^a_i
//phi(i)为p_i^a_i的欧拉函数
struct factor {
    int prime[20], expo[20], pk[20], tot;
    void factor_integer(int n) {
        tot = 0;
        for (int i = 2; i * i <= n;   i) if (n % i == 0) {
            prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
            do   expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
              tot;
        }
        if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot  ] = n;
    }
    int phi(int id) const {
        return pk[id] / prime[id] * (prime[id] - 1);
    }
} mods, _p;

//求x对pk[i](M)的逆元
//x^-1=1*x^-1=x^phiM*x^-1=x^(phiM-1) mod M
int p_inverse(int x, int id) {
    assert(x % mods.prime[id] != 0);
    return pow(x, mods.phi(id) - 1, mods.pk[id]);
}
//求ax by=gcd(a,b)
void exgcd(int a, int b, int &x, int &y) {
    if (!b) x = 1, y = 0;
    else exgcd(b, a % b, y, x), y -= a / b * x;
}
//求xret mody=gcd(x, mod), if(ret>=2e31 && (mod&1)) ret
int inverse(int x, int mod) {
    assert(__gcd(x, mod) == 1);
    int ret, tmp;
    exgcd(x, mod, ret, tmp), ret %= mod;
    return ret   ((ret >> 31) & mod);
}

vector<int> sol[20];

void solve_2(int id, int k) {
    int mod = 1 << k;
    if (k == 0) { sol[id].emplace_back(0); return; }
    else {
        solve_2(id, k - 1); vector<int> t;
        for (int s : sol[id]) {
            if (!((pow(s, A) ^ B) & (mod - 1)))
                t.emplace_back(s);
            if (!((pow(s | (1 << (k - 1)), A) ^ B) & (mod - 1)))
                t.emplace_back(s | (1 << (k - 1)));
        }
        swap(sol[id], t);
    }
}

// g^x = B (mod M) => g^iL = B*g^j (mod M) : iL - j
int BSGS(int B, int g, int mod) {
    unordered_map<int, int> map;
    int L = ceil(sqrt(mod)), t = 1;
    for (int i = 1; i <= L;   i) {
        t = (LL) t * g % mod;
        map[(LL) B * t % mod] = i;
    }
    int now = 1;
    for (int i = 1; i <= L;   i) {
        now = (LL) now * t % mod;
        if (map.count(now)) return i * L - map[now];
    }
    assert(0);
    return -1;//不会走到这，防warning 

}

int find_primitive_root(int id) {
    int phi = mods.phi(id); _p.factor_integer(phi);
    auto check = [&] (int g) {
        for (int i = 0; i < _p.tot;   i)
            if (pow(g, phi / _p.prime[i], mods.pk[id]) == 1)
                return 0;
        return 1;
    };
    for (int g = 2; g < mods.pk[id];   g) if (check(g)) return g;
    assert(0);
    return -1;//不会走到这，防warning
}

//  ax = b (mod M)
void division(int id, int a, int b, int mod) {
    int M = mod, g = __gcd(__gcd(a, b), mod);
    a /= g, b /= g, mod /= g;
    if (__gcd(a, mod) > 1) return;
    int t = (LL) b * inverse(a, mod) % mod;
    for (; t < M; t  = mod) sol[id].emplace_back(t);
}

void solve_p(int id, int B = ::B) {
    int p = mods.prime[id], e = mods.expo[id], mod = mods.pk[id];
    if (B % mod == 0) {
        int q = pow(p, (e   A - 1) / A);
        for (int t = 0; t < mods.pk[id]; t  = q)
            sol[id].emplace_back(t);
    } else if (B % p != 0) {
        int phi = mods.phi(id);
        int g = find_primitive_root(id), z = BSGS(B, g, mod);
        division(id, A, z, phi);
        for (int &x : sol[id]) x = pow(g, x, mod);
    } else {
        int q = 0; while (B % p == 0) B /= p,   q;
        int pq = pow(p, q);
        if (q % A != 0) return;
        mods.expo[id] -= q, mods.pk[id] /= pq;
        solve_p(id, B);
        mods.expo[id]  = q, mods.pk[id] *= pq;
        if (!sol[id].size()) return;

        int s = pow(p, q - q / A);
        int t = pow(p, q / A);
        int u = pow(p, e - q);

        vector<int> res;
        for (int y : sol[id]) {
            for (int i = 0; i < s;   i)
                res.emplace_back((i * u   y) * t);
        }
        swap(sol[id], res);
    }
}

vector<int> allans;
void dfs(int dep, int ans, int mod) {
    if (dep == mods.tot) {allans.emplace_back(ans); return;}
    int p = mods.pk[dep], k = p_inverse(mod % p, dep);
    for (int a : sol[dep]) {
        int nxt = (LL) (a - ans % p   p) * k % p * mod   ans;
        dfs(dep   1, nxt, mod * p);
    }
}
//x^A=B mod mod
//m可以不为素数，输入顺序为A, mod, B
//输出第一行为不同解的个数，保证所有解的个数和小于等于1e6
//第二行为解，从小到大排序
void solve() {
    cin >> A >> mod >> B, mods.factor_integer(mod);
    allans.clear();
    for (int i = mods.tot - 1; ~i; --i) {
        sol[i].clear();
        mods.prime[i] == 2 ? solve_2(i, mods.expo[i]) : solve_p(i);
        if (!sol[i].size()) {return cout << 0 << '\n', void(0);}
    }
    dfs(0, 0, 1), sort(allans.begin(), allans.end());
    cout << allans.size() << '\n';
    for (int i : allans) cout << i << ' '; cout << '\n';
}

int main() {
    //ios::sync_with_stdio(0), cin.tie(0);

    int tc; cin >> tc; while (tc--) solve();
    return 0;
}

卢卡斯定理

$m$为素数时，用于求解大组合数取模的问题。

模数为$p$，定理内容：$C_n^m\bmod p=C_{\lfloor{\frac n p}\rfloor}^{\lfloor{\frac m p}\rfloor}\times C_{n\bmod p}^{m\bmod p}\bmod p$

其中：$m\bmod p,\ n\bmod p\ <p$，可直接计算，

$C_{\lfloor{\frac n p}\rfloor}^{\lfloor{\frac m p}\rfloor}$继续使用卢卡斯定理求解，$m=0$时结束，返回$1$。

复杂度：$O(f(p) g(n)\log n)$，$f(p)$为预处理的复杂度，$g(n)$为单次求组合数的复杂度。

ll Lucas(ll n, ll m, ll p) {
  if(m==0) return 1;
  return (C(n%p, m%p, p)*Lucas(n/p, m/p, p))%p;
}

exLucas定理

$m$不为素数。

完全分解：$n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}$

结合中国剩余定理，列出$i$个方程$C_n^m\equiv a_i\pmod{p_i^{q_i}}$并求解。

然而，$p_i^{q_i}$也不一定为质数，以下为求解$C_n^m\bmod p^t$：

BSGS

BSGS（baby-step gaint-step），大步小步算法。求解离散对数，可以在$O\sqrt p$内求解$a^x\equiv b\pmod p$。注意区分N次剩余，这里求解指数，而N次剩余求解底数，且这里是$p$非$m$，要求$\gcd(a,p)=1$。得到$[0, p)$中所有满足的$x$。

实际上，就算$p$不是质数，只要$\gcd(a,m)=1$就可以用BSGS求解。下面我就写奇素数$p$了。

$\because a^0\equiv1\pmod p$且$a^{\varphi(p)}\equiv1\pmod p\Rightarrow a^{p-1}\equiv 1\pmod p$，出现循环，则$[0,p-2]$为一个简化剩余系，其间出现所有（p-1个不同的）$b$。

因此，只需在$[0, p-2]$中枚举$x$，暴力即可得到解。复杂度为$p$。

改变枚举方式，将一轮枚举变为两轮枚举：设$[0, p-2]$内的两个数$x_1\geq x_2$，令$a^{x_1 x_2}\equiv b\pmod p$，则$a^{x_1}\equiv b\times a^{-x_2}\pmod p$。则原问题等价于是否存在$(x_1, x_2)$对，使得上式成立。先枚举$x_2$存map（或者hash一下），再枚举$x_1$看左边值是否出现过，若出现过，则为一组解，枚举次数总和仍为$p$。

进一步约束$(x_1,x_2)$的相对关系，减少枚举量：创建约束$x=x_1\lceil\sqrt p\rceil-x_2$，当$x_2\geq\sqrt p$，$x_2-=\sqrt p,\ x_1 =1$，使得$[0,p-2]$的$x$可在$[0,\lceil\sqrt p\rceil]$的$x_1,x_2$中表示出来。同上，依次枚举$x_2, x_1$，获得答案。时间复杂度$2\sqrt p$，若用map就加个$\log(\sqrt p)$。

当$p$换为$m$非质数时，仍要满足$\gcd(a,m)=1$。$x$在$[0,l]$为一个循环节，计算所得$b$对应$m$的简化剩余系，故每个$b$与$m$互质。在以后的循环节，其它解依次加上$l$。以下模板为求在第一个循环节中的解。

kuangbin板子：

//hash
#define MOD 76543
int hs[MOD], head[MOD], nex[MOD], id[MOD], top;
void insert(int x,int y)
{
    int k = x%MOD;
    hs[top] = x, id[top] = y, nex[top] = head[k], head[k] = top  ;
}
int find(int x)
{
    int k = x%MOD;
    for(int i = head[k]; i != -1; i = nex[i])
        if(hs[i] == x)
            return id[i];
    return -1;
}

//a^x=b mod n, n可以不为素数，只要(a,n)=1即可
//无解返回-1，有解返回一个最小（在[0,l-1]中的）整数解
//l为a关于模m的阶，之后是循环，要全部解，依次之后加上l即可
int BSGS(int a,int b,int n)
{
    memset(head,-1,sizeof(head));
    top = 1;
    if(b == 1)return 0;
    int m = sqrt(n*1.0), j;
    ll x = 1, p = 1;
    for(int i = 0; i < m;   i, p = p*a%n)insert(p*b%n,i);
    for(ll i = m; ;i  = m)
    {
        if( (j = find(x = x*p%n)) != -1 )return i-j;
        if(i > n)break;
    }
    return -1;
}
int main(){
int a, b, p;
while(scanf("%d%d%d", &p, &a, &b)!=EOF){
int res=BSGS(a, b, p);
if(res==-1)printf("no solution\n");
else printf("%d\n", res);
}
return 0;
}

扩展BSGS

解决BSGS中$\gcd(a,m)\not=1$的情况。

解方程$a^x\equiv b\pmod m$，其中$d=\gcd(a,m),\ d\not=1$。

若$d\not\mid b$，则除$x=0,\ b=1$外无其他$(x,b)$对满足。

若$d\mid b$，则$a^x\equiv b\pmod m$，令$k$为使得$d_2=\gcd(a^k,m)$最大前提下的最小整数，即不断让左侧乘$a$，直到$d_2$不再增大，则$a^{x-k}\times \frac{a^k}{d_2}\equiv\frac{b}{d_2}\pmod{\frac m{d_2}}$，

问题变为$a^{x-k}\equiv\frac b{d_2}\times(\frac{a^k}{d_2})^{-1}\pmod {\frac m{d_2}}$，$\gcd(a,\frac m{d_2})=1$。

若$d_2\not\mid b$，显然只有$x=0,b=1$一对解；

对$d_2\mid p$，只需在$[0,k]$枚举$x$，后面部分直接BSGS，注意求得$x$要加上$k$。

以下采用将方程化为$a^{x-k}\equiv\frac b{d_2}\times(\frac{a^k}{d_2})^{-1}\pmod {\frac m{d_2}}$，$\gcd(a,\frac m{d_2})=1$的形式后，直接使用BSGS。所以，指数$k$最多$31$，加号连接所以前面部分忽略不计，直接考虑BSGS复杂度为$k 2\sqrt{\frac m{d_2}}$，不计$k$，则复杂度为$2\sqrt{\frac m{d_2}}$。也是只求出最小循环节的解，要别的继续加$\text{ord}_{\frac m{d_2}}a$。

kuangbin板子：

//hash用的，计算中不管数组大小
const int  MAXN= 99991;
struct LINK{
    ll data;
    ll j;
    ll next;
}HASH_LINK[1000000];
ll ad, head[MAXN];

ll Gcd(ll a, ll b){
return b ? Gcd(b, a % b) : a;
}

ll Ext_Gcd(ll a, ll b, ll &x, ll &y){
    if(!b){
       x = 1; y = 0;
       return a;
    }
    ll r = Ext_Gcd(b, a % b, x, y);
    ll t = x; x = y; y = t - a / b * y;
    return r;
}

ll POWER(ll a, ll b, ll c){
    ll ans = 1;
    while(b){
       if(b & 1) ans = ans * a % c;
       a = a * a % c;
       b >>= 1;
    }
    return ans;
}

void init(){
    memset(head, -1, sizeof(head));
    ad = 0;
}

ll Hash(ll a){
    return a % MAXN;
}

void INSERT_HASH(ll i, ll buf){
    ll hs = Hash(buf), tail;
    for(tail = head[hs]; ~tail; tail = HASH_LINK[tail]. next)
       if(buf == HASH_LINK[tail]. data) return;
    HASH_LINK[ad]. data = buf;
    HASH_LINK[ad]. j    = i;
    HASH_LINK[ad]. next = head[hs];
    head[hs] = ad   ;
}

ll exBSGS(ll a, ll b, ll c){
    ll i, buf, m, temp, g, D, x, y, n = 0;
    for(i = 0, buf = 1; i < 100; i   , buf = buf * a % c)
       if(buf == b) return i;
    D = 1;
    while((g = Gcd(a, c)) != 1){
       if(b % g) return -1; // g | b 不满足，则说明无解
       b /= g;
       c /= g;
       D = D * a / g % c;
          n;
    }
    init();
    m = ceil(sqrt((long double) c));
    for(i = 0, buf = 1; i <= m; buf = buf * a % c, i   ) INSERT_HASH(i, buf);
    for(i = 0, temp = POWER(a, m, c), buf = D; i <= m; i   , buf = temp * buf % c){
       Ext_Gcd(buf, c, x, y);
       x = ((x * b) % c   c) % c;
       for(ll tail = head[Hash(x)]; ~tail; tail = HASH_LINK[tail].next)
           if(HASH_LINK[tail]. data == x) return HASH_LINK[tail].j   n   i * m;
    }
    return -1;
}
//a^x=b mod m
int main(){
ll a, b, m;
while(scanf("%lld%lld%lld", &a, &m, &b)!=EOF && (a || b || m)){
ll res=exBSGS(a, b, m);
if(res==-1) printf("No Solution\n");
else printf("%lld\n", res);
}

}

原根

$\gcd(a,m)=1$，使$a^l\equiv1\pmod m$成立的最小$l$，当然$l>0$，称为$a$关于模$m$的阶，记作$\text{ord}_ma$。

作用：$\text{ord}_ma=l$，$\because\ a^0\equiv1\pmod m,\ a^l\equiv1\pmod m$，模数重复了，之后等价于模数依次乘$a$，故模数变化相同。因此，$a^0,a^1,a^2,\dots,a^{l-1}$是对模$m$的第一个循环节，其中没有相同的余数，此后一直循环。

即：若$g$是$m$的原根，则$g^1,g^2,\dots ,g^{\varphi(m)}$在模$m$意义下两两不同，恰好为$1$到$p-1$间与$m$互质的数，共$\varphi(m)$个。

$\text{ord}_ma=l$，则$\text{ord}_ma^t=\frac l{(t,l)}$，显然找最小正整数$x$满足$l\mid tx$，故$x=\frac l{(t,l)}$。

完全分解$m=p_1^{a_1}p_2^{a_2}\dots p_k^{a_k}$，显然$\text{ord}_ma=[\text{ord}_{p_1}^{a_1},\text{ord}_{p_2}^{a_2},\dots,\text{ord}_{p_k}^{a_k}]$

特殊的，当$m$为素数$p$：$\text{ord}_pa=l$，显然$l\mid\varphi(p)$，则对$p$阶为$l$的所有$a$本质不同$a\mod p$有$\varphi(l)$个。eg：$\text{ord}_{17}2=8$，$\varphi(8)=4$，当$a=2,8,9,15$时，有$l=8$。

特殊的：

$(a, m)=1$，当$\text{ord}_ma=\varphi(m)$时，称$a$为$m$的一个原根。

$a$为$m$的一个原根充要条件为$a,a^2,\dots ,a^{\varphi(m)}$构成模$m$的一个既约剩余系。因为$\gcd(a, m)=1$，所以不会出现$a^x\equiv b\pmod p$使$\gcd(b, m)!=1$，又因为$a^0\equiv a^{\varphi(m)-1}\equiv1\pmod m$，形成循环，所以，当$(a,m)=1$且$\text{ord}_ma=\varphi(m)$一定形成一个简化剩余系。
（简化剩余系也称既约剩余系或缩系，是$m$的完全剩余系中与$m$互素的数构成的子集）

判断原根存在：存在原根的$m$一定为：$2,4,p^a,2p^a$的形式，$p$为奇素数，$a$为正整数。(充要条件)

求一个原根：$(g,m)=1$，设$p_1,p_2,\dots,p_k$为$\varphi(m)$的所有不同素因数，则$g$
为$m$原根的充要条件为：对所有素因数都有$g^{\frac {\varphi(m)}{p_i}}\not\equiv 1\pmod m$。

只需枚举最小原根，在$O(m^{0.25})$级别可枚举到。

求所有原根：先求出一个，则集合$S=\{g^s\mid 1\leq s\leq\varphi(m),\gcd(s,\varphi(m))=1\}$，为$m$的全部原根。即：若$m$存在原根则一定有$\varphi(\varphi(m))$个。

自己写的模板，虽然过了题，但是很慢，要是tle就提前打素数表，分解因式时大于2个判断，并提前退出，其他不变。

#define ll long long

ll qpm(ll a, ll b, ll m) {
a%=m;
ll res=1;
while(b>0) {
if(b&1) res=res*a%m;
a=a*a%m;
b>>=1;
}
return res;
}
//来自N次剩余
struct factor {
int m;
    int prime[20], expo[20], pk[20], tot;
    void factor_integer(int n) {
    m=n;
        tot = 0;
        for (int i = 2; i * i <= n;   i) if (n % i == 0) {
            prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
            do   expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
              tot;
        }
        if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot  ] = n;
    }
    int phi() const {//与N次剩余此处修该
    int res=m;
    for(int i=0; i<tot; i  )
        res=res/prime[i] * (prime[i] - 1);
        return res;
    }
}f, fp;
bool check(int a, int m){//O(m^0.25)
for(int i=0; i<fp.tot; i  ){
if(qpm(a, f.phi()/fp.prime[i], m)==(ll)1)
return 0;
}
return 1;
}

//O(m^0.25)找到最小原根，再O(m)实验全部根
int rt(int m){
f.factor_integer(m);
if(m==2 || m==4 || (f.tot==1 && f.prime[0]!=2) || (f.tot==2
&& f.pk[0]==2));
else return -1;//不是2,4,p^2,2p^a，p是奇素数的形式
fp.factor_integer(f.phi());
for(int i=1; i<m; i  ){
if(__gcd(i, m)!=1) continue;
if(check(i, m)) return i;
}
return -1;//不会到这一步
}
vector<int> v;
void art(int m){
v.clear();
int mrt=rt(m);
if(mrt==-1){
v.push_back(-1);
return;
}
for(int i=1; i<=f.phi(); i  ){
if(__gcd(i, f.phi())!=1) continue;
v.push_back(qpm(mrt, i, m));
}
return;
}
void p(){
if(v[0]==-1){
printf("0\n\n");
return;
}
int si=v.size();
printf("%d\n", si);
sort(v.begin(), v.end());
for(int i=0; i<si; i  ){
printf("%d ", v[i]);
}
printf("\n");
return;
}
int main(){
int m, T;
cin>>T;
while(T--){
scanf("%d", &m);
art(m);
p();
}
return 0;
}

杜教筛

用于求一些数论函数（积性函数）的前缀和，即对$f$，求$S(n)=\sum_{i=0}^nf(i)$。

想要快于$O(n)$的快速计算，要构造一个$S(n)$关于$S(\lfloor\frac ni\rfloor)$的递推式。

结论是：$g(1)S(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{i=2}^ng(i)S(\lfloor\frac ni\rfloor)$。

后面部分用数论分块可快速计算，只需要找到一个使前面部分快速计算的$g$即可。

常见数论函数前缀和对应的$g$函数：

（1）$f=\mu,\ g=u(\equiv1),\ I=\mu\ast u,$

$\therefore S(n)=1-\sum_{i=2}^nS(\lfloor\frac ni\rfloor)$

（2）$f=\varphi,\ g=u,\ e=\varphi\ast u,$

$\therefore S(n)=\frac 12n(n 1)-\sum_{i=2}^nS(\lfloor\frac ni\rfloor)$

乘法逆元

线性同余方程$ax\equiv1\pmod b$，则称$x$为$a\mod b$的逆元，记作$a^{-1}$。

费马小定理：$p$为质数且$gcd(a,p)=1\Rightarrow a^{-1}=a^{p-2}$，快速幂即可

扩展欧几里得：n

用于求$ax by=gcd(a, b)$的一组解。

两边同除以$gcd(a, b)$后，等价于$gcd(a, b)=1$，求$ax by=1$的一组解。

$a=q_1b r_1,\quad r_1=a\%b$

$b=q_2r_1 r_2,\quad r_2=b\%r_1$

$r_1=q_3r_2 r_3,\quad r_3=r_1\%r_2$

$\dots$

$r_{n-1}=q_{n 1}r_n,\quad r_{n 1}=0$，结束。

$gcd(a,b)=r_n$；

则有$Q_ka-P_kb=(-1)^{k-1}r_k,\ k=1,2,3,\dots,n$，

其中$P_0=1,\ P_1=q_1,\ P_k=q_kP_{k-1} P_{k-2}$

$Q_0=0,\ Q_1=1,\ Q_k=q_kQ_{k-1} Q_{k-2},\quad 2<=k<=n$

计算：

我的想法：

使用元组$(x,\ y)$表示$ax by$的值，即：$a=(1,0),\ b=(0,1)$。有：

$r_1=a-q_1b=(1,-q_1)$

$r_2=b-q_2r_1=(0, 1)-q_2(1,-q_1)=(-q_2，q_1q_2 1)$

$r_3=r_1-q_3r_2=(1,-q_1)-q_3(-q_2,q_1q_2 1)$

$\dots$

$r_n=r_{n-2}-q_nr_{n-1}=(x_n,y_n)$

$r_{n 1}=0=r_{n-1}-q_{n 1}r_n=(x_{n 1},y_{n 1})$

$\because r_n=gcd(a,b)$

$\therefore r_n=(x_n,y_n)$为$ax by=gcd(a,b)$的一个解。

exgcd：

$ax_1 by_1=gcd(a, b)$

$bx_2 (a\mod b)y_2=gcd(b, a\mod b)=gcd(a, b)$

$ax_1 by_1=bx_2 (a\mod b)y_2$

$a\mod 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$

$x_1=y_2,\quad y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2$

ll Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x=1, y=0;
return a;
}
ll d=Exgcd(b, a%b, x, y);
ll t=x;
x=y;
y=t-(a/b)*y;
return d;
}
ll Inv(ll a, ll m){
ll x, y;
ll d=exgcd(a, m, x, y);
if(d==1) return (x%m m)%m;
    return -1;
}

线性递推：1到n

$O(n)$

求$a$在膜$m$下的逆元：

$1^{-1}\equiv1\pmod m$

$m=qa r,\quad qa r\equiv0\pmod m$

两边同乘以$a^{-1},\ r^{-1}$：

$qr^{-1} a^{-1}\equiv0\pmod m$

$a^{-1}\equiv-qr^{-1}\pmod m$

$a^{-1}\equiv-(\frac{m}{a})\times(m\bmod a)^{-1}\pmod m$

避免负号出现，原式加上$m$：

$a^{-1}\equiv(m-\frac{m}{a})\times(m\bmod a)^{-1}\pmod m$

inv[1]=1;
for(int i=2; i<=n; i  ) inv[i]=(ll)(M-M/i)*inv[M%i]%M;

线性递推：1!到n!

$O(n \log p)$

可用于求任意$n$个数的逆元，其中$1\leq a_i<p$，因为倒着递推，不必计算$1$到$\min(a_i)-1$的逆元。

使用费马小定理或扩展欧几里得计算最后一个数$a_t$的前缀积的逆元，

乘以最后一个数后得到$a_{t-1}$的前缀积的逆元，以此类推。

void init(){
fact[0]=1;
for(int i=1; i<maxn; i  ){
fact[i]=fact[i-1]*i%mod;
}
inv[maxn-1]=qpm(fact[maxn-1], mod-2);
for (int i=maxn-2; i>=0; i--){
inv[i]=inv[i 1]*(i 1)%mod;
}
}

线性递推：给定随机n个数

$1\leq a_i\leq p$的$n$个数存在$a_1$到$a_n$中，类似$1!$到$n!$，先求前缀积$s_n$，求$s_n$逆元$invs_n$，线性递推每个$s_i$逆元$invs_i$，最后一步$inv_i=invs_i\times s_{i-1} \mod p$。

//a[maxn]存给定数，inv[maxn]存对应逆元，n为a的大小，从1开始存
//辅助 s[maxn], invs[maxn], qpm，p=M
void init() {
s[0]=1;
for (int i=1; i<=n; i  ) s[i]=s[i-1]*a[i]%p;
invs[n]=qpm(s[n], p-2);
for (int i=n; i>=1; i--) invs[i-1]=invs[i]*a[i]%p;
for (int i=1; i<=n; i  ) inv[i]=invs[i]*s[i-1]%p;
}

常见数论函数

1）u(n), e(n), I(n)

$u(n)\equiv1,\ n\geq1$

$e(n)=n,\ n\geq1$

$I(n)=\begin{cases}1,\quad n=1\\0,\quad n>1\end{cases}$

2）d(n) 约数个数

$n$的所有正除数的个数，含1和$n$本身：$d(n)=\sum_{d\mid n}1$

完全分解：$n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}$，$d(1)=1$

$d(n)=\prod_{i=1}^s(k_i 1)$

可乘函数可使用线性筛：

$p$存放质数，用d表示约数个数，num表示最小质因子出现次数。考虑$n>1$。

$d(p)=2, num(p)=1$

设$d=\frac n p$，$p$为$n$的最小质因子：

$p\mid d:\quad num(n)=num(d) 1,\quad d(n)=\frac{d(d)\times (num(n) 1)}{num(d) 1}$

考虑$d,p$互质，因为为可乘函数：

$p\nmid d:\quad num(n)=1,\quad d(n)=d(d)\times d(p)$

//数组 d, num, p
//bool v   int tot, n
void pre(){
d[1]=1, tot=0;
for(int i=2; i<=n; i  ){
if(!v[i]) p[  tot]=i, d[i]=2, num[i]=1;
for(int j=1; j<=tot && i<=n/p[j]; j  ){
v[p[j]*i]=1;
if(i%p[j]==0){
num[i*p[j]]=num[i] 1;
d[i*p[j]]=d[i]/(num[i] 1)*(num[i*p[j]] 1);
break;
}
num[i*p[j]]=1;
d[i*p[j]]=d[i]*d[p[j]];
}
}
}

3）Omega(n) 可重复素因子个数

$n$的全部素因子的个数$\Omega(n)$：

$\Omega(1)=0$

$\Omega(n)=\sum_{i=1}^sa_i$

int b[maxn];
int cnt[maxn];
void pre() {
for(int i=2; i<maxn; i  ) {
if(b[i]!=0) continue;
for(int j=i; j<maxn; j =i) {
if(b[j]==0) b[j]=i;
}
}
cnt[2]=1;
for(int i=3; i<maxn; i  )
cnt[i]=cnt[i/b[i]] 1;
return;
}

4）omega(n) 不可重复素因子个数

$n$的不同素因子个数$\omega(n)$：

$\omega(1)=0,\quad\omega(n)=s$

bool vis[maxn];
int cnt[maxn];
void pre() {
for(int i=2; i<maxn; i  ) {
if(vis[i]!=0) continue;
cnt[i]=1;
for(int j=i i; j<maxn; j =i) {
cnt[j]  , vis[j]=1;
}
}
return;
}

5）sigma_lambda(n) 约数幂和

$n$的正除数的幂和函数$\sigma_\lambda(n),\ \lambda$为实数：

$\sigma_\lambda(n)=\sum_{d\mid n}d^\lambda$

特殊的，约数和：$\lambda=1$时，$\sigma(n)=\sum_{d\mid n}d$

完全分解：$n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}$，$d(1)=1$

$\sigma(n)=\prod_{i=1}^s(1 p_i p_i^2 \dots p_i^{k_i})=\prod_{i=1}^s\sum_{j=0}^{k_i}p_i^j$

又因：$\sum_{j=0}^{k_i}p_i^j=\frac{p_i^{k_i 1}-1}{p_i-1}$

$\sigma(n)=\prod_{i=1}^s(1 p_i p_i^2 \dots p_i^{k_i})=\prod_{i=1}^s\frac{p_i^{k_i 1}-1}{p_i-1}$

以下计算均采用前者推导。

可乘函数可使用线性筛：

$p$存放质数，$f_i$表示$i$的约数和，$g_i$表示$i$的最小质因子$p$的$\sum_{j=0}^{k}p^j$。考虑$n>1$。

$g(p)=p 1,\quad f(p)=i 1$

设$d=\frac n p$，$p$为$n$的最小质因子：

$p\mid d:\quad g(n)=g(d)\times p 1,\quad f(n)=\frac{f(d)\times g(n)}{g(d)}$

考虑$d,p$互质，因为为可乘函数：

$p\nmid d:\quad g(n)=p 1,\quad f(n)=f(d)\times f(p)$

void pre(){
g[1]=f[1]=1, tot=0;
for(int i=2; i<=n; i  ){
if(!v[i]) p[  tot]=i, g[i]=i 1, f[i]=i 1;
for(int j=1; j<=tot && i<=n/p[j]; j  ){
v[p[j]*i]=1;
if(i%p[j]==0){
g[i*p[j]]=g[i]*p[j] 1;
f[i*p[j]]=f[i]/g[i]*g[i*p[j]];
break;
}
            g[i*p[j]]=p[j] 1;
f[i*p[j]]=f[i]*f[p[j]];
}
}
}

注：求$\sigma(n)$是否为奇数，根据第一个公式，$p_i=2$时，$k_i$为任意数可，对其他$p_i,\ k_i$必须为偶数，所以$n$为平方数或$n/2$为平方数为充要条件。求$1$到$n$中，$\sigma(i)$为奇数的个数，$ans=sqrt(n) sqrt(n/2)$即可。

6）varphi(n) 互质个数

欧拉函数$\varphi(n)$：小于等于$n$，和$n$互质的个数。

1与任何正整数互质，$\varphi(1)=1$

求解：完全分解：$n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s}$

推导表达式，由容斥原理：$\varphi(n)=n (-\frac{n}{p_1}-\frac{n}{p_2}-\dots-\frac{n}{p_s}) (\frac{n}{p_1p_2} \frac{n}{p_1p_3} \dots \frac{n}{p_{s-1}p_s})\dots\varphi(n)=\sum_{S\subseteq\{p_1,p_2\dots p_s\}}(-1)^{\mid S\mid}\frac{n}{\prod_{p_i\in S} p_i}=n\prod_{i=1}^s(1-\frac{1}{p_i})$

$\varphi(n)=n\times\prod_{i=1}^{s}\frac{p_i-1}{p_i}$

int phi(int n) {
int m=(int)sqrt(n 0.5);
int ans=n;
for (int i=2; i<=m; 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;
}

性质：

（1）为可乘函数，$gcd(a,b)=1$：$\varphi(a\times b)=\varphi(a)\times\varphi(b)$。

（2）特殊的，$n$为奇数：$\varphi(2n)=\varphi(n)$，$\because 2\times\frac{2-1}{2}=1$

（3）$n=\sum_{d\mid n}\varphi(d)$，含$1$和$n$，证明：

法一：要证明$e(n)=n\ (n\geq1)$为$\varphi(n)$的莫比乌斯变换，

只需证明$\varphi(n)$为$e(n)$的反莫比乌斯变换即可。

$e (n)$的反莫比乌斯变换：$\varphi(n)=\sum_{d\mid n}\mu(d)\frac n d$

$\varphi(n)$从定义式正经变换：$\varphi(n)=\sum_{1\leq d\leq n, (d,n)=1}1=\sum_{1\leq d\leq n}\sum_{l\mid (d, n)}\mu(l)=\sum_{l\mid n}\mu(l)\sum_{1\leq d\leq n,l\mid d}1=\sum_{l\mid n}\mu(l)\frac n l$

可见上下两式相等，故$\varphi(n)$为$e(n)$的反莫比乌斯变换。

法二：$gcd(a,n)=d$，按$d$对$a$进行分类，令$n=dk$则$gcd(\frac a d, k)=1$。固定的$d$对应$a$的个数为

$\varphi(k)=\varphi(\frac n d)$，所以$n=\sum_{d\mid n} \varphi(\frac n d)=\sum_{d\mid n} \varphi(d)$。

（4）特殊的$n=p^k$：$\varphi(n)=p^k-p^{k-1}$

$=(p-1)\times p^{k-1}=n\times\frac{p-1}{p}$，由定义可得。

可用于求$\sum_{i=1}^{n}\sum_{j=i 1}^n gcd(i, j)$：

显然$ans[n]=ans[n-1] \sum_{i=1}^{n-1}gcd(i,n)$，仅需证明此处$\sum_{i=1}^{n-1}gcd(i,n)=\sum_{d\mid n}\varphi(\frac n d)\times d$：

按$d$对$a$进行分类，令$n=dk$则$gcd(\frac a d, k)=1$。固定的$d$对应$a$的个数为

$\varphi(k)=\varphi(\frac n d)$，和为$\varphi(\frac n d)\times d$。$f(n)=\sum_{d\mid n}\varphi(\frac n d)\times d=(\varphi\ast e)(n)$，解题关键是$f(n)$的表达式也可用线性筛，所以不会超时。因为此题$d<n$而不是$d\leq n$，故内层循环从$i*2$开始，而不是从$i$开始。$(\varphi\ast e)(n)$筛法：

for(int i=1; i<=n; i  ){
for(int j=i i; j<=n; j =i){
f[j] =(ll)i*(ll)phi[j/i];
}
}

筛法：

void phi_table(int n){
memset(phi, 0, sizeof(phi));
phi[1]=1;
for(int i=2; i<=n; i  )
if(!phi[i])
for(int j=i; j<=n; j =i){
if(!phi[j]) phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}

7）mu(n)

莫比乌斯函数$\mu(n)$：

$\mu(n)=\begin{cases}1,\quad &n=1\\(-1)^s,\quad &n=p_1p_2\dots p_s,\ p_1<p_2<\dots<p_s\\0,\quad &其他\end{cases}$

容斥原理的系数啊：利如，给定集合S，给出q个查询x，求S中与x不互质的数的个数。

注意题目中出现的其他说法：存在$k(k>1)$使$k^2\mod x$则$f(x)=0$，否则$f(x)=1$。

线性筛：

//mu p v
//tot
void pre() {
mu[1]=1;
for(int i=2; i<=maxn-13; i  ) {//注意这里的上界
if(!v[i]) mu[i]=-1, p[  tot]=i;
for(int j=1; j<=tot && i*p[j]<=maxn-13; j  ) {//改上界
v[i*p[j]]=1;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
}

8）Lambda(n)

Mangoldt函数$\Lambda(n)$：

$\Lambda(n)=\begin{cases}\log p,\quad &n=p^k,\ k\geq1\\0\quad &其他\end{cases}$

9）lambda(n)

Liouville函数$\lambda(n)$：

$\lambda(n)=(-1)^{\Omega(n)}$

可乘函数（积性函数）

$f(mn)=f(m)\times f(n),\quad\gcd(m,n)=1$

欧拉函数$\varphi(n)$，莫比乌斯函数$\mu(n)$，正因子和$\sigma(n)$，正因子数$d(n)$。

均可采用线性筛法。

可乘函数的卷积仍是可乘函数，可用线性筛：

求$f=g\ast h$提前处理$g,h$后，线性筛法如下：

for(int i=1; i<=n; i  ){
for(int j=i i; j<=n; j =i){
f[j] =g[i]*h[j/i];
}
}

完全（绝对）可乘函数

$f(mn)=f(m)f(n)$

$u(n),\ e(n),\ I(n)$

Dirichlet 乘积

数论分块

引理1：$\forall a,b,c\in\mathbb{Z},\lfloor\frac{a}{bc}\rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor$

引理2：$\forall n\in N,\vert\{\lfloor\frac{n}{d}\rfloor\mid d\in N\}\vert\leq\lfloor2\sqrt{n}\rfloor$，$\vert V\vert$表示集合$V$的元素个数

解决类似计算$\sum_{i=1}^n \lfloor\frac ki\rfloor$的问题，其中$n,\ k$为常量：

对每个$i(i\leq k)$要找到相应最大$j(i\leq j\leq k)$，使得$\lfloor\frac ki\rfloor=\lfloor\frac kj\rfloor$，即可按$j$分块计算，显然，$j=\lfloor\frac k{\lfloor\frac ki\rfloor}\rfloor$。所以，想办法让求和表达式中出现$\lfloor\frac ki\rfloor$，且其他部分的求和可合并后$O(1)$直接计算即可。

例：求$\sum_{i=1}^n (k\mod i)$

$=\sum_{i=1}^n k-i\lfloor\frac ki\rfloor$

ll n, k;
int main() {
scanf("%lld%lld", &n, &k);
ll ans=n*k;
for(ll l=1, r; l<=n; l=r 1){
if(k/l!=0) r=min(k/(k/l), n);
else r=n;
ans-=(k/l)*(l r)*(r-l 1)/2;
//(k/l)为题干中k/i，后部分为区间求和
}
printf("%lld\n", ans);
return 0;
}

标准计算(只计算到i=n，而不是无穷大)，若要乘以某个函数，提前处理其前缀和。（二维见莫比乌斯反演），若想计算到无穷大，令$n==k$即可。计算$\sum_{i=1}^n\lfloor\frac ki\rfloor$，若$i$固定从$s$开始，而不是$1$，把$l$初始值设为$s$即可：

ll fk(ll n, ll k){//n为i从1到n，k为k/i
    ll ans=0;
for(ll l=1, r; l<=n; l=r 1){
if(k/l!=0) r=min(k/(k/l), n);
else r=n;//l大于k时
ans =k/l*(r-l 1);//此处写带k/i的表达式
        // ans =(s[r]-s[l-1])*(k/l);
}
    return ans;
}

常见卷积

设$f(n),\ h(n)$是两个数论函数，

$h(n)=\sum_{d\mid n}f(d)g(\frac{n}{d})$为$f(n),\ g(n)$的Dirichlet 乘积或卷积，记作：

$h=f\ast g,\quad (f\ast g)(n)=\sum_{ab=n}f(a)g(b),$

满足结合律和交换律。

$f\ast g\ast h=\sum_{abc=n}f(a)g(b)h(c)$

$I(n)=\begin{cases}1,\quad n=1\\0,\quad n>1\end{cases}$为卷积单位元素。

常见卷积：

$I=\mu\ast u=\sum_{d\mid n}\mu(d)$

$d=u\ast u=\sum_{d\mid n}1$

$\sigma=e\ast u=\sum_{d\mid n}d$

$e=\varphi\ast u=\sum_{d\mid n}\varphi(d)$

$\varphi=\mu\ast e=\sum_{d\mid n}d\cdot\mu(\frac nd)$

互为Dirichlet逆：$I=\mu\ast u$

若$F=f\ast u$，则称$F$为$f$的莫比乌斯变换，即：

$F(n)=\sum_{d\mid n}f(d)$，等于$f\ast \mu^{-1}$，即再卷积一个莫比乌斯函数可变为原函数。

还原：$f=F\ast\mu,\ f(n)=\sum_{d\mid n}F(d)\mu(\frac n d)$

常见F为f的莫比乌斯变换有：

$d(n),\quad u(n)$，正除数个数$=\sum_{d\mid n}1$

$I(n),\quad \mu(n)$，$[n=1]=\sum_{d\mid n}\mu(d)$

$e(n),\quad \varphi(n)$，$n=\sum_{d\mid n}\varphi(d)$

$\sigma(n),\quad e(n)$，正除数和$=\sum_{d\mid n}d$

$\log n,\quad\Lambda(n)$，$\log n=\sum_{d\mid n}\Lambda(d)$

莫比乌斯反演

带莫比乌斯函数希望计算不是$O(n)$，先分解质因数，dfs质因数选与不选，虽然可以卡，但总的来说1e8可以开个根。

很难求出原函数的值，而容易求出其倍数和或约数和 ，则可通过莫比乌斯反演简化运算，求得原函数的值。

$F=f\ast u\Rightarrow F\ast\mu=f$

$F(n)=\sum_{d\mid n}f(d)$

$f(n)=\sum_{d\mid n}F(d)\mu(\frac{n}{d})$

因为$[n=1]=\sum_{d\mid n}\mu(d)$，将$n$换为想计数的变量，即可求得该变量等于1的数量，例：

$[\gcd(i,j)=1]\iff\sum_{d\mid\gcd(i,j)}\mu(d)$

数论分块的二维形式，求$\sum \lfloor\frac ni\rfloor\lfloor\frac mi\rfloor$，$i,j$均到无穷大：

int sol(int n, int m) {
int ans=0;
for (int l=1, r; l<=min(n, m); l=r 1) {
r=min(n/(n/l), m/(m/l));
ans =(mu[r]-mu[l-1])*(n/l)*(m/l);
}
return ans;
}

若$i,j$不是到无穷大时，参考数论分块代码。

一些变换：

$\sum_{1\leq i\leq n\\\gcd(i,n)=1} \frac 1i=\sum_{d\mid n}\mu(d)h_{\frac nd}\frac 1d$，其中$h$为调和级数，可直接线性预处理。

FT

快速傅里叶变换

FFT，求多项式卷积。$O(n\log n)$

若求正常定义的卷积$h=f\ast g,\ (f\ast g)(n)=\sum_{a b=n}f(a)g(b)$，值为$h$的$n$次项系数。

输入$f,\ h$从低次到高次多项式的系数表达，求$g=f*h$。$f,\ h$做$DFT$转换为点标记的表示方法$F,\ H$，对应点相乘后，得到$G=F\cdot H$的离散$H$方式的各点，再对$G$做$IFT$转换为系数表达式。

蝴蝶变换：把位置$i$的元素直接移到$i$进行$\log_2len$次递归后的位置，为$i$的二进制数的翻转，eg：$11010\Rightarrow01011$，而取消递归。

$IFT$：把多项式$g$的离散傅里叶变换$G$结果作为另一个多项式的系数，取单位根的倒数即$\omega_n^i\ \{i=0,-1,\dots,-(n-1)\}$作为$x$代入$G$，得到的$G(x)$再除以$n$，就是$g$的各项系数。

注：数组大小乘2。

struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator (const Complex &b) const {
return Complex(x   b.x, y   b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y   y * b.x);
}
};

//  蝴蝶变换 len必须为2的幂，i与二进制反转位置交换
void change(Complex y[], int len) {
int k;
for (int i = 1, j = len / 2; i < len - 1; i  ) {
if (i < j) swap(y[i], y[j]);
k = len / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k) j  = k;
}
}
//  len必须为2的幂
//  on=-1逆运算是带入共轭复数再除n
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j  = h) {
Complex w(1, 0);
for (int k = j; k < j   h / 2; k  ) {
Complex u = y[k];
Complex t = w * y[k   h / 2];
y[k] = u   t;
y[k   h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; i  ) {
y[i].x /= len;
}
}
}

Complex x1[maxn<<2], x2[maxn<<2];
//len为最高次 1，x[i]表示x^i对应的系数
int main() {
int len1, len2, len=1;
cin>>len1>>len2;
len1  , len2  ;
double t;
while(len<len1*2 || len<len2*2) len<<=1;
for(int i=0; i<len1; i  ){
scanf("%lf", &t);
x1[i]=Complex(t, 0);
}
for(int i=len1; i<len; i  ) x1[i]=Complex(0, 0);
for(int i=0; i<len2; i  ){
scanf("%lf", &t);
x2[i]=Complex(t, 0);
}
for(int i=len2; i<len; i  ) x2[i]=Complex(0, 0);
fft(x1, len, 1);
fft(x2, len, 1);
for(int i=0; i<len; i  ) x1[i]=x1[i]*x2[i];
fft(x1, len, -1);
int m=len1 len2-1;
for(int i=0; i<m; i  )
printf("%d ", (int)(x1[i].x 0.5));
printf("\n");
return 0;
}

快速沃尔什变换

$C_k=\sum_{i|j=k}A_i\times B_j$

FWT，将多项式卷积中的$*$换为其他位运算符。

不取模的话，开long long，把逆元改为除2即可。

void FWT_or(int *a, int N, int opt)
{
    for(int i=1;i<N;i<<=1)
        for(int p=i<<1,j=0;j<N;j =p)
            for(int k=0;k<i;  k)
                if(opt==1)a[i j k]=(a[j k] a[i j k])%MOD;
                else a[i j k]=(a[i j k] MOD-a[j k])%MOD;
}
void FWT_and(int *a, int N, int opt)
{
    for(int i=1;i<N;i<<=1)
        for(int p=i<<1,j=0;j<N;j =p)
            for(int k=0;k<i;  k)
                if(opt==1)a[j k]=(a[j k] a[i j k])%MOD;
                else a[j k]=(a[j k] MOD-a[i j k])%MOD;
}
void FWT_xor(int *a, int N, int opt)
{
    for(int i=1;i<N;i<<=1)
        for(int p=i<<1,j=0;j<N;j =p)
            for(int k=0;k<i;  k)
            {
                int X=a[j k],Y=a[i j k];
                a[j k]=(X Y)%MOD;a[i j k]=(X MOD-Y)%MOD;
                if(opt==-1)a[j k]=1ll*a[j k]*inv2%MOD,a[i j k]=1ll*a[i j k]*inv2%MOD;
            }
}

快速幂

ll inv(ll a, ll b) {
ll res=1;
while(b>0) {
if(b&1) res=res*a%M;
a=a*a%M;
b>>=1;
}
return res;
}

快速乘

求a>1e9 && b>1e9但模M时，可防爆ll：

ll mul(ll a,ll b, ll mod){
bool f=0;
if(a<0) f^=1, a=-a;
if(b<0) f^=1, b=-b;
    a%=mod, b%=mod;
    ll ans=0;
    while(b){
    if(b&1) ans=(ans a)%mod;
    b>>=1;
    a=(a a)%mod;
}
if(f) ans=-ans;
    return ans;
}

常见数列

斐波那契数

$F_0=0$，$F_1=1$，$F_n=F_{n-1} F_{n-2}$

前几项：0，1，1，2，3，5，8，13，21，34，55，89

性质：1.$F_{(n-1)}F_{(n 1)}-F_n^2=(-1)^n$

2.$F_{(n k)}=F_kF_{(n 1)} F_{k-1}F_n$ $\Rightarrow$ $F_{2n}=F_n(F_{n 1} F_{n-1})$

3.$\Rightarrow$

4.$gcd$性质：$(F_m,F_n)=F(m,n)$

卡特兰数

$H_0=1$，$H_1=1$

前几项：1，1，2，5，14，42，132

描述： 1.将$n$分为$n_1,n_2$，将$n_1$分为$n_{11},n_{12}$，...分法数，$H_n=\begin{cases} \sum_{i=1}^{n}H_{i-1}H_{n-i},n\geq2,n\in{N_ }\\1,n=0,1 \end{cases}$

2.$num(a)==num(b)==n$，当对任意$p$有：$\sum_{i=1}^{p}a\geq\sum_{i=1}^{p}b$的排列数，

$H_n=\frac{C_{2n}^n}{n 1}=C_{2n}^n-C_{2n}^{n-1}$，$a_1$与$a_2$无区别，或认为$a_i$间顺序已定

一般使用递推式：$H_n=\frac{H_{n-1}\times(4n-2)}{n 1}$

Harmonic(调和级数)

$h_n=\frac 11 \frac 12 \frac 13 \dots \frac 1n$

存储：

$n<2^{31}$，且多次查询。不能开那么大的数组，正常求法是分段打表，仅存100整数倍函数值，再递推。可以整数取逆元，也可算小数。

不太正常的做法，当$n$较大时，只能算小数：

欧拉常数$\gamma=\lim_{n\to \infty}[(\sum_{k=1}^n\frac 1k)-\ln(n)]=\int_1^\infty(\frac 1{\lfloor x\rfloor}-\frac 1x)dx$

$h_n$为发散无穷级数，但$S_n=(\sum_{k=1}^n\frac 1k)-\ln(n)$存在极限$\gamma$。

$\gamma= 0.57721566490153286060651209 \dots$

$\therefore h_n \approx \gamma \ln(n)$

因为$\ln(n) \gamma$永远在$\sum_{k=1}^n \frac 1k$下端，所以多加一项余项不亏，能使差的绝对值减小：

$\therefore h_n \approx\gamma \ln(n) \frac 1{2n}$

常见常数

e

$e=1 \sum \frac 1{i!}= 2.718281828459\dots$

pi

$\arctan=x-\frac{x^3}3 \frac{x^5}5-\dots$，带入$x=1$得到

莱布尼兹级数：$\frac{\pi}4=\sum_{i=0}^\infty\frac{(-1)^k}{2k 1}$

收敛很慢，基本上没用。

快速收敛公式：

拉马努金公式：$\frac1\pi=\frac{2\sqrt2}{99^2}\sum_{k=0}^\infty\frac{(4k!)(1103 26390k)}{(k!)^4396^{4k}}$

每计算一项得到8个有效数字。其变体，现有最快的

Chudnovsky公式：$\frac1\pi=\frac1{53360\sqrt{640320}}\sum_{k=0}^\infty(-1)^k\frac{(6k!)}{(k!)^3(3k)!}\times\frac{13591409 545140134k}{640320^{3k}}$

每计算一项得到14个有效数字。

特殊的：

BBP公式：$\pi=\sum_{k=0}^\infty[\frac1{16^k}(\frac4{8k 1}-\frac2{8k 4}-\frac1{8k 5}-\frac1{8k 6})]$计算一项可得到1.25个有效数字，虽然慢，但可在$O(n)$计算$\pi$的特定位(16进制)。分成四部分计算，两边同乘以$16^{n-1}$整数部分不要，小数部分相加，最后$ans=res-(int)res$，若小于零则加1。计算完前$n-1$项后，可适当添加一些常数加项防止没进位，加22位都过题了。

泰勒公式

$f(x)=\sum_{i=0}^\infty\frac {f^i(x_0)}{i!}(x-x_0)^i$

麦克劳林：令$x_0=0,\ f(x)=\sum_{i=0}^\infty\frac {f^i(0)}{i!}x^i$

常见：

$e^x=\sum_{i=0}^\infty \frac 1{i!}x^i$1

$\ln(1 x)=\sum_{i=0}^\infty \frac{(-1)^i}{i 1}x^{i 1}$2

$\frac 1{1-x}=\sum_{i=0}^\infty x^i$3

$\frac 1{e^n-1}=\sum_{i=1}^\infty e^{-in}$3

组合数学

排列组合

$C_n^k C_n^{k-1}=C_{n 1}^k$

$(x_1 x_2 \cdots x_t)^n=\sum(_{n_1n_2\cdots n_t}^{n})x_1^{n_1}x_2^{n_2}\cdots x_t^{n_t}$，其中含满足$n_1 n_2 \cdots n_t=n$的所有非负整数

其它数学

Cantor

排列组合的字典序，从0开始，结果经常加1。

康托展开预处理后，可在$o(1)$查询较小表格、字符串等映射。

int cantor(string s){
int ans=0;
for(int i=0; i<9; i  ){
int k=0;
for(int j=i 1; j<9; j  )
if(s[i]>s[j]) k  ;
ans =fac[9-i-1]*k;
}
return ans;

}
string decantor(int ans){
int num=0;
string s;
set<int> st;
for(int i=0; i<9; i  ) st.insert(i);
while(num<9){
auto it=st.begin();
for(int i=0; i<ans/fac[9-num-1]; i  ) it  ;
s =*it '0';
st.erase(it);
ans%=fac[9-num-1];
num  ;
}
return s;
}

数位dp

例题模板

//不要连续的62和4
int dp[20][2], num[20];

dfs(int len, int sta, int lim) {
if(len<1) return 1;
if(!lim && dp[len][sta]!=-1) return dp[len][sta];
int ans=0, mi=lim ? num[len] : 9;
for(int i=0; i<=mi; i  ) {
if(i==4 || (sta && i==2)) continue;
ans =dfs(len-1, i==6, lim && i==mi);
}
if(!lim) dp[len][sta]=ans;
return ans;
}

int jj(int x) {
int cnt=0;
while(x) {
num[  cnt]=x;
x/=10;
}
return dfs(cnt, 0, 1);
}

int main() {
memset(dp, -1, sizeof(dp));
int n, m;
while(cin>>n>>m && n && m) {
cout<<jj(m)-jj(n-1)<<endl;
}
return 0;
}

计算基础

$1^2 2^2 3^2 \cdots n^2=\frac{n\times(n 1)\times(2n 1)}{6}$

$1^3 2^3 3^3 \cdots n^3=(1 2 3 \cdots n)^2=(\frac{n\times(n 1)}{2})^2$

求区间$[i, j]$的两两组合乘积的和：

即$[1,3]=1 2 3 6$：

法一：$dp[i][j 1]=dp[i][j] a[j]*(sum[j-1]-sum[i-1])$

$dp[i 1][j]=dp[i][j]-a[i]*(sum[j]-sum[i])$

法二：$2w[i][j]=(sum[j] sum[i-1])^2-(squaresum[j]-squaresum[i-1])$

求$n^k$的最高$m$位：

令$10^w=n^k$，$w$为小数，$x=(ll)w,\ y=w-x$，则$x=100\dots0$，$1$与$n^k$的最高位同位，$10^y=\frac{n^k}{10^x}$，相当于右移$x$位，最后留一位在小数点前。所以，$(ll)10^y$位$n^k$的最高位，$(ll)10^{y m-1}$为$n^k$的最高$m$位。

关键在于把$k$次方在$\log10$中变成乘$k$：w=$\log10(n^k)=k\times\log10(n)$

$10^y$使用$pow(10, y)$即可。

注：c 中$\log(n)$是以$e$为底，而不是$10$，$\log10(n)$才是以十为底。且结果可为小数，此处利用可为小数，关键在于我不用关注如何计算$10$的小数次方。输出变(ll)前一定要加eps，不然会减一。

求$lcm(i,j)=n,\ i\leq j$的$(i,j)$对数：

完全分解：$n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}$。

仅看$p_i$，$(i,j)$组合数有$2*(k_i 1)-1$，最后再加上一次$i==j$的情况，所以$ans=(\prod_{i=1}^s(2\times k_i 1) 1)/2$

求$n!$中$p$的个数：

即$n!$中可有多少个$p$相乘。

完全分解：$n!=\prod_{p\leq n}p^{\sum[\frac n{p^r}]}$。证明见数论笔记本。

对于素因子$p$的指数：$h=[\frac n p] [\frac n{p^2}] \dots=\sum_{r=1}^{\infty}[\frac n {p^r}]$

因为单次计算快，可用于二分答案。

ll check(ll m){
ll ans=0;
while(m){
ans =m/5;
m/=5;
}
return ans;
}

求$ax bz=c$且$x\geq0,b\geq0$的$(x,y)$对数，$\gcd(a,b)=1$：

对固定$a,b$，变动$c$，外层循环枚举$x$，内层循环枚举$y$，给$vis[ax by]$赋值1，即求出$a*b$(即$lcm(a,b)$内的)次数。对于$c>a*b$，在查询$c$时，$cnt=c/(ab) vis [c\�]$。

求$ax (a 1)y (a 2)z==b$，$x,y,z$为非负整数的$(x,y,z)$对数：

只选$a,a 1$等于$b$时$ans $，且$y>=2$时可将$a(x 1)=x (x 2)$，所以贡献为$y/2 1$,同理计算$a 1,a 2$等于$b$时，贡献相同。

其中，将$2(a 1)=a (a 2)$的变化，不会产生重复的$(x,y,z)$，为什么呢，因为相当于按把$x,z$全部尽可能移向$y$分类，此时的$(x,z)$不同，则最终不同。

其中有重复计算$(0,0)$，若存在$(a 1)\mid b$可计算后减去重复计算的$\frac{\frac b{a 1}}2 1$。或者直接把$x,y$的初始值设为$0,1$或$1,0$。

//s[i]存可变动a使ax (a 1)y (a 2)z=i的(x,y,z)对数。
//若要计算特定a，去掉最外层循环即可
void pre2() {
for(ll x=1; x<=n/3-1; x  ) {//枚举a
ll d=n-x*3-3;
for(ll i=0; x*i<=d; i  ) {
for(ll j=0; x*i (x 1)*j<=d; j  ) {
s[x*i (x 1)*j 3*x 3] =j/2 1;
}
}
for(ll i=1; (x 2)*i<=d; i  ) {
for(ll j=0; (x 2)*i (x 1)*j<=d; j  ) {
s[(x 2)*i (x 1)*j 3*x 3] =j/2 1;
}
}
}
}

求有限$\pm2^i$表示数：

有$2^i$若干个，能否用$\sum\pm2^i$表示$n$。直接指数从大到小贪心。

for(int i=25; i>=0; i--){//题目i最大25
while(cnt[i]--){
if(w>=0) w-=a[i];
else w =a[i];
}
}
if(w) cout<<"No\n";
else cout<<"Yes\n";

四边形不等式

若二元函数满足当$a\le b\le c\le d,\ a\le b\le c\le d$,有$w(a,d) w(b,c)\ge w(a,c) w(b,d)w(a,d) w(b,c)\ge w(a,c) w(b,d)$，则称二元函数满足四变形不等式。

判定：若$a<b$，$w(a,b 1) w(a 1,b)\ge w(a,b) w(a 1,b 1)a<b,w(a,b 1) w(a 1,b)\ge w(a,b) w(a 1,b 1)，则二元函数w$满足四变形不等式。

应用：方程$f_i=min_{0\le j<i}\{f_j w(j,i)\}$中$w$满足四边形不等式，则决策点单调递增。

多项式

拉格朗日插值

$f(x)=\sum_{i=1}^ny_i\prod_{i\neq j}\frac{x-x_j}{x_i-x_j},\ O(n^2)$

特殊的，当$x_i$为等差数列时，可以预处理分母部分为阶乘逆元，复杂度降为$O(nlog(n))$

高精度

平平给的模板

typedef long long LL;

struct BigInteger {
static const int BASE = 100000000;
static const int WIDTH = 8;
vector<LL> s;

BigInteger(long long num = 0) {
*this = num;
}
BigInteger operator = (long long num) {
s.clear();
do {
s.push_back(num % BASE);
num /= BASE;
} while(num > 0);
return *this;
}
BigInteger operator = (const string & str) {
s.clear();
int x, len = (str.size() - 1) / WIDTH   1;
for(int i = 0 ; i < len ; i   ) {
int ed = str.size() - i * WIDTH;
int sta = max(0,ed - WIDTH);
sscanf(str.substr(sta,ed - sta).c_str(),"%d",&x);
s.push_back(x);
}
strip();
return *this;
}
bool operator < (const BigInteger & b) const {
if(s.size() != b.s.size()) return s.size() < b.s.size();
for(int i = s.size() - 1 ; i >= 0 ; i-- )
if(s[i] != b.s[i]) return s[i] < b.s[i];
return false;
}
bool operator > (const BigInteger & b) const {
return b < *this;
}
bool operator <= (const BigInteger & b) const {
return !(b < *this);
}
bool operator == (const BigInteger & b) const {
return !(b < *this) && !(*this < b);
}
bool operator != (const BigInteger & b) const {
return !(*this == b);
}
BigInteger operator * (const BigInteger & b) const {
///printf("*1\n");
BigInteger c;
c.s.resize(s.size()   b.s.size()   1,0);
for(int i = 0 ; i < (int)s.size() ; i   )
for(int j = 0 ; j < (int)b.s.size() ; j   )
c.s[i   j]  = s[i] * b.s[j];
for(int i = 0 ; i < (int)s.size()   (int)b.s.size() ; i   ) if(c.s[i] > BASE) {
c.s[i   1]  = c.s[i] / BASE;
c.s[i] %= BASE;
}
c.strip();
///printf("*2\n");
return c;
}
int mult(const BigInteger & d,const BigInteger & b) const {
int L = 0,R = BASE - 1;
while(L < R) {
int m = (L   R) / 2;
if((2 * m) != L   R)   m;
BigInteger p = b * m;
if(p == d) return m;
else if(p < d) L = m;
else R = m - 1;
}
return L;
}
BigInteger operator / (const BigInteger & b) {
///printf("\1\n");
if(b == 0) return -1;
BigInteger c,d;
c.s.resize(s.size());
for(int i = s.size() - 1 ; i >= 0 ; i-- ) {
d.s.insert(d.s.begin(),s[i]);
d.strip();
int cnt = mult(d,b);
d = d - b * cnt;
c.s[i] = cnt;
}
///printf("\2\n");
c.strip();
return c;
}
BigInteger operator - (const BigInteger & b) {
BigInteger c;
c.s.resize(s.size(),0);
for(int i = 0 ; i < (int)c.s.size() ; i   )
c.s[i] = s[i] - (i < (int)b.s.size() ? b.s[i] : 0);
int i = c.s.size() - 1;
for( ; i >= 0 ; i-- ) if(c.s[i]) break;
for(i-- ; i >= 0 ; i-- ) if(c.s[i] < 0) {
int j = i   1,k = 0;
while(j < (int)c.s.size() && !c.s[j]) {
j  ;
k  ;
}
c.s[j]--;
c.s[i]  = BASE;
for( ; k > 0 ; k-- ) c.s[i   k] = BASE - 1;
}
c.strip();
return c;
}
void strip() {//去掉首位0
while(!s.empty() && !s.back()) s.pop_back();
if(s.empty()) s.push_back(0);
return ;
}
BigInteger operator   (const BigInteger & b) const {
BigInteger c;
c.s.clear();
for(int i = 0 ,g = 0 ; ; i   ) {
if(g == 0 && i >= (int)s.size() && i >= (int)b.s.size()) break;
int x = g;
if(i < (int)s.size()) x  = s[i];
if(i < (int)b.s.size()) x  = b.s[i];
c.s.push_back(x % BASE);
g = x / BASE;
}
return c;
}
};
istream& operator >> (istream &in,BigInteger & x) {
string s;
if(!(in >> s)) return in;
x = s;
return in;
}
ostream& operator << (ostream &out ,const BigInteger & x) {
out << x.s.back();
for(int i = x.s.size() - 2 ; i >= 0 ; i-- ) {
char buf[20];
sprintf(buf,"lld",x.s[i]);
for(int j = 0 ; j < (int)strlen(buf) ; j   ) out << buf[j];
}
return out;
}
const int N = 300;
BigInteger f[N],p[N][N],k;
int n,m,m_comple,ans[N];
int vis[N],sub[N],comple[N];
BigInteger frac(int x);
void dfs(BigInteger accum,int sub_i,int taken);
BigInteger P(int x,int y);
int main() {
BigInteger a = BigInteger((long long)1e18);
BigInteger b = BigInteger((long long)1e18);
a=a*b;
cout<<a;
return 0;
}
BigInteger frac(int x) {
if(x <= 1) return 1;
if(f[x] != 0) return f[x];
f[x] = frac(x - 1) * x;
/// printf("f[%d] = ",x);
///cout << f[x] << endl;
return f[x];
}
void dfs(BigInteger accum,int sub_i,int taken) {
if(taken == n) {
printf("%d",ans[0]);
for(int i = 1 ; i < n ; i   ) printf(" %d",ans[i]);
printf("\n");
return ;
}
BigInteger pp;
int flag = 1,i = 0;
while(true) {
if(vis[comple[i]]) {
  i;
continue;
}
if(i >= m_comple) {
ans[taken] = sub[sub_i];
dfs(accum,sub_i   1,taken   1);
return ;
}
if(sub_i < m && flag && sub[sub_i] < comple[i]) {
///printf("11111\n");
flag = 0;
pp = P(n - taken - 1,m - sub_i - 1);
///cout << k << " 1 " << accum << " 1 " << pp << endl;
///printf("11111\n");
if(k <= accum   pp) {
ans[taken] = sub[sub_i];
dfs(accum,sub_i   1,taken   1);
return ;
} else accum = accum   pp;
} else {
///printf("222222\n");
pp = P(n - taken - 1,m - sub_i);
///cout << k << " 2 " << accum << " 2 " << pp << endl;
if(k <= accum   pp) {
ans[taken] = comple[i];
vis[comple[i]] = 1;
dfs(accum,sub_i,taken   1);
return ;
} else {
accum = accum   pp;
  i;
}
}
}
return ;
}
BigInteger P(int x,int y) {
if(x <= 0) x = 1;
if(y <= 0) y = 1;
if(p[x][y] != 0) return p[x][y];
///printf("x = %d;y = %d\n",x,y);
///cout << frac(x) << "--------" << frac(y);
return p[x][y] = frac(x) / frac(y);
}

来源：https://www.icode9.com/content-4-875651.html

bqlp的数学模板笔记

数学

数论

素数

素性探测

Fermat 素性测试

卡迈克尔数

二次探测定理

实现

miller-rabin Pollard_rho(kuangbin模板)

短小互推的

欧拉定理

扩展欧拉定理

费马小定理

裴蜀定理

类欧几里得算法

线性同余方程

中国剩余定理

二次剩余

N次剩余

卢卡斯定理

exLucas定理

BSGS

扩展BSGS

原根

杜教筛

乘法逆元

扩展欧几里得：n

线性递推：1到n

线性递推：1!到n!

线性递推：给定随机n个数

常见数论函数

1）u(n), e(n), I(n)​

2）d(n)​ 约数个数

3）Omega(n)​ 可重复素因子个数

4）omega(n)​ 不可重复素因子个数

5）sigma_lambda(n)​ 约数幂和

6）varphi(n) 互质个数

7）mu(n)​

8）Lambda(n)​

9）lambda(n)​

可乘函数（积性函数）

完全（绝对）可乘函数

Dirichlet 乘积

数论分块

常见卷积

莫比乌斯反演

FT

快速傅里叶变换

快速沃尔什变换

快速幂

快速乘

常见数列

斐波那契数

卡特兰数

Harmonic(调和级数)

常见常数

e​

pi

泰勒公式

组合数学

排列组合

其它数学

Cantor

数位dp

计算基础

求区间\([i, j]\)的两两组合乘积的和：

求\(n^k\)的最高\(m\)位：

求\(lcm(i,j)=n,\ i\leq j\)的\((i,j)\)对数：

求\(n!\)中\(p\)的个数：

求\(ax bz=c\)且\(x\geq0,b\geq0\)的\((x,y)\)对数，\(\gcd(a,b)=1\)：

求\(ax (a 1)y (a 2)z==b\)，\(x,y,z\)为非负整数的\((x,y,z)\)对数：

求有限\(\pm2^i\)表示数：

四边形不等式

多项式

拉格朗日插值

高精度

相关推荐

1）u(n), e(n), I(n)

2）d(n) 约数个数

3）Omega(n) 可重复素因子个数

4）omega(n) 不可重复素因子个数

5）sigma_lambda(n) 约数幂和

7）mu(n)

8）Lambda(n)

9）lambda(n)

e