hdu4686Arc of Dream(推公式 或者 矩阵乘法)

题意：

 

a ₀ = A0 
a _i = a _i-1*AX+AY 
b ₀ = B0 
b _i = b _i-1*BX+BY

 

给定0<=n<=10^18，A0,Ax,Ay，B0,Bx,By。求和。

分析：

【推公式】线性递推式，比较显然得可以求出通项公式，然作乘法求和。就是一系列等比数列和了。注意为1的时候需要特判。麻烦死。

ai = Ax^i * A0 + Ay*(Ax^i-1 + Ax^i-2 + ...+ 1)同理bi

然后化简求之

【矩阵乘法】线性递推方程，很容易构造出变换矩阵。做矩阵乘法。

si 为前i项和：

si = si-1 + Ax*Bx *ai-1 * bi-1 + Ax*By*ai-1 + Bx*Ay*bi-1 + Ay*By

递推式包含5项，构造5*5的矩阵即可。

 

Mark::注意n为0的时候，输出0！！！！特判。跪在这里了！

 

推公式版（略麻烦）：

  1 #include<stdio.h>
  2 #include<math.h>
  3 #include<string.h>
  4 #define MOD 1000000007
  5 #define LL long long
  6 
  7 LL A0,Ax,Ay;
  8 LL B0,Bx,By;
  9 LL n;
 10 
 11 
 12 LL power(LL a, LL b)
 13 {
 14     LL res = 1;
 15     while (b)
 16     {
 17         if (b & 1)
 18         {
 19             res = (res * a) % MOD;
 20         }
 21         a = (a * a) % MOD;
 22         b = b / 2 ;
 23     }
 24     return res;
 25 }
 26 LL get(LL x)
 27 {
 28     if (x == 0)
 29         return 1;
 30     return power(x,MOD-2);
 31 }
 32 LL ggg(LL n)
 33 {
 34     LL a = n % MOD;
 35     LL b = (n+1) % MOD;
 36     LL c = (2*n + 1) % MOD;
 37     LL res = (((a * b) % MOD * c) % MOD * get(6))%MOD;
 38     return res;
 39 }
 40 LL hhh(LL n)
 41 {
 42     LL a = n % MOD;
 43     LL b = (n + 1) % MOD;
 44     LL res = ((a*b)% MOD * get(2)) % MOD;
 45     return res;
 46 }
 47 LL calc(LL x,LL n)
 48 {
 49     if (x == 1)
 50         return n;
 51     LL res = power(x,n+1);
 52     res = ((res - x) %  MOD + MOD) % MOD;
 53     res = (res * get(x-1)) % MOD;
 54     return res;
 55 }
 56 
 57 LL www(LL x,LL n)
 58 {
 59     LL sss = (n*power(x, n + 1) ) % MOD;
 60     sss = ((sss - calc(x,n)) % MOD + MOD) % MOD;
 61     sss = (sss * get(x - 1)) % MOD;
 62     return sss;
 63 }
 64 LL solve(LL n)
 65 {
 66     if (n == 0)
 67         return (A0 * B0) % MOD;
 68     LL ans = 0;
 69     LL res = 0;
 70     LL ABx = (Ax * Bx) % MOD;
 71     LL AB0 = (A0 * B0) % MOD;
 72     LL ABy = (Ay * By) % MOD;
 73 
 74 
 75     LL ttt;
 76     LL sss;
 77     LL tmp = calc(ABx,n);
 78     res  = (res + AB0 * tmp) % MOD;
 79 
 80 
 81     ans = (ans + res) % MOD;
 82     // printf("%I64d    %I64d\n",res,ans);
 83     res = 0;
 84 
 85     if (Ax == 1)
 86     {
 87         if (Bx == 1)
 88         {
 89             sss = (ABy * ggg(n)) % MOD;
 90             res = sss;
 91         }
 92         else
 93         {
 94             sss = www(Bx,n);
 95             sss = ((sss - hhh(n)) % MOD + MOD) % MOD;
 96 
 97             res = (ABy * get(Bx - 1)) % MOD;
 98             res = (res * sss) % MOD;
 99         }
100     }
101     else
102     {
103         if (Bx == 1)
104         {
105             sss = www(Ax,n);
106             sss = ((sss - hhh(n)) % MOD + MOD) % MOD;
107 
108             res = (ABy * get(Ax - 1)) % MOD;
109             res = (res * sss) % MOD;
110         }
111         else
112         {
113             sss = (ABy * get(((Ax-1)*(Bx-1))% MOD) ) % MOD;
114 
115             // printf("%I64d\n",sss);
116 
117             tmp = calc(ABx,n);
118             res = (res + tmp*sss) % MOD;
119 
120             tmp = calc(Ax,n);
121             res = ((res - tmp*sss) % MOD + MOD) % MOD;
122 
123             tmp = calc(Bx,n);
124             res = ((res - tmp*sss) % MOD + MOD) % MOD;
125 
126             tmp = n % MOD ;
127             res = (res + sss * tmp) % MOD;
128         }
129     }
130 
131 
132 
133     ans = (ans + res) % MOD;
134   //    printf("changshu = %I64d        %I64d\n",res,ans);
135     res = 0;
136 
137     //a0
138     if (Bx != 1)
139     {
140         ttt = (A0 * By) % MOD;
141         ttt = (ttt * get(Bx - 1)) % MOD;
142 
143         tmp = calc(ABx,n);
144         res = (res + tmp * ttt) % MOD;
145 
146         tmp = calc(Ax,n);
147         res = ((res - tmp * ttt) % MOD + MOD) % MOD;
148 
149     }
150     else
151     {
152         if (Ax == 1)
153         {
154             ttt = (A0 * By) % MOD;
155             res = (ttt * hhh(n)) % MOD;
156         }
157         else
158         {
159             ttt = (A0 * By) % MOD;
160             sss = www(Ax,n);
161             res = (ttt * sss) % MOD;
162         }
163     }
164 
165     ans = (ans + res) % MOD;
166   //  printf("now = %I64d %I64d\n",res,ans);
167 
168 
169     //b0
170     res = 0;
171     if (Ax != 1)
172     {
173         ttt = (B0 * Ay) % MOD;
174         ttt = (ttt * get(Ax - 1)) % MOD;
175 
176         tmp = calc(ABx,n);
177         res = (res + tmp * ttt) % MOD;
178 
179         tmp = calc(Bx,n);
180         res = ((res - tmp * ttt) % MOD + MOD) % MOD;
181     }
182     else
183     {
184         if (Bx == 1)
185         {
186             ttt = (B0 * Ay) % MOD;
187             res = (ttt * hhh(n)) % MOD;
188         }
189         else
190         {
191             ttt = (B0 * Ay) % MOD;
192 
193             sss = www(Bx,n);
194 
195             res = (ttt * sss) % MOD;
196         }
197     }
198 
199 
200     ans = (ans + res) % MOD;
201   //  printf("cur = %I64d %I64d\n",res,ans);
202     res = 0;
203 
204     ans = (ans + A0 * B0) % MOD;
205     ans =  ( (ans + MOD) % MOD + MOD ) % MOD ;
206     return ans;
207 }
208 int main()
209 {
210   //  printf("%I64d\n",ggg(3));
211   //  printf("%I64d\n",www(2,4));
212     while (scanf("%I64d",&n)==1)
213     {
214         scanf("%I64d%I64d%I64d",&A0,&Ax,&Ay);
215         scanf("%I64d%I64d%I64d",&B0,&Bx,&By);
216         A0 %= MOD;
217         Ax %= MOD;
218         Ay %= MOD;
219         B0 %= MOD;
220         Bx %= MOD;
221         By %= MOD;
222         printf("%I64d\n",solve(n-1));//n等于0的时候没有特判，但是返回了正确值。
223     }
224     return 0;
225 }

