POJ 1845 (约数和+二分等比数列求和)

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

题目大意：A^B的所有约数和，mod 9901.

解题思路：

①整数唯一分解定理：

一个整数A一定能被分成：A=(P₁^K₁)*(P₂^K₂)*(P₃^K₃).....*(P_n^K_n)的形式。其中P_n为素数。

如2004=(2²)*3*167。

那么2004^x=(2^2x)*(3^x)*(167^x)。

②约数和公式

对于一个已经被分解的整数A=(P₁^K₁)*(P₂^K₂)*(P₃^K₃).....*(P_n^K_n),

有约数和S=(1+P₁²+P₁³+.....P₁^k₁)*.....(1+P_n²+P_n³+.....P_n^k_n)。

(1+P1²+P1³+.....P1^k1)是一个等比数列，化简为(P1^k1+1 -1)/(P₁-1)，由于有除法同余式，很容易想到乘法逆元。

但是这题和HDU 1452不同，对于逆元表达式ax=1 mod n，乘法逆元存在的条件是gcd(a,n)=1,即a,n互质，但是这题的gcd(P₁-1,9901)≠1, 所以不能用乘法逆元求解。

所以有必要对等比数列求和公式改一改：

（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);

这样，在对A质因数分解后，对于每一个质因数，累乘sum(质因数，次数)%mod即可，注意sum计算的时候都要mod防止溢出。

注意一下A的范围，A=0或A=1时无法分解质因数，所以特判结果分别是0和1。

#include "cstdio"
#include "map"
using namespace std;
#define LL long long
#define mod 9901
map<LL,LL> prime_factor(LL n)
{
    map<LL,LL> res;
    for(LL i=2;i*i<=n;i++)
        while(n%i==0) {++res[i];n/=i;}
    if(n!=1) res[n]=1;
    return res;
}
LL pow(LL a,LL n)
{
    LL base=a,ret=1;
    while(n)
    {
        if(n&1) ret=(ret*base)%mod;
        base=(base*base)%mod;
        n>>=1;
    }
    return ret%mod;
}
LL sum(LL p,LL n)
{
    if(n==0) return 1;
    if(n&1)  return ((1+pow(p,(n>>1)+1))*sum(p,n>>1))%mod;
    else     return ((1+pow(p,(n>>1)+1))*sum(p,(n-1)>>1)+pow(p,n>>1))%mod;
}
int main()
{
    //freopen("in.txt","r",stdin);
    LL a,b,res=1;
    scanf("%I64d%I64d",&a,&b);
    if(a==0) {printf("0\n");return 0;}
    map<LL,LL> fac=prime_factor(a);
    for(map<LL,LL>::iterator i=fac.begin();i!=fac.end();i++)
    {
        LL tmp=sum(i->first,i->second*b)%mod;
        res=(tmp*res)%mod;
    }
    printf("%I64d\n",res);
}

POJ 1845 (约数和+二分等比数列求和)

原文：http://www.cnblogs.com/neopenx/p/4094705.html