*模板--数学

 

加减乘除

//compare比较函数：相等返回0，大于返回1，小于返回-1
int compare(string str1,string str2)
{
    if(str1.length()>str2.length()) return 1;
    else if(str1.length()<str2.length())  return -1;
    else return str1.compare(str2);
}
//高精度加法
//只能是两个正数相加
string add(string str1,string str2)//高精度加法
{
    string str;

    int len1=str1.length();
    int len2=str2.length();
    //前面补0，弄成长度相同
    if(len1<len2)
    {
        for(int i=1;i<=len2-len1;i++)
           str1="0"+str1;
    }
    else
    {
        for(int i=1;i<=len1-len2;i++)
           str2="0"+str2;
    }
    len1=str1.length();
    int cf=0;
    int temp;
    for(int i=len1-1;i>=0;i--)
    {
        temp=str1[i]-'0'+str2[i]-'0'+cf;
        cf=temp/10;
        temp%=10;
        str=char(temp+'0')+str;
    }
    if(cf!=0)  str=char(cf+'0')+str;
    return str;
}
//高精度减法
//只能是两个正数相减，而且要大减小
string sub(string str1,string str2)//高精度减法
{
    string str;
    int tmp=str1.length()-str2.length();
    int cf=0;
    for(int i=str2.length()-1;i>=0;i--)
    {
        if(str1[tmp+i]<str2[i]+cf)
        {
            str=char(str1[tmp+i]-str2[i]-cf+'0'+10)+str;
            cf=1;
        }
        else
        {
            str=char(str1[tmp+i]-str2[i]-cf+'0')+str;
            cf=0;
        }
    }
    for(int i=tmp-1;i>=0;i--)
    {
        if(str1[i]-cf>='0')
        {
            str=char(str1[i]-cf)+str;
            cf=0;
        }
        else
        {
            str=char(str1[i]-cf+10)+str;
            cf=1;
        }
    }
    str.erase(0,str.find_first_not_of('0'));//去除结果中多余的前导0
    return str;
}
//高精度乘法
//只能是两个正数相乘
string mul(string str1,string str2)
{
    string str;
    int len1=str1.length();
    int len2=str2.length();
    string tempstr;
    for(int i=len2-1;i>=0;i--)
    {
        tempstr="";
        int temp=str2[i]-'0';
        int t=0;
        int cf=0;
        if(temp!=0)
        {
            for(int j=1;j<=len2-1-i;j++)
              tempstr+="0";
            for(int j=len1-1;j>=0;j--)
            {
                t=(temp*(str1[j]-'0')+cf)%10;
                cf=(temp*(str1[j]-'0')+cf)/10;
                tempstr=char(t+'0')+tempstr;
            }
            if(cf!=0) tempstr=char(cf+'0')+tempstr;
        }
        str=add(str,tempstr);
    }
    str.erase(0,str.find_first_not_of('0'));
    return str;
}

//高精度除法
//两个正数相除，商为quotient,余数为residue
//需要高精度减法和乘法
void div(string str1,string str2,string &quotient,string &residue)
{
    quotient=residue="";//清空
    if(str2=="0")//判断除数是否为0
    {
        quotient=residue="ERROR";
        return;
    }
    if(str1=="0")//判断被除数是否为0
    {
        quotient=residue="0";
        return;
    }
    int res=compare(str1,str2);
    if(res<0)
    {
        quotient="0";
        residue=str1;
        return;
    }
    else if(res==0)
    {
        quotient="1";
        residue="0";
        return;
    }
    else
    {
        int len1=str1.length();
        int len2=str2.length();
        string tempstr;
        tempstr.append(str1,0,len2-1);
        for(int i=len2-1;i<len1;i++)
        {
            tempstr=tempstr+str1[i];
            tempstr.erase(0,tempstr.find_first_not_of('0'));
            if(tempstr.empty())
              tempstr="0";
            for(char ch='9';ch>='0';ch--)//试商
            {
                string str,tmp;
                str=str+ch;
                tmp=mul(str2,str);
                if(compare(tmp,tempstr)<=0)//试商成功
                {
                    quotient=quotient+ch;
                    tempstr=sub(tempstr,tmp);
                    break;
                }
            }
        }
        residue=tempstr;
    }
    quotient.erase(0,quotient.find_first_not_of('0'));
    if(quotient.empty()) quotient="0";
}

---

 扩展欧几里得

// ax+by=gcd(a,b)
// d为a和b的最大公因数
void exgcd(LL a,LL b,LL &d,LL &x,LL &y)  
{
    if(!b)
    {
        d=a;x=1;y=0;
    }
    else {
        exgcd(b,a%b,d,y,x);
        y-=x*(a/b);
    }
}