hdu4686(推公式)

 

矩阵乘法：

  1 #include<stdio.h>
  2 #include<math.h>
  3 #include<string.h>
  4 #define LL long long
  5 #define MOD 1000000007
  6 struct matrix{
  7     LL mat[5][5];
  8     void init(){
  9         memset(mat,0,sizeof(mat));
 10     }
 11     void mod()
 12     {
 13         for (int i=0;i<5;i++)
 14             for (int j=0;j<5;j++)
 15             {
 16                 mat[i][j] = (mat[i][j] % MOD + MOD) % MOD;
 17             }
 18         return ;
 19     }
 20     void EE()
 21     {
 22         init();
 23         for (int i=0;i<5;i++)
 24             mat[i][i] = 1;
 25         return ;
 26     }
 27 
 28 }M,E;
 29 
 30 LL A0,Ax,Ay;
 31 LL B0,Bx,By;
 32 LL n;
 33 
 34 LL S[5];
 35 
 36 void init()
 37 {
 38     A0 %= MOD;
 39     Ax %= MOD;
 40     Ay %= MOD;
 41     B0 %= MOD;
 42     Bx %= MOD;
 43     By %= MOD;
 44 
 45     M.init();
 46     M.mat[0][0] = 1;
 47     M.mat[0][1] = M.mat[1][1] = Ax * Bx;
 48     M.mat[0][2] = M.mat[1][2] = Ax * By;
 49     M.mat[0][3] = M.mat[1][3] = Ay * Bx;
 50     M.mat[0][4] = M.mat[1][4] = Ay * By;
 51 
 52     M.mat[2][2] = Ax;
 53     M.mat[2][4] = Ay;
 54     M.mat[3][3] = Bx;
 55     M.mat[3][4] = By;
 56     M.mat[4][4] = 1;
 57     M.mod();
 58 
 59     S[0] = ((A0 * B0) % MOD + MOD) % MOD;
 60     S[1] = ((A0 * B0) % MOD + MOD) % MOD;
 61     S[2] = A0;
 62     S[3] = B0;
 63     S[4] = 1;
 64     return ;
 65 }
 66 matrix mul(matrix a,matrix b)
 67 {
 68     matrix c;
 69     c.init();
 70 
 71     for (int i=0;i<5;i++)
 72         for (int j=0;j<5;j++)
 73         {
 74             c.mat[i][j] = 0;
 75             for (int k=0;k<5;k++)
 76             {
 77                 c.mat[i][j] = ((c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % MOD + MOD) % MOD;
 78             }
 79         }
 80     c.mod();
 81     return c;
 82 }
 83 matrix power(matrix mm,LL n)
 84 {
 85     E.EE();
 86     while (n)
 87     {
 88         if (n %2 == 1)
 89         {
 90             E = mul(E,mm);
 91         }
 92         mm = mul(mm,mm);
 93         n /= 2;
 94     }
 95     return E;
 96 }
 97 LL solve(LL n)
 98 {
 99     if (n == 0)
100         return ((A0 * B0) % MOD + MOD) % MOD;
101     M = power(M,n);
102     LL res = 0;
103     for (int i=0;i<5;i++)
104     {
105       //  printf("%I64d %I64d\n",M.mat[0][i],S[i]);
106         res = ((res + M.mat[0][i]*S[i]) % MOD + MOD) % MOD;
107     }
108    // puts("");
109     return (res % MOD + MOD) % MOD;
110 }
111 int main()
112 {
113     while (scanf("%I64d",&n)==1){
114         scanf("%I64d%I64d%I64d",&A0,&Ax,&Ay);
115         scanf("%I64d%I64d%I64d",&B0,&Bx,&By);
116 
117         if (n == 0){//1的时候需要特判！！！
118             printf("0\n");
119             continue;
120         }
121         init();
122 
123         LL ans = solve(n - 1);
124         ans = (ans % MOD + MOD) %  MOD;
125        // while (ans < 0 || ans >= MOD);
126         printf("%I64d\n",ans);
127     }
128     return 0;
129 }

hdu4686

 

 

转载于:https://www.cnblogs.com/wangsouc/articles/3272045.html