JZOJ_5845. 西西算数

Description

西西是可爱的孩子,他总是喜欢算数,每天都在算啊算啊......这天,西西学会了如何进行指数运算,
于是他算出了√5=2.2360679774997896964091736687313
然后他又算出了 (3 + √5)^2 = 27.4164079... (3 + √5)^5 = 3935.73982.... 西西顿时对(3 + √5)^n产生了浓厚的兴趣。但对于N=2000000000时，西西一下子晕了，算了好久也没算出来，他只好求助于你。为了你的计算方便，你只需要告诉西西(3 + √5)^n 的整数部分的最后三位就行了（不够的在前面补0，如果有前导0也一起输出）

Input

第一行一个T（T<=100），表示有T组数据。
每组数据一个N（2<=N<=2000000000）

Output

T行, 每组数据输出：Case #X: Y
X表示是第几组数据,Y是一个三位的数字串表示对应的答案.

Solution

我们设 $a_n=(3+sqrt{5})^n+(3-sqrt{5})^n$

可以知道 $a_n$ 是一个整数，在这不给予证明。

因为 $(3-sqrt{5})^n$ 的值域为(0,1)，所以 $a_n-1$ 即为答案。

怎么求出 $a_n$ 呢？

如上等式可以变为 $a_n=xcdot a_{n-1}+ycdot a_{n-2}$

我们要做的就是求出系数x和y，

把式子 $a_{n-1}=(3+sqrt{5})^{n-1}+(3-sqrt{5})^{n-1}$ 两边同时乘上 $3+sqrt{5}$

化简得： $(3+sqrt{5})a_{n-1}=(3+sqrt{5})^n-4(3-sqrt{5})^{n-2}$

同理可得： $(3-sqrt{5})a_{n-1}=(3-sqrt{5})^n-4(3+sqrt{5})^{n-2}$

两个式子相加后，运用最初的式子可得： $6a_{n-1}=a_n-4a_{n-2}$

移项得： $a_n=6a_{n-1}-4a_{n-2}$ ,系数x和y的值就知道了。

因为n有2*10^9这么大，for一遍都会超时，我们考虑矩阵乘法。

转移矩阵看代码，不会矩阵乘法的问度娘。

代码

 1 type
 2   arr=array [1..2,1..2] of longint;
 3 var
 4   n,t,ans,nm:longint;
 5   a,c:arr;
 6 procedure mi(a,b:arr);
 7 var
 8   i,j,k:longint;
 9 begin
10   fillchar(c,sizeof(c),0);
11   for i:=1 to 2 do
12     for j:=1 to 2 do
13       for k:=1 to 2 do
14         c[i,k]:=(c[i,k]+a[i,j]*b[j,k]) mod 1000;
15 end;
16 
17 procedure main(n:longint);
18 begin
19   if n<=1 then exit;
20   main(n div 2);
21   mi(c,c);
22   if n mod 2=1 then mi(c,a);
23 end;
24 
25 begin
26   assign(input,'Count.in');
27   assign(output,'Count.out');
28   reset(input);
29   rewrite(output);
30   readln(t); nm:=0;
31   while t>0 do
32     begin
33       inc(nm);
34       readln(n);
35       a[1,1]:=0; a[1,2]:=-4;
36       a[2,1]:=1; a[2,2]:=6;
37       c:=a;
38       main(n);
39       ans:=(c[1,1]+c[2,2]+1999) mod 1000;
40       write('Case #',nm,': ');
41       if ans<10 then write('00') else
42         if ans<100 then write('0');
43       writeln(ans);
44       dec(t);
45     end;
46   close(input);
47   close(output);
48 end.