poj 1845 Sumdiv (算术基本定理求一个数因子和)

求一个数的所有因子和可以用算术基本定理，下面是它的两个重要应用：
        （1）一个大于1的正整数N，如果它的标准分解式为： N=(P1^a1)*(P2^a2)......(Pn^an)
            那么它的正因数个数为（1+a1)(1+a2).....(1+an)。
        （2） 它的全体正因数之和为d(N)=(1+p1+...p1^an)(1+p2+...p2^a2)...(1+pn+...+pn^an)
            和求一个数正因数个数的方法类似.
    可以先打表出sqrt(n)以内的所有素数（当然也可以不打表）,因为n的素因数中最多只有一个大于sqrt(n)(以前题目中有证明过),
    所以可以最后处理它。对于每一项都是一个等比数列，求和很容易。


    下面以poj 1845 Sumdiv为例。
    求n^m%p,p为素数9901.
    1)因为n和m很大，所以可以处理出n的所有素因数,然后^m即每个因数指数*m。此题可以不用预处理所有素数，直接分解素因数,对于这道题时间复杂度可接受;
    2)用到快速幂运算。

    3）在求等比数列和取模的时候两种方法：1. 用等比数列求和公式：S=a1*（q^n-1）/(q-1)， 要用到(q-1)模9901的逆元，可以用欧拉定理或扩展欧几里得来求，但因为 mod 为素数,所以a的逆元为a^(mod-2)%mod,用快速幂来求。 2.递归形式的二分

    对于逆元补充：费马小定理(Fermat Theory)是数论中的一个重要定理，其内容为： 假如p是质数，且Gcd(a,p)=1，那么 a(p-1) ≡1（mod p）    即：假如a是整数，p是质数，且a,p互质(即两者只有一个公约数1)，那么a的(p-1)次方除以p的余数恒等于1。

  //逆元法求解

// 逆元解法，有些小细节。。 ac: 16ms
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<cmath>
#define maxn 10000
#define ll long long
using namespace std;
const ll MOD=9901;
ll prime[maxn],num,flag[maxn+5];
ll pow_mod(ll a,ll b,ll p){ //快速幂
    ll ret=1;
    while(b){
        if(b&1) ret=(ret*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return ret;
}
void get_prime(){          //筛素数
    memset(flag,0,sizeof(flag));
    num=0;
    for(ll i=2;i<maxn;i++){
        if(!flag[i]) prime[++num]=i;
        for(ll j=1;j<=num&&i*prime[j]<maxn;j++){
            flag[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}

ll get_sum(ll k,ll n){   //计算 (1+k+k^2+...+k^n)%mod ,等比数列和为(k^(n+1)-1)/(k-1);
    if(n==0) return 1;
    ll ret=1;
    ret=ret*(pow_mod(k,n+1,MOD)-1)%MOD;
    ret=ret*pow_mod(k-1,MOD-2,MOD)%MOD;       // (k-1)的逆元,要注意(k-1)%9901 != 0
    return ret;
}
int main(){
    ll n,m,k,a;
    get_prime();
    while(~scanf("%lld%lld",&n,&m)){
        ll ans=1; if(!n&&m) ans=0;     //0^0=1 , 0^1..=0;
        //素因数分解
        for(int i=1;i<=num&&prime[i]*prime[i]<=n;i++){
            if(n<prime[i]) break;
            if(n%prime[i]==0){
                a=0;               //指数
                while(n%prime[i]==0){
                    n/=prime[i];
                    a++;
                }
                a*=m;
                k=prime[i];
                ans=ans*get_sum(k,a)%MOD;
            }
        }
        //如果还有一个大于sqrt(n)的素因数,它的指数一定为1*m=m
        if(n>1&&(n-1)%MOD==0) ans=ans*(pow_mod(n,m+1,MOD*(n-1))/(n-1))%MOD;   //如果(n-1)%mod==0,不能求逆元了
        else if(n>1) ans=ans*get_sum(n,m)%MOD;
        printf("%lld\n",(ans+MOD)%MOD);                //刚开始没加mod wa了好几遍，因为这种方法ans可能为负
    }
    return 0;
}

 //二分法求解

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<cmath>
#define maxn 10000
#define ll long long
using namespace std;
const ll MOD=9901;
ll prime[maxn],num,flag[maxn+5];
ll pow_mod(ll a,ll b,ll p){ //快速幂
    ll ret=1;
    while(b){
        if(b&1) ret=(ret*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return ret;
}
void get_prime(){          //筛素数
    memset(flag,0,sizeof(flag));
    num=0;
    for(ll i=2;i<maxn;i++){
        if(!flag[i]) prime[++num]=i;
        for(ll j=1;j<=num&&i*prime[j]<maxn;j++Accepted){
            flag[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}
ll get_sum(ll k,ll n){    //可以用递归形式的二分来计算 (1+k+k^2+...+k^n)%mod
    if(n==0) return 1;
    if(n&1)
        return ((1+pow_mod(k,n/2+1,MOD))*get_sum(k,n/2))%MOD;
    else
        return ((1+pow_mod(k,n/2+1,MOD))*get_sum(k,n/2-1)+pow_mod(k,n/2,MOD))%MOD;
}
int main(){
    ll n,m,k,a;
    get_prime();
    while(~scanf("%lld%lld",&n,&m)){
        ll ans=1; if(!n&&m) ans=0;     //0^0=1 , 0^1..=0;
        for(int i=1;i<=num&&prime[i]*prime[i]<=n;i++){
            if(n<prime[i]) break;
            if(n%prime[i]==0){
                a=0;               //指数
                while(n%prime[i]==0){
                    n/=prime[i];
                    a++;
                }
                a*=m;
                k=prime[i];
                ans=ans*get_sum(k,a)%MOD;
            }
        }
        if(n>1) ans=ans*get_sum(n,m)%MOD;
        printf("%lld\n",ans);
    }
    return 0;
}