POJ 1845 Sumdiv

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 10515		Accepted: 2477

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8. 
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.

15 modulo 9901 is 15 (that should be output). 

Source

Romania OI 2002

求 A^B 的约数和 mod 9901

把A因式分解成 p1^a1*p2^a2*..*pk^ak
然后用约数求和函数

这里关键问题是 有可能出现 gcd(pi-1，9901)!=1
这样的话就不能求逆元了、、、、 

如果没有逆元
那么 对于  1+p1+p1^2...p^n

作以下处理 
n%2==1
 原式=(1+p+p^2+..p^(n/2))(1+p^(n/2+1))
n%2==0
n=n-1;
 原式=(1+p+p^2+..p^(n/2))(1+p^(n/2+1))+p^(n+1);

二分递归求解

我忘记考虑 a是9901的倍数情况了，那样 tp会为负，Wa的我好痛苦、、、

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
using namespace std;
int x,y;
#define Mod 9901
void EGcd(int a,int b)
{
    if(b==0)
    {
        x=1;
        y=0;
        return ;
    }
    EGcd(b,a%b);
    int x0=y;
      y=x-a/b*y;
      x=x0;
}
int gcd(int a,int b)
{
    int r;
    while(r=a%b){a=b;b=r;}
    return b;
}
#define N 10002
bool h[N];
int pm[N>>1];
void Getprime()
{
    int cnt=1;
    int i,j;
    pm[0]=2;
    for(i=2;i<N;i+=2)
     h[i]=1;
    for(i=3;i<N;i+=2)
     if(!h[i])
     {
         pm[cnt++]=i;
         for(j=i*i;j<N;j+=i)
          h[j]=1;
     }
}
int Pw(int a,int b)
{
    int t;
    for(t=1;b>0;b=(b>>1),a=(a%Mod)*(a%Mod)%Mod)
        if(b&1) t=(t*a)%Mod;
    return t;
}
int bn(int a,int n)
{
    if(n==2)
    {
        return (1+a+(a%Mod)*(a%Mod)%Mod)%Mod;
    }
    if(n==1)
    {
        return (1+a)%Mod;
    }
    if(n%2)
    {
        return bn(a,n/2)*(1+Pw(a,n/2+1)%Mod)%Mod;
    }
       n--;
        return bn(a,n/2)*(1+Pw(a,n/2+1)%Mod)%Mod+Pw(a,n+1);
}
int rc[12][2];
int main()
{
    Getprime();
    int a,b;
    int tp;
    while(scanf("%d %d",&a,&b)!=EOF)
    {
        if(a==0){printf("0\n");continue;}
        if(b==0||a==1){printf("1\n");continue;}
        tp=1;
        int m=(int)sqrt(a*1.0);
        int i,cnt=0;
        for(i=0;pm[i]<=m;i++)
         if(a%pm[i]==0)
         {
             rc[cnt][0]=pm[i];
             rc[cnt][1]=0;
             while(a%pm[i]==0){a/=pm[i];rc[cnt][1]++;}
             cnt++;
         }
         if(a>1) { rc[cnt][0]=a;rc[cnt++][1]=1;}

         for(i=0;i<cnt;i++)
         {

             if(gcd(rc[i][0]-1,Mod)==1)
             {
                 tp=tp*(Pw(rc[i][0],rc[i][1]*b+1)-1)%Mod;
                 while(tp<0)tp+=9901;//忘记加这个、、、Wa好久
                 EGcd(rc[i][0]-1,9901);
                 while(x<0)x+=9901;
                 tp=tp*(x%9901)%9901;
             }
             else
             { 
                 tp=tp*bn(rc[i][0],rc[i][1]*b)%9901;
             }
         }
        printf("%d\n",tp);
    }


    return 0;
}