逆元

一、扩展欧几里得求逆元

　　逆元的定义即是对于a，满足a*x≡1（mod b），就把x称作a在模b意义下的逆元，而扩展欧几里得可以就线性同余方程，这里显然适用。

二、快速幂求逆元

　　这里用到了费马小定理和欧拉定理。

　　费马小定理：对于素数p，有a^p-1≡1（mod p）

　　证明：首先根据同余的性质，我们可知若a≡b（mod m），则a*c≡b*c（mod m）

　　　　∴我们对于素数p的剩余系：1,2,3,4...p-1中的两个数x，y

　　　　　有a*x≡a*y（mod p）必定不成立，因为p为素数，所以可以约去a，得到x=y，这显然不成立

　　　　∴1,2,3,4...p-1和a,2a,3a,4a...（p-1）a形成映射

　　　　∴（p-1）！≡（p-1）！a^p-1（mod p）

　　　　∴约去（p-1）！，得a^p-1≡1（mod p）

　　欧拉定理：对于任意互质正整数a，n，满足a^Φ（n）≡1（mod n）

　　证明：与费马小定理类似，只是将剩余系改为小于n与n互质的数。

三、线性递推（递归）

　　对于i在模p意义下的逆元，我们可以用假设p=k*i+r（r<i），因此显然：k*i+r≡0（mod p）

　　而将等式两边同时乘上i^-1和r^-1得到：k*r^-1+i^-1≡0（mod p），移项后为i^-1≡-k*r^-1（mod p）

　　因此i^-1≡-（p div i）*（p mod i）^-1（mod p）

　　所以递推式为inv[i]=-（p/i）*inv[p%i]

四、阶乘逆元

　　显然对于阶乘n！的逆元为，所以阶乘逆元的递推公式即为inv[i]=inv[i+1]*（i+1）

例1:Sumdiv

　　求出a^b的所有约数的和模9901的值。

　　这道题实际上就是让我们求（1+p₁+p₁²+...+p₁^r₁）（1+p₂+p₂²+...+p₂^r₁）...（1+p_n+p_n²+...+p_n^r_n）

　　而对于每个素数的分解我们可以直接通过枚举因数得到，但我们还需要快速求每个因数的和。

　　方法一：对于每个求和，我们根据r的奇偶性分类，可以不断递归求解，递归式如下：

　　　　　　①若n为奇数，cal（p，n）=cal（p，n/2）×（1+p^n/2+1）

　　　　　　②若n为偶数，cal（p，n）=cal（p，n/2-1）×（1+p^n/2+1）+p^n/2

　　方法二：直接运用等比数列求和公式

　　代码

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long ll;
const ll N=1000000; 
const ll mod=9901;

ll prime[N],r[N]; 
ll qpow(ll a,ll b)
{
    ll res=1;
    for(;b;b>>=1)
    {
        if(b&1)res=res*a%mod;
        a=a*a%mod;
    }
    return res;
}
ll cal(ll p,ll n)
{ 
    if(!n)return 1;
    if(n&1)
        return cal(p,n>>1)*(1+qpow(p,(n>>1)+1))%mod;
    else return (cal(p,(n>>1)-1)*(1+qpow(p,(n>>1)+1))%mod+qpow(p,n>>1))%mod;
}

int main()
{
    ll a,b,cnt=0;
    scanf("%lld%lld",&a,&b);
    for(ll i=2;i<=sqrt(a);i++)
        if(!(a%i))
        {
            prime[++cnt]=i;
            while(!(a%i))r[cnt]++,a/=i;
        }
    if(a!=1)prime[++cnt]=a,r[cnt]=1;
    ll ans=1;
    for(ll i=1;i<=cnt;i++)
        ans=ans*cal(prime[i],r[i]*b)%mod;
    printf("%lld",ans);
}