NOI2012 随机数生成器

http://www.lydsy.com/JudgeOnline/problem.php?id=2875

【问题描述】
栋栋最近迷上了随机算法，而随机数生成是随机算法的基础。栋栋准备使用线性同余法（Linear Congruential Method）来生成一个随机数列，这种方法需要设置四个非负整数参数ｍ, a, c, X0，按照下面的公式生成出一系列随机数<Xn>：
Xn+1= (a Xn + c) mod m
其中 mod m 表示前面的数除以m的余数。从这个式子可以看出，这个序列的下一个数总是由上一个数生成的。
用种方法生成的序列具有随机序列的性质，因此这种方法被广泛地使用，包括常用的C++ 和Pascal 的产生随机数的库函数使用的也是这种方法。
知道这样产生的序列具有良好的随机性，不过心急的他仍然想尽快知道Xn是多少。由于栋栋需要的随机数是 0, 1,…, g − 1 之间的，他需要将Xn除以g。取余得到他想要的数，即Xn mod g，你只需要告诉栋栋他想要的数Xn mod g 是多少就可以了。

【时间限制】 1000ms

【内存限制】 512MB

【输入格式】
　　输入文件random.in 中包含 6 个用空格分割的整数 m，a，c，X0，n和g ，其中a，c，X0 是非负整数，m，n，g 是正整数。

【输出格式】
　　输入到文件random.out 中，输出一个数，即Xn mod g。

【样例输入】
11 8 7 1 5 3

【样例输出】
2

【样例说明】
<Xn>前几项依次是：

　　k   0 1 2 3 4 5
　　Xk 1 4 6 0 7 8

因此答案为 X5 mod g = 8 mod 3 = 2。

【解题报告】
　　纯裸可以得50分。
　　被几位大神虐了之后我终于懂了满分的写法。只看Xn+1=(a*Xn+c) 是一个递推数列，于是可以用矩阵求解。写出来是这样：

|Xn 1| * |a 0| = |a*Xn+c 1|
|c 1|

　　而a*Xn+c mod m 之后就是Xn+1。我们先不看mod m，有了这个矩阵，知道X0就可以求出Xn，因为中间2*2的矩阵没有发生变化，所以将中间的矩阵乘以n次幂即可得解。
　　再看数据范围，10^18，一次一次乘肯定会T，于是采用矩阵快速幂。最后还有一个问题，就是两个接近10^18的数相乘肯定会爆掉，所以可以将所有的乘法写成喵老师假期讲的快速加法。这样做的好处是可以边加边模，这样就不会爆掉了。
　　二掉的几个地方：
　　1、矩阵快速幂的初值是单位矩阵，调了半天才发现……
　　2、快速加法里每一个加的地方都要模的（其实貌似就两个地方），开始漏了一个，挂了三个点……

【心得】
　　被虐成洪特了。。。。

type
    matrix=record
        d:array[1..5,1..5]of int64;
        sx,sy:longint;                        //size
    end;
var
    ms,mm,me:matrix;
    i,j:longint;
    n,a,c,x0,s,g,m:int64;
function q_plus(t1,t2:int64):int64;
var
    t:int64;
begin
    q_plus:=0;
    t:=t1;
    while (t2<>0) do begin
        if odd(t2) then q_plus:=(q_plus+t)mod m;
        t2:=t2 shr 1;
        t:=(t+t)mod m;
    end;
end;
function m_multi(t1,t2:matrix):matrix;
var
    x,y,z:longint;
    t:int64;
begin
    for x:=1 to t1.sx do
        for y:=1 to t2.sy do begin
            t:=0;
            for z:=1 to t1.sy do
                inc(t,q_plus(t1.d[x,z],t2.d[z,y]));
            m_multi.d[x,y]:=t;
        end;
    m_multi.sx:=t1.sx;
    m_multi.sy:=t2.sy;
end;
function q_power(t1:matrix;t2:int64):matrix;
var
    t:matrix;
begin
    t:=t1;
    q_power.sx:=2;
    q_power.sy:=2;
    q_power.d[1,1]:=1;
    q_power.d[1,2]:=0;
    q_power.d[2,1]:=0;
    q_power.d[2,2]:=1;
    while (t2<>0) do begin
        if odd(t2) then q_power:=m_multi(q_power,t);
        t2:=t2 shr 1;
        t:=m_multi(t,t);
    end;
end;
begin
    assign(input,'random.in');
    reset(input);
    assign(output,'random.out');
    rewrite(output);
    
    read(m,a,c,x0,n,g);
    ms.sx:=1;
    ms.sy:=2;
    ms.d[1,1]:=x0;
    ms.d[1,2]:=1;
    mm.sx:=2;
    mm.sy:=2;
    mm.d[1,1]:=a;
    mm.d[1,2]:=0;
    mm.d[2,1]:=c;
    mm.d[2,2]:=1;

    if (m=1)or(g=1)then writeln('0') else begin
        mm:=q_power(mm,n);
        me:=m_multi(ms,mm);
        writeln((me.d[1,1] mod m)mod g);
    end;

    close(input);
    close(output);
end.