首页 > 其他 > 详细

乘法逆元(转)

时间:2016-04-20 20:06:05      阅读:173      评论:0      收藏:0      [点我收藏+]

转自我敬爱的大神:ACdreamers(http://blog.csdn.net/acdreamers/article/details/8220787)

今天我们来探讨逆元在ACM-ICPC竞赛中的应用,逆元是一个很重要的概念,必须学会使用它。

 

对于正整数技术分享技术分享,如果有技术分享,那么把这个同余方程中技术分享的最小正整数解叫做技术分享技术分享的逆元。

 

逆元一般用扩展欧几里得算法来求得,如果技术分享为素数,那么还可以根据费马小定理得到逆元为技术分享

 

推导过程如下

                            技术分享

 

求现在来看一个逆元最常见问题,求如下表达式的值(已知技术分享

 

         技术分享  

 

当然这个经典的问题有很多方法,最常见的就是扩展欧几里得,如果技术分享是素数,还可以用费马小定理。

 

但是你会发现费马小定理和扩展欧几里得算法求逆元是有局限性的,它们都会要求技术分享技术分享互素。实际上我们还有一

种通用的求逆元方法,适合所有情况。公式如下

 

          技术分享

 

现在我们来证明它,已知技术分享,证明步骤如下

 

          技术分享

 

接下来来实战一下,看几个关于逆元的题目。

 

题目:http://poj.org/problem?id=1845

 

题意:给定两个正整数技术分享技术分享,求技术分享的所有因子和对9901取余后的值。

 

分析:很容易知道,先把技术分享分解得到技术分享,那么得到技术分享,那么技术分享

     的所有因子和的表达式如下

 

    技术分享

 

所以我们有两种做法。第一种做法是二分求等比数列之和。

 

代码:

  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4.   
  5. using namespace std;  
  6. typedef long long LL;  
  7. const int N = 10005;  
  8. const int MOD = 9901;  
  9.   
  10. bool prime[N];  
  11. int p[N];  
  12. int cnt;  
  13.   
  14. void isprime()  
  15. {  
  16.     cnt = 0;  
  17.     memset(prime,true,sizeof(prime));  
  18.     for(int i=2; i<N; i++)  
  19.     {  
  20.         if(prime[i])  
  21.         {  
  22.             p[cnt++] = i;  
  23.             for(int j=i+i; j<N; j+=i)  
  24.                 prime[j] = false;  
  25.         }  
  26.     }  
  27. }  
  28.   
  29. LL power(LL a,LL b)  
  30. {  
  31.     LL ans = 1;  
  32.     a %= MOD;  
  33.     while(b)  
  34.     {  
  35.         if(b & 1)  
  36.         {  
  37.             ans = ans * a % MOD;  
  38.             b--;  
  39.         }  
  40.         b >>= 1;  
  41.         a = a * a % MOD;  
  42.     }  
  43.     return ans;  
  44. }  
  45.   
  46. LL sum(LL a,LL n)  
  47. {  
  48.     if(n == 0) return 1;  
  49.     LL t = sum(a,(n-1)/2);  
  50.     if(n & 1)  
  51.     {  
  52.         LL cur = power(a,(n+1)/2);  
  53.         t = (t + t % MOD * cur % MOD) % MOD;  
  54.     }  
  55.     else  
  56.     {  
  57.         LL cur = power(a,(n+1)/2);  
  58.         t = (t + t % MOD * cur % MOD) % MOD;  
  59.         t = (t + power(a,n)) % MOD;  
  60.     }  
  61.     return t;  
  62. }  
  63.   
  64. void Solve(LL A,LL B)  
  65. {  
  66.     LL ans = 1;  
  67.     for(int i=0; p[i]*p[i] <= A; i++)  
  68.     {  
  69.         if(A % p[i] == 0)  
  70.         {  
  71.             int num = 0;  
  72.             while(A % p[i] == 0)  
  73.             {  
  74.                 num++;  
  75.                 A /= p[i];  
  76.             }  
  77.             ans *= sum(p[i],num*B) % MOD;  
  78.             ans %= MOD;  
  79.         }  
  80.     }  
  81.     if(A > 1)  
  82.     {  
  83.         ans *= sum(A,B) % MOD;  
  84.         ans %= MOD;  
  85.     }  
  86.     cout<<ans<<endl;  
  87. }  
  88.   
  89. int main()  
  90. {  
  91.     LL A,B;  
  92.     isprime();  
  93.     while(cin>>A>>B)  
  94.         Solve(A,B);  
  95.     return 0;  
  96. }  

 

第二种方法就是用等比数列求和公式,但是要用逆元。用如下公式即可

 

                     技术分享

 

因为技术分享可能会很大,超过int范围,所以在快速幂时要二分乘法。

 

代码:

  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4.   
  5. using namespace std;  
  6. typedef long long LL;  
  7. const int N = 10005;  
  8. const int MOD = 9901;  
  9.   
  10. bool prime[N];  
  11. int p[N];  
  12. int cnt;  
  13.   
  14. void isprime()  
  15. {  
  16.     cnt = 0;  
  17.     memset(prime,true,sizeof(prime));  
  18.     for(int i=2; i<N; i++)  
  19.     {  
  20.         if(prime[i])  
  21.         {  
  22.             p[cnt++] = i;  
  23.             for(int j=i+i; j<N; j+=i)  
  24.                 prime[j] = false;  
  25.         }  
  26.     }  
  27. }  
  28.   
  29. LL multi(LL a,LL b,LL m)  
  30. {  
  31.     LL ans = 0;  
  32.     a %= m;  
  33.     while(b)  
  34.     {  
  35.         if(b & 1)  
  36.         {  
  37.             ans = (ans + a) % m;  
  38.             b--;  
  39.         }  
  40.         b >>= 1;  
  41.         a = (a + a) % m;  
  42.     }  
  43.     return ans;  
  44. }  
  45.   
  46. LL quick_mod(LL a,LL b,LL m)  
  47. {  
  48.     LL ans = 1;  
  49.     a %= m;  
  50.     while(b)  
  51.     {  
  52.         if(b & 1)  
  53.         {  
  54.             ans = multi(ans,a,m);  
  55.             b--;  
  56.         }  
  57.         b >>= 1;  
  58.         a = multi(a,a,m);  
  59.     }  
  60.     return ans;  
  61. }  
  62.   
  63. void Solve(LL A,LL B)  
  64. {  
  65.     LL ans = 1;  
  66.     for(int i=0; p[i]*p[i] <= A; i++)  
  67.     {  
  68.         if(A % p[i] == 0)  
  69.         {  
  70.             int num = 0;  
  71.             while(A % p[i] == 0)  
  72.             {  
  73.                 num++;  
  74.                 A /= p[i];  
  75.             }  
  76.             LL M = (p[i] - 1) * MOD;  
  77.             ans *= (quick_mod(p[i],num*B+1,M) + M - 1) / (p[i] - 1);  
  78.             ans %= MOD;  
  79.         }  
  80.     }  
  81.     if(A > 1)  
  82.     {  
  83.         LL M = MOD * (A - 1);  
  84.         ans *= (quick_mod(A,B+1,M) + M - 1) / (A - 1);  
  85.         ans %= MOD;  
  86.     }  
  87.     cout<<ans<<endl;  
  88. }  
  89.   
  90. int main()  
  91. {  
  92.     LL A,B;  
  93.     isprime();  
  94.     while(cin>>A>>B)  
  95.         Solve(A,B);  
  96.     return 0;  
  97. }  


 

其实有些题需要用到技术分享技术分享的所有逆元,这里技术分享为奇质数。那么如果用快速幂求时间复杂度为技术分享

如果对于一个1000000级别的素数技术分享,这样做的时间复杂度是很高了。实际上有技术分享的算法,有一个递推式如下

 

                   技术分享

 

它的推导过程如下,设技术分享,那么

 

       技术分享

 

对上式两边同时除技术分享,进一步得到

 

       技术分享

 

再把技术分享技术分享替换掉,最终得到

 

       技术分享

 

初始化技术分享,这样就可以通过递推法求出技术分享模奇素数技术分享的所有逆元了。

 

另外技术分享技术分享的所有逆元值对应技术分享中所有的数,比如技术分享,那么技术分享对应的逆元是技术分享

 

 

题目:http://www.lydsy.com/JudgeOnline/problem.php?id=2186

 

题意:技术分享技术分享互质的数的个数,其中技术分享

 

分析:因为技术分享,所以技术分享,我们很容易知道如下结论

     对于两个正整数技术分享技术分享,如果技术分享技术分享的倍数,那么技术分享中与技术分享互素的数的个数为技术分享

 

     本结论是很好证明的,因为技术分享中与技术分享互素的个数为技术分享,又知道技术分享,所以

     结论成立。那么对于本题,答案就是

 

     技术分享

 

      其中技术分享为小于等于技术分享的所有素数,先筛选出来即可。由于最终答案对一个质数取模,所以要用逆元,这里

      求逆元就有技巧了,用刚刚介绍的递推法预处理,否则会TLE的。

 

代码:

  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4. #include <bitset>  
  5.   
  6. using namespace std;  
  7. typedef long long LL;  
  8. const int N = 10000005;  
  9.   
  10. bitset<N> prime;  
  11.   
  12. void isprime()  
  13. {  
  14.     prime.set();  
  15.     for(int i=2; i<N; i++)  
  16.     {  
  17.         if(prime[i])  
  18.         {  
  19.             for(int j=i+i; j<N; j+=i)  
  20.                 prime[j] = false;  
  21.         }  
  22.     }  
  23. }  
  24.   
  25. LL ans1[N],ans2[N];  
  26. LL inv[N];  
  27.   
  28. int main()  
  29. {  
  30.     isprime();  
  31.     int MOD,m,n,T;  
  32.     scanf("%d%d",&T,&MOD);  
  33.     ans1[0] = 1;  
  34.     for(int i=1; i<N; i++)  
  35.         ans1[i] = ans1[i-1] * i % MOD;  
  36.     inv[1] = 1;  
  37.     for(int i=2;i<N;i++)  
  38.     {  
  39.         if(i >= MOD) break;  
  40.         inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;  
  41.     }  
  42.     ans2[1] = 1;  
  43.     for(int i=2; i<N; i++)  
  44.     {  
  45.         if(prime[i])  
  46.         {  
  47.             ans2[i] = ans2[i-1] * (i - 1) % MOD;  
  48.             ans2[i] = ans2[i] * inv[i % MOD] % MOD;  
  49.         }  
  50.         else  
  51.         {  
  52.             ans2[i] = ans2[i-1];  
  53.         }  
  54.     }  
  55.     while(T--)  
  56.     {  
  57.         scanf("%d%d",&n,&m);  
  58.         LL ans = ans1[n] * ans2[m] % MOD;  
  59.         printf("%lld\n",ans);  
  60.     }  
  61.     return 0;  
  62. }  


 

接下来还有一个关于逆元的有意思的题目,描述如下

 

     技术分享

 

证明:

 

     技术分享

 

     其中技术分享

 

     所以只需要证明技术分享,而我们知道技术分享技术分享的逆元对应全部

     技术分享中的所有数,既是单射也是满射。

 

     所以进一步得到

 

      技术分享

 

      证明完毕!

乘法逆元(转)

原文:http://www.cnblogs.com/topW2W/p/5413879.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!