题意: 求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 ",solve(A,B)); } return 0; }
也可以用二分来求 1+p+...+p^m
1.当m=2*k时:
p+...+p^2k=(1+p^k)(p+p^2+...+p^k)
2.当m=2*k+1时:
p+...+p^2k+p^(2k+1)=(1+p^k)(p+p^2+...+p^k)+p^(2k+1)
#include <iostream> #include <stdio.h> #include <string.h> #include <algorithm> using namespace std; const int maxn=50000010; int cnt,prilen; long long A,B; int shu[10010]; //存储质因数 int num[10010]; //存储相应的质因数个数 /** 用于处理n的质因数 */ void eular(int n) { memset(num,0,sizeof(num)); int i; cnt=0; for (i=2; i*i<=n; i++) { if (n%i==0) { cnt++; n/=i; shu[cnt]=i; num[cnt]++; while (n%i==0) { n/=i; num[cnt]++; } } } if(n>1){ cnt++; shu[cnt]=n; num[cnt]++; } } //这里传入参数的a和b值要定义为long long,否则答案不对 long long quickPow(long long a,long long b) { long long ans=1; while(b>0) { if(b&1) { ans=(ans*a)%9901; } a=(a*a)%9901; b=b>>1; } return ans; } /** 计算p+...+p^n; */ long long cal(long long p,long long n){ long long ans=1; if(n==0) return 0; if(n&1){ ans=((cal(p,n-1)%9901)+(quickPow(p,n)%9901))%9901; } else{ ans=(((1+quickPow(p,n/2))%9901)*(cal(p,n/2)%9901))%9901; } return ans%9901; } long long sum(){ long long ans=1; for(int i=1;i<=cnt;i++){ ans=(ans*(1+cal(shu[i],num[i]*B))%9901)%9901; } return ans%9901; } int main() { scanf("%lld%lld",&A,&B); eular(A); printf("%I64d ",sum()); return 0; }