Cod sursa(job #916550)

Utilizator romircea2010FMI Trifan Mircea Mihai romircea2010 Data 16 martie 2013 17:20:47
Problema Suma divizorilor Scor 100
Compilator cpp Status done
Runda Arhiva de probleme Marime 2.88 kb
/**
    suma divizorilor modulo 9901 lui A^B
    1 <= A, B <= 50000000

    SUMA = (1 + p1 + p1^2 + ... + p1^d1) * (1 + p2 + p2^2 + ... + p2^d2) * ... * (1 + pk + pk^2 + ... + pk^dk) =
    = produs de progresii geometrice cu primul termen 1 si ratia p1, p2, ... ,pk =
    (p1 ^ (d1+1) - 1)/(p1 - 1) * (p2 ^ (d2+1) - 1)/(p2 - 1) * ... * (pk ^ (dk+1) - 1)/(pk - 1)

    descompun pe A in factori primi si la fiecare exponent inmultesc B si aplic formula de mai sus.
*/

#include <iostream>
#include <fstream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define MOD 9901

using namespace std;

int A, B, sol, na, radA;
bool viz[8000];
int ciur[1500], nciur;

inline void CreareCiur()
{
    radA = (int)sqrt(A*1.0);
    ciur[++nciur] = 2;
    int i, j, doii, radradA;
    radradA = (int)sqrt(radA*1.0);
    for (i=3; i<=radA; i+=2)
        if (viz[i] == false)
        {
            ciur[++nciur] = i;
            if (i <= radradA)
                for (j=i*i, doii = i+i; j<=radA; j+=doii)
                    viz[j] = true;
        }
}

inline int Putere (int x, int y)
{
    /// x^y - 1
    int put = 1;
    while (y)
    {
        if (y&1)
        {
            put = (put * x%MOD)%MOD;
            y--;
        }
        y >>= 1;
        x = ((x%MOD)*(x%MOD))%MOD;
    }
    put--;
    if (put<0)
        put+=MOD;
    return put;
}

inline int InversModular(int x)
{
    /// x^(MOD-2)%MOD
    int y, put;
    put = 1;
    y = MOD-2;
    while (y)
    {
        if (y&1)
        {
            put = (put*(x%MOD))%MOD;
            y--;
        }
        y >>= 1;
        x = ((x%MOD)*(x%MOD))%MOD;
    }
    return put;
}

inline void Solve()
{
    int i, div, put;
    i = 1;
    sol = 1;
    while (i <= nciur && ciur[i]*ciur[i] <= A)
    {
        div = ciur[i];
        i++;
        if (A % div == 0)
        {
            put = 0;
            while (A % div == 0)
            {
                put++;
                A /= div;
            }
            put = put*B;
            sol = (((sol * Putere(div, put+1))%MOD) * InversModular(div-1))%MOD;
        }
    }
    if (A!=1)
    {
        if (A%MOD == 1)
        {
            /**
                daca A%MOD == 1
                sol = sol * (B+1) % MOD deoarece:
                A^(B+1) - 1 = (A-1) * (A^B + A^(B-1) + ... + A + 1)
                care se simplifica cu (A-1) de la numitor si obtinem
                A^B + A^(B-1) + ... + A + 1 care are B+1 termeni si care
                %MOD = B+1;
            */
            sol = sol*((B+1)%MOD)%MOD;
        }
        else
            sol = (((sol * Putere(A, B+1))%MOD) * InversModular(A-1))%MOD;
    }
}

int main()
{
    ifstream f("sumdiv.in");
    f>>A>>B;
    f.close();

    CreareCiur();
    Solve();

    ofstream g("sumdiv.out");
    g<<sol<<"\n";
    g.close();
    return 0;
}