POJ 1845- Sumdiv

[mathjax]
题目链接：POJ1845

题意:求\(A^B\)的所有约数之和\(\mod 9901\), \(1 \le A,B \le 5*10^7 \).

把A分解质因数,可以得到约数和为:
\((1+p_1+p_1^2+…+p_1^{B*c_1})*(1+p_2+p_2^2+…+p_2^{B*c_2})*…*(1+p_n+p_n^2+…+p_n^{B*c_n})\)

其中每一项都是等比数列求和:\( \frac{p_i^{B*c_i+1}-1}{p_i-1}\),如果分母逆元存在,则用乘以逆元代替除法运算,如果不存在,有\(p_i\mod 9901=1\),则:
\( 1+p_i+p_i^2+…+p_i^{B*c_i}\equiv 1+1+1^2+…+1^{B*c_i}\equiv B*c_i+1 \pmod {9901} \).

code：

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <string>
#include <vector>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;

int c[20],p[20];
int m;
void divide(int a){
    m=0;
    for(int i=2;i<=sqrt(a);++i){
      if(a%i==0){
        p[++m]=i;
        c[m]=0;
        while(a%i==0){
          a/=i;
          c[m]++;
        }
      }
    }
    if(a>1){
      p[++m]=a;
      c[m]=1;
    }
}

int quickpow(int mod,int a,int n){
  int res=1;
  while(n){
    if(n&1) res=(ll)res*a%mod;
    a=(ll)a*a%mod;
    n>>=1;
  }
  return res;
}

int inv(int mod,int x){
  return quickpow(mod,x,mod-2);
}

int main(){
  int a,b;
  const int mod=9901;
  cin>>a>>b;
  divide(a);
  int ans=1;
  for(int i=1;i<=m;++i){
    int x=(quickpow(mod,p[i],(ll)b*c[i]+1)-1+mod)%mod;
    int y=(p[i]-1)%mod;
    if(y){
      ans=(ll)ans*(x*inv(mod,y))%mod;
    }
    else{
      ans=ans*((ll)b*c[i]+1)%mod;  
    }
  }
  cout<<ans<<endl;
  return 0;
}