poj 1845--Sumdiv

本题要求A^B的所有约数和对9901取模的值。

1. 求约数：当A=0,1或者B=0时，A^B不为合数，其取模值易知。当A为合数时，根据唯一因子分解定理：A = p1^e1*p2^e2*...*pr^er，其中pi为质数，ei为其指数。易知A^B = p1^(e1*B)*p2^(e2*B)*...*pr^(er*B)。所以我们可以用类似筛选法的方法选出A的所有质数因子。
2. 约数和公式：易知A^B的所有因子均可表示为：p1^k1*p2^k2*...*pr^kr，其中0 <= kr <= er*B。这些因子加起来得到：sum = (1+p1+...+p1^(e1*B))*...*(1+pr+...+pr^(er*B))。
3. 等比公式：易知可以使用等比公式对sum每一个pi的因子求和对9901取模，sum1 = (p1^(e1*B+1)-1)/(p1-1)，我们设：sum1 = k*9901+r，其中k>=0，0 <= r < 9901，根据除法定理k以及模值r都是唯一的。所以将上式变形得到，(p1^(e1*B+1)-1) = k*9901*(p1-1)+r*(p1-1).
4. 对式子两边求模9901*(p1-1)，得到(p1^(e1*B+1)-1)%(9901*(p1-1))=r*(p1-1)，所以我们就可以计算出左边的值再除以(p1-1)就是模值r。但是因为p1*9901最大为10^11，在我们利用反复平方法求左边的模值时，有可能超出64位导致结果错误。
5. 在分解素数因子时，立即对素数因子p对9901做求模运算，这时要注意使用等比公式，当取模以后的p1 = 0时，显然sum1%9901 = 1，当p1 = 1时，sum1%9901 = e1*B+1。

#include<cstdio>
#include<cstring>
using namespace std;
#define maxN 2000

int A,B;
int p[maxN],e[maxN];
int factorNum;

int generateFactor(int x)
{
    int i;
    factorNum = 0;
    memset(e,0,sizeof(e));
    for(i = 2;i*i <= x;i++)
    {
        if((i%2||i == 2)&&(x%i == 0))   //只有素数因子才进入统计因子个数，且通过判定的i肯定是素数
        {
            p[factorNum++] = i%9901;   //首先就对素数因子进行取模
            while(x%i == 0)
            {
                e[factorNum-1]++;
                x /= i;
            }
        }
    }
    if(x != 1)   //有一个大于等于sqrt(x)的因子，此种因子不可能有两个以上
    {
        p[factorNum++] = x%9901;
        e[factorNum-1] = 1;
    }
    return 0;
}

int modular_exponentiation(long long p,int e,long long n)
{
    int b[32];
    long long d = 1;
    int maxDigit = 0;
    while(e)
    {
        b[maxDigit++] = e&1;   //求出e中的每个二进制位
        e >>= 1;
    }
    while(maxDigit--)
    {
        d = (d*d)%n;       //若p不是取模以后的素数因子，n可能超过32位，d*d可能溢出超过64位
        if(b[maxDigit])
            d = (d*p)%n;
    }
    return (d-1)/(p-1);
}

int module()
{
    int ans,tmpAns;
    int i;
    int n;
    ans = 1;
    for(i = 0;i < factorNum;i++)
    {
        if(p[i] == 0)          //素数因子为9901*k，其等比数列对9901取模等于1
            tmpAns = 1;
        else if(p[i] == 1)     //素数因子为9901*k+1,其等比数列对9901取模等于e[i]*B+1
            tmpAns = (e[i]*B+1)%9901;
        else
        {
        n = (p[i]-1)*9901;
        tmpAns = modular_exponentiation(p[i],e[i]*B+1,n);
        }
        ans = (ans*tmpAns)%9901;
    }
    return ans;
}

int main()
{
    int ans;
    while(~scanf("%d%d",&A,&B))
    {
        if(A == 0) ans = 0;
        else if(A == 1 || B == 0) ans = 1;
        else
        {
            generateFactor(A);
            ans = module();
        }
        printf("%d\n",ans);
    }
    return 0;
}