poj1845 逆元，快速模幂

题目大意：

给定两个正整数和，求的所有因子和对9901取余后的值。

分析：

很容易知道，先把分解得到，那么得到，那么

     的所有因子和的表达式如下

 

    

因为要取模且存在除法，所以要用到逆元。

对于正整数和，如果有，那么把这个同余方程中的最小正整数解叫做模的逆元。

 

逆元一般用扩展欧几里得算法来求得，如果为素数，那么还可以根据费马小定理得到逆元为。

 

推导过程如下

                            

 

求现在来看一个逆元最常见问题，求如下表达式的值（已知）

 

           

 

当然这个经典的问题有很多方法，最常见的就是扩展欧几里得，如果是素数，还可以用费马小定理。

 

但是你会发现费马小定理和扩展欧几里得算法求逆元是有局限性的，它们都会要求与互素。实际上我们还有一

种通用的求逆元方法，适合所有情况。公式如下

 

          

 

现在我们来证明它，已知，证明步骤如下

 

          

等比数列求和公式，用如下公式即可

 

                     

 

因为可能会很大，超过int范围，所以在快速幂时要二分乘法。

接下来就是代码了：（一定要讲乘法二分，还有取模后一定要加模）

#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <queue>
#include <stack>
#include <vector>
#include <map>
#include <set>
#include <iomanip>

using namespace std;
#pragma comment(linker, "/STACK:102400000,102400000")
#define maxn 10005
#define MOD 1000000007
#define mem(a , b) memset(a , b , sizeof(a))
#define LL long long
#define ULL long long
const long long INF=0x3fffffff;
bool prime[maxn];
int p[maxn];

void is_prime()
{
    mem(prime , true);
    int id = 0;
    prime[0] = prime[1] = false;
    for(int i = 2 ; i < maxn ; i ++)
    {
        if(prime[i])
        {
            p[id++] = i;
            int j = 2 * i;
            while(j <= maxn) prime[j] = false , j += i;
        }
    }
}

LL mulit(LL a , LL  b , LL mod)
{
    a %= mod;
    b %= mod;
    LL res = 0;
    while(b)
    {
        if(b & 1)
        {
            res = (res + a) % mod;
        }
        b >>= 1;
        a = (a + a) % mod;
    }
    return res % mod;
}

LL quick_mod(LL a , LL b , LL mod)
{
    LL res = 1;
    while(b)
    {
        if(b & 1)
        {
            res = mulit(res , a , mod);
        }
        a = mulit(a , a , mod);
        b >>= 1;
    }
    return res;
}
//ofstream  ofile;
int main()
{
    LL a , b;
    is_prime();
    while(scanf("%lld %lld",&a , &b) != EOF )
    {
        LL ans = 1;
        for(int i = 0 ; p[i] * p[i] <= a ; i++)
        {
            if(a % p[i] == 0)
            {
                int num = 0;
                while(a % p[i] == 0)
                {
                    num ++;
                    a /= p[i];
                }
                ans *= (quick_mod(p[i] , num * b+1 , 9901 *(p[i]-1) ) + (9901 *(p[i]-1))- 1) / (p[i] - 1);
                ans %= 9901;
            }
        }
        if(a > 1)
        {
            ans *= (quick_mod(a ,  b+1 , 9901 *(a-1) ) + 9901 *(a-1) - 1) / (a - 1);
            ans %= 9901;
        }
        cout << ans << endl;
    }
    return 0;
}