poj 1845 约数和问题

最新推荐文章于 2019-08-26 10:30:00 发布

zxc106

最新推荐文章于 2019-08-26 10:30:00 发布

阅读量430

点赞数

分类专栏：数论

本文链接：https://blog.csdn.net/u013983192/article/details/37756823

版权

数论专栏收录该内容

11 篇文章 0 订阅

订阅专栏

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

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 20007
using 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;
}