P1593 因子和

P1593 因子和
新算法：
#define ni 逆元
先质因数分解，
（1+p1^1+p1^2...p1^x）*(1+p2^1+p2^2...p2^x)
然后套等比数列公式就可以了。

#include<iostream>
#include<cstdio>
#include<queue>
#include<algorithm>
#include<cmath>
#include<ctime>
#include<set>
#define mod 9901
#include<cstring>
#define inf long long_MAX
#define For(i,a,b) for(register long long i=a;i<=b;i++)
#define p(a) putchar(a)
#define g() getchar()
//by war
//2017.10.27
using namespace std;
long long num,b;
struct node
{
    long long cnt;
    long long p;
}a[10000];//40M
long long prime[10000];//40M
long long cnt;
long long tot;
bool vis[10000];
long long ans;

void in(long long &x)
{
    long long y=1;
    char c=g();x=0;
    while(c<'0'||c>'9')
    {
    if(c=='-')
    y=-1;
    c=g();
    }
    while(c<='9'&&c>='0')x=x*10+c-'0',c=g();
    x*=y;
}
void o(long long x)
{
    if(x<0)
    {
        p('-');
        x=-x;
    }
    if(x>9)o(x/10);
    p(x%10+'0');
}

void Euler(long long x)
{
    For(i,2,x)
    {
        if(!vis[i])prime[++cnt]=i;
        for(register long long j=1;j<=cnt&&prime[j]*i<=x;j++)
        {
            vis[prime[j]*i]=true;
            if(i%prime[j]==0)
            break;
        }
    }
}

void resolve(long long x)
{
    For(i,1,cnt)
    {
        if(x%prime[i]==0)
        {
            a[++tot].p=prime[i];
            while(x%prime[i]==0)
            {
            a[tot].cnt++;
            x/=prime[i];
            }
        }
    }
    if(x>1)
    {
    a[++tot].p=x;
    a[tot].cnt++;    
    }
}

long long ksm(long long a,long long b)
{
    if(b==0)
    return 1;
    while(b%2==0)
    {
        a=(a*a)%mod;
        b>>=1;
    }
    long long r=1;
    while(b>0)
    {
        if(b%2==1)
        r=(r*a)%mod;
        a=(a*a)%mod;
        b>>=1;
    }
    return r%mod;
}

void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    long long temp=x;
    x=y;
    y=temp-(a/b)*y;
}

long long series(long long q,long long n)
{
    long long fz=ksm(q,n)-1;
    long long x,y,b;
    exgcd(q-1,mod,x,y);
    long long ni=x;
    ni=(ni%mod+mod)%mod;
    return (fz*ni%mod+mod)%mod;
}

int main()
{
    in(num),in(b);
    Euler(sqrt(num));
    resolve(num);
    ans=1;
    For(i,1,tot)
    ans=ans*series(a[i].p,a[i].cnt*b+1)%mod;
    o(ans%mod);
    return 0;
}