题目链接： http://poj.org/problem?id=1845

定义： 满足a*k≡1 (mod p)的k值就是a关于p的乘法逆元。

为什么要有乘法逆元呢？

当我们要求（a/b） mod p的值，且a很大，无法直接求得a/b的值时，我们就要用到乘法逆元。 我们可以通过求b关于p的乘法逆元k，将a乘上k再模p，

即(a*k) mod p。其结果与(a/b) mod p等价。

n=p1^a1*p2^a2*p3^a3****Ps^as,那么

s(n)=[（p1^a1+1 -1)/(p1-1)]*[（p2^a2+1 -1)/(p2-1)]*[（p3^a3+1 -1)/(p3-1)]***[（ps^as+1 -1)/(ps-1)];（因子和）

又因为s(n)%mod等于每一个部分取模，所以可以逐步求解，如求（p1^a1+1  -1)/(p1-1)%mod,在这里就要运用除法取模所以要用到乘法逆元的概念，

即(a/b) %p= ( a *b^(-1)%p) ，又因为(a^b) % p = ((a % p)^b) % p ，

所以（p1^a1+1  -1)/(p1-1)%mod==((（p1%mod)^a1+1 -1)%mod*(p1-1)^-1)%mod;

当然存在逆元的前提是gcd(a,p)==1;

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define N  500010
#define mod 9901
typedef __int64 ll;
using namespace std;
ll a,b,X,Y;
ll ans[N],num[N],top;
ll pow(ll x,ll k)
{
    ll t=1;
    while(k)
    {
        if(k&1) t=((t%mod)*(x%mod))%mod;
        k>>=1;
        x=((x%mod)*(x%mod))%mod;
    }
    return t;
}
void extend(__int64 A,__int64 B,__int64 &x1, __int64 &y1)
{
    if(B==0)
    {
        x1=1;
        y1=0;
        return ;
    }
    extend(B,A%B,x1,y1);
    ll t=x1;
    x1=y1;
    y1=t-(A/B)*y1;
    return ;
}
void solve()
{
    ll sum=1,A,xx;
    for(int i=0; i<top; i++)
    {
        if(ans[i]%mod==0) continue;//关键的两个判断，关系到求逆元。 如果ans[i]%mod=0,那么有等级公式可以看出，原式小于0，所以也只能利用原式求，结果为1
        if(ans[i]%mod==1)//即mod|(ans[i]-1),因为ans[i]>=2,所以ans[i]不可能等于1,这是gcd(ans[i]-1,mod)==mod,不存在逆元，无法利用扩展欧几里得求逆元
        {                                  //这时为（1+ans[i]^1+ans[i]^2+.....+ans[i]^num[i])%mod=(num[i]+1)%mod;
            sum=(sum*(num[i]+1))%mod;
            continue;
        }
        A=pow(ans[i],num[i]+1);
        A=(A-1)%mod;
        extend(ans[i]-1,mod,X,Y);//因为ans[i]为素数，ans[i]-1为偶数，所以ans[i]-1与9901互质
        xx=(X%mod+mod)%mod;
        A=((A%mod)*(xx%mod))%mod;
        sum=(sum*A)%mod;
    }
    printf("%I64d\n",sum);
}
int main()
{
    while(scanf("%I64d%I64d",&a,&b)!=EOF)
    {
        if(a==0)
        {
            printf("0\n");
            continue;
        }
        else if(a==1||b==0)
        {
            printf("1\n");
            continue;
        }
        ll t=a;
        top=0;
        memset(num,0,sizeof(num));
        for(int i=2; i*i<=a; i++)
        {
            if(t%i==0)
            {
                num[top]++;
                ans[top]=i;
                t/=i;
                while(t%i==0)
                {
                    num[top]++;
                    t/=i;
                }
                top++;
            }
        }
        if(t>1)
        {
            num[top]++;
            ans[top++]=t;
        }
        for(int i=0; i<top; i++)
        {
            num[i]*=b;
        }
        solve();
    }
    return 0;
}
  • 1.
  • 2.
  • 3.
  • 4.
  • 5.
  • 6.
  • 7.
  • 8.
  • 9.
  • 10.
  • 11.
  • 12.
  • 13.
  • 14.
  • 15.
  • 16.
  • 17.
  • 18.
  • 19.
  • 20.
  • 21.
  • 22.
  • 23.
  • 24.
  • 25.
  • 26.
  • 27.
  • 28.
  • 29.
  • 30.
  • 31.
  • 32.
  • 33.
  • 34.
  • 35.
  • 36.
  • 37.
  • 38.
  • 39.
  • 40.
  • 41.
  • 42.
  • 43.
  • 44.
  • 45.
  • 46.
  • 47.
  • 48.
  • 49.
  • 50.
  • 51.
  • 52.
  • 53.
  • 54.
  • 55.
  • 56.
  • 57.
  • 58.
  • 59.
  • 60.
  • 61.
  • 62.
  • 63.
  • 64.
  • 65.
  • 66.
  • 67.
  • 68.
  • 69.
  • 70.
  • 71.
  • 72.
  • 73.
  • 74.
  • 75.
  • 76.
  • 77.
  • 78.
  • 79.
  • 80.
  • 81.
  • 82.
  • 83.
  • 84.
  • 85.
  • 86.
  • 87.
  • 88.
  • 89.
  • 90.
  • 91.
  • 92.
  • 93.
  • 94.
  • 95.
  • 96.
  • 97.
  • 98.
  • 99.
  • 100.
  • 101.