【摘记】
数论题目。ACM。
细节需考虑完整。
【题目描述】
给定N和P，以及F(1)=A,F(2)=B,F(i)=F(i-1)*F(i-2) (i>=3)。求F(n) mod P。1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000
【解题分析】
令S(i)为斐波那契第I项，易得F(i)=A^S(i-2)*B^S(i-1)。暴力显然不能做，由于N比较大，所以考虑优化。首先我们遇到的一个问题是求斐波那契第I项，这可以用矩阵加速做到O(logn)的时间（不考虑精度），可以如果真的要把高精度数算出来，时空复杂度都是明显不够的。我想到了欧拉函数，但是欧拉函数要满足底数和模数互质，这个题目却不一定满足，无法用，于是一直想不出来。
Sol1
P不大，用找循环节的方法在O(P)时间内找出循环长度（易证必定<=P），这样我们就可以把指数缩小到<=P，然后快速幂求解。
对于a和b指数小于1000000，为防止特殊情况，直接做；
有可能a^x mod p=0，需特判
Sol2
公式a^b mod p=a^(b mod φ(p)+φ(p))
（考虑到b<=φ(p) 时，可能a^b mod p>0 但a^(b mod φ(p)+φ(p)) =0，所以b<=φ(p) 需特判）
你懂的，不解释。
共同需要注意的问题：
N很小直接算，避免矩阵加速时出错；
//找循环节
type mat=array[1..2,1..2] of int64;
var t,n,i,j,l,aa,bb,zzz,qq,a2,b2:longint;xx,yy,zz,a,b,p,z,ansa,ansb:int64;
yuan,tot:mat;
hash:array[0..1000000] of longint;
function mul(a,b:mat):mat; {矩阵乘法}
var i,j,k:longint;
begin
fillchar(mul,sizeof(mul),0);
for k:=1 to 2 do
for i:=1 to 2 do
for j:=1 to 2 do
mul[i,j]:=(mul[i,j]+a[i,k]*b[k,j] mod Z) mod Z;
end;
function mul2(a,b:mat):mat;
var i,j,k:longint;
begin
fillchar(mul2,sizeof(mul2),0);
for k:=1 to 2 do
for i:=1 to 2 do
for j:=1 to 1 do
mul2[i,j]:=(mul2[i,j]+a[i,k]*b[k,j] mod Z) mod Z;
end;
function get(i:int64):mat; {快速幂}
begin
if i=1 then get:=yuan
else
if i=2 then get:=mul(yuan,yuan)
else
if odd(i) then begin get:=get(i div 2);exit(mul(mul(get,get),yuan));end
else begin get:=get(i div 2);exit(mul(get,get));end;
end;
function get_2(i,s:int64):int64;
begin
if i=0 then exit(1)
else
if odd(i) then exit(sqr(get_2(i div 2,s)) mod P*s mod P)
else exit(sqr(get_2(i div 2,s)) mod P);
end;
function getlong(s:int64):longint; {找循环节长度}
var i,j:longint;tmp:int64;
begin
fillchar(hash,sizeof(hash),0);
tmp:=1;
for i:=1 to maxlongint do
begin
tmp:=(tmp*s) mod P;
if tmp=0 then begin getlong:=-1;exit;end;
if hash[tmp]<>0 then begin getlong:=i-hash[tmp];exit;end;
hash[tmp]:=i;
//writeln(s,' ',tmp);
end;
end;
begin
assign(input,'bruteforce.in');
assign(output,'bruteforce.out');
reset(input);rewrite(output);
readln(t);
for l:=1 to t do
begin
readln(a,b,p,n);
write('Case #',l,': ');
if n=1 then writeln(a mod P)
else
if n=2 then writeln(b mod P)
else
if n=3 then writeln(a*b mod P)
else
if n=4 then writeln(a*b mod P*b mod P)
else
if n=5 then writeln(a*a mod P*b mod P*b mod P*b mod P) {小数据直接输出}
else
begin
qq:=3;xx:=1;yy:=1;
repeat
inc(qq);
zz:=yy;
yy:=xx+yy;
xx:=zz;
until (qq=n)or(xx>1000000);
if xx<=1000000 then
//对于a和b指数小于1000000，为防止特殊情况，直接做
begin
writeln(get_2(xx,a)*get_2(yy,b) mod P);
end
else
begin
aa:=getlong(a);
bb:=getlong(b);
if aa=-1 then ansa:=0
else
begin
yuan[1,1]:=1;yuan[1,2]:=1;yuan[2,1]:=1;yuan[2,2]:=0;
Z:=aa;
tot:=get(n-4);
xx:=(tot[1,1]+tot[1,2]) mod Z;
xx:=xx+(1000000 div Z+1)*Z; {因为循环可能不是从第一个开始的，所以要找后面某一个}
ansa:=get_2(xx,a);
end;
if bb=-1 then ansb:=0
else
begin
Z:=bb;
tot:=get(n-3);
xx:=(tot[1,1]+tot[1,2]) mod Z;
xx:=xx+(1000000 div Z+1)*Z;
ansb:=get_2(xx,b);
end;
//writeln(ansa,' ',ansb);
writeln(ansa*ansb mod P);
end;
end;
end;
close(input);close(output);
end.
//套公式
type mat=array[1..2,1..2] of int64;
var t,n,i,j,l,aa,bb,zzz,qq,a2,b2:longint;q,ph,a,b,p,z,ansa,ansb:int64;
yuan,tot:mat;
hash:array[0..1000000] of longint;
function mul(a,b:mat):mat; {矩阵乘法}
var i,j,k:longint;
begin
fillchar(mul,sizeof(mul),0);
for k:=1 to 2 do
for i:=1 to 2 do
for j:=1 to 2 do
mul[i,j]:=(mul[i,j]+a[i,k]*b[k,j] mod Z) mod Z;
end;
function get(i:int64):mat; {快速幂}
begin
if i=1 then get:=yuan
else
if i=2 then get:=mul(yuan,yuan)
else
if odd(i) then begin get:=get(i div 2);exit(mul(mul(get,get),yuan));end
else begin get:=get(i div 2);exit(mul(get,get));end;
end;
function get_2(i,s:int64):int64;
begin
if i=0 then exit(1)
else
if odd(i) then exit(sqr(get_2(i div 2,s)) mod P*s mod P)
else exit(sqr(get_2(i div 2,s)) mod P);
end;
begin
assign(input,'bruteforce.in');
assign(output,'bruteforce.out');
reset(input);rewrite(output);
readln(t);
for l:=1 to t do
begin
readln(a,b,p,n);
write('Case #',l,': ');
if n=1 then writeln(a mod P)
else
if n=2 then writeln(b mod P)
else
if n=3 then writeln(a*b mod P)
else
if n=4 then writeln(a*b mod P*b mod P)
else
if n=5 then writeln(a*a mod P*b mod P*b mod P*b mod P) {小数据直接输出}
else
begin
ph:=p;
q:=p;
for i:=2 to trunc(sqrt(p)) do
if q mod i=0 then
begin
ph:=ph*(i-1) div i;
while q mod i=0 do q:=q div i;
end;
if q>1 then ph:=ph*(q-1) div q;
Z:=ph;
yuan[1,1]:=1;yuan[1,2]:=1;yuan[2,1]:=1;yuan[2,2]:=0;
tot:=get(n-4);
if tot[1,1]+tot[1,2]<=ph then ansa:=get_2(tot[1,1]+tot[1,2],a)
else ansa:=get_2((tot[1,1]+tot[1,2]) mod ph+ph,a);
tot:=get(n-3);
if tot[1,1]+tot[1,2]<=ph then ansb:=get_2(tot[1,1]+tot[1,2],b)
else ansb:=get_2((tot[1,1]+tot[1,2]) mod ph+ph,b);
writeln(ansa*ansb mod P);
end;
end;
close(input);close(output);
end.