小M的因子和
时间限制：1000 ms | 内存限制：65535 KB
难度：2
- 描述
-
小M在上课时有些得意忘形，老师想出道题目难住他。小M听说是求因子和，还是非常得意，但是看完题目是求A的B次方的因子和，有些手足无措了，你能解决这个问题吗？
- 输入
- 有多组测试样例
每行两个数 A ,B ，（1≤A,B≤10^9） - 输出
- 输出A的B次方的因子和，并对9901取余。
- 样例输入
-
2 3
- 样例输出
-
15
- 上传者
Sumdiv
| Time Limit: 1000MS | Memory Limit: 30000K | |
| Total Submissions: 13547 | Accepted: 3321 |
Description
Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).
Input
The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.
Output
The only line of the output will contain S modulo 9901.
Sample Input
2 3
Sample Output
15
Hint
2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).
Source
都是一样的。
1 /** 2 这道题，类似于hdu happy 2004. 3 这题是一个通法。 4 A^B%P, 5 拆分A变成素数，因为素数满足因子之和 s(x*y)=s(x)*s(y); 6 更重要的是，素数的x^n的因子个数是可以求出来的。就是 7 和(x^n) = 1 + x + x^2 +x^3 +......+x^n; 8 9 这样的话，我们就可以轻松解决这样的一个问题了。 10 提供两种思路。 11 1 + x + x^2 +x^3 +......+x^n，直接求它对%p的值。运用快速幂也可以的。 12 这就是第一种方法，也是下面的ac方法。 13 14 第二种方法：1 + x + x^2 +x^3 +......+x^n= 等比数列前n+1和。 15 很据 S(p^X)=1+p+p^2+...+p^X = (p^(X+1)-1)/(p-1)； 16 这样就等于求这个式子了。好的，怎么求呢？ 17 18 p^(X+1)-1 这个应该没有问题，快速幂取模 19 关键是1/(p-1); 这个不能直接取模。转化为求乘法的逆元。 20 21 乘法的逆元？？恩。 22 一开始，我就是这样做的，后来想用费马小定理，联想到了一道题C(n,m)的求法时候 23 也出现过 n!/(m!*(n-m)!) 对于费马小定理 a%p == a^p-1%p; 那么这样的话，我就能 24 转化一下，对于 1/(p-1) ,转化为 (p-1)^-1 ==> (p-1)^-1 % mod = (p-1)^mod-2 %mod; 25 26 好像这样是对的，是的。 27 费马小定理的前提是什么？ mod是一个素数，这个满足了。 28 还有一个条件gcd(mod,p-1)==1 这个就不一定了.当p为 mod的倍数+1而且是素数的时候。 29 就很感慨的发现，p-1就是mod的倍数。 30 那么费马小定理的路，就不好走了。 31 32 那我用扩展欧几里得的算法来求逆元。我看到很多人的解题思路可能都是这个吧。 33 其实，我依然有一个疑问。 34 对于(p^(X+1)-1)/(p-1)，显然我能对其转化 p%Euler(mod) == t 35 ==> (t ^(x+1)-1)/(t-1); 但是如果 t ^(x+1)%p 为1的时候，这个值就为0了。 36 37 38 例子 A B P 39 59407 1 9901 40 **/ 41 42 #include<iostream> 43 #include<stdio.h> 44 #include<cstring> 45 #include<cstdlib> 46 #include<vector> 47 using namespace std; 48 typedef long long LL; 49 50 const LL p = 9901; 51 LL prime[1000],len; 52 LL num[1000]; 53 LL dp[100],dlen; 54 void Euler(LL n) 55 { 56 LL i,k; 57 len=0; 58 for(i=2; i*i<=n;i++) 59 { 60 if(n%i==0){ 61 k=0; 62 while(n%i==0){ 63 n=n/i; 64 k++; 65 } 66 prime[++len]=i%p; 67 num[len]=k; 68 } 69 } 70 if(n!=1){ 71 prime[++len]=n%p; 72 num[len]=1; 73 } 74 } 75 LL sum_mod(LL a,LL n) 76 { 77 LL ans=0; 78 n=n%p; 79 while(n) 80 { 81 if(n&1) ans=(ans+a)%p; 82 n=n>>1; 83 a=(a+a)%p; 84 } 85 return ans; 86 } 87 LL solve(LL a,LL n) 88 { 89 LL p1=a,p2=a,ans,i; 90 dlen=0; 91 while(n) 92 { 93 dp[++dlen]=(n&1); 94 n=n>>1; 95 } 96 ans=dlen-1; 97 for(i=ans;i>=1;i--) 98 { 99 p1=sum_mod(p1,p2+1); 100 p2=sum_mod(p2,p2); 101 if(dp[i]==1) 102 { 103 p2=sum_mod(p2,a); 104 p1=(p1+p2)%p; 105 } 106 } 107 return (p1+1)%p; 108 } 109 int main() 110 { 111 LL n,m,i; 112 while(scanf("%lld%lld",&n,&m)>0) 113 { 114 if(n==0){ 115 printf("0 "); 116 continue; 117 } 118 else if(m==0) 119 { 120 printf("1 "); 121 continue; 122 } 123 Euler(n); 124 LL hxl=1; 125 for(i=1;i<=len;i++) 126 { 127 hxl=(hxl*solve(prime[i],num[i]*m))%p; 128 } 129 printf("%lld ",hxl); 130 } 131 return 0; 132 }