POJ 1845 Sumdiv 【二分 || 逆元】

任意门：http://poj.org/problem?id=1845、

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 30268		Accepted: 7447

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

 

 

题目概括：

给出 A 和 B ，求 A^B 所有因子之和，答案 mod 9901；

 

解题思路：

可对 A 先进行素因子分解 A = p1^a1 * p2^a2 * p3^a3 * ...... * pn^an

即 A^B = p1^a1*B * p2^a2*B * p3^a3*B * ...... * pn^an*B

那么 A^B 所有因子之和 = (1 + P1¹ + P1² + P1³ + ... + P1^a1*B) *  (1 + P2¹ + P2² + P2³ + ... + P2^a2*B) * ......* (1 + Pn¹ + Pn² + Pn³ + ... + Pn^an*B) ;

 

对于  (1 + P¹ + P² + P³ + ... + P^a1*B) 我们可以

①二分求和 （47 ms)

AC code：

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define INF 0x3f3f3f3f
#define LL long long
using namespace std;
const LL MOD = 9901;
const int MAXN = 1e4+10;
LL p[MAXN];
bool isprime[MAXN];
int cnt;
LL A, B;

void init()                               //打表预处理素因子
{
    cnt = 1;
    memset(isprime, 1, sizeof(isprime));
    for(LL i = 2; i < MAXN; i++){
        if(isprime[i]){
            p[cnt++] = i;
            for(int j = 1; j < cnt && p[j]*i < MAXN; j++)
                isprime[p[j]*i] = false;
        }
    }
}

LL qpow(LL x, LL n)
{
    LL res = 1LL;
    while(n){
        if(n&1) res = (res%MOD*x%MOD)%MOD;
        x = (x%MOD*x%MOD)%MOD;
        n>>=1LL;
    }
    return res;
}

LL sum(LL x, LL n)
{
    if(n == 0) return 1;
    LL res = sum(x, (n-1)/2);
    if(n&1){
        res = (res + res%MOD*qpow(x, (n+1)/2)%MOD)%MOD;
    }
    else{
        res = (res + res%MOD*qpow(x, (n+1)/2)%MOD)%MOD;
        res = (res + qpow(x, n))%MOD;
    }
    return res;
}

int main()
{
    LL ans = 1;
    scanf("%lld %lld", &A, &B);
    init();
    //printf("%d
", cnt);
    for(LL i = 1; p[i]*p[i] <= A && i < cnt; i++){          //素因子分解
        //printf("%lld
 ", p[i]);
        if(A%p[i] == 0){
            LL num = 0;
            while(A%p[i] == 0){
                num++;
                A/=p[i];
            }
            ans = ((ans%MOD*sum(p[i], num*B)%MOD)+MOD)%MOD;
        }
    }
    if(A > 1){
       // puts("zjy");
        ans = (ans%MOD*sum(A, B)%MOD+MOD)%MOD;
    }

    printf("%lld
", ans);

    return 0;
}

②应用等比数列求和公式 ( 0ms )

因为要保证求模时的精度，所以要求逆元。

这里用的方法是一般情况都适用的 ： A/B mod p = A mod (p*B) / B;

但是考虑乘法会爆int，所以自定义一个二分乘法。

AC code:

// 逆元 + 二分乘法
#include <cstdio>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>
#define INF 0x3f3f3f3f
#define LL long long
using namespace std;
const int MAXN = 1e4+10;
const LL mod  = 9901;
int p[MAXN], cnt;
bool isp[MAXN];
LL A, B;

void prime()
{
    memset(isp, 1, sizeof(isp));
    isp[0] = isp[1] = false;
    for(int i = 2; i < MAXN; i++){
        if(isp[i]){
            p[cnt++] = i;
            for(int k = 0; k < cnt && i*p[k] < MAXN; k++){
                isp[i*p[k]] = false;
            }
        }
    }
}

LL mutil(LL x, LL y, LL MOD)
{
    LL res = 0;
    while(y){
        if(y&1) res = (res+x)%MOD;
        x = (x+x)%MOD;
        y>>=1LL;
    }
    return res;
}

LL q_pow(LL x, LL n, LL MOD)
{
    LL res = 1LL;
    while(n){
        if(n&1) res = mutil(res, x, MOD)%MOD;
        x = mutil(x, x, MOD)%MOD;
        n>>=1LL;
    }
    return res;
}

int main()
{
    LL num = 0;
    LL ans = 1;
    scanf("%lld %lld", &A, &B);
    prime();
    for(int i = 0; p[i]*p[i] <= A; i++){
        if(A%p[i] == 0){
            num = 0;
            while(A%p[i] == 0){
                A/=p[i];
                num++;
            }
            LL m = mod*(p[i]-1);
            ans *= (q_pow(p[i], num*B+1, m) + m-1)/(p[i] - 1);
            ans = ans%mod;
            //ans += ((q_pow(p[i], num*B, mod*(p[i]-1))-p[i])/(p[i]-1))%mod;
        }
    }

    if(A != 1){
//        ans += ((q_pow(A, B, mod*(A-1))-A)/(A-1))%mod;
        LL m = mod*(A-1);
        ans*=(q_pow(A, B+1, m)+m-1)/(A-1);
        ans = ans%mod;
    }

    printf("%lld
", ans);

    return 0;
}