pku 1845 Sumdiv 数论

http://poj.org/problem?id=1845

这道题真是纠结死人了。

首先要明白整数分解素因子的唯一性贴个连接

http://ganlu.name/blog/archives/%E8%AE%B0%E4%B8%A4%E4%B8%AA%E8%AF%81%E6%98%8E%E8%B4%B9%E9%A9%AC%E5%AE%9A%E7%90%86%E4%B8%8E%E6%95%B4%E6%95%B0%E5%88%86%E8%A7%A3%E7%B4%A0%E5%9B%A0%E5%AD%90%E7%9A%84%E5%94%AF%E4%B8%80%E6%80%A7

其次还要明白

A可以唯一分解成p1^a1*p2^a2*...*pn^an,
A^B=p1^(a1*B)*p2^(a2*B)*...*pn^(an*B);

约数和公式: 

将A^B 分解成素因数形式：A^B=(p1^k1)*(p2^k2)*(p3^k3)………

那么A^B所有因子之和就是 S=(1+p1+p1^2+p1^3+…..p1^k1)*(1+p2+p2^2+p2^3+…..p2^k2)*(1+p3+…)*…………..这里是为什么呢？因为A^B%1 == 0 A^B%p1 == 0

A^B%p1^2 == 0 依次类推所以1,p1,p1^2......p2,p2^2,......都是A^B的因子

 

计算1+p+p^2+...p^n可以利用二分进行加速

 

当n为奇数时,例如n=5

则1+p+p^2+p^3+p^4+p^5=(1+p+p^2)+p^3*(1+p+p^2)=(1+p+p^2)*(1+p^3)

可以发现1+p+p^2+...p^n=(1+p+p^2+...+p^(k + 1)/2 - 1)*(1 + p^(k + 1)/2)

 

当n为偶数时，例如n=4

则1+p+p^2+p^3+p^4=(1+p)+p^3*(1+p)+p^2=(1+p)*(1+p^3)+p^2

可以发现1+p+p^2+...p^n=(1+p+p^2+...+p^(k/2-1))*(1+p^(k/2+1))+p^(k/2)

 

View Code

#include <cstdio>
#include <cstring>
#include <iostream>
#define maxn 9901
using namespace std;

int fac[maxn],num[maxn];
int len;
//计算p^k,快速幂
__int64 Pow(__int64 p,__int64 k)
{
    __int64 sum = 1,tmp = p;
    while (k)
    {
        if (k%2 == 1) sum = (sum*tmp)%maxn;
        tmp = (tmp*tmp)%maxn;
        k = k/2;
    }
    return sum;
}
//计算1 + p + p^2 + p^3 +.....+p^k，推公式
__int64 getnum(__int64 p,__int64 k)
{
    if (k == 0) return 1;
    if (k%2 == 0)
    return ((getnum(p,k/2 - 1)%maxn)*((1 + Pow(p,k/2 + 1))%maxn) + Pow(p,k/2)%maxn)%maxn;
    else
    return ((getnum(p,(k + 1)/2 - 1)%maxn)*((1 + Pow(p,(k + 1)/2))%maxn))%maxn;
}
int main()
{
    int i;
    __int64 a,b;
    scanf("%I64d%I64d",&a,&b);
    memset(num,0,sizeof(num));
    len = 0;
    
    //分解质因子,注意这里i*i<a的理解，i处理完后如果i*i都大于a了，i加1后a肯定小于(i + 1)^2所以只剩下一个减完后的a了
    for (i = 2; i*i <= a; ++i)
    {
        printf(">>>>%d\n",i);
        if (a%i == 0)
        {
            printf("<><><>%d\n",i);
            fac[len] = i;
            while (a%i == 0)
            {
                num[len]++;
                a /= i;
            }
            len++;
        }
    }
    if (a > 1)
    {
        fac[len] = a;
        num[len++] = 1;
    }
    //计算约数和
    __int64 ans = 1;
    for (i = 0; i < len; ++i)
    {
        __int64 tmp = getnum(fac[i],b*num[i]);
        ans = (ans*tmp)%maxn;
    }
    printf("%I64d\n",ans);
    return 0;
}