码迷,mamicode.com
首页 > 其他好文 > 详细

Sumdiv POJ - 1845

时间:2018-12-17 23:57:23      阅读:288      评论:0      收藏:0      [点我收藏+]

标签:不用   span   偶数   put   const   lap   lan   sub   约数和   

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.OutputThe 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).

 

题意:求AB的所有约数的和  % MOD (9901)    题意中有点问题,我们知道0是没有约数的,我觉得A、B应该都是>0的

思路:我们可以把A分解质因数(p1c1  *  p2c2  *  .... * pncnB

约数和:(1 + p1 + p12 + ... + p1B*c1)* (1 + p2 + p22 + ... + p2B*c2)* .... * (1 + pn + pn2 + ... + pnB*cn)    ( 排列组合问题)

这样我们可以看出这是多个等比数列乘积,可以用等比数列求和公式  (a1 *(1-qn))/(1-q),我们注意到这里有除法,但是同余模定理是对于加减乘的,那么我们可以利用费马小定理,

求出(1-q)的逆元,然后把除变成乘逆元

坑点:应为 9901 这个质数较小,很容易找到一个数x,(x-1)% MOD == 0 ,就说明这个数是没有逆元的(例217823),那么对于这种情况,我们不能用逆元算,你会发现这种情况下,

pn % MOD == 1  ((p-1)% MOD == 0)),这样(1 + pn + pn2 + ... + pnB*cn) == (1 % MOD + pn %MOD + pn2 %MOD + ... + pnB*cn %MOD) == B*cn+1

技术分享图片
 1 #include<iostream>
 2 #include<cstdio>
 3 #include<math.h>
 4 using namespace std;
 5 
 6 const int maxn = 1e4;
 7 const int mod = 9901;
 8 int a,b;
 9 int p[maxn];
10 int c[maxn];
11 int calc(int x)
12 {
13     int m= 0;
14     int up = sqrt(x);
15     for(int i=2;i<=up;i++)
16     {
17         if(x % i == 0)p[++m] = i,c[m] = 0;
18         while(x % i == 0)x/= i,c[m]++;
19     }
20     if(x > 1)p[++m] = x,c[m] = 1;
21     return m;
22 }
23 
24 typedef long long ll;
25 
26 ll qpow(ll a,ll b)
27 {
28     ll ans = 1;
29     ll base = a;
30     while(b)
31     {
32         if(b&1)ans = (ans * base)%mod;
33         base = (base * base)%mod;
34         b >>= 1;
35     }
36     return ans;
37 }
38 int main()
39 {
40     scanf("%d%d",&a,&b);
41     int n = calc(a);
42     ll ans = 1;
43     for(int i=1;i<=n;i++)
44     {
45         if((1-p[i])%mod == 0)
46         {
47             ans = (ans * (b * c[i]+ 1))%mod;
48             continue;
49         }
50         ll Ni = qpow(1-p[i],mod-2);
51         ll tmp = 1-qpow(p[i],c[i]*b+1);
52         ans = (ans * (tmp*Ni%mod+mod)%mod)%mod;
53     }
54     printf("%lld\n",ans);
55 }
View Code

 

还有一种写法,就是不用公式计算等比数列和,这样就避免了逆元的问题

sum(p,c) = (1 + p + p2 + ... + pk

(1)c为奇数,sum(p,k)= sum(p,(k-1 )/2)*(1+p(k+1)/2)   

sum(p,c) = (1 + p + p2 + ... + p(k-1)/2)+ (p(k+1)/2 + ... + pk)           (c为奇数,加上0次幂,变成偶数,刚好可以分成两个等长的数列)

(2)c为偶数,sum(p,k)= sum(p,k/2-1)*(1+pk/2)+ p

sum(p,c) = (1 + p + p2 + ... + pk/2-1)+ (pk/2 + ... + pk-1)+ pk                         (c+1是奇数)

技术分享图片
 1 #include<iostream>
 2 #include<cstdio>
 3 #include<math.h>
 4 using namespace std;
 5 
 6 const int maxn = 1e4;
 7 const int mod = 9901;
 8 int a,b;
 9 int p[maxn];
10 int c[maxn];
11 int calc(int x)
12 {
13     int m= 0;
14     for(int i=2;i*i<=x;i++)
15     {
16         if(x % i == 0)p[++m] = i,c[m] = 0;
17         while(x % i == 0)x/= i,c[m]++;
18     }
19     if(x > 1)p[++m] = x,c[m] = 1;
20     return m;
21 }
22 
23 typedef long long ll;
24 ll qpow(ll a,ll b)
25 {
26     ll ans = 1;
27     ll base = a;
28     while(b)
29     {
30         if(b&1)ans = (ans * base)%mod;
31         base = (base * base)%mod;
32         b >>= 1;
33     }
34     return ans;
35 }
36 ll sum(ll p,ll c)
37 {
38     if(c == 0)return 1;
39     if(c&1)return ((1+qpow(p,(c+1)/2))%mod*(sum(p,(c-1)/2)%mod))%mod;
40     else return ((1+qpow(p,c/2))%mod*(sum(p,c/2-1))%mod+qpow(p,c))%mod;
41 }
42 
43 int main()
44 {
45     scanf("%d%d",&a,&b);
46     int n = calc(a);
47     ll ans = 1;
48     for(int i=1;i<=n;i++)
49     {
50         ans = (ans * sum(p[i],c[i]*b))%mod;
51     }
52     printf("%lld\n",ans);
53 }
View Code

 

Sumdiv POJ - 1845

标签:不用   span   偶数   put   const   lap   lan   sub   约数和   

原文地址:https://www.cnblogs.com/iwannabe/p/10134388.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!