POJ 1845 Sumdiv

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

题意：给出A,B两个整数(0 <= A,B <= 50000000)，输出A^B(A的B次，不是异或）的所有因子和模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).

思路：首先，一个数的因子和是他的所有质因数的最大次方的因此和的乘积。（例如12的因子和是1+2+3+4+6+12=1+2+3+2*2+2*3+2*2*3=（1+2+2*2）*(1+3)=28）。然后，算出所有的质因数和质因数的个数，然后，算出（1+a+a*a+..），我使用等比数列求和公式算的。

等比数列求和公式，我按首项为一来算，1*（1-q^n）/(1-q);

之前每步都要%9901，因为等比数列求和公式有除法，涉及到求逆元运算，用扩展欧几里和定理算出逆元。

扩展欧几里和定理：若c%gcd(a,b)==0,则a*x+b*y=c,xy有解。

求逆元：使a/b%p==a*c%p; 则b*c%p==1。-->b*c=1+k*p;-->b*c-k*p=1;然后用扩展欧几里何定里求出c，c便是逆元。

还有别的方法求逆元，待补充（心情好时来补充。）

用逆元算的话，要考虑除数是9901的倍数的问题，因为除数是a-1,a%9901==1,这种情况直接展开，答案是（1+a+a*a+a*a*a）%9901，共n+1项,答案是(n+1);

代码：

#include<cstdio>
#include<cmath>
#include<cstring>
#include<queue>
#include<stack>
#include<cstdlib>
#include<iomanip>
#include<string>
#include<vector>
#include<map>
#include<string>
#include<iostream>
#include<algorithm>
using namespace std;
#define INF 0x3f3f3f3f
typedef long long ll;
#define Max(a,b) (a>b)?a:b
#define lowbit(x) x&(-x)
ll pow1(ll a,ll b)//二分快速幂
{
    ll ans=1,k=a;
    while(b)
    {
        if(b%2)
            ans*=k,ans%=9901;
        k*=k;
        k%=9901;
        b/=2;
    }
    return ans;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;y=0;
        return a;
    }
    ll r=exgcd(b,a%b,x,y);
    ll t=x;x=y;y=t-a/b*y;
    return r;
}
ll work(ll a,ll b,ll x,ll y)
{
    if(a>1&&(a-1)%9901==0)
    {
        return (b+1)%9901;
    }
    b++;
    b%=9900;//这边我用费马小定理优化了一下，如果a,p互质，则（a^(p-1)）%p=1；
    b=pow1(a,b);
     b=(b-1+9901)%9901;
    a=exgcd(a-1,9901,x,y);//x就是a-1的逆元。
    x%=9901;
    x=(x+9901)%9901;
    b*=x;
    b%=9901;
    return b;
}
int main()
{
    ll a,b,x,y,sum=1,xx=0,n;
    scanf("%lld%lld",&a,&b);
    //求质因数及个数
    n=a;
    for(int i=2;i*i<=a;i++)
    {
        xx=0;
        while(n%i==0)
        {
            n/=i;
            xx++;
        }
        if(xx)
        sum*=work(i,xx*b,x,y);//求（1+i*i+i*i*i+...)(xx项）
        sum%=9901;
    }
    if(n!=1)//避免因数是质数并且该因数平方>a的情况，例如28*31
    {
        sum*=work(n,b,x,y);
        sum%=9901;
    }
    if(a==0)
        sum=0;
    printf("%lld\n",sum);
}

还可以用矩阵快速幂来求（1+a+a*a+..），还有用二分递归来求的，待补充。（心情好时再写。）