POJ 1845 Sumdiv（逆元）

题目链接：Sumdiv

题意：给定两个自然数A，B，定义S为A^B所有的自然因子的和，求出S mod 9901的值。

题解：了解下以下知识点

 

1.整数的唯一分解定理

任意正整数都有且只有唯一的方式写出其质因子的乘积表达式

$A={p_1}^{k_1}*{p_2}^{k_2}*{p_3}^{k_3}*...*{p_n}^{k_n}$

2.整数因数个数

$B=(k_1+1)*(k_2+1)*(k_3+1)...*(k_n+1)$

3.整数因数总和

$S=(1+p_1+p_1^2+p_1^3+...+{p_1}^{k_1})*(1+p_2+p_2^2+p_2^3+...+{p_2}^{k_2})*...(1+p_n+p_n^2+p_n^3+...+{p_n}^{k_n})$

明显地：先分解整数A，A的整数因数总和能表示，A^B只不过是pi变成pi^B。

通过等比数列求和公式：$dfrac{{p_1}^{k_1b+1}-1}{p_1-1}=1+p_1+p_1^2+p_1^3+...+{p_1}^{k_1b}$

费马小定理：

$a^{p-1} ≡1 (mod p)$

两边同除以a

$a^{p-2} ≡inv(a) (mod p)$

$inv(a) = a^{p-2} (mod p)$

快速幂求一下即可：

ll fast_mod(ll x,ll y,ll mod){
    ll res=1;
    while(y){
        if(y&1) res=fast_mul(res,x,mod);
        x=fast_mul(x,x,mod);
        y>>=1;
    }
    return res;
}

直接使用快速幂可能爆long long，结合快速乘防止爆long long：

ll fast_mul(ll x,ll y,ll mod){
    ll res=0;
    while(y){
        if(y&1) res=(res+x)%mod;
        x=(x+x)%mod;
        y>>=1;
    }
    return res;
}

最后一个问题，存在逆元的情况是$ax ≡1(mod p)$

如果a是p的倍数时候左边为0，明显不等于右边。可以使用这个公式

$ans=dfrac{a}{b}modm=amod(mb)/b$

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long ll;
const int N=1e7;
const ll MOD=9901;
ll p[N],k[N];
ll a,b,c=1;

ll fast_mul(ll x,ll y,ll mod){
    ll res=0;
    while(y){
        if(y&1) res=(res+x)%mod;
        x=(x+x)%mod;
        y>>=1;
    }
    return res;
}

ll fast_mod(ll x,ll y,ll mod){
    ll res=1;
    while(y){
        if(y&1) res=fast_mul(res,x,mod);
        x=fast_mul(x,x,mod);
        y>>=1;
    }
    return res;
}

void solve(){
    ll res=1;
    for(ll i=1;i<c;i++){
        ll M=MOD*(p[i]-1);
        res=res*(fast_mod(p[i],k[i]*b+1,M)+M-1)/(p[i]-1)%MOD;
    }
    printf("%lld
",(res+MOD)%MOD);
}

int main(){
    while(scanf("%lld%lld",&a,&b)!=EOF){
        c=1;
        for(ll i=2;i*i<=a;i++){
            if(a%i==0){
                ll cnt=0;
                p[c]=i;
                while(a%i==0) a/=i,cnt++;
                k[c]=cnt;
                c++;
            }
        }
        if(a>1) {p[c]=a;k[c]=1;c++;}
        solve();
    }
    return 0;
}