POJ1845:Sumdiv

点击打开题目链接

Sumdiv

Time Limit: 1000MS Memory Limit: 30000KTotal Submissions: 11902 Accepted: 2859

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

Source

Romania OI 2002

=====================================题目大意=====================================

计算A^B的所有正因子之和取模9901。

=====================================算法分析=====================================

一、根据质因数分解定理：

    每个大于1的整数都能分解成质因数乘积的形式，并且若将质因数由小到大排列在一起，相同因数的积写成幂的形式，则

   这种分解方法是唯一的。

    比如整数90可表示为：90 = (2^1) * (3^2) * (5^1)。

    从而得知A可表示为：A = (prim1^n1) * (prim2^n2) * ...... * （primK^nK)。

    进而得知A^B可表示为：A^B = (prim1^(n1*B)) * (prim2^(n2*B)) * ...... * (primK^(nK*B))。                      

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

二、根据约数和公式：

    

    一个数的所有约数之和等于每个质因数从0次幂加到其最高次幂之和的乘积。

    比如整数90的约数和为：( 2^0 + 2^1 ) * ( 3^0 + 3^1 + 3^2 ) * ( 5^0 + 5^1 )。

    从而得知A^B的所有约数之和为：Sum = ( 1 + Prim1^1 + Prim1^2 + ...... + Prim1^(n1*B) ) 

                                     * ( 1 + Prim2^1 + Prim2^2 + ...... + Prim2^(n2*B) )

                                     ......

                                     * ( 1 + PrimK^1 + PrimK^2 + ...... + PrimK^(nK*B) )。

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

     

三、根据同余模公式：

    

   ( a + b ) % c = ( ( a % c ) + ( b % c ) ) % c。

   ( a - b ) % c = ( ( a % c ) - ( b % c ) ) % c。

   ( a * b ) % c = ( ( a % c ) * ( b % c ) ) % c。

    可以推广得知：如果将一个整数分解为若干个“加法/减法/乘法”因子，则其取模等于所有“加法/减法/乘法”取模后

   的“和，差，积”。

    比如本题所求即为：Sum % MOD = ( ( ( 1 + Prim1^1 + Prim1^2 + ...... + Prim1^(n1*B) ) % MOD )

                                  * ( ( 1 + Prim2^1 + Prim2^2 + ..,,,. + Prim2^(n2*B) ) % MOD )

                                  ......

                                  * ( ( 1 + PrimK^1 + PrimK^2 + ...... + PrimK^(nK*B) ) % MOD ) ) % MOD。

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

四、分析得出本题主要问题：

    求一个首项为1的等比数列的和模9901（以下令Q=等比数列公比，N=等比数列项数，S=等比数列的和，MOD=9901）。

     

    1、朴素算法：逐项相加

        

       在本题的数据大小下如此纯暴力的算法必然会超时！

         

    1、改进算法：提取公因式

       对一个多项式提取公因式（当然是说化简的那种）是减少计算量的一个有效途径，并且被提取公因式的项越多减少的

      计算量也就越多。

       而对于上述的等比数列显然不存在从所有项提取公因式使得计算量减少的方法。 

       那么就应该想到组合等比数列的项而后从组合项中提取公因式，首先自然是考虑两两组合了，而这显然是需要对N进行奇

       偶性讨论的：

       a、当N为偶数时，令Mid=N/2则有：S = 1 + Q^1 + Q^2 + ...... + Q^(N-1)

                                        = ( 1 + Q^Mid ) + ( Q^1 + Q^(Mid+1) ) + ...... + ( Q^(Mid-1) + Q^(Mid+Mid-1) )

                                        = ( 1 + Q^Mid )* ( 1 + Q^1 + Q^2 + ...... + Q^(Mid-1))。

          通过组合等比数列的第K项与第K+Mid-1项并提取公因式(1+Q^Mid)后，等比数列的和S分解成为了黄色因子与蓝色因子的乘积。

          对于黄色因子，只需求一个整数幂即可，用快速幂即可高效解决。

          至于蓝色因子，它不还是一个首项为1的等比数列么（而且还是原等比数列的前半部分）？那么它可以通过此算法递归解决。

       b、当N为奇数时，用以上算法求其N-1项的和，用快速幂求其末项，最后相加不就行了？

       c、综上可得S(N)的计算方法为：S(N) =( N % 2 == 0) ? ( S(N/2) * ( Q^(N/2) + 1 ) ) : ( S(N-1) + Q^(N-1) )。

       d、关于所求为S(N)%MOD这一问题请参考第三点的“推广”解决。

    2、朴素算法：等比数列求和公式（以下令X=求和公式分子Q^N-1，Y=求和公式分母Q-1）

 

       由于只能求得X模MOD的结果而不可能求得X的精确值（因为X可能极大)

       而且不存在同余模公式( X / Y ) % MOD = ( ( X % MOD ) / ( Y % MOD ) ) % MOD。

       因此直接使用等比数列求和公式不可行。

   

    2、改进算法：乘法逆元

       无法直接使用等比数列求和公式的原因是无法拆分除法取模(X/Y)%MOD，而乘法逆元正是解决除法取模的有效方法之一(其

       他还有没有我不知道)。

       假设可以找到满足等式(Y*Z)%MOD=1的Z。

       那么就可以得到 ( X / Y ) % MOD = ( ( X / Y ) % MOD ) % MOD

                                      = ( ( ( X / Y ) % MOD ) * 1 ) % MOD

                                      = ( ( ( X / Y ) % MOD ) * ( ( Y * Z ) % MOD ) ) % MOD

                                      = ( ( X / Y ) * ( Y * Z ) ) % MOD

                                      = ( X * Z ) % MOD

                                      = ( ( X % MOD ) * ( Z * MOD ) ) % MOD

       这里的X%MOD（X是求和公式分子Q^N-1）可以通过快速幂取模求得，那么当Z已知时，问题就迎刃而解了！

       而在数学上，称满足等式(Y*Z)%MOD=1的Y与Z即互为取模MOD的乘法逆元。

       重要性质：当Y与MOD互素时Y关于取模MOD的乘法逆元有唯一解，否则无解。

                 特别的，当MOD为素数时，由于1到MOD-1的任何数都与MOD互素，因此1到MOD-1之间的任何数都存在唯一的关于取模

                 MOD的乘法逆元。

       求解算法：扩展欧几里德公式。

       相关资料：百度百科——乘法逆元 

               ChinaUnix博客——乘法逆元的作用

               网易博客——扩展欧几里德算法求乘法逆元

       注意事项：在求Y关于模MOD的乘法逆元时要求Y与MOD必须互素，而本题的MOD=9901为素数，那么应该确保Y介于1至MOD-1之间。

                 为此可让等比数列公比模MOD（根据第三点的“推广”可知这样并不会影响等比数列的和模MOD的结果），此时Q便

                 介于0与MOD-1之间，而若Q为0(不构成等比数列)应直接返回0，若Q为1（公比为1不可应用求和公式）应直接返回N。

                 最后Q便介于2与MOD-1之间，那么Y就介于1与MOD-2之间，根据乘法逆元的性质可知Y关于模MOD的乘法逆元此时有

                 唯一解。

