SDOI2010古代猪文 总结

猪王国的文明源远流长，博大精深。
　　iPig在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为N。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁逐渐进行了简化，去掉了一些不常用的字。
　　iPig打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。
　　iPig觉得只要符合文献，每一种能整除N的k都是有可能的。他打算考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig预计，如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。
　　现在他想知道猪王国研究古代文字的代价是多少。由于iPig觉得这个数字可能是天文数字，所以你只需要告诉他答案除以999911659的余数就可以了。
　　Input

　　输入文件ancient.in有且仅有一行：两个数N、G，用一个空格分开。

Output

　　输出文件ancient.out有且仅有一行：一个数，表示答案除以999911659的余数。

Sample Input

4 2

Sample Output

2048

Data Constraint

　　10%的数据中，1 <= N <= 50；
　　20%的数据中，1 <= N <= 1000；
　　40%的数据中，1 <= N <= 100000；
　　100%的数据中，1 <= G <= 1000000000，1 <= N <= 1000000000。

简化版题意：实际上是让你求G^M mod p,其中M为∑ C（i,N）(i|N)，p=999911659
说实话我刚开始以为是所有的i加起来，算了半天不知道样例怎么算的，一脸懵逼。
后来看见“所有可能的k的所有情况数”才知道是组合数，知道了也没卵用，草草打了个40分的暴力。后来看题解才知道这题不是数学大触做不出来，弃疗QAQ。。

首先p是质数，G^M mod p直接上费马小定理。。化成G^(M mod p-1)
现在就是要知道M mod p-1，所以可以用lucas定理，此时发现p-1不是质数，p-1=999911658，嗯~o(￣▽￣)o，有点奇怪。质因数分解一下：
999911658=2*3*4679*35617。可以做了，lucas套完之后出来四个一次同余方程，可以用CRT解啦啦啦啦。

代码不长，但是很恶心。。

const p=999911659;
var
        t:Array[1..4]of longint=(2,3,4679,35617);
        fac:array[0..4,0..40000]of longint;
        m:array[0..4]of longint;
        i,j,k,n,g,tmp:longint;
function pow(a,b,x:int64):int64;
var
        ans:int64;
begin
        ans:=1;
        while b>0 do
        begin
                if b and 1=1 then ans:=(ans*a)mod x;
                b:=b>>1;
                a:=(a*a)mod x;
        end;
        exit(ans);
end;
function c(a,b,x:int64):int64;
begin
        if b>a then exit(0);
        exit(fac[x,a]*pow(fac[x,a-b]*fac[x,b],t[x]-2,t[x])mod t[x]);
end;
function lucas(a,b,x:int64):int64;
begin
        if b=0 then exit(1);
        exit( c(a mod t[x],b mod t[x],x)*lucas(a div t[x],b div t[x],x)mod t[x] );
end;
procedure gcd(a,b:int64;var x,y:int64);
var
        t:longint;
begin
        if b=0 then
        begin
                x:=1;
                y:=0;
                exit;
        end;
        gcd(b,a mod b,x,y);
        t:=x;
        x:=y;
        y:=t-(a div b)*y;
end;
function chinese:int64;
var
        n,d,x0,y0,ans:int64;
        i:longint;
begin
        n:=1;
        ans:=0;

        for i:=1 to 4 do
        begin
                n:=n*t[i];
                //inc(i);
        end;

        for i:=1 to 4 do
        begin
                d:=n div t[i];
                gcd(d,t[i],x0,y0);
                ans:=(ans+d*x0*m[i])mod n;
        end;
        while ans<=0 do ans:=ans+n;
        exit(ans);
end;
begin
        {assign(input,'Yakumo.in');
        assign(output,'Yakumo.out');
        reset(input);
        rewrite(output);}
        readln(n,g);
        if g=p then
        begin
                writeln(0);
                exit;
        end;
        g:=g mod p;
        for i:=1 to 4 do
        begin
                fac[i,0]:=1;
                for j:=1 to t[i] do
                fac[i,j]:=(fac[i,j-1]*j)mod t[i];
        end;
        i:=1;
        while i*i<=n do
        begin
                if n mod i=0 then
                begin
                        tmp:=n div i;
                        for j:=1 to 4 do
                        begin
                                if tmp<>i then
                                m[j]:=(m[j]+lucas(n,i,j))mod t[j];
                                //else 因为是所有情况所以不能用else
                                m[j]:=(m[j]+lucas(n,tmp,j))mod t[j];
                        end;
                end;
                inc(i);
        end;
        writeln(pow(g,chinese,p));
        {close(input);
        close(output);}
end.

4ms，超快。。