POJ 1845 Sumdiv【唯一分解定理+快速幂+分治法(不用逆元)】

题目传送门： POJ 1845 Sumdiv

思路

题意很简单，求x的y次方这个数的因子和对9901取模的值。

首先我们知道，要求一个数的因子和，先要把它分解质因数，得到它的所有素因子（唯一分解定理如下）

假设 x 素因子分解后是p1^c1+p2^c2+…+pn^cn
则 x^y 素因子分解后是p1^c1*y+p2^c2*y+…+pn^cn*y=(1+p1¹+…+p1^c1*y)*…*(1+pn¹+…+pn^cn*y)

那么现在的关键就在于如何快速求等比数列的和对p取模的值，比如说求 s=(a⁰+a¹+a²+…+aⁿ)%p

尝试用逆元求等比数列和，总是WA。因为逆元使用有限制，如果两个数不互质，就没有逆元，可能出现这种情况。

逆元搞了半天，也不会特判，后来看别人题解，发现可以用分治的方法快速求等比数列的和。

假设我现在要求 s=a⁰+a¹+a²+…+aⁿ

我们可以分奇偶讨论，为了便于理解，先举两个数字比较小的例子：

n为奇数，比如n=3时，序列有4项
s=a⁰+a¹+a²+a³
=(a⁰+a¹)+(a²+a³)【把序列分成两半】
=(a⁰+a¹)+a²*(a⁰+a¹)
=(1+a²)*(a⁰+a¹)

n为偶数，比如n=4时，序列有5项
s=a⁰+a¹+a²+a³+a⁴
=(a⁰+a¹)+(a³+a⁴)+a²【把最中间的数单独拿出后，再把序列分成两半】
=(a⁰+a¹)+a³*(a⁰+a¹)+a²
=(1+a³)*(a⁰+a¹)+a²

那么，找下标的规律，容易得到：
n为奇数的一般情况下，s=(1+a^n/2+1)*(a⁰+a¹+…+a^n/2)
n为偶数的一般情况下，s=(1+a^n/2+1)*(a⁰+a¹+…+a^n/2-1)+a^n/2

其中，1+a^n/2+1 可以用快速幂直接求出，a⁰+a¹+…+a^n/2 或者 a⁰+a¹+…+a^n/2-1 可以继续递归求解，因为它们又是等比数列了。

所以，我们可以得到分治求等比数列和的模板（注意求和之后要取模）

ll sum(ll a,ll n)//分治(二分)求等比数列和  a^0+a^1+a^2+...+a^n
{
    if(n==0)return 1;
    if(n&1)return (1+qpow(a,n/2+1))*sum(a,n/2)%p;//qpow是快速幂
    else return ((1+qpow(a,n/2+1))*sum(a,n/2-1)%p+qpow(a,n/2)%p);
}

这个比用逆元求等比数列和取模好用多了，不需要考虑特判的情况，这个在任何情况下都能完美求出等比数列和取模的值。

AC代码

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int p=9901;
ll x,y;
ll qpow(ll a,ll b)//快速幂
{
    ll s=1;
    while(b)
    {
        if(b&1)s=s*a%p;
        b=b/2;
        a=a*a%p;
    }
    return s;
}
ll sum(ll a,ll n)//分治(二分)求等比数列和  a^0+a^1+a^2+...+a^n
{
    if(n==0)return 1;
    if(n&1)return (1+qpow(a,n/2+1))*sum(a,n/2)%p;
    else return ((1+qpow(a,n/2+1))*sum(a,n/2-1)%p+qpow(a,n/2)%p);
}
ll get_ans()
{
    ll ans=1;
    for(ll i=2;i*i<=x;i++)
    {
        ll c=0;
        if(x%i==0)//i是当前找到的素因子
        {
            while(x%i==0)
                x=x/i,c++;//c是素因子的幂次
        }
        ans=ans*sum(i,c*y)%p;
    }
    if(x>1)ans=ans*sum(x,y)%p;//如果最后x还未被除尽到1，则再乘以(x^0+x^1+...+x^y)
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    cin>>x>>y;
    printf("%lld\n",get_ans());
    return 0;
}