HDU1852 Beijing 2008（快速幂+特殊公式）

As we all know, the next Olympic Games will be held in Beijing in 2008. So the year 2008 seems a little special somehow. You are looking forward to it, too, aren't you? Unfortunately there still are months to go. Take it easy. Luckily you meet me. I have a problem for you to solve. Enjoy your time. 

Now given a positive integer N, get the sum S of all positive integer divisors of 2008 ^N. Oh no, the result may be much larger than you can think. But it is OK to determine the rest of the division of S by K. The result is kept as M. 

Pay attention! M is not the answer we want. If you can get 2008 ^M, that will be wonderful. If it is larger than K, leave it modulo K to the output. See the example for N = 1，K = 10000: The positive integer divisors of 20081 are 1、2、4、8、251、502、1004、2008，S = 3780, M = 3780, 2008 ^M % K = 5776. 

InputThe input consists of several test cases. Each test case contains a line with two integers N and K (1 ≤ N ≤ 10000000, 500 ≤ K ≤ 10000). N = K = 0 ends the input file and should not be processed. 
OutputFor each test case, in a separate line, please output the result. 
Sample Input

1  10000
0  0

Sample Output

5776

题意：求M=2008^n的所有因子和%K；输出2008^M%K。

思路：后面部分就是普通快速幂了；问题关键在于解决前面部分。

         因子和：由唯一分解，2008=(2^3)*(251^1）。所以2008^n=(2^3n)*(251^n)，则因子有（3n+1）*(n+1)个；sum=(1+2+4+...2^3n)*(1+251+..+251*n)。

                        对其合并,得ans=(2^(3n+1)-1)*(251^(n+1)-1)/250。（等比数列）

         解决分母250：ans%K，而ans有分母250，且gcd（250，K）不一定就为1，即不一定互素，所以不能用逆元解决。

         公式：(a/b)%mod=a%(b*mod)/b%mod 

（emmm，注意代码里两个等比数列求和都是取余250*k，尽管前面一个没有分母250.）

对比：这个公式使用于分母比较小的情况，而且不要求b,Mod互素，条件是a|b。而一般求组合数(n!)/(m!*(n-m)!)%Mod，由于分母较大，所以采用求逆元求解的方法。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;
#define ll long long
ll q_pow(ll a,ll x,ll Mod)
{
    ll res=1;
    while(x){
        if(x&1LL) res=res*a%Mod;
        a=a*a%Mod;
        x>>=1;
    }return  res;
}
int main()
{
    ll n,k,sum;
    while(~scanf("%lld%lld",&n,&k)){
        if(n==0&&k==0) return 0;
        sum=(q_pow(2,3*n+1,250*k)-1)*(q_pow(251,n+1,250*k)-1)%(250*k)/250;
        printf("%lld
",q_pow(2008,sum,k));
    } return 0;
}