=======================================代码=======================================

一、提取公因式

#include<stdio.h> const int MOD=9901;int A,B; //快速幂取模MOD:(X^N)%MODint QuickPowerMOD(int X,int N)    {int Base=X,Res=1;while(N){if(N&1) { Res=(Res*Base)%MOD; }Base=(Base*Base)%MOD;N>>=1;}return Res;}//扩展欧几里德算法求解Num关于模MOD的乘法逆元:(?*Num)%MOD=1int ExtendedEuclid(int Num)        {int X[3]={1,0,MOD};int Y[3]={0,1,Num};while(1){if(Y[2]==0) {    return 0; }if(Y[2]==1) { return Y[1]; }int Quot=X[2]/Y[2];int Tmp0=X[0]-Quot*Y[0];  X[0]=Y[0];  Y[0]=Tmp0;  int Tmp1=X[1]-Quot*Y[1];  X[1]=Y[1];  Y[1]=Tmp1; int Tmp2=X[2]-Quot*Y[2];  X[2]=Y[2];  Y[2]=Tmp2;}}//首项为1项数为N公比为Q的等比数列的和取模MODint GeomeSeriesSumMOD(int Q,int N){//公比首先取模MODQ%=MOD;//公比取模MOD为0的特殊处理if(Q==0) { return 1; }//公比取模MOD为1的特殊处理if(Q==1) { return N; }//等比数列求和公式分子取模MODint NumeratorMOD=QuickPowerMOD(Q,N)-1;    //扩展欧几里得算法求解(Q-1)关于取模MOD的乘法逆元int InverseElement=ExtendedEuclid(Q-1);//对得出的为负数的InverseElement需处理转化为正数while(InverseElement<0) { InverseElement+=MOD; }//返回计算结果return (NumeratorMOD*InverseElement)%MOD;}int main(){while(scanf("%I64d%I64d",&A,&B)==2){ if(A==0)        { printf("0\n"); continue; }    if(A==1||B==0)  { printf("1\n"); continue; }int SumMOD=1,Prim=2;while(Prim*Prim<=A){int Cnt=0;while(A%Prim==0) { ++Cnt;  A/=Prim; }if(Cnt!=0) { SumMOD=(SumMOD*GeomeSeriesSumMOD(Prim,Cnt*B+1))%MOD; }Prim=(Prim==2?Prim+1:Prim+2);}//特殊处理当A自身为素数时的情况if(A!=1) { SumMOD=(SumMOD*GeomeSeriesSumMOD(A,B+1))%MOD; }printf("%d\n",SumMOD);}return 0;}

二、乘法逆元

#include<stdio.h> const int MOD=9901;int A,B; //快速幂取模MOD:(X^N)%MODint QuickPowerMOD(int X,int N)    {//令X取模MOD可确保之后的计算数据不超出Int32的表示范围int Base=X%MOD,Res=1;while(N){if(N&1) { Res=(Res*Base)%MOD; }Base=(Base*Base)%MOD;N>>=1;}return Res;}//首项为1项数为N公比为Q的等比数列的和取模MODint GeomeSeriesSumMOD(int Q,int N){if(N==1){return 1;}if(N&1) {int Factor1=GeomeSeriesSumMOD(Q,N-1);int Factor2=QuickPowerMOD(Q,N-1);return (Factor1+Factor2)%MOD;}else{int Factor1=GeomeSeriesSumMOD(Q,N/2);int Factor2=QuickPowerMOD(Q,N/2)+1;return (Factor1*Factor2)%MOD; }}int main(){while(scanf("%I64d%I64d",&A,&B)==2){ if(A==0)        { printf("0\n"); continue; }    if(A==1||B==0)  { printf("1\n"); continue; }int SumMOD=1,Prim=2;while(Prim*Prim<=A){int Cnt=0;while(A%Prim==0) { ++Cnt;  A/=Prim; }if(Cnt!=0) { SumMOD=(SumMOD*GeomeSeriesSumMOD(Prim,Cnt*B+1))%MOD; }Prim=(Prim==2?Prim+1:Prim+2);}//特殊处理当A自身为素数时的情况if(A!=1) { SumMOD=(SumMOD*GeomeSeriesSumMOD(A,B+1))%MOD; }printf("%d\n",SumMOD);}return 0;}