POJ1845 Sumdiv

题意

求 $A^{B}$ 的约数和。对 $9901$ 取模。

其中 $A,B\le 5\times 10^{7}$ 。

题解

将 $A$ 分解：

$$A=\prod_{i=1}^n p_i^{a_i}$$

那么约数和为：

$$\prod_{i=1}^n \sum_{j=0}^{a_i\times B} {p_i}^j $$

所以可以先预处理出 $\le \sqrt{A}$ 的质数，然后对 $A$ 进行分解。

剩下的东西就是对很多等比数列求和，然后相乘了。求和直接用公式：

$$S_n=\dfrac{a_1-a_n\times q}{1-q}$$

很显然需要对除法取模，可以用取模的某个性质。我也不知道这个叫什么，反正是有这个性质。

对于

$$\dfrac{x}{y} \ mod \ p$$

可以化为

$$\dfrac{x \ mod \ (y\times p)}{y} \tag{1}$$

然后就可以求了。快速幂中的模数也需要改成 $y\times p$ 这种形式。

我就很愉快的这么做了，然后很愉快的得了 $95$ 分。

然后发现当 $A$ 是质数时，$1-q$ 就会很大，这就意味着转化为 $(1)$ 式中 $y\times p$ 会爆 int ，然后过程中就会爆 long long，所以我就用 __int128 水过了。

#include<bits/stdc++.h>
#define ll __int128
#define ha 9901

using namespace std;

inline ll read()
{
    char ch=getchar();
    ll f=1,x=0;
    while (ch<'0' || ch>'9')
    {
        if (ch=='-') f=-1;
        ch=getchar();
    }
    while (ch>='0' && ch<='9')
    {
        x=x*10+ch-'0';
        ch=getchar();
    }
    return f*x;
}

bool ss[7100];
ll zs[2005],cnt,A,B,c[2005];

inline void get_prime()
{
    memset(ss,1,sizeof(ss));
    ss[0]=ss[1]=0;
    for (ll i=2;i*i<=A;i++)
    {
        if (ss[i]) zs[++cnt]=i;
        for (ll j=1;j<=cnt && zs[j]*i*zs[j]*i<=A;j++)
        {
            ss[zs[j]*i]=0;
            if (i%zs[j]==0) break;
        }
    }
}

inline ll get_c()
{
    ll n=A;
    for (ll i=1;i<=cnt;i++)
    {
        while (n%zs[i]==0)
        {
            n/=zs[i];
            c[i]++;
        }
    }
    return (n==1) ? 0 : n;
}

inline ll POW(ll x,ll y,ll mod)
{
    ll ans=1;
    while (y)
    {
        if (y&1) ans=ans*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ans;
}

void print(ll x)
{
    if (x==0) return;
    print(x/10);
    putchar(x%10+'0');
}

signed main()
{
    A=read(); B=read();
    get_prime();
    if (ll x=get_c()) c[++cnt]=1,zs[cnt]=x; //left
    ll ans=1;
    for (ll i=1;i<=cnt;i++)
    {
        if (!c[i]) continue;
        //a1=1,an=zs[i]^(B*c[i])
        ll mod=(zs[i]-1)*ha;
        ll an=POW(zs[i]%mod,B*c[i],mod);
        ll now=( (an*zs[i]-1)%mod )/(zs[i]-1); //(anq-a1)/q-1;
        ans=ans*now%ha;
    }
    print(ans);
    return 0;
}