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;
}