LightOj 1054 Efficient Pseudo Code

/*
*   Author: johnsondu
*   time: 2013-4-25
*   meanning: 求出约数和n ^ m的约数和,快速模取幂，扩展欧几里德，素数筛选
*   problem: LightOj 1054
*   url: http://lightoj.com/volume_showproblem.php?problem=1054
*
*/

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
using namespace std ;

#define max(x, y) (x > y ? x : y)
#define min(x, y) (x < y ? x : y)
#define LL long long
#define M 1000005
#define N 1005
const int Mod = 1000000007LL ;
bool ip[M] ;
int p[M], pl ;
LL n, m, pnum, ans ;

struct Node
{
    int num ;
    int prime ;
}node[N] ;

void init ()
{
    for (int i = 2; i < M; i ++)
        ip[i] = true ;
    pl = 0 ;
    for (int i = 2; i < M; i ++)
        if (ip[i])
        {
            p[pl ++] = i ;
            for (int j = 2 * i; j < M; j += i)
                ip[j] = false ;
        }
}

void divide (int t)
{
    pnum = 0 ;
    for (int i = 0; i < N; i ++)
        node[i].num = 0 ;
    for (int i = 0;i < pl; i ++)
    {
        if (t % p[i] == 0)
        {
            node[pnum].prime = p[i] ;
            while (t % p[i] == 0)
            {
                node[pnum].num ++ ;
                t /= p[i] ;
            }
            pnum ++ ;
        }
        if (p[i]*p[i] > t)
            break ;
    }
    if (t != 1)
    {
        node[pnum].prime = t ;
        node[pnum].num ++ ;
        pnum ++ ;
    }
}

void exgcd (LL a, LL b, LL &d, LL &x, LL &y)
{
    if (b == 0)
    {
        x = 1 ;
        y = 0 ;
        d = a ;
        return ;
    }
    exgcd (b, a%b, d, x, y) ;
    LL tmp = x ;
    x = y ;
    y = tmp - a/b * y ;
}

LL powerMod (LL a, LL b)
{
    LL res = 1 ;
    while (b)
    {
        if (b & 1)
        {
            res = ((res % Mod) * a) % Mod ;
            b -- ;
            continue ;
        }
        a = ((a % Mod) * a) % Mod ;
        b >>= 1 ;
    }
    return res ;
}

void solve ()
{
    ans = 1 ;
    for (int i = 0; i < pnum; i ++)
    {
        LL u = node[i].prime, v = node[i].num * m ;
        ans = (ans * (powerMod (u, v+1) - 1) % Mod + Mod) % Mod ;
        LL x, y, d ;
        exgcd (u-1, Mod, d, x, y) ;
        ans = (ans * x % Mod + Mod) % Mod ;
    }
    printf ("%lld\n", ans) ;
}

int main ()
{
    //freopen ("data.txt", "r", stdin) ;
    init () ;
    int tcase, cs = 1 ;
    scanf ("%d", &tcase) ;
    while (tcase --)
    {
        scanf ("%lld%lld", &n, &m) ;
        printf ("Case %d: ", cs ++) ;
        divide (n) ;

        solve () ;
    }
    return 0 ;
}