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

标签：数论   算法   poj   算术基本定理   

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

  //逆元法求解

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

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

标签：数论   算法   poj   算术基本定理   

原文地址：http://blog.csdn.net/crazy__c/article/details/44904913