一本通1633【例 3】Sumdiv

1633：【例 3】Sumdiv

时间限制: 1000 ms 内存限制: 524288 KB

【题目描述】

原题来自：Romania OI 2002

求 ABAB 的所有约数之和 mod9901。

【输入】

输入两个整数 A,B。

【输出】

输出答案 mod9901。

【输入样例】

2 3

【输出样例】

【提示】

样例说明

2³=8，8 的所有约数为 1,2,4,8，1+2+4+8=15,15mod9901=15，因此输出 15。

数据范围与提示：

对于全部数据，0≤A,B≤5×10⁷。

sol：写了半天几乎吐血终于水过这道板子题

首先有个很显然的东西

对于一个数

如果他是p1^a1*p1^a2*p3^a3*~~*pn^an

那么他的因数和就是（p1⁰+p1¹+p1²+...+p1^a1）*（p2⁰+p2¹+...+p2^a2）*...*（pn¹+pn²+...+pn^an）

于是这道题就是要把所有质因数的幂次的乘积，先来推一下等比公式

令 S = p⁰+p¹+p²+...+p^a，（1）

则p*S = p¹+p²+...+p^a+p^a+1，（2）

（2）-（1）可以得到 S*(p-1) = p^a+1-p⁰，所以S= （p^a+1-p⁰）/ (p-1) ，前面一项显然Ksm就可以解决，考虑后一项

需要用到p-1的逆元，如何求呢？？

定义一个数a（就求a的逆元）

对于任意一个数 X/a ≡ X*Inva (%Mod)

易知 a*Inva ≡ 1 (%Mod)

随便推导下

原式: a*Inva = 1 (%Mod)
--> a*Inva-k*Mod = 1 (%Mod)
--> a*Inva+P*Mod = 1 (%Mod) (类ax+by=c的形式)

配上巨丑无比的代码

/*

原式: a*Inva = 1 (%Mod)

  --> a*Inva-k*Mod = 1 (%Mod)

  --> a*Inva+P*Mod = 1 (%Mod) (类ax+by=c的形式)

*/

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

inline ll read()

{

    ll s=;

    bool f=;

    char ch=' ';

    while(!isdigit(ch))

    {

        f|=(ch=='-'); ch=getchar();

    }

    while(isdigit(ch))

    {

        s=(s<<)+(s<<)+(ch^); ch=getchar();

    }

    return (f)?(-s):(s);

}

#define R(x) x=read()

inline void write(ll x)

{

    if(x<)

    {

        putchar('-'); x=-x;

    }

    if(x<)

    {

        putchar(x+''); return;

    }

    write(x/);

    putchar((x%)+'');

    return;

}

#define W(x) write(x),putchar(' ')

#define Wl(x) write(x),putchar('\n')

const ll Mod=;

ll A,B;

int cnt=;

struct Prime

{

    ll Shuz,Xis;

}Prim[];

inline ll Ksm(ll x,ll y)

{

    ll ans=;

    while(y)

    {

        if(y&)

        {

            ans=ans*x%Mod;

        }

        x=x*x%Mod;

        y>>=;

    }

    return ans;

}

ll gcd(ll x,ll y)

{

    return (!y)?(x):(gcd(y,x%y));

}

inline void Exgcd(ll a,ll b,ll &X,ll &Y)

{

    if(b==)

    {

        X=;

        Y=;

        return;

    }

    Exgcd(b,a%b,X,Y);

    ll XX=X,YY=Y;

    X=YY;

    Y=XX-a/b*YY;

    return;

}

//a*Inva+Mod*k = 1 (%Mod)

inline ll Inv(ll Num)

{

    ll a,b,c,r,X,Y;

    a=Num;

    b=Mod;

    c=;

    r=gcd(a,b);

    Exgcd(a,b,X=,Y=);

    X=X*r/c;

    ll tmp=b/r;

    X=(X>=)?(X%tmp):(X%tmp+tmp);

//    printf("X=%lld Test=%lld\n",X,Num*X%Mod);

    return X;

}

inline void Solve(ll Num)

{

    int i;

    for(i=;i<=sqrt(Num);i++)

    {

        if(Num%i==)

        {

            Prim[++cnt].Shuz=i;

            Prim[cnt].Xis=;

            while(Num%i==)

            {

                Prim[cnt].Xis++;

                Num/=i;

            }

            Prim[cnt].Xis*=B;

        }

    }

    if(Num>)

    {

        Prim[++cnt]=(Prime){Num,B};

    }

    return;

}

int main()

{

    int i;

    ll ans=1ll;

    R(A); R(B);

    Solve(A);

    for(i=;i<=cnt;i++)

    {

//        printf("%lld %lld\n",Prim[i].Shuz,Prim[i].Xis);

        ll tmp=(Ksm(Prim[i].Shuz,Prim[i].Xis+)-+Mod)%Mod;

        tmp=tmp*Inv(Prim[i].Shuz-)%Mod;

        ans=ans*tmp%Mod;

    }

    Wl(ans);

    return ;

}

/*

input

2 3

output

15

input

10 6

output

5187

input

45279441 39876543

output

6060

input

9 8

output

5660

*/

秒客网