poj 1845 Sumdiv （同余定理，快速幂取余）

链接：poj 1845

题意：求A^B的所有因子的和对9901取余后的值

如：2^3=8,8的因子有 1,2,4,8,所有和为15，取余后也是15

应用定理主要有三个：

（1）整数的唯一分解定理：

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

A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)   其中pi均为素数

（2）约数和公式：

对于已经分解的整数A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)

有A的所有因子之和为

S = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2)

* (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)

（3）同余定理

(a+b)%m=(a%m+b%m)%m

(a*b)%m=(a%m*b%m)%m

A^B的所有约数之和为：

sum = [1+p1+p1^2+...+p1^(k1*B)] * [1+p2+p2^2+...+p2^(k2*B)] *...

* [1+pn+pn^2+...+pn^(kn*B)]

必须得先求A的素因子及其个数，再求和

但是A,B都到5000万了，直接求和肯定会超时，我们再来看

用递归二分求等比数列1+pi+pi^2+pi^3+...+pi^n：

（1）若n为奇数,一共有偶数项，则：

1 + p + p^2 + p^3 +...+ p^n

= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))

= (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

（2）若n为偶数,一共有奇数项,则:

1 + p + p^2 + p^3 +...+ p^n

= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)

= (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

由以上式子可知，前半部分恰好就是原式的一半，可以递归二分求解

其中p^k暴力求解也会超时，可以用快速幂取余

#include<stdio.h>
#include<string.h>
#include<math.h>
#define M 10000
#define N 9901
int p[M],cnt[M];
__int64 PowMod(__int64 a,__int64 b)  //a^b%N快速幂取余
{
    __int64 ans=1;
    a%=N;
    while(b){
        if(b%2)
            ans=ans*a%N;
        b/=2;
        a=a*a%N;
    }
    return ans;
}
__int64 sum(int p,int n)  //递归二分求和
{
    if(n==0)
        return 1;
    if(n%2)
        return sum(p,n/2)*(1+PowMod(p,n/2+1))%N;  //奇数时
    return (sum(p,n/2-1)*(1+PowMod(p,n/2+1))%N+PowMod(p,n/2))%N;  //偶数时
}
int main()
{
    int a,b,i,k;
    __int64 ans;
    while(scanf("%d%d",&a,&b)!=EOF){
        memset(cnt,0,sizeof(cnt));
        k=1;
        for(i=2;i*i<=a;i++){   //计算a的素因子和个数
            if(a%i==0){
                p[k]=i;
                while(a%i==0){
                    cnt[k]++;
                    a/=i;
                }
                k++;
            }
        }
        if(a!=1){   //记住最后判断一下
            p[k]=a;
            cnt[k++]++;
        }
        ans=1;
        for(i=1;i<k;i++)
            ans=ans*(sum(p[i],cnt[i]*b)%N)%N;
        printf("%I64d\n",ans);
    }
    return 0;
}