POJ 1845 Sumdiv（质因数分解+快速幂+二分法求等比数列的和）

题目大意：求 A ^ B 的所有约数的和 模 9901 的值。

解题思路：

1、对A进行素因子分解得：

     A  =  p1^a1  *  p2^a2  *  p3^a3  * ...  *  pn^an.
     故 A ^ B  =   p1^(a1*B)  *  p2^(a2*B)  * ...  *  pn^(an*B);

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

      sum  =  [ 1 + p1 + p1^2  + ... + p1^(a1*B) ]  *  [  1 + p2 + p2^2 + ... + p2^(a2*B) ]  * ... 

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

3、等比求和公式 (p^(n+1)-1)/(p-1)，由于这里有1/(p-1)，需要求(p-1)关于9901的逆元。

      但是，存在逆元需要 gcd(a,m)  = 1，此题不能保证这一点。所以我们用二分+快速求幂。

4、二分法对等比数列求和：

      求 Sn = 1 + p + p^2 + p^3 + ... + p^n，

      若n == 0 ， Sn = 1;

      若n为奇数，

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

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

            + ... + p^(n/2) )  

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

         = ( 1 + p^(n/2+1）) * Sn/2

     若n为偶数，同理  Sn = ( 1 + p^(n/2+1）) * S(n-1)/2    +  p^n/2   

     按照这样递推下去，二分求 

     求：1 + p + p^2 + ... + p^n    返回上式结果  函数Pow(p,n)返回p^n的值

__int64 Sum(__int64 p,__int64 n)
{
    if(n==0)
        return 1;
    if(n&1)
        return ((1+Pow(p,n/2+1))*Sum(p,n/2));
    else 
        return ((1+Pow(p,n/2+1))*Sum(p,(n-1)/2)+Pow(p,n/2));
}

代码如下：

#include <cstdio>
#include <cstring>
#include <cmath>
#define  ll long long
#define  M  9901
#define  N  300000
bool  is[N];
ll   prm[N];

struct div{
     ll  val;
     ll  cnt;
}p[300];

void  getPrm()
{
      int i,j,k=0;
      int s,e=(int)(sqrt(N+0.0)+1);
      memset(is,1,sizeof(is));
      prm[k++]=2;
      for(i=4;i<N;i+=2) is[i]=0;
      for(i=3;i<e;i+=2)
      {
              if(is[i])
              {
                     prm[k++]=i;
                     for(s=i*2,j=i*i;j<N;j+=s)
                              is[j]=0;
              }
      }
      for(;i<N;i+=2)
             if(is[i])    prm[k++]=i;

}

int  division(ll a,ll b)
{
     int k=0;
     for(int i=0;;i++)
     {
             if(prm[i] > sqrt(a*1.0) + 1)  break;
             if(a%prm[i] == 0)
             {
                      while(a%prm[i] == 0)
                      {
                               a = a / prm[i];
                               p[k].cnt++;
                      }
                      p[k].cnt *= b;
                      p[k++].val = prm[i];
            }
     }
     if(a != 1)
     {
            p[k].val = a;
            p[k++].cnt = b;
     }
     return k;
}

ll   Pow(ll a, ll b)
{
     ll ret = 1;
     while(b)
     {
            if(b & 1) ret = ret * a % M;
            a = a * a % M;
            b /= 2;
     }
     return ret;
}

ll   getSum(ll val, ll cnt)
{
     if(cnt == 0) return 1;
     if(cnt & 1)
            return (1 + Pow(val,cnt/2+1)) * getSum(val,cnt/2) % M;
     else
            return ((1 + Pow(val,cnt/2+1)) * getSum(val,cnt/2-1) + Pow(val,cnt/2)) % M;
}

int  main()
{
     ll a,b;
     getPrm();
     while(~scanf("%I64d %I64d",&a,&b))
     {
             memset(p,0,sizeof(p));
             int k = division(a,b);
             ll  ret=1;
             for(int i = 0; i < k; i++)
             {
                    ret = ret * getSum(p[i].val,p[i].cnt) % M;
             }
             printf("%I64d\n",ret);

     }
     return 0;
}