nylg 小M的因子和

  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\n");
116             continue;
117         }
118         else if(m==0)
119         {
120             printf("1\n");
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\n",hxl);
130     }
131     return 0;
132 }