// ax+by=gcd(a,b)
LL exgcd(LL a,LL &x, LL b, LL &y)  
{
    LL x1=0,y1=1,x0=1,y0=0;
    LL r=(a%b+b)%b;
    LL q=(a-r)/b;
    x=0;y=1;
    while(r)
    {
        x=x0-q*x1;
        y=y0-q*y1;
        x0=x1;y0=y1;
        x1=x;y1=y;
        a=b;b=r;
        r=a%b;q=(a-r)/b;
    }
    return b;  //返回a和b的最大公因数
}

如果得到ax ≡ c (mod b)的某一特解X，那么令r = b/gcd(a, b)，可知x在[0, r-1]上有唯一解

所以用x = (X % r + r) % r就可以求出最小非负整数解x

---

求解n个方程x%m[i]=b[i]，返回x

void exgcd(int a,int b,int &d,int &x,int &y){
    if(!b){d=a;x=1;y=0;}
    else {exgcd(b,a%b,d,y,x);y-=x*(a/b);}
}
//  求解n个方程x%m[i]=b[i]，返回x
int solve(int* b,int* m,int n){
    int m1=m[0],b1=b[0];;
    for(int i=1;i<n;i++){
        int m2=m[i],b2=b[i];
        int g,x,y;
        exgcd(m1,m2,g,x,y);
        int c=b2-b1;
        if(c%g) return -1;  //无解
        x=x*(c/g);
        int r=m2/g;
        x=(x%r+r)%r;
        b1+=m1*x;
        m1*=m2/g;
        b1%=m1;
    }
    return (b1%m1+m1-1)%m1+1;  //返回正整数解
}

---

中国剩余定理

求解x%m[i]=b[i]，其中m[i]两两互质

//  求解x%m[i]=b[i]，其中m[i]两两互质
int china(int* b,int* m,int n){
    int M=1,ans=0;
    for(int i=0;i<n;i++) M*=m[i];
    for(int i=0;i<n;i++){
        int N=M/m[i];
        int g,x,y;
        exgcd(N,m[i],g,x,y);
        ans=(ans+b[i]*N*x)%M;
    }
    return (ans+M)%M;
}

---

BSGS算法

求解A^x≡B(mod C)，C是素数。

#define ll long long
const int mod=99991;
ll head[mod],nex[mod];
ll cnt;
ll hs[mod],id[mod];
void init(){
    memset(head,-1,sizeof(head));
    cnt=0;
}
void insert_(ll x,ll i){
    hs[cnt]=x;id[cnt]=i;
    ll u=x%mod;
    nex[cnt]=head[u];
    head[u]=cnt++;
}
ll find_(ll x){
    int u=x%mod;
    for(int i=head[u];~i;i=nex[i]){
        if(hs[i]==x) return id[i];
    }
    return -1;
}
//求解A^x≡B(mod C),返回x，若无解返回-1
ll bsgs(ll a,ll b,ll c){
    if(b==1) return 0;
    ll m=ceil(sqrt(c*1.0));
    ll q=1;
    for(ll i=0;i<=m;i++){
        if(i==0) insert_(b%c,i);
        else{
            q=q*a%c;
            insert_(q*b%c,i);
        }
    }
    ll temp=1;
    for(ll i=1;i<=m;i++){
        temp=temp*q%c;
        ll j=find_(temp);
        if(j!=-1) return i*m-j;
    }
    return -1;
}

扩展BSGS算法

求解A^x≡B(mod C)。

#define ll long long
const int mod=99991;
ll head[mod],nex[mod];
ll cnt;
ll hs[mod],id[mod];
void init(){
    memset(head,-1,sizeof(head));
    cnt=0;
}
void insert_(ll x,ll i){
    hs[cnt]=x;id[cnt]=i;
    ll u=x%mod;
    nex[cnt]=head[u];
    head[u]=cnt++;
}
ll find_(ll x){
    int u=x%mod;
    for(int i=head[u];~i;i=nex[i]){
        if(hs[i]==x) return id[i];
    }
    return -1;
}
ll gcd(ll a,ll b){
    ll r=a%b;
    while(r){a=b;b=r;r=a%b;}
    return b;
}
//求解a^x≡b(mod c)，返回x，无解返回-1 
ll BSGS(ll a,ll b,ll c){
    a%=c;b%=c;
    if(b==1) return 0;
    ll g=1,d=1,k=0;
    while((g=gcd(a,c))!=1){
        if(b%g) return -1;
        b/=g;c/=g;
        k++;
        d=d*(a/g)%c;
        if(b==d) return k;
    }
    ll m=ceil(sqrt(c*1.0)),q=1;
    for(ll i=0;i<=m;i++){
        if(i==0) insert_(b,i);
        else {
            q=q*a%c;
            insert_(q*b%c,i);
        }
    }
    for(ll i=1;i<=m;i++){
        d=d*q%c;
        ll j=find_(d);
        if(j!=-1) return i*m-j+k;
    }
    return -1;
}

