[poj1845] Sumdiv(数论,质因数分解,约数和,快速幂)

Sumdiv

Time Limit: 1000MS Memory Limit: 30000K

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


解题思路

题意:求 AB A B 的所有约数和模9901

首先我们需要知道约数和公式

由唯一分解定理,对 AZ ∀ A ∈ Z ,可以将 A A 写作:A=p1k1p2k2pnkn,其中 pi p i 是质数,那么 A A 的约数和为i=1nj=0kipij

这道题中,可以先把 A A 分解质因数,假设A=p1k1p2k2pnkn,那么 AB=pBk11pBk22pBknn A B = p 1 B k 1 p 2 B k 2 ⋯ p n B k n ,所以答案就是 j=0Bkipjimod9901 ∏ ∑ j = 0 B k i p i j mod 9901
所以现在的问题变成了:
1. 如何分解质因数
A A 不断除2直至除尽,然后不断除3直至除尽……由于是从小到大不断相
除,可以保证A不会除到一个合数(因为 A A 一定已经被这个合数的最小质因子除尽了),除的同时记录A的因子及其出现次数即可

int fac = 2;
while(fac * fac <= A){
    while(A % fac == 0){
        if(p[p[0]] != fac)  p[++p[0]] = fac;
        k[p[0]]++;
        A /= fac;
    }
    if(fac == 2)    fac++;
    else    fac += 2;
}
if(A != 1)  p[++p[0]] = A, k[p[0]] = 1;

2.如何快速地求出 j=0npjmod9901 ∑ j = 0 n p j mod 9901 ,即 (p0+p1++pn)mod9901 ( p 0 + p 1 + ⋯ + p n ) mod 9901

  • 法一:
    • 如果 n n 是奇数,
      (p0+p1+...+pn1)+(pn+pn+1+...+p2n1)=(p0+p1+p2+...+pn1)(1+pn)
    • 如果 n n 是偶数,
      (p0+p1+...+pn1)+pn+(pn+1+...+p2n)=pn+(p0+p1+p2+...+pn1)(1+pn+1)
    • 求和部分递归即可, pn p n pn+1 p n + 1 快速幂即可
      这样就可以 O(log2N) O ( l o g 2 N ) 地计算了
  • 法二:
    当然,也可以直接用等比数列求和公式, (p0+p1++pn)=1pn+11p ( p 0 + p 1 + ⋯ + p n ) = 1 − p n + 1 1 − p ,但是由于要取模,所以注意乘逆元,另外,又因为 9901 9901 是一个小质数,用快速幂时还需要特判

Code

#include<cstdio>
#include<cstring>

using namespace std;
typedef long long LL;

const int MOD = 9901;
int A, B, p[10001], k[10001], ans;

LL fPow(LL bs, LL idx){
    LL res = 1;
    bs %= MOD;
    while(idx){
        if(idx & 1) (res *= bs) %= MOD;
        (bs *= bs) %= MOD;
        idx >>= 1;
    }
    return res;
}

LL cal(int p, LL k){
    if(k == 0)  return 1ll;
    if(k & 1)   return cal(p, k/2) * (1ll + fPow(p, k/2+1)) % MOD;
    else    return (fPow(p, k/2) + cal(p, k/2-1) * (1ll + fPow(p, k/2+1)) % MOD) % MOD;
}

int main(){
    while(scanf("%d%d", &A, &B) != EOF){
        ans = 1;
        p[0] = 0;
        int fac = 2;
        while(fac * fac <= A){
            while(A % fac == 0){
                if(p[p[0]] != fac)  p[++p[0]] = fac;
                k[p[0]]++;
                A /= fac;
            }
            if(fac == 2)    fac++;
            else    fac += 2;
        }
        if(A != 1)  p[++p[0]] = A, k[p[0]] = 1;
        for(int i = 1; i <= p[0]; i++)
            (ans *= cal(p[i], 1ll * k[i] * B)) %= MOD;
        printf("%d\n", ans);
    }
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值