hdu 1452 Happy 2004(数论(模P乘法逆元+快速幂取模))

题目：http://acm.hdu.edu.cn/showproblem.php?pid=1452
大意：给出整数n，求出2004^n的约数和模29的结果
2004=2^2*3*167。
f[1]=2^2,f[2]=3,f[3]=167三者两两互质。2004所有因子的情况：(1+2+2^2)*(1+3)*(1+167)。所以2004的因子个数和就是3*2*2=12。并且，因为三者两两互质，所以因子和就是s=s1*s2*s3=（1+2+4）*（1+3）*（1+167）=1+2+3+4+6+12+167+334+501+668+1002+2004=4704。这也就是积性函数的性质：gcd(a,b)==1 && s(a,b)==s(a)*s(b)满足这种条件的s叫做积性函数。s(2004^n)=s(2^(2n)*3^n*167^n)=s(2^（2n）)*s(3^n)*s(167^n)。【互质的数同次方（非0）也是互质的】。
si=s(a[i]^n)=1+a[i]+a[i]^2+……+a[i]^n=a1*(1-q^(n+1))/(1-q)=(1-a[i]^(n+1))/(1-a[i])=(a[i]^(n+1)-1)/(a[i]-1)。
那么s(2^(2n))=1+2+2^2+……+2^2n=(2^(2n+1)-1)/(2-1)=2^(2n+1)-1.
s(3^n)=(3^(n+1)-1)/2
s(167^n)=(167^(n+1)-1)/166 因为167%29=22. 所以s(167^n)=（22^(n+1)-1)/22.
相关知识点：
"除以一个数模q等于乘以一个数模q的逆元"。这和"除以一个数等于乘以一个数的倒数"及其相似。
例如：
4关于1模7的乘法逆元为多少？
4X≡1 mod 7
这个方程等价于求一个X和K，满足
4X=7K+1
其中X和K都是整数。
若ax≡1 mod f, 则称a关于1模f的乘法逆元为x。也可表示为ax≡1(mod f)。
当a与f互素时，a关于模f的乘法逆元有唯一解。如果不互素，则无解。如果f为素数，则从1到f-1的任意数都与f互素，即在1到f-1之间都恰好有一个关于模f的乘法逆元。
例如，求5关于模14的乘法逆元：【关于1】
14=5*2+4
5=4+1
说明5与14互素，存在5关于14的乘法逆元。
1=5-4=5-(14-5*2)=5*3-14
因此，5关于模14的乘法逆元为3。【关于1】 即：5X≡1 mod 14 --》5X+14Y=1 -->X=3，Y=-1。
程序：
欧几里德算法求乘法逆元：
Extended Euclid (d，f) //算法求d关于模f的乘法逆元d-1 ，即 d* d-1 mod f = 1
1 (X1，X2，X3) := (1，0，f)； (Y1，Y2，Y3) := (0，1，d)
2 if (Y3=0) then return d-1 = null //无逆元
3 if (Y3=1) then return d-1 = Y2 //Y2为逆元
4 Q := X3 div Y3 //整除
5 (T1，T2，T3) := (X1 - Q*Y1，X2 - Q*Y2，X3 - Q*Y3)
6 (X1，X2，X3) := (Y1，Y2，Y3)
7 (Y1，Y2，Y3) := (T1，T2，T3)
8 goto 2
证明：

详见代码例子：

#include<cstdio>#include<iostream>using namespace std;//ax+by=1int Extend_Eulid(int b,int a)//(b,a)~(x,mod)~(b,a){ //初始化：X=(1,0,a); Y=(0,1,b) Y2对应y，也是最终的结果。    int x1,x2,x3,y1,y2,y3 ;    x1=1,x2=0,x3=a,y1=0,y2=1,y3=b ;    printf("%d*%d+%d*%d=%d\n",x1,a,x2,b,x3);    printf("%d*%d+%d*%d=%d\n",y1,a,y2,b,y3);    while(y3 && y3!=1)    {        int q=x3/y3 ;        int t1,t2,t3 ;        cout<<q<<endl;        t1=x1-q*y1,t2=x2-q*y2,t3=x3-q*y3 ;        x1=y1,x2=y2,x3=y3 ;        y1=t1,y2=t2,y3=t3 ;        printf("%d*%d+%d*%d=%d\n",y1,a,y2,b,y3);    }    if(!y3)return -1 ;    return y2 ;}/*两个式子不断相减，确保ax+by=c的右值c不断变小直至趋近1.最后得到的y2即是b关于1模a的乘法逆元。*/int main(){    int x,mod;    while(cin>>x>>mod){        int ans=Extend_Eulid(x,mod);        printf("%d关于1模%d的乘法逆元是%d\n",x,mod,ans);    }    return 0;}

乘法逆元模板：

int Extend_Eulid(int b,int a)
{ 
    int x1,x2,x3,y1,y2,y3 ;
    x1=1,x2=0,x3=a,y1=0,y2=1,y3=b ;
    while(y3 && y3!=1)
    {
        int q=x3/y3 ;
        int t1,t2,t3 ;
        t1=x1-q*y1,t2=x2-q*y2,t3=x3-q*y3 ;
        x1=y1,x2=y2,x3=y3 ;
        y1=t1,y2=t2,y3=t3 ;
    }
    if(!y3)return -1 ;
    return y2 ;
}
另外，还需要快速幂取模：
int supower(int x,int p){
    int ans=1;
    while(p){
        if(p&1)ans=ans*x%29;
        x=x*x%29;
        p=p>>1;
    }
    return ans;
}
代码：

#include <iostream>#include<cstdio>#include<cstring>using namespace std;int Extend_Eulid(int d,int f){    int x1,x2,x3,y1,y2,y3 ;    x1=1,x2=0,x3=f,y1=0,y2=1,y3=d ;    while(y3 && y3!=1)    {        int q=x3/y3 ;        int t1,t2,t3 ;        t1=x1-q*y1,t2=x2-q*y2,t3=x3-q*y3 ;        x1=y1,x2=y2,x3=y3 ;        y1=t1,y2=t2,y3=t3 ;    }    if(!y3)return -1 ;    return y2 ;}int supower(int x,int p){    int ans=1;    while(p){        if(p&1)ans=ans*x%29;        x=x*x%29;        p=p>>1;    }    return ans;}/*a=s(2^(2n))=(2^(2n+1)-1)b=s(3^n)=(3^(n+1)-1)/2c=s(22^n)=(22^(n+1)-1)/21*/int main(){    int qb=Extend_Eulid(2,29),qc=Extend_Eulid(21,29);    //freopen("cin.txt","r",stdin);    int n;    while(cin>>n&&n){        int a=supower(2,2*n+1)-1;        int b=(supower(3,n+1)-1)*qb%29;        int c=(supower(22,n+1)-1)*qc%29;        printf("%d\n",a*b*c%29);    }    return 0;}