poj 1845 约数和问题

这是一道约数和问题，就是数据特别大的约数和问题。所以就需要求等比数列的时候能快一些。

A的B次幂，首先把A分解成质因子乘积的形式，然后每个质因子的质数乘以B，就变成了典型的求约数和了。

我这里不说什么了，推荐一个博客，讲的很详细：

http://blog.csdn.net/lyy289065406/article/details/6648539

约数和公式：

对于已经分解的整数A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)

有A的所有因子之和为

    S = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)

 

这道题学到的东西就是二分求等比数列和：

用递归二分求等比数列1+pi+pi^2+pi^3+...+pi^n：

（1）若n为奇数,一共有偶数项，则：
      1 + p + p^2 + p^3 +...+ p^n

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

（2）若n为偶数,一共有奇数项,则:
      1 + p + p^2 + p^3 +...+ p^n

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

   上式红色加粗的前半部分恰好就是原式的一半，依然递归求解

 

至于分解质因子肯定都会，我习惯打表，就先求出了根号5千万以内的素数，然后再去分解。

代码如下：

#include<iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>#define MAX 20007using namespace std;int u[MAX],su[MAX];int p[MAX],n[MAX];int len;int num=1;long long A,B,ans=1;void init()  //素数表{    int i,j;    memset(u,1,sizeof(u));    for(i=2;i<MAX;i++)    {        if(u[i]) su[num++]=i;        for(j=1;j<num;j++)        {            if(i*su[j]>=MAX) break;            u[i*su[j]]=0;            if(i%su[j]==0) break;        }    }}void div() //将数A分解质因子，p[i]为第i质因子，n[i]为p[i]个数。{    int i,j=1;    memset(n,0,sizeof(n));    for(i=1;i<num&&A!=1;i++)    {        if(A%su[i]==0)        {            p[j]=su[i];            while(A%su[i]==0)            {                A/=su[i];                n[j]++;            }            j++;        }    }    if(A!=1) {p[j]=A;n[j]++;j++;}    len=j;}long long power(long long p,long long k) //p的k次幂 mod 9901，快速幂法{    int i,j;    long long ret=1;    while(k)    {        if(k&1) ret=(ret*p)%9901;        p=(p*p)%9901;        k>>=1;    }    return ret;}long long work(int p,long long k) //二分递归计算等比数列的和{    if(k==0) return 1;    if(k%2)        return (work(p,k/2)*(1+power(p,k/2+1)))%9901;    else return (work(p,k/2-1)*(1+power(p,k/2+1))+power(p,k/2))%9901;}int main(){    int i;    init();    //scanf("%lld%lld",&A,&B);    while(scanf("%lld%lld",&A,&B)!=EOF)    {        if(!A) {cout<<"1"<<endl;continue;}        div();        ans=1;        for(i=1;i<len;i++)        {            ans*=work(p[i],n[i]*B);            ans%=9901;        }        cout<<ans<<endl;    }    return 0;}