1. 程式人生 > >[洛谷P1593][POJ-1845Sumdiv] 因子和

[洛谷P1593][POJ-1845Sumdiv] 因子和

題目連結:

洛谷

POJ

題意:

輸入兩個正整數a和b,求a^b的因子和。結果太大,只要輸出它對9901的餘數。0≤a,b≤50000000

思路:根據唯一分解定理,a^b=(b1^(p1*b))*(b2^(p2*b))*...*(bn^(pn*b)),那麼a^b的因子和就是

(b1^0+b1^1+...+b1^(p1*b))*...*(bn^0+bn^1+...+bn^(pn*b)),可用等比數列求和來做,分母是(bi-1),當bi-1是9901的倍數的時候,求的逆元是0,不符合我們的要求,所以對於(bi-1)%9901=0的情況特判一下,因為bi%9901=1,所以根據唯一分解定理,

(1+bi^1+bi^2+...+bi^pi)%9901=pi+1。

程式碼:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define R register
#define ll long long int
#define ull unsigned long long
using namespace std;
const int N=7100;
ll a,b,c,mod=9901,ans=1,num,prime[N],notprime[N];
ll su[N],yue[N],top;
inline ll ksm(R ll x,R ll p){
    R ll tot
=1; while(p){ if(p&1){ tot=(tot*x)%mod; } x=(x*x)%mod; p>>=1; } return tot%mod; } void pre(){ notprime[1]=1; for(R int i=2;i<=sqrt(c);i++){ if(notprime[i]==0) prime[++num]=i; for(R int j=1
;(j<=num)&&(i*prime[j]<=sqrt(c));j++){ notprime[i*prime[j]]=1; if(i%prime[j]==0)break; } } } int main(){ scanf("%lld%lld",&a,&b); c=a; pre(); for(R int i=1;i<=num;i++){ R int k=prime[i]; if(c%k==0){ su[++top]=k; yue[top]=0; while(c%k==0){ c/=k; yue[top]++; } } } if(c>1){ su[++top]=c; yue[top]=1; } for(R int i=1;i<=top;++i){ if((su[i]-1)%mod==0) ans*=(yue[i]+1); else ans=(ans*((ksm(su[i],yue[i]*b+1)-1)%mod)*(ksm(su[i]-1,mod-2)%mod))%mod; } printf("%lld",(ans+mod)%mod); return 0; }