poj1845--Sumdiv(数论篇3--真滴是数论啊。。。。)

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 14953		Accepted: 3680

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).

Source

Romania OI 2002

给出n和m，求出n^m的所有的约数和。

首先 n^m = p1^k1 + p2^k2 + ....pn^kn ; 那么约数和 = ( 1 + p1 + p1^2 .. p1^k1 ) * ( 1 + p2 + p2^2 ... p2^k2 ) *()...(1 + pn + pn^2 + ...pn^kn) ;

得到约数和的计算公式后，就可以向那个公式转化，n的素数组成（p1,p2,p3。。）和n^m是相同的，所以先计算的n的p1->k1, p2->k2，将k1,k2乘以m就得到了n^m的组合。

求1 + p + p^2 ... + p^n 有两种方式

1.递归求解。

2.公式  = a1*(1 - p^(n+1))/(1-p) = ( p^(n+1) - 1 ) / (p - 1 ) ;

减法取模，如果结果小于0，加上mod ；

要对公式进行取模的话，要用到除法取模。也有两种方法。a / b  % mod ;

1.转化 ( a % (mod * b) ) / b % mod ;

2.逆元，a / b % mod = a * b^( mod-2 ) % mod ;

逆元要求 b不能是mod的倍数。

 

 

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
#define MOD 9901
#define LL __int64
LL vis[1100000] , pri[1100000] , cnt ;
LL a[1100000] , num[1100000] ;
void get_pri()
{
    LL i , j ;
    memset(vis,0,sizeof(vis)) ;
    cnt = 0 ;
    for(i = 2 ; i < 1010000 ; i++)
    {
        if( !vis[i] )
            pri[cnt++] = i ;
        for(j = 0 ; j < cnt ; j++)
        {
            if( pri[j] * i > 1010000 ) break ;
            vis[ pri[j]*i ] = 1 ;
            if(i%pri[j]==0) break ;
        }
    }
    return ;
}
LL pow(LL a,LL k,LL mod)
{
    if( k == 0 ) return 1 ;
    LL b = pow(a,k/2,mod) ;
    b = b*b % mod ;
    if( k%2 ) b = b*a%mod ;
    return b ;
}
int main()
{
    get_pri() ;
    LL n , m , i , j ;
    LL ans , temp ;
    while( scanf("%I64d %I64d", &n, &m) != EOF )
    {
        ans = 1 ;
        memset(num,0,sizeof(num)) ;
        j = 0 ;
        for(i = 0 ; i < cnt ; i++)
        {
            while( n % pri[i] == 0 )
            {
                n /= pri[i] ;
                a[j] = pri[i] ;
                num[j]++ ;
            }
            if( num[j] )
                j++ ;
            if( n == 1 ) break ;
        }
        if( n > 1 )
        {
            a[j] = n ;
            num[j++] = 1 ;
        }
        for(i = 0 ; i < j ; i++)
        {
            if( (a[i]-1)%MOD == 0 )
            {
                temp = pow(a[i],num[i]*m+1,MOD*(a[i]-1)) - 1 ;
                if( temp < 0 ) temp += MOD*(a[i]-1) ;
                ans = ans * ( temp / (a[i]-1) % MOD ) % MOD ;
            }
            else
            {
                temp = pow(a[i],num[i]*m+1,MOD) - 1 ;
                if( temp < 0 ) temp += MOD ;
                ans = ans * (temp) % MOD * pow(a[i]-1,MOD-2,MOD) % MOD ;

            }
        }
        printf("%I64d\n", ans) ;
    }
    return 0 ;
}