hdu1845(a^b的因子和%p)

题目链接：http://poj.org/problem?id=1845

思路：

１．整数唯一分解定理：

      任意正整数都有且只有一种方式写出其素因子的乘积表达式。

      a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)   其中pi均为素数

２．约数和公式：

　　对于已经分解的整数a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)

　　有a的所有因子之和为

    　S` = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)

那么 a^b 的所有因子和为

　　S = (1+p1+p1^2+p1^3+...p1^(k1*b)) * (1+p2+p2^2+p2^3+….p2^(k2*b)) * (1+p3+ p3^3+…+ p3^(k3*b)) * .... * (1+pn+pn^2+pn^3+...pn^(kn*b))

对于数列 1, p, p^2, p^3 ...  p^n % mod，其中  mod 为质数，打个表可以发现该数列是一个循环数列，其中存在一个循环节为 1, p, p^1, ... p^(mod-2)．其实这点在费马小定理中是有体现的，a^(mod-1) = 1 (% mod)．

那么对于求 cnt = 1 + p + p^2 + ... + p^n % mod，可以令

　　cc1 = (n + 1) / (mod - 1)

　　cc2 = (n + 1) % (mod - 1)

　　cnt1 = 1 + p + p^2 + ... + p^(mod - 2)

　　cnt2 = 1 + p + p^2 + ... + p^(cc2 - 1)

那么 cnt = cc1 * cnt1 + cnt2 % mod

对于求 a^b，可以先将 a 质因分解，得到 S 的表达式，对于 S 表达式，只需要按照上面的方法求出其中每个乘数项即可．时间复杂度为 O(loga * mod), 本题中 mod = 9901, 时间上是允许的．

代码：

#include <iostream>
#include <stdio.h>
#include <map>
#define ll long long
using namespace std;

const int mod = 9901;
const int MAXN = 1e5 + 10;
ll prime[MAXN], indx = 0;
map<int, int> num;

int get(ll a, ll b){//计算sigma(a^i),其中0<=i<=b
    ll sol1 = 0, sol2 = -1, cnt = 1;
    ll cc1 = (b + 1) / (mod - 1);
    ll cc2 = (b + 1) % (mod - 1);
    for(int i = 0; i < mod - 1; i++){
        sol1 = (sol1 + cnt) % mod;
        if(sol2 == -1 && i + 1 == cc2) sol2 = sol1;
        cnt = (cnt * a) % mod;
    }
    return (sol1 * cc1 % mod + sol2) % mod;
}

int main(void){
    ll a, b, sol = 1;
    cin >> a >> b;
    for(int i = 2; i * i <= a; i++){
        if(a % i == 0){
            prime[indx] = i;
            while(a % i == 0){
                num[i]++;
                a /= i;
            }
            indx++;
        }
    }
    if(a > 1) prime[indx++] = a, num[a]++;
    for(int i = 0; i < indx; i++){
        sol = (sol * get(prime[i], b * num[prime[i]])) % mod;
    }
    cout << sol << endl;
    return 0;
}

但是当 mod 比较大时这个方法就不行了，mod 比较大时可以用下面这个代码，贴的 Kuangbin 模板

代码：

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <algorithm>
#include <string.h>
#define ll long long
using namespace std;

//******************************************
//素数筛选和合数分解
const int MOD = 9901;
const int MAXN=10000;
int prime[MAXN+1];

void getPrime(void){
    memset(prime, 0, sizeof(prime));
    for(int i = 2; i <= MAXN; i++)
    {
        if(!prime[i]) prime[++prime[0]] = i;
        for(int j=1; j<= prime[0]&&prime[j] <= MAXN/i; j++)
        {
            prime[prime[j]*i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

ll factor[100][2];
int fatCnt;

int getFactors(ll x){
    fatCnt = 0;
    ll tmp = x;
    for(int i=1; prime[i] <= tmp/prime[i]; i++){
        factor[fatCnt][1] = 0;
        if(tmp % prime[i] == 0){
            factor[fatCnt][0] = prime[i];
            while(tmp % prime[i] == 0){
                factor[fatCnt][1]++;
                tmp /= prime[i];
            }
            fatCnt++;
        }
    }
    if(tmp!=1)
    {
        factor[fatCnt][0]=tmp;
        factor[fatCnt++][1]=1;
    }
    return fatCnt;
}

//******************************************
ll pow_m(ll a, ll n)//快速模幂运算
{
    ll res = 1;
    ll tmp = a % MOD;
    while(n){
        if(n & 1){ 
            res *= tmp;
            res%=MOD;
        }
        n >>= 1;
        tmp *= tmp;
        tmp %= MOD;
    }
    return res;
}

ll sum(ll p, ll n){//计算1+p+p^2+...+p^n
    if(p == 0)return 0;
    if(n == 0)return 1;
    if(n & 1){//奇数
        return ((1 + pow_m(p, n/2 + 1)) % MOD * sum(p, n / 2) % MOD) % MOD;
    }else return ((1 + pow_m(p, n / 2 + 1)) % MOD * sum(p, n / 2 - 1) + pow_m(p, n / 2) % MOD) % MOD;

}

int main(void){
    int A, B;
    getPrime();
    while(scanf("%d%d", &A, &B) != EOF){
        getFactors(A);
        ll ans = 1;
        for(int i = 0; i < fatCnt; i++){
            ans *= (sum(factor[i][0], B * factor[i][1]) % MOD);
            ans %= MOD;
        }
        printf("%lld
",ans);
    }
    return 0;
}