数学模板

dream_qwq · 2022-11-21 13:03:34 · 算法·理论

这是一篇关于*数学模板*的文章这些是其他模板分类的链接	数据结构	字符串	杂项	图论

高精度计算

P1932 A+B A-B A*B A/B A%B Problem


class high{
int l,a[N];
bool f;

int cmp(high x){
    if(l>x.l) return 1;
    if(l<x.l) return -1;
    for(int i=l; i; i--){
        if(a[i]>x.a[i]) return 1;
        if(a[i]<x.a[i]) return -1;
    }
    return 0;
}

void cpy(high x, int y){
    l=x.l+y-1;
    for(int i=1; i<=x.l; i++) a[i+y-1]=x.a[i];
}

public:

high(){
    l=f=0, memset(a,0,sizeof a);
}

high(int x){
    l=f=0, memset(a,0,sizeof a);
    while(x) a[++l]=x%10, x/=10;
}

operator int(){
    int s=0;
    for(int i=l; i; i--) s=s*10+a[i];
    return s;
}

void read(){
    string s;
    getline(cin,s);
    while(s[s.size()-1]<'0' || s[s.size()-1]>'9') s.resize(s.size()-1);

    if(s=="0") return;
    l=s.size();
    for(int i=0; i<l; i++) a[l-i]=s[i]-'0';
}

high operator +(high x){
    x.l=max(l,x.l);
    for(int i=1; i<=x.l; i++) x.a[i]+=a[i], x.a[i+1]+=x.a[i]/10, x.a[i]%=10;
    if(x.a[x.l+1]) x.l++;
    return x;
}

high operator -(high x){
    high y=*this;
    if(y.cmp(x)<0) swap(x,y), y.f=1;

    y.l=max(y.l,x.l);
    for(int i=1; i<=y.l; i++){
        if(y.a[i]<x.a[i]) y.a[i]+=10, y.a[i+1]--;
        y.a[i]-=x.a[i];
    }
    while(y.l && !y.a[y.l]) y.l--;
    return y;
}

high operator *(int x){
    high y=*this;
    y.l=l+log10(x)+1;
    for(int i=1; i<=l; i++) y.a[i]*=x;
    for(int i=1; i<=y.l; i++) y.a[i+1]+=y.a[i]/10, y.a[i]%=10;
    if(!y.a[y.l]) y.l--;
    return y;
}

high operator *(high x){
    high y;
    y.l=l+x.l;
    for(int i=1; i<=l; i++){
        for(int j=1; j<=x.l; j++) y.a[i+j-1]+=a[i]*x.a[j], y.a[i+j]+=y.a[i+j-1]/10, y.a[i+j-1]%=10;
    }
    if(!y.a[y.l]) y.l--;
    return y;
}

pair<high,int> div(int x){
    if(l<(int)log10(x)+1) return {0,(int)(*this)};

    long long s=0;
    high y=*this;
    y.l=l-(int)log10(x);
    for(int i=l; i; i--) s=s*10+a[i], y.a[i]=s/x, s%=x;
    if(!y.a[y.l]) y.l--;
    return {y,s};
}

high operator /(int x){
    return this->div(x).first;
}

high operator %(int x){
    return this->div(x).second;
}

pair<high,high> div(high x){
    if(this->cmp(x)<0) return {(high)0,*this};

    high y=*this, z;
    z.l=l-x.l+1;
    for(int i=z.l; i; i--){
        high t;
        t.cpy(x,i);
        while(y.cmp(t)>=0) z.a[i]++, y=y-t;
    }
    if(!z.a[z.l]) z.l--;
    return {z,y};
}

high operator /(high x){
    return this->div(x).first;
}

high operator %(high x){
    return this->div(x).second;
}

void print(){
    if(!l){
        puts("0");
        return;
    }

    if(f) putchar('-');
    for(int i=l; i; i--) printf("%d",a[i]);
    putchar('\n');
}
};


