HDU 1452 Happy 2004

题目意思：求2004^x的所有正因数的和对29求余

解析：

我们用s(x)表示x的因子和：

2的因子为1,2，s(2)=3;

3的因子为1,3，s(3)=4;

6的因子为1,2,3,6,s(6)=12;

可以发现：s(6)=s(2)*s(3)=3*4=12;

4的因子为1,2,4,，s(4)=7;

5的因子为1,5,s(5)=6;

20的因子为1,2,4,5,10,20，s(20)=42；

可以发现：s(20)=s(4)*s(5)=7*6=42;

s(25)=1+5+25=31

再看s(50)=1+2+5+10+25+50=93=3*31=s(2)*s(25)

这在数论中称为积性函数。

数论中积性函数的定义：对于正整数n的一个算术函数 f(n)，若f(1)=1，且当a,b互质时f(ab)=f(a)f(b)，在数论上就称它为积性函数。若对于某积性函数 f(n) ，就算a, b不互质，也有f(ab)=f(a)f(b)，则称它为完全积性的。

因为2004=2^2*3^1*167,因此s(2004^x)=s(2^(2x)*3^x*167^x)=s(2^(2x))*s(3^x)*s(167^x)

令a=s(2^(2x)),b=s(3^x),c=s(167^x),由于同余性质：167=22(mod29)，167可以用22代替，c=s(22^x).

同时，如果p是一个素数，则s(p^n)=p^0+p^1+P^2+......p^n，可以发现这是一个等比数列求和，因此s(p^n)=p^0(1-p^(n+1))/(1-p)=(p^(n+1)-1)/(p-1)......①

由①式得：a=s(2^(2x))=(2^(2x+1)-1)/1

b=s(3^x)=(2^(x+1)-1)/2

c=s(22^x)=(22^(x+1)-1)/21

而b与c不能直接计算，对于取模运算有如下规则：1. (a*b) %p= ( a%p) *(b%p)

2. (a/b) %p= ( a *b^(-1)%p)

b^(-1)(mod p)是b的逆元，2的逆元素为15,因为2*15%29=1%29,21的逆元素为18，因为21*18%29=1%29

在计算(a/b)%Mod时，往往需要先计算b%Mod的逆元p（b有逆元的条件是gcd(b,Mod)==1，显然素数肯定有逆元），然后由(a*p)%Mod得结果c。这里b的逆元p满足(b*p)%Mod=1

这里做简单的证明：(a/b)%Mod=c;    (b*p)%Mod=1;    ==》   (a/b)*(b*p) %Mod=c;    ==》    (a*p)%Mod=c;

求逆元可以用扩展欧几里得实现，假设已知b,mod,用扩展欧几里得求出bx+mody=1的一组解，两边同时%mod可得bx%mod=1%mod，因此x即为b%mod的逆元

综上所诉，题目中的a=(2^(2x+1)-1)%29

b=(2^(x+1)-1)*15%29

c=(22^(x+1)-1)/*18%29

ans=a*b*c%29

abc三项直接用快速幂求解即可，代码如下

#include <cstdio>
#include <algorithm>
#include <iostream>
#include <cstring>
#define LL long long
using namespace std;
LL pow_mod(LL a,LL n)
{
    LL ans=1;
    a=a%29;
    while(n)
    {
        if(n&1)ans=(ans*a)%29;
        a=a*a%29;
        n>>=1;
    }
    return ans;
}
int main()
{
    LL x,a,b,c;
    while(scanf("%I64d",&x)&&x)
    {
        a=(pow_mod(2,2*x+1)-1)%29;
        b=(pow_mod(3,x+1)-1)*15%29;
        c=(pow_mod(22,x+1)-1)*18%29;
        printf("%I64d\n",(a*b)%29*c%29);
    }
}