---

欧拉函数（求与n互质的数有多少个）

//返回n的欧拉函数值
int phi(int n){
    int m=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到n的欧拉函数值
void get_phi(int n,int  *phi){
    for(int i=2;i<=n;i++) phi[i]=0;
    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);
        }
    }
}

---

高斯消元

const int inf=0x3f3f3f3f;
//const int mod=2;
const int maxn=25;
int a[maxn][maxn];  //系数
int x[maxn],free_x[maxn];
int n,m;
/*
int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
}
int lcm(int a,int b){
    return a/gcd(a,b)*b;
}
void exgcd(int a,int b,int &d,int &x,int &y){
    if(!b){d=a;x=1;y=0;}
    else {exgcd(b,a%b,d,y,x); y-=x*(a/b);}
}
*/

//n个方程m个未知量
int gauss(int n,int m){
    int r,c;
    int num=0;  //自由变元
    for(r=0,c=0;r<n&&c<m;c++){
        int max_r=r;
        for(int i=r+1;i<n;i++) if(abs(a[i][c]) > abs(a[max_r][c])) max_r=i;
        if(max_r!=r) for(int j=c;j<=m;j++) swap(a[r][j],a[max_r][j]);
        if(!a[r][c]) {free_x[num++]=c;continue;}  //
        for(int i=r+1;i<n;i++) if(a[i][c]){
            /*
             int d=lcm(abs(a[i][c]),abs(a[r][c]));
            int t1=d/a[i][c],t2=d/a[r][c];
            for(int j=c;j<=m;j++) a[i][j]=((a[i][j]*t1-a[r][j]*t2)%mod+mod)%mod;
            */
            for(int j=c;j<=m;j++) a[i][j]^=a[r][j];
        }
        r++;
    }
    for(int i=r;i<n;i++) if(a[i][m]) return -1;  //无解
    return num; //返回自由变元数
    //求解
    /*
    for(int i=r-1;i>=0;i--){
        x[i]=a[i][m];
        for(int j=i+1;j<m;j++){
            x[i]=((x[i]-a[i][j]*x[j])%mod+mod)%mod;
        }
        int x1,y1,d;
        exgcd(a[i][i],mod,d,x1,y1);
        x1=((x1%mod)+mod)%mod;
        x[i]=x[i]*x1%mod;
    }
    */
    
    //枚举自由变元（开关问题只能取0或1）
    /*
    int sta=1<<(m-r);
    int res=inf;
    for(int i=0;i<sta;i++){
        int cnt=0;
        int id=i;
        for(int j=0;j<num;j++){
            x[free_x[j]]=(id&1);
            if(x[free_x[j]]) cnt++;
            id>>=1;
        }
        for(int j=r-1;j>=0;j--){
            int temp=a[j][m];
            for(int k=j+i;k<m;k++){
                if(a[j][k]) temp^=x[k];
            }
            x[j]=temp;
            if(x[j]) cnt++;
        }
        if(cnt<res) res=cnt;
    }
    return res; 
    */
}

const int maxn=75;
const int mod=131;  //模数
int a[maxn][maxn];
int x[maxn];

int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
}
int lcm(int a,int b){
    return a/gcd(a,b)*b;
}
void exgcd(int a,int b,int &d,int &x,int &y){
    if(!b){d=a;x=1;y=0;}
    else {exgcd(b,a%b,d,y,x); y-=x*(a/b);}
}

int gauss(int n,int m){
    int r,c;
    for(r=0,c=0;r<n&&c<m;c++){
        int max_r=r;
        for(int i=r+1;i<n;i++) if(abs(a[i][c]) > abs(a[max_r][c])) max_r=i;
        if(max_r!=r) for(int j=c;j<=m;j++) swap(a[r][j],a[max_r][j]);
        if(!a[r][c]) continue;
        for(int i=r+1;i<n;i++) if(a[i][c]){
            int d=lcm(abs(a[i][c]),abs(a[r][c]));
            int t1=d/a[i][c],t2=d/a[r][c];
          //  if(a[i][c]*a[r][c]<0) t2=-t2;
            for(int j=c;j<=m;j++) a[i][j]=((a[i][j]*t1-a[r][j]*t2)%mod+mod)%mod;
        }
        r++;
    }
    for(int i=r;i<n;i++) if(a[i][m]) return -1; //无解
    for(int i=r-1;i>=0;i--){
        x[i]=a[i][m];
        for(int j=i+1;j<m;j++){
            x[i]=((x[i]-a[i][j]*x[j])%mod+mod)%mod;
        }
        int x1,y1,d;
        exgcd(a[i][i],mod,d,x1,y1);
        x1=((x1%mod)+mod)%mod;
        x[i]=x[i]*x1%mod;
    }
    if(r<m) return m-r;  //自由变元
    return 1;//唯一解
}