# 快速幂
- [P1226 【模板】快速幂||取余运算](https://www.luogu.com.cn/problem/P1226)
```cpp
int power(long long a, int b, int p){
    long long s=1;
    for(int i=b; i; i>>=1){
        if(i&1) s=s*a%p;
        a=a*a%p;
    }
    return s;
}

应用

模意义下大整数乘法

快速乘

Laoj P2196 64位整数乘法


long long mul(long long a, long long b, long long p){
return ((unsigned long long)a*b-(unsigned long long)((long double)a/p*b+0.5)*p+p)%p;
}


# 数论
## 最大公约数
- [B3634 最大公约数和最小公倍数](https://www.luogu.com.cn/problem/B3634)
```cpp
int gcd(int x, int y){
    for(int i=x%y; i; i=x%y) x=y, y=i;
    return y;
}

最小公倍数

B3634 最大公约数和最小公倍数


long long lcm(long long x, long long y){
return x/gcd(x,y)*y;
}


## 欧拉函数
- [Laoj P1468 [模板]欧拉函数](http://luanoj.cn:9090/oj/Problem_Show.asp?id=1468)
```cpp
int euler(int x){
    int s=x;
    for(int i=2; i*i<=x; i++){
        if(!(x%i)){
            s=s/i*(i-1);
            while(!(x%i)) x/=i;
        }
    }
    if(x>1) s=s/x*(x-1);
    return s;
}

筛法

素数筛法

P3383 【模板】线性筛素数
```
int p[P];
```

void sieve(int n){ int c=0; bitset<N> u;

for(int i=2; i<=n; i++){
    if(!u[i]) p[++c]=i;

    for(int j=1; j<=c && i*p[j]<=n; j++){
        u[i*p[j]]=1;
        if(!(i%p[j])) break;
    }
}

}


#### 非前缀区间素数筛法
- [P1835 素数密度](https://www.luogu.com.cn/problem/P1835)
```cpp
int p[P];

void sieve(int l, int r){
    int c=0;
    bitset<N> u;

    for(int i=2; (long long)i*i<=r; i++){
        if(!u[i]) p[++c]=i;

        for(int j=1; j<=c && i*p[j]<=sqrt(r); j++){
            u[i*p[j]]=1;
            if(!(i%p[j])) break;
        }
    }

    u.reset();
    for(int i=1; i<=c; i++){
        for(int j=max(2,(l-1)/p[i]+1); j<=r/p[i]; j++) u[j*p[i]-l]=1;
    }
}

筛法求欧拉函数

Laoj P1460 [模板]筛法求欧拉函数
```
int p[P],e[N];
```

void sieve(int n){ int c=0; bitset<N> u;

e[1]=1;
for(int i=2; i<=n; i++){
    if(!u[i]) p[++c]=i, e[i]=i-1;

    for(int j=1; j<=c && i*p[j]<=n; j++){
        u[i*p[j]]=1;
        if(i%p[j]) e[i*p[j]]=e[i]*e[p[j]];
        else{
            e[i*p[j]]=e[i]*p[j];
            break;
        }
    }
}

}


## 乘法逆元
### 扩展欧几里得法
- [P1082 [NOIP2012 提高组] 同余方程](https://www.luogu.com.cn/problem/P1082)
```cpp
long long x,y;

long long exgcd(int a, int b){
    if(!b){
        x=1, y=0;
        return a;
    }

    long long r=exgcd(b,a%b), t=x;
    x=y, y=t-a/b*y;
    return r;
}

int inv(int a, int p){
    return exgcd((a%p+p)%p,p)==1 ? (x%p+p)%p : -1;
}

线性求逆元

P3811 【模板】乘法逆元
```
long long inv[N];
```

int main(){ int n,p; scanf("%d%d",&n,&p);

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

}


## 线性同余方程
- [Laoj P1467 [模板]线性同余方程](http://luanoj.cn:9090/oj/Problem_Show.asp?id=1467)
```cpp
int sol(int a, int b, int p){
    if(!a) return b%p ? -1 : 0;

    int d=exgcd(a,p);
    if(b%d) return -1;

    b/=d, p/=d;
    return (x*b%p+p)%p;
}

中国剩余定理

P1495 【模板】中国剩余定理（CRT）/ 曹冲养猪
```
int a[N],m[N];
```

long long crt(){ long long p=1; for(int i=1; i<=n; i++) p*=m[i];

long long s=0;
for(int i=1; i<=n; i++){
    long long t=p/m[i];
    s=(s+a[i]*t%p*inv(t,m[i])%p)%p;
}
return s;

}


### 扩展：模数不互质的情况
- [P4777 【模板】扩展中国剩余定理（EXCRT）](https://www.luogu.com.cn/problem/P4777)
```cpp
long long a[N],m[N];

long long excrt(){
    for(int i=2; i<=n; i++){
        if(a[1]>a[i]) swap(a[1],a[i]), swap(m[1],m[i]);

        long long d=gcd(m[1],m[i]), p=lcm(m[1],m[i]);
        if((a[i]-a[1])%d) return -1;

        exgcd(m[1]/d,m[i]/d);
        a[1]=(a[1]+mul(mul((a[i]-a[1])/d,(x%p+p)%p,p),m[1],p)%p+p)%p, m[1]=p;
    }
    return a[1];
}

卢卡斯定理

P3807 【模板】卢卡斯定理/Lucas 定理


int com(int n, int m, int p){
return n<m ? 0 : fac[n]%p*inv(fac[m]%p,p)%p*inv(fac[n-m]%p,p)%p;
}

int lucas(int n, int m, int p){ return n<m ? 0 : n ? com(n%p,m%p,p)*lucas(n/p,m/p,p)%p : 1; }


### 扩展卢卡斯定理
- [P4720 【模板】扩展卢卡斯定理/exLucas](https://www.luogu.com.cn/problem/P4720)
```cpp
long long a[L],m[L];

long long fac(long long n, long long p, long long pk){
    if(!n) return 1;

    long long s=1;
    for(int i=1; i<pk; i++){
        if(i%p) s=s*i%pk;
    }
    s=power(s,n/pk,pk);

    for(int i=1; i<=n%pk; i++){
        if(i%p) s=s*i%pk;
    }
    return s*fac(n/p,p,pk)%pk;
}

long long com(long long n, long long m, long long p, long long pk){
    if(n<m) return 0;

    long long s=0, t1=fac(n,p,pk), t2=fac(m,p,pk), t3=fac(n-m,p,pk);
    for(long long i=n; i; i/=p) s+=i/p;
    for(long long i=m; i; i/=p) s-=i/p;
    for(long long i=n-m; i; i/=p) s-=i/p;
    return t1*inv(t2,pk)%pk*inv(t3,pk)%pk*power(p,s,pk)%pk;
}

long long exlucas(long long n, long long m, long long p){
    int c=0;
    for(int i=2; i*i<=p; i++){
        long long t=1;
        while(!(p%i)) p/=i, t*=i;
        if(t>1) a[++c]=com(n,m,i,t), ::m[c]=t;
    }
    if(p>1) a[++c]=com(n,m,p,p), ::m[c]=p;
    return crt(c);
}

离散对数

P2485 [SDOI2011]计算器


int bsgs(int a, int b, int p){
if(!(a%p)) return b%p==1 ? 0 : b%p ? -1 : 1;

int n=ceil(sqrt(p));
long long s=b, t=power(a,n,p);
map<int,int> m;
for(int i=0; i<=n; i++) m[s]=i, s=s*a%p;

s=1;
for(int i=1; i<=n; i++){
    if(m.count(s=s*t%p)) return i*n-m[s];
}
return -1;
}


# 线性代数
## 矩阵
- [B2105 矩阵乘法](https://www.luogu.com.cn/problem/B2105)
- [P3390 【模板】矩阵快速幂](https://www.luogu.com.cn/problem/P3390)
```cpp
struct mat{
    int n,m;
    long long a[N][N];

    mat(){
        n=m=0, memset(a,0,sizeof a);
    }

    mat(int x, int y){
        n=x, m=y, memset(a,0,sizeof a);
    }

    void read(){
        for(int i=1; i<=n; i++){
            for(int j=1; j<=m; j++) scanf("%lld",&a[i][j]);
        }
    }

    mat operator *(mat x){
        mat y={n,x.m};
        for(int i=1; i<=n; i++){
            for(int j=1; j<=x.m; j++){
                for(int k=1; k<=m; k++) y.a[i][j]=(y.a[i][j]+a[i][k]*x.a[k][j])%MOD;
            }
        }
        return y;
    }

    mat pow(long long x){
        mat s={n,n}, a=*this;
        for(int i=1; i<=n; i++) s.a[i][i]=1;

        for(long long i=x; i; i>>=1){
            if(i&1) s=s*a;
            a=a*a;
        }
        return s;
    }

    void print(){
        for(int i=1; i<=n; i++){
            for(int j=1; j<=m; j++) printf("%lld ",a[i][j]);
            putchar('\n');
        }
    }
};

高斯消元

P3389 【模板】高斯消元法
```
const double EPS=1e-8;
```

struct mat{ bool cmp(int p, int x, int y){ if(fabs(fabs(a[x][p])-fabs(a[y][p]))>EPS) return fabs(a[x][p])>fabs(a[y][p]); for(int i=p+1; i<=n; i++){ if(fabs(fabs(a[x][i])-fabs(a[y][i]))>EPS) return fabs(a[x][i])<fabs(a[y][i]); } return 0; }

int gauss(){
    for(int i=1; i<=n; i++){
        int t=i;
        for(int j=i+1; j<=n; j++){
            if(cmp(i,j,t)) t=j;
        }
        if(t!=i) swap(a[t],a[i]);
        if(fabs(a[i][i])<EPS) continue;

        for(int j=m; j>=i; j--) a[i][j]/=a[i][i];
        a[i][i]=1;
        for(int j=i+1; j<=n; j++){
            for(int k=m; k>=i; k--) a[j][k]-=a[i][k]*a[j][i];
        }
    }

    for(int i=n; i; i--){
        if(fabs(a[i][i])<EPS) continue;
        for(int j=1; j<i; j++) a[j][m]-=a[j][i]*a[i][m], a[j][i]=0;
    }

    int f=1;
    for(int i=1; i<=n; i++){
        if(fabs(a[i][i])<EPS){
            if(a[i][m]>EPS) return -1;
            else f=0;
        }
    }
    return f;
}

};