Cod sursa(job #1937863)

Utilizator GoogalAbabei Daniel Googal Data 24 martie 2017 12:46:37
Problema Suma divizorilor Scor 40
Compilator cpp Status done
Runda Arhiva de probleme Marime 1.65 kb
#include <fstream>
#define MOD 9901
#define NMAX 1000005

using namespace std;

ifstream fin("sumdiv.in");
ofstream fout("sumdiv.out");

long long a,b,prim[NMAX],k;
bool ciur[NMAX];

void ciur_er()
{
    long long i,j;

    for(i=2;i<=NMAX;i++)
    {
        if(!(ciur[i]))
        {
            prim[++k]=i;
            for(j=i+i;j<=NMAX;j+=i)
                ciur[j]=true;
        }
    }
}

long long ridic_log(long long n,long long p)
{
    long long a=n,i,sol=1;
    for(i=0;(1<<i)<=p;i++)
    {
        if(((1<<i)&p)>0)
            sol=(sol*a)%MOD;
        a=(a*a)%MOD;
    }

    return sol;
}

void inv_mod(long long &x, long long &y, long long a, long long b)
{
    if(!b)
    {
        x=1;
        y=0;
    }

    else
    {
        inv_mod(x,y,b,a%b);
        swap(x,y);
        y-=x*(a/b);
    }
}

long long turn_it(long long a)
{
    long long inv=0,ins;
    inv_mod(inv,ins,a,MOD);
    if(inv<=0)
        inv=MOD+inv%MOD;
    return inv;
}

int main()
{
    long long j,p,p1,p2,r=1;
    ciur_er();
    fin>>a>>b;
    fin.close();
    for(j=1;j<=k && prim[j]*prim[j]<=a;j++)
    {
        if(a%prim[j])
            continue;

        p=0;
        while(!(a%prim[j]))
        {
            p++;
            a/=prim[j];
        }
        p*=b;

        p1=(ridic_log(prim[j],p+1)-1)%MOD;
        p2=turn_it(prim[j]-1)%MOD;

        r=(((r*p1)%MOD)*p2)%MOD;
    }

    if(a>1)
    {
        long long x=a,i;
        r=1;
        for(i=1;i<=b;i++)
        {
            r=(r+x)%MOD;
            x=(x*a)%MOD;
        }
    }
    fout<<r;
    fout.close();
    return 0;
}