const int inf=0x3f3f3f3f;
const int maxn=25;
int a[maxn][maxn];
int x[maxn],free_x[maxn];
int n,m;

int gauss(int n,int m){
    int r,c;
    int num=0;
    for(r=0,c=0;r<n&&c<m;c++){
        int max_r=r;
        for(int i=r+1;i<n;i++) if(abs(a[i][c]) > abs(a[max_r][c])) max_r=i;
        if(max_r!=r) for(int j=c;j<=m;j++) swap(a[r][j],a[max_r][j]);
        if(!a[r][c]) {free_x[num++]=c;continue;}  //
        for(int i=r+1;i<n;i++) if(a[i][c]){
            for(int j=c;j<=m;j++) a[i][j]^=a[r][j];
        }
        r++;
    }
    for(int i=r;i<n;i++) if(a[i][m]) return -1;
    int sta=1<<(m-r);  //方案数
    int res=inf;
    //枚举自由变元
    for(int i=0;i<sta;i++){
        int cnt=0;
        int id=i;
        for(int j=0;j<m-r;j++){
            x[free_x[j]]=(id&1);
            if(x[free_x[j]]) cnt++;
            id>>=1;
        }
        for(int j=r-1;j>=0;j--){
            int temp=a[j][m];
            for(int k=j+1;k<m;k++){
                if(a[j][k]) temp^=x[k];
            }
            x[j]=temp;
            if(x[j]) cnt++;
        }
        if(cnt<res) res=cnt;
    }
    return res;
}

线性基

pow(1, 62) ==  4 611 686 018 427 387 904  

struct LiBase{
    ll a[63],p[63];
    int cnt;  // 重建后，该线性基的范围是[1,1<<(cnt)-1]
    //初始化
    void init(){
        cnt=0;
        memset(a,0,sizeof(a));
        memset(p,0,sizeof(p));
    }
    //插入
    bool insert_(ll x){
        for(int i=62;i>=0&&x;i--)if(x&(1LL<<i)){
            if(!a[i]){
                a[i]=x;
                break;
            }else x^=a[i];
        }
        return x>0;
    }
    //最大
    ll query_max(){
        ll ans=0;
        for(int i=62;i>=0;i--){
            if((ans^a[i])>ans) ans^=a[i];
        }
        return ans;
    }
    //最小
    ll query_min(){
        for(int i=0;i<=62;i++) if(a[i]) return a[i];
        return 0;
    }
    //重建
    void rebuid(){
        for(int i=62;i>=0;i--)
            for(int j=i-1;j>=0;j--) if(a[i]&(1LL<<j)) a[i]^=a[j];
        for(int i=0;i<=62;i++) if(a[i]) p[cnt++]=a[i];
    }
    //第k小
    ll kth(ll k){
        ll ans=0;
        if(k>=(1LL<<cnt)) return -1;
        for(int i=cnt-1;i>=0;i--){
            if(k&(1LL<<i)) ans^=p[i];
        }
        return ans;
    }
};
/*
合并
LiBase merge_(LiBase &x,LiBase &y){
    for(int i=60;i>=0;i--)if(y.a[i]) x.insert_(y.a[i]);
    return x;
}
*/

康拓展开（逆）

int  fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320}; //i的阶乘为fac[i]  
// 康托展开-> 表示数字a是 a的全排列中从小到大排，排第几  
// n表示1~n个数  a数组表示数字。  
int kangtuo(int n, char a[]){  
    int sum = 0;  
    for(int i = 0; i < n; i++){  
        int t = 0;  
        for(int j = i + 1; j < n; j++)  
            if(a[i] > a[j])  t++;  
        sum += t * fac[n - i - 1];  
    }
    return sum+1;  
}

//康托展开的逆运算,{1...n}的全排列，中的第k个数为s[]
void reverse_kangtuo(int n, int k, char s[]){
    int vis[8]={0};
    s[k] = 0;
    --k;
    for(int i = 0; i < n; i++){
        int t = k / fac[n - i - 1];
        int j;
        for (j = 1; j <= n; j++){
            if (!vis[j]){
                if (t == 0) break;
                --t;
            }
        }
        s[i] = '0' + j;
        vis[j] = 1;
        k %= fac[n - i - 1];
    }
}