POJ 1845 Sumdiv （求某个数的所有正因子的和）

题意： 求A^B的所有正因子的和，最后模9901的结果。

思路：

若对一个数n进行素数分解，n=p1^a1*p2^a2*p3^a3*...*pk^ak
那么n的所有正因子之和sum=(1+p1+...+p1^a1)*(1+p2+...+p2^a2)*...*(1+pk+...+pk^ak)
然后可以用等比数列求和公式(pk^(ak+1)-1)/(pk-1)求每项的和，再累乘。
用等比数列求1+pk+...+pk^ak时候要注意几点：

1.这里有除法，所以模的时候要将除以分母转化成乘以分母的逆元
a = (b/c) ==> a%m = b*c^(m-2)%m ( m为素数 )
证明：
b = a * c
根据费马小定理 a^(p-1)= 1 %p;(p是素数且a不能整除p)
所以 c^(m-1)%m=1%m
因此 a % m = a*1%m = a * c^(m-1)%m = a*c*c^(m-2)%m = b*c^(m-2)%m;

2.等比求和公式要注意，分母pk-1不能为0。
当pk%mod=1的时候， (1+pk+...+pk^ak)%mod=ak+1

3.当pk%mod=0的时候，(1+pk+...+pk^ak)=1，可以直接pass即可

 

还有从discuss里看到别人求和的做法：
1.可以首先计算rem=p^(cB+1)%(MOD*(p-1))
2.然后计算rem=(rem-1+MOD*(p-1))/(p-1)
3.最后计算rem%MOD
这里MOD*(p-1)可能会超过32位，所以在计算p^(cB+1)%(MOD*(p-1))时候可能乘法会溢出....careful....
证明：
令t*(p-1)=(p^0+p^1+p^2+...+p^a)*(p-1)=p^(a+1)-1 (mod m*(p-1))
因为gcd(p-1,m*(p-1))=p-1且 (p-1)|p^(a+1)-1 ，所以t=(p^(a+1)-1)/(p-1) (mod m)
因此我们可以先求p^(a+1)-1 (mod m*(p-1))，再把这个值除以(p-1)后取mod m

#include <iostream>
#include <string.h>
#include <algorithm>
#include <stdio.h>
#include <math.h>

using namespace std;
long long A,B;
const int mod=9901;
const int maxn=100005;
bool isprime[maxn];
int prime[maxn];
int cnt=0;

void init(){
    memset(isprime,true,sizeof(isprime));
    for(int i=2;i<maxn;i++){
        if(isprime[i]){
            prime[cnt++]=i;
            for(int j=2*i;j<maxn;j+=i)
                isprime[j]=false;
        }
    }
}
long long quickPow(long long a,long long b){
    long long ret=1;
    while(b){
        if(b&1)
            ret=ret*a%mod;
        a=a*a%mod;
        b=b>>1;
    }
    return ret;
}
long long solve(long long a,long long b){
    int c;
    long long ans=1;
    for(int i=0;i<cnt;i++){
        c=0;
        if(a%prime[i]==0){
            while(a%prime[i]==0){
                c++;
                a=a/prime[i];
            }
            if(prime[i]%mod==0)
                continue;
            else if(prime[i]%mod==1){
                //(1+p+p^2+...+p^cb),此时求和不能用等比，因为分母(p-1)!=0。模9901后的和为c+1
                ans=(ans*(c*b+1))%mod;  //注意：是c*b+1,一开始写成了c+1。
            }
            else{
                long long ret=(quickPow((long long)prime[i],(long long)(c*b+1))-1+mod)%mod;
                ans=(ans*(ret*quickPow((long long)(prime[i]-1),(long long)(mod-2))%mod))%mod;
            }
        }
    }
    //不要忘记了最后的a若大于1，则为素因子
    if(a>1){
        if(a%mod==0)
            return ans;
        else if(a%mod==1){
            ans=(ans*(b+1))%mod;  //注意：由于这里c=1，所以是b+1,一开始写成了2，也忘记乘以b了。
        }
        else{
            long long ret=(quickPow(a,(long long)(b+1))-1+mod)%mod;
            ans=(ans*(ret*quickPow(a-1,(long long)(mod-2))%mod))%mod;
        }
    }
    return ans;
}
int main()
{
    init();
    while(scanf("%I64d%I64d",&A,&B)!=EOF){
        printf("%I64d\n",solve(A,B));
    }
    return 0;
}

View Code