扩展欧几里得，解线性同余方程 逆元 poj1845

定理：对于任意整数a,b存在一堆整数x,y，满足ax+by=gcd(a,b)

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

当d可以整除c时，一般方程ax+by=c的一组特解求法：

　　1.求ax+by=d的特解x0,y0

　　2.ax+by=c的特解为（c/d）x0,（c/d）y0

上述方程的通解：（c/d）x0+k(b/d) ,(c/d)y0-k(a/d)

乘法逆元有自然数倒数的类似性质

乘法逆元：b,m互质，并且b整除a，则存在x，有a/b = a*x(mod m)，即a/b模m的结果和a*x模m的结果是相同的，这个x称为b的模m的乘法逆元，记作b^(-1) (mod m)

　　可得b*b^(-1) = 1（mod m)

　　那么当m是质数时，根据费马小定理，有b^(m-1)=1(mod m),那么b的逆元就是b^(m-2)

　　如果只是保证b,m互质，那么解同余方程b*x=1(mod m)可以求出x

所以当遇到除法取模运算时，可以先求出逆元，转换成乘法取模运算

/*
如果单独是个A，那么就可以分解质因数后用公式求约数个数
那么B个A相乘，其约数个数就是mul{1+p^1+p^2...+p^B*ci} 
结果是比数列求和后再相乘，每项等比数列的结果是 
(pi^(B*ci+1)-1)/(pi-1) mod9901,
1.pi-1不是9901的倍数，（pi-1）^(9901-2)就是逆元 
2.pi-1是9901的倍数，逆元不存在，但是pi mod 9901=1。。。 

先把A分解质因数，再等比数列求和(快速幂+逆元)， 
*/
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long  
#define mod 9901

int m,p[200],c[200];
void divide(int n){
    m=0;
    for(int i=2;i*i<=n;i++)
        if(n%i==0){
            p[++m]=i,c[m]=0;
            while(n%i==0) n/=i,c[m]++; 
        }
    if(n>1) p[++m]=n,c[m]=1;
}
ll pow(ll a,ll b){
    ll res=1;
    while(b){
        if(b&1) res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}

int main(){
    ll a,b,ans=1;
    scanf("%lld%lld",&a,&b);
    divide(a);//分解质因数
    for(int i=1;i<=m;i++){
        if((p[i]-1)%mod==0){
            ans=ans*(b*c[i]+1)%mod;
            continue;
        }
        //求分子和分母逆元 
        ll x=pow(p[i],b*c[i]+1)%mod;
        x=(x-1+mod)%mod;
        ll y=pow(p[i]-1,mod-2)%mod;
        ans=ans*x%mod*y%mod; 
    }
    printf("%lld
",ans);
}

求解同余方程：a*x=b(mod m)等价于a*x-b是m的倍数，等价于a*x+m*y=b，当gcd(a,m)|b时，有解

按照拓展欧几里得算法，可解得特解x=x0*b/gcd(a,m)就是原线性同余方程的一个解

　　通解为所有模m/gcd(a,m)与x同余的整数

求解同余方程：noip2012：a*x=1(mod b)的最小整数解 

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

ll a,b,x,y;
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 z=x; x=y,y=z-y*(a/b);
    return d; 
} 
int main(){
    cin >> a >> b;
    exgcd(a,b,x,y);//x可能是负数 
    cout << (x%b+b)%b<<endl;    
}