1 #ifndef matrix_h
2 #define matrix_h
3 #include <vector>
4 #include <time.h>
5 #include <math.h>
6 #include <iomanip>
7 #include <random>
8 using std::vector;
9 using std::ios;
10
11 /*
12 Normal distribution random number
13 About File
14 */
15
16 #define abs(T) (T>0?T:-T)
17
18 typedef enum Direction{
19 LT, RT, LB, RB
20 }Dir;
21
22
23 template <typename T>
24 class Vector;
25
26 template <typename T>
27 class Matrix {
28 public:
29 int x, y;
30 vector<vector<T> > M;
31 Matrix(int a,int b) {
32 this->Set(a,b);
33 }
34 Matrix() {
35 this->x = this->y=0;
36 }
37 Matrix(const Matrix &b) {
38 this->x = b.x;
39 this->y = b.y;
40 this->M = b.M;
41 }
42 void operator +=(const Matrix a) {
43 if (this->x != a.x || this->y != a.y) {
44 return;
45 }
46 for (int i = 0;i < this ->x;i++) {
47 for (int j = 0;j < this->y;j++) {
48 this->M[i][j] += a.M[i][j];
49 }
50 }
51 }
52 void operator -=(const Matrix a) {
53 if (this->x != a.x || this->y != a.y) {
54 return;
55 }
56 for (int i = 0;i < this->x;i++) {
57 for (int j = 0;j < this->y;j++) {
58 this->M[i][j] -= a.M[i][j];
59 }
60 }
61 }
62
63 void Plus(const Matrix a,const Matrix b) {
64 if (a.x != b.x || a.y != b.y) {
65 return;
66 }
67 this->x = a.x;
68 this->y = a.y;
69 this->M.clear();
70 this->Set(a.x, a.y);
71 for (int i = 0;i < a.x;i++) {
72 for (int j = 0;j < a.y;j++) {
73 this->M[i][j] = a.M[i][j] + b.M[i][j];
74 }
75 }
76 }
77 void Sub(const Matrix a, const Matrix b) {
78 if (a.x != b.x || a.y != b.y) {
79 return;
80 }
81 this->x = a.x;
82 this->y = a.y;
83 this->M.clear();
84 this->Set(a.x, a.y);
85 for (int i = 0;i < a.x;i++) {
86 for (int j = 0;j < a.y;j++) {
87 this->M[i][j] = a.M[i][j] - b.M[i][j];
88 }
89 }
90 }
91 void Plus(const Matrix a, const Matrix b,int p,int q) {
92 if (p == 0 && q == 0) {
93 if (a.x != b.x || a.y != b.y) {
94 return;
95 }
96 this->x = a.x;
97 this->y = a.y;
98 this->M.clear();
99 this->Set(a.x, a.y);
100 for (int i = 0;i < a.x;i++) {
101 for (int j = 0;j < a.y;j++) {
102 this->M[i][j] = a.M[i][j] + b.M[i][j];
103 }
104 }
105 }else if (p == 1 && q == 1) {
106 if (a.x != b.x || a.y != b.y) {
107 return;
108 }
109 this->x = a.y;
110 this->y = a.x;
111 this->M.clear();
112 this->Set(a.y, a.x);
113 for (int i = 0;i < a.y;i++) {
114 for (int j = 0;j < a.x;j++) {
115 this->M[i][j] = a.M[j][i] + b.M[j][i];
116 }
117 }
118 }else if (p == 0 && q == 1) {
119 if (a.x != b.y || a.y != b.x) {
120 return;
121 }
122 this->x = a.x;
123 this->y = a.y;
124 this->M.clear();
125 this->Set(a.x, a.y);
126 for (int i = 0;i < a.x;i++) {
127 for (int j = 0;j < a.y;j++) {
128 this->M[i][j] = a.M[i][j] + b.M[j][i];
129 }
130 }
131 }else if (p == 1 && q == 0) {
132 if (a.x != b.y || a.y != b.x) {
133 return;
134 }
135 this->x = a.y;
136 this->y = a.x;
137 this->M.clear();
138 this->Set(a.y, a.x);
139 for (int i = 0;i < a.y;i++) {
140 for (int j = 0;j < a.x;j++) {
141 this->M[i][j] = a.M[j][i] + b.M[i][j];
142 }
143 }
144 }
145 }
146 void Mult(const Matrix<T> a, const Matrix<T> b) {
147 if (a.y != b.x) {
148 return;
149 }
150 this->x = a.x;
151 this->y = b.y;
152 this->M.clear();
153 this->Set(a.x, b.y);
154 for (int i = 0;i < a.x;i++) {
155 for (int j = 0;j < b.y;j++) {
156 T tmp = 0;
157 for (int k = 0;k < a.y;k++) {
158 tmp += a.M[i][k] * b.M[k][j];
159 }
160 this->M[i][j] = tmp;
161 }
162 }
163 }
164
165 void operator *=(const Matrix<T> b) {
166 if (this->y != b.x) {
167 return;
168 }
169 Matrix<T> result(this->x,b.y);
170 for (int i = 0;i < this->x;i++) {
171 for (int j = 0;j < b.y;j++) {
172 T tmp = 0;
173 for (int k = 0;k < this->y;k++) {
174 tmp += this->M[i][k] * b.M[k][j];
175 }
176 result.M[i][j] = tmp;
177 }
178 }
179 this->x = result.x;
180 this->y = result.y;
181 this->M = result.M;
182 }
183 void Negative() {
184 for (int i = 0;i < this->x;i++) {
185 for (int j = 0;j < this->y;j++) {
186 this->M[i][j] = -this->M[i][j];
187 }
188 }
189 }
190 void Zero() {
191 for (int i = 0;i < this->x;i++) {
192 for (int j = 0;j < this->y;j++) {
193 this->M[i][j] =0;
194 }
195 }
196 }
197 void DelError(T eps) {
198 for (int i = 0;i < this->x;i++) {
199 for (int j = 0;j < this->y;j++) {
200 if(abs(this->M[i][j])<eps)
201 this->M[i][j] = 0;
202 }
203 }
204 }
205 void Mult(Matrix &b,int landa,int p) {
206 int i, j;
207 if (p == 0) {
208 this->Set(b.x,b.y);
209 this->M = b.M;
210 for (i = 0;i < this->x;i++) {
211 for (j = 0;j < this->y;j ++ ) {
212 this->M[i][j] = this->M[i][j] * landa;
213 }
214 }
215 }else if (p == 1) {
216 this->Set(b.y, b.x);
217 for (i = 0;i < this->x;i++) {
218 for (j = 0;j < this->y;j++) {
219 this->M[i][j] = b.M[j][i] * landa;
220 }
221 }
222 }
223 }
224 void DiagPlus(const Vector<T> &a) {
225 if (this->x != this->y || this->x != a.x) return;
226 for (int i = 0;i < this->x;i++) {
227 this->M[i][i] += a.V[i];
228 }
229 }
230 void LineSum(Vector<T> &a) {
231 a.Set(this->y);
232 for (int i = 0;i < this->y;i++) {
233 T sum = 0;
234 for (int j = 0;j < this->x;j++) {
235 sum += this->M[j][i];
236 }
237 a.V[i] = sum;
238 }
239 }
240 void operator >> (Vector<T> &a) {
241 a.Set(this->x*this->y);
242 int pos = 0;
243 for (int i = 0;i < this->x;i++) {
244 for (int j = 0;j < this->y;j++) {
245 a.V[pos++]=this->M[i][j];
246 }
247 }
248 }
249 void operator <<(const Vector<T> &a) {
250 this->Zero();
251 int pos = 0;
252 while (pos < a.x&&pos<this->x*this->y) {
253 htis->M[pos/this->y][pox%this->y] = a[pos];
254 pos++;
255 }
256 }
257 void ColumSum(Vector<T> &a) {
258 a.Set(this->x);
259 for (int i = 0;i < this->x;i++) {
260 T sum = 0;
261 for (int j = 0;j < this->y;j++) {
262 sum += this->M[i][j];
263 }
264 a.V[i] = sum;
265 }
266 }
267 void DiagSub(const Vector<T> &a) {
268 if (this->x != this->y || this->x != a.x) return;
269 for (int i = 0;i < this->x;i++) {
270 this->M[i][i] -= a.V[i];
271 }
272 }
273 void Mult(const Vector<T> &a,const Vector<T> &b) {
274 this->Set(a.x,b.x);
275 for (int i = 0;i < this->x;i++) {
276 for (int j = 0;j < this->y;j++) {
277 this->M[i][j] = a.V[i] * b.V[j];
278 }
279 }
280 }
281 void Mult(const Matrix a, const Matrix b, int p, int q) {
282 if (p == 0 && q == 0) {
283 if (a.y != b.x) {
284 return;
285 }
286 this->x = a.x;
287 this->y = b.y;
288 this->M.clear();
289 this->Set(a.x, b.y);
290 for (int i = 0;i < a.x;i++) {
291 for (int j = 0;j < b.y;j++) {
292 T tmp = 0;
293 for (int k = 0;k < a.y;k++) {
294 tmp += a.M[i][k] * b.M[k][j];
295 }
296 this->M[i][j] = tmp;
297 }
298 }
299 }
300 else if (p == 1 && q == 1) {
301 if (a.x != b.y) {
302 return;
303 }
304 this->x = a.y;
305 this->y = b.x;
306 this->M.clear();
307 this->Set(a.y, b.x);
308 for (int i = 0;i < a.y;i++) {
309 for (int j = 0;j < b.x;j++) {
310 T tmp = 0;
311 for (int k = 0;k < a.x;k++) {
312 tmp += a.M[k][i] * b.M[j][k];
313 }
314 this->M[i][j] = tmp;
315 }
316 }
317 }
318 else if (p == 0 && q == 1) {
319 if (a.y != b.y) {
320 return;
321 }
322 this->x = a.x;
323 this->y = b.x;
324 this->M.clear();
325 this->Set(a.x, b.x);
326 for (int i = 0;i < a.x;i++) {
327 for (int j = 0;j < b.x;j++) {
328 T tmp = 0;
329 for (int k = 0;k < a.y;k++) {
330 tmp += a.M[i][k] * b.M[j][k];
331 }
332 this->M[i][j] = tmp;
333 }
334 }
335 }
336 else if (p == 1 && q == 0) {
337 if (a.x != b.x) {
338 return;
339 }
340 this->x = a.y;
341 this->y = b.y;
342 this->M.clear();
343 this->Set(a.y, b.y);
344 for (int i = 0;i < a.y;i++) {
345 for (int j = 0;j < b.y;j++) {
346 T tmp = 0;
347 for (int k = 0;k < a.x;k++) {
348 tmp += a.M[k][i] * b.M[k][j];
349 }
350 this->M[i][j] = tmp;
351 }
352 }
353 }
354 }
355 void Sub(const Matrix a, const Matrix b, int p, int q) {
356 if (p == 0 && q == 0) {
357 if (a.x != b.x || a.y != b.y) {
358 return;
359 }
360 this->x = a.x;
361 this->y = a.y;
362 this->M.clear();
363 this->Set(a.x, a.y);
364 for (int i = 0;i < a.x;i++) {
365 for (int j = 0;j < a.y;j++) {
366 this->M[i][j] = a.M[i][j] - b.M[i][j];
367 }
368 }
369 }
370 else if (p == 1 && q == 1) {
371 if (a.x != b.x || a.y != b.y) {
372 return;
373 }
374 this->x = a.y;
375 this->y = a.x;
376 this->M.clear();
377 this->Set(a.y, a.x);
378 for (int i = 0;i < a.y;i++) {
379 for (int j = 0;j < a.x;j++) {
380 this->M[i][j] = a.M[j][i] - b.M[j][i];
381 }
382 }
383 }
384 else if (p == 0 && q == 1) {
385 if (a.x != b.y || a.y != b.x) {
386 return;
387 }
388 this->x = a.x;
389 this->y = a.y;
390 this->M.clear();
391 this->Set(a.x, a.y);
392 for (int i = 0;i < a.x;i++) {
393 for (int j = 0;j < a.y;j++) {
394 this->M[i][j] = a.M[i][j] - b.M[j][i];
395 }
396 }
397 }
398 else if (p == 1 && q == 0) {
399 if (a.x != b.y || a.y != b.x) {
400 return;
401 }
402 this->x = a.y;
403 this->y = a.x;
404 this->M.clear();
405 this->Set(a.y, a.x);
406 for (int i = 0;i < a.y;i++) {
407 for (int j = 0;j < a.x;j++) {
408 this->M[i][j] = a.M[j][i] - b.M[i][j];
409 }
410 }
411 }
412 }
413 void Init(vector<T> arr) {
414 if (arr.size() != x*y) {
415 return;
416 }
417 int pos = 0;
418 for (int i = 0;i < this->x;i++) {
419 for (int j = 0;j < this->y;j++) {
420 this->M[i][j] = arr[pos++];
421 }
422 }
423 }
424 void Unit(int n) {
425 this->Set(n,n);
426 for (int i = 0;i < this->x;i++) {
427 this->M[i][i] = 1;
428 }
429 }
430 void Random(T m,T n) {
431 int i, j;
432 if (m == n) {
433 for (i = 0;i < this->x;i++) {
434 for (j = 0;j < this->y;j++) {
435 this->M[i][j] = m;
436 }
437 }
438 return;
439 }
440 if (m > n) {
441 m = m+n;
442 n = m-n;
443 m = m-n;
444 }
445 time_t t=time(NULL);
446 srand(t*rand());
447 for (i = 0;i < this->x;i++) {
448 for (j = 0;j < this->y;j++) {
449 double t = rand() / (RAND_MAX + 0.0);
450 this->M[i][j] = (T)(t*(n-m)+m);
451 }
452 }
453 }
454 void RandomInt(int n) {
455 int i, j;
456 time_t t = time(NULL);
457 srand(t * rand());
458 for (i = 0;i < this->x;i++) {
459 for (j = 0;j < this->y;j++) {
460 double t = rand() / (RAND_MAX + 0.0);
461 this->M[i][j] = (int)(t*2*n-n);
462 }
463 }
464 }
465 void TimePlus(Matrix b,int landa) {
466 if (this->x != b.x || this->y != b.y) {
467 return;
468 }
469 for (int i = 0;i < b.x;i++) {
470 for (int j = 0;j < b.y;j++) {
471 this->M[i][j] += b.M[i][j] * landa;
472 }
473 }
474 }
475 void Turn() {
476 Matrix tmp(this->y,this->x);
477 for (int i = 0;i < this->y;i++) {
478 for (int j = 0;j < this->x;j++) {
479 tmp.M[i][j] = this->M[j][i];
480 }
481 }
482 this->x = tmp.x;
483 this->y = tmp.y;
484 this->M = tmp.M;
485 }
486 void Turn(const Matrix<T> b) {
487 this->Set(b.y,b.x);
488 for (int i = 0;i < b.y;i++) {
489 for (int j = 0;j < b.x;j++) {
490 this->M[i][j] = b.M[j][i];
491 }
492 }
493 }
494 void Print() {
495 cout.setf(ios::left);
496 for (int i = 0;i < this->x;i++) {
497 for (int j = 0;j < this->y;j++) {
498 cout << setw(12) <<this->M[i][j] << " ";
499 }
500 cout << endl;
501 }
502 cout << endl;
503 }
504 void Set(int a, int b) {
505 this->x = a;
506 this->y = b;
507 this->M.clear();
508 for (int i = 0;i < a;i++) {
509 vector<T> tmp;
510 for (int j = 0;j < b;j++) {
511 tmp.push_back(0);
512 }
513 this->M.push_back(tmp);
514 }
515 }
516 void ExchangeLine(int a,int b) {
517 if (a >= this->x || b >= this->x||a==b) {
518 return;
519 }
520 for (int i = 0;i < this->M[a].size();i++) {
521 this->M[a][i] = this->M[a][i] + this->M[b][i];
522 this->M[b][i] = this->M[a][i] - this->M[b][i];
523 this->M[a][i] = this->M[a][i] - this->M[b][i];
524 }
525 }
526 void LineTimes(int k,int landa) {
527 if (k >= this->x) {
528 return;
529 }
530 for (int i = 0;i < this->M[k].size();i++) {
531 this->M[k][i] *= landa;
532 }
533 }
534 void ColumTimes(int k, int landa) {
535 if (k >= this->y) {
536 return;
537 }
538 for (int i = 0;i < this->M.size();i++) {
539 this->M[i][k] *= landa;
540 }
541 }
542 void LineTimesPlus(int a,int b,int landa) {
543 if (a >= this->x || b >= this->x) {
544 return;
545 }
546 for (int i = 0;i < this->M[a].size();i++) {
547 this->M[a][i] += this->M[b][i]*landa;
548 }
549 }
550 void ColumTimesPlus(int a, int b, int landa) {
551 if (a >= this->y || b >= this->y) {
552 return;
553 }
554 for (int i = 0;i < this->M.size();i++) {
555 this->M[i][a] += this->M[i][b] * landa;
556 }
557 }
558 void AddLine() {
559 this->x++;
560 vector<T> tmp(this->y);
561 this->M.push_back(tmp);
562 }
563 void AddLine(const Vector<T> &a) {
564 if (this->y != a.x) return;
565 this->x++;
566 this->M.push_back(a.V);
567 }
568 void AddColum() {
569 this->y++;
570 for (int i = 0;i < this->x;i++) {
571 this->M[i].push_back(0);
572 }
573 }
574 void AddColum(const Vector<T> &a) {
575 if (this->x != a.x) return;
576 this->y++;
577 for (int i = 0;i < this->x;i++) {
578 this->M[i].push_back(a.V[i]);
579 }
580 }
581 void ExchangeColum(int a, int b) {
582 if (a >= this->y || b >= this->y || a == b) {
583 return;
584 }
585 for (int i = 0;i < this->M.size();i++) {
586 this->M[i][a] = this->M[i][a] + this->M[i][b];
587 this->M[i][b] = this->M[i][a] - this->M[i][b];
588 this->M[i][a] = this->M[i][a] - this->M[i][b];
589 }
590 }
591 void InsLine(int k) {
592 this->AddLine();
593 this->ExchangeLine(k,this->x-1);
594 }
595 void InsLine(const Vector<T> &a,int k) {
596 this->AddLine(a);
597 this->ExchangeLine(k, this->x - 1);
598 }
599 void InsColum(int k) {
600 this->AddColum();
601 this->ExchangeColum(k,this->y-1);
602 }
603 void InsColum(const Vector<T> &a,int k) {
604 this->AddColum(a);
605 this->ExchangeColum(k, this->y - 1);
606 }
607 void SetLine(const Vector<T> &a, int k) {
608 if (k >= this->x||this->y!=a.x) {
609 return;
610 }
611 this->M[k] = a.V;
612 }
613 void SetColum(const Vector<T> &a, int k) {
614 if (k >= this->y || this->x != a.x) {
615 return;
616 }
617 for (int i = 0;i < a.x;i++) {
618 this->M[k][i] = a.V[i];
619 }
620 }
621 void DelLine(int k) {
622 if (k >= this->x) return;
623 this->M.erase(this->M.begin()+k);
624 this->x--;
625 }
626 const int LineSize() {
627 return this->x;
628 }
629 void GetBlock(Matrix<T> a,vector<int> p,vector<int> q) {
630 this->Set(p.size(), q.size());
631 for (int i = 0;i < p.size();i++) {
632 for (int j = 0;j < q.size();j++) {
633 this->M[i][j] = a.M[p[i]][q[j]];
634 }
635 }
636 }
637 void GetLineBlock(const Matrix<T> a,int p,int q) {
638 if (p > q || p >= a.x || q >= a.x) return;
639 this->Set(q-p+1,a.y);
640 for (int i = 0;i < this->x;i++) {
641 for (int j = 0;j < this->y;j++) {
642 this->M[i][j] = a.M[p+i][j];
643 }
644 }
645 }
646 void GetColumBlock(const Matrix<T> a,int p,int q) {
647 if (p > q || p >= a.y || q >= a.y) return;
648 this->Set(a.x,q-p+1);
649 for (int i = 0;i < this->x;i++) {
650 for (int j = 0;j < this->y;j++) {
651 this->M[i][j] = a.M[i][p+j];
652 }
653 }
654 }
655 void FillIn(Matrix<T> &a,int p,int q) {
656 if (this->x + p > a.x || this->y + q > a.y) return;
657 for (int i = 0;i < this->x;i++) {
658 for (int j = 0;j < this->y;j++) {
659 a.M[p+i][q+j] = this->M[i][j];
660 }
661 }
662 }
663 void RightLink(const Matrix<T> a,const Matrix<T> b) {
664 if (a.x != b.x) return;
665 this->Set(a.x,a.y+b.y);
666 int i, j;
667 for (i = 0;i < this->x;i++) {
668 for (j = 0;j < a.y;j++) {
669 this->M[i][j] = a.M[i][j];
670 }
671 for (j = 0;j < b.y;j++) {
672 this->M[i][j+a.y] = b.M[i][j];
673 }
674 }
675 }
676 void DownLink(const Matrix<T> a, const Matrix<T> b) {
677 if (a.y != b.y) return;
678 this->Set(a.x+b.x, a.y);
679 int i, j;
680 for (i = 0;i < this->y;i++) {
681 for (j = 0;j < a.x;j++) {
682 this->M[j][i] = a.M[j][i];
683 }
684 for (j = 0;j < b.x;j++) {
685 this->M[j+a.x][i] = b.M[j][i];
686 }
687 }
688 }
689 void RightLink(const Matrix<T> a) {
690 if (this->x != a.x) return;
691 int i, j;
692 Matrix <T> tmp;
693 tmp.Set(this->x,this->y+a.y);
694 for (i = 0;i < this->x;i++) {
695 for (j = 0;j < this->y;j++) {
696 tmp.M[i][j] = this->M[i][j];
697 }
698 for (int j = 0;j < a.y;j++) {
699 tmp.M[i][tmp.y - a.y + j] = a.M[i][j];
700 }
701 }
702 this->x = tmp.x;
703 this->y = tmp.y;
704 this->M = tmp.M;
705 }
706 void DownLink(const Matrix<T> a) {
707 if (this->y != a.y) return;
708 int i, j;
709 Matrix <T> tmp;
710 tmp.Set(this->x+a.x, this->y);
711 for (i = 0;i < this->y;i++) {
712 for (j = 0;j < this->x;j++) {
713 tmp.M[j][i] = this->M[j][i];
714 }
715 for (int j = 0;j < a.x;j++) {
716 tmp.M[j+tmp.x-a.x][i] = a.M[j][i];
717 }
718 }
719 this->x = tmp.x;
720 this->y = tmp.y;
721 this->M = tmp.M;
722 }
723 void operator >> (vector<T> &data) {
724 data.clear();
725 for (int i = 0;i < this->x;i++) {
726 for (int j = 0;j < this->y;j++) {
727 data.push_back(this->M[i][j]);
728 }
729 }
730 }
731 void operator <<(const vector<T> data) {
732 this->Zero();
733 int pos = 0;
734 for (int i = 0;i < this->x;i++) {
735 for (int j = 0;j < this->y;j++) {
736 if (pos >= data.size()) return;
737 this->M[i][j] = data[pos++];
738 }
739 }
740 }
741 void DiagPlus(T landa) {
742 if (this->x != this->y) return;
743 for (int i = 0;i < this->x;i++) {
744 this->M[i][i] += landa;
745 }
746 }
747 T Trace() {
748 if (this->x != this->y) return 0;
749 T result = 0;
750 for (int i = 0;i < this->x;i++) {
751 result += this->M[i][i];
752 }
753 return result;
754 }
755 T AbsMax(int &a,int &b) {
756 T result=0;
757 for (int i = 0;i < this->x;i++) {
758 for (int j = 0;j < this->y;j++) {
759 if (abs(this->M[i][j]) > result) {
760 result = abs(this->M[i][j]);
761 a = i;
762 b = j;
763 }
764 }
765 }
766 return result;
767 }
768 T AbsMin(int &a, int &b) {
769 T result = 1e9;
770 for (int i = 0;i < this->x;i++) {
771 for (int j = 0;j < this->y;j++) {
772 if (abs(this->M[i][j]) < result) {
773 result = abs(this->M[i][j]);
774 a = i;
775 b = j;
776 }
777 }
778 }
779 return result;
780 }
781 T Max(int &a, int &b) {
782 T result = -1e9;
783 for (int i = 0;i < this->x;i++) {
784 for (int j = 0;j < this->y;j++) {
785 if (this->M[i][j] > result) {
786 result = this->M[i][j];
787 a = i;
788 b = j;
789 }
790 }
791 }
792 return result;
793 }
794 void CoVar(Matrix<T> a) {
795 int i, j;
796 this->Set(a.x,a.x);
797 vector<T> avgC(a.x);
798 for (i = 0;i < a.x;i++) {
799 T tmp=0;
800 for (j = 0;j < a.y;j++) {
801 tmp += a.M[i][j];
802 }
803 avgC[i] = tmp/a.y;
804 }
805 for (i = 0;i < a.x;i++) {
806 for (j = 0;j < a.y;j++) {
807 a.M[i][j] -= avgC[i];
808 }
809 }
810 for (i = 0;i < a.x;i++) {
811 for (j = 0;j < a.x;j++) {
812 T t=0;
813 for (int k = 0;k < a.y;k++) {
814 t += a.M[i][k] * a.M[j][k];
815 }
816 this->M[i][j] = t/(this->x-1);
817 }
818 }
819 }
820 T Min(int &a, int &b) {
821 T result = 1e9;
822 for (int i = 0;i < this->x;i++) {
823 for (int j = 0;j < this->y;j++) {
824 if (this->M[i][j] < result) {
825 result = this->M[i][j];
826 a = i;
827 b = j;
828 }
829 }
830 }
831 return result;
832 }
833 T Mean() {
834 T result = 0;
835 for (int i = 0;i < this->x;i++) {
836 for(int j=0;j<this->y;j++)
837 result += this->M[i][j];
838 }
839 return result / (this->x*this->y);
840 }
841 void operator +=(T landa) {
842 for (int i = 0;i < this->x;i++) {
843 for (int j = 0;j < this->y;j++) {
844 this->M[i][j] += landa;
845 }
846 }
847 }
848 void FillIn(Matrix<T> &a,Dir d) {
849 int i, j;
850 if (this->x > a.x || this->y > a.y) return;
851 switch (d) {
852 case LT:
853 for (i = 0;i < this->x;i++) {
854 for (j = 0;j < this->y;j++) {
855 a.M[i][j] = this->M[i][j];
856 }
857 }
858 break;
859 case RT:
860 for (i = 0;i < this->x;i++) {
861 for (j = 0;j < this->y;j++) {
862 a.M[i][j+a.y-this->y] = this->M[i][j];
863 }
864 }
865 break;
866 case LB:
867 for (i = 0;i < this->x;i++) {
868 for (j = 0;j < this->y;j++) {
869 a.M[i+a.x-this->x][j] = this->M[i][j];
870 }
871 }
872 break;
873 case RB:
874 for (i = 0;i < this->x;i++) {
875 for (j = 0;j < this->y;j++) {
876 a.M[i + a.x - this->x][j + a.y - this->y] = this->M[i][j];
877 }
878 }
879 break;
880 }
881 }
882 void GetBlock(const Matrix<T> a,int p,int q,Dir d) {
883 if (p >= a.x || q >= a.y) return;
884 this->Set(p,q);
885 switch (d) {
886 case LT:
887 for (int i = 0;i < p;i++) {
888 for (int j = 0;j < q;j++) {
889 this->M[i][j] = a.M[i][j];
890 }
891 }
892 break;
893 case RT:
894 for (int i = 0;i < p;i++) {
895 for (int j = a.y-q;j < a.y;j++) {
896 this->M[i][j-a.y+q] = a.M[i][j];
897 }
898 }
899 break;
900 case LB:
901 for (int i = a.x-p;i < a.x;i++) {
902 for (int j = 0;j <q;j++) {
903 this->M[i-a.x+p][j] = a.M[i][j];
904 }
905 }
906 break;
907 case RB:
908 for (int i = a.x-p;i < a.x;i++) {
909 for (int j = a.y - q;j < a.y;j++) {
910 this->M[i-a.x+p][j-a.y+q] = a.M[i][j];
911 }
912 }
913 break;
914 default:
915 return;
916 }
917
918
919 }
920 const int ColumSize() {
921 return this->y;
922 }
923 void DelColum(int k) {
924 if (k >= this->y) return;
925 for (int i = 0;i < this->x;i++) {
926 this->M[i].erase(this->M[i].begin()+k);
927 }
928 this->y--;
929 }
930 void Set(int a, int b,T value) {
931 if (a >= this->x || b >= this->y) {
932 return;
933 }
934 this->M[a][b] = value;
935 }
936 T Get(int a,int b) {
937 if (a >= this->x || b >= this->y) {
938 return 0;
939 }
940 return this->M[a][b];
941 }
942 Matrix<T> operator -(Matrix<T> b) {
943 Matrix<T> result;
944 if (this->x != b.x || this->y != b.y) return result;
945 result.Set(this->x, this->y);
946 for (int i = 0;i < this->x;i++) {
947 for (int j = 0;j < this->y;j++) {
948 result.M[i][j] = this->M[i][j] - b.M[i][j];
949 }
950 }
951 return result;
952 }
953 Matrix<T> operator +(Matrix<T> b) {
954 Matrix<T> result;
955 if (this->x != b.x || this->y != b.y) return result;
956 result.Set(this->x, this->y);
957 for (int i = 0;i < this->x;i++) {
958 for (int j = 0;j < this->y;j++) {
959 result.M[i][j] = this->M[i][j] + b.M[i][j];
960 }
961 }
962 return result;
963 }
964
965 };
966 template <typename T>
967 class Vector {
968 public:
969 vector<T> V;
970 int x;
971 Vector(const Vector<T> &a){
972 this->x = a.x;
973 this->V = a.V;
974 }
975 Vector() {
976 this->x = 0;
977 this->V.clear();
978 }
979 Vector(int n) {
980 this->Set(n);
981 }
982 void Set(int n) {
983 V.clear();
984 for (int i = 0;i < n;i++) {
985 this->V.push_back(0);
986 }
987 this->x = n;
988 }
989 void Print() {
990 int i;
991 for (i = 0;i < this->x;i++) {
992 cout << setw(12) << this->V[i] << " ";
993 if ((i + 1) % 5 == 0) cout << endl;
994 }
995 if(i%5!=0)
996 cout << endl;
997 cout << endl;
998 }
999 void Init(vector<T> data) {
1000 this->V.clear();
1001 this->Set(data.size());
1002 for (int i = 0;i < data.size();i++) {
1003 this->V[i] = data[i];
1004 }
1005 }
1006 void Unit(int n,int pos) {
1007 this->Set(n);
1008 this->V[pos] = 1;
1009 }
1010 void Const(T landa) {
1011 for (int i = 0;i < this->x;i++) {
1012 this->V[i] = landa;
1013 }
1014 }
1015 void Cut(T a, T b, int num) {
1016 this->Set(num+1);
1017 T interval = (b-a) / num;
1018 for (int i = 0;i < this->x;i++) {
1019 this->V[i] = i*interval;
1020 }
1021 }
1022 void Cal(T f(T x),const Vector<T> &xs) {
1023 this->Set(xs.x);
1024 for (int i = 0;i < this->x;i++) {
1025 this->V[i] = f(xs.V[i]);
1026 }
1027 }
1028 void Random(T a,T b) {
1029 if (b < a) return;
1030 time_t t = time(NULL);
1031 srand(t*rand());
1032 for (int i = 0;i < this->x;i++) {
1033 double t = rand() / (RAND_MAX + 0.0);
1034 this->V[i] = (T)(t*(b - a) + a);
1035 }
1036 }
1037 void RandomInt(int n) {
1038 int i;
1039 time_t t = time(NULL);
1040 srand(t * rand());
1041 for (i = 0;i < this->x;i++) {
1042 double t = rand() / (RAND_MAX + 0.0);
1043 this->V[i] = (int)(t * 2 * n - n);
1044 }
1045 }
1046 double GaussRand(double mu,double sigma) {
1047 const double epsilon = std::numeric_limits<double>::min();
1048 const double two_pi = 2.0*3.14159265358979323846;
1049 static double z0, z1;
1050 static bool generate;
1051 generate = !generate;
1052 if (!generate)
1053 return z1 * sigma + mu;
1054 double u1, u2;
1055 do
1056 {
1057 u1 = rand() * (1.0 / RAND_MAX);
1058 u2 = rand() * (1.0 / RAND_MAX);
1059 } while (u1 <= epsilon);
1060 z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
1061 z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
1062 return z0 * sigma + mu;
1063 }
1064 double GaussRand() {
1065 return this->RandomGauss(0,1);
1066 }
1067 //Generate Gauss random number
1068 void RandomGauss(double mu, double sigma) {
1069 for (int i = 0;i < this->x;i++) {
1070 double number = GaussRand(mu,sigma);
1071 this->V[i] = number;
1072 }
1073 }
1074 //Generate Gauss random number(0,1)
1075 void RandomGauss() {
1076 for (int i = 0;i < this->x;i++) {
1077 double number = GaussRand();
1078 this->V[i] = number;
1079 }
1080 }
1081 //Generate Gamma random number
1082 void RandomGamma(double alpha, double lambda){
1083 std::default_random_engine generator;
1084 std::tr1::gamma_distribution<double> distribution(alpha, lambda);
1085 for (int i = 0;i < this->x;i++) {
1086 double number = distribution(generator);
1087 this->V[i] = number;
1088 }
1089 }
1090 //Generate Exponential random number
1091 double ExponentialRand(double lambda) {
1092 double pV = 0.0;
1093 while (true)
1094 {
1095 pV = (double)rand() / (double)RAND_MAX;
1096 if (pV != 1)
1097 {
1098 break;
1099 }
1100 }
1101 pV = (-1.0 / lambda)*log(1 - pV);
1102 return pV;
1103 }
1104 void RandomExponential(double alpha) {
1105 for (int i = 0;i < this->x;i++) {
1106 double number = ExponentialRand(alpha);
1107 this->V[i] = number;
1108 }
1109 }
1110 double randomUniform() {
1111 return rand()*0.1 / RAND_MAX;
1112 }
1113 double BetaRand(double alpha, double beta){
1114 /*Johnk's beta generator*/
1115 double u, v;
1116 double x, y;
1117 do{
1118 u = randomUniform();
1119 v = randomUniform();
1120 x = pow(u, 1 / alpha);
1121 y = pow(v, 1 / beta);
1122 } while (x + y>1);
1123 return x / (x + y);
1124 }
1125 //Generate Beta random number
1126 void RandomBeta(double alpha, double beta) {
1127 for (int i = 0;i < this->x;i++) {
1128 double number = BetaRand(alpha,beta);
1129 this->V[i] = number;
1130 }
1131 }
1132 void Frequency(const Vector<T> &a) {
1133 if (a.x == 0) return;
1134 this->Set(10);
1135 int j;
1136 T min = a.V[0], max = a.V[0];
1137 for (j = 1;j < a.x;j++) {
1138 if (a.V[j] > max) {
1139 max = a.V[j];
1140 }
1141 if (a.V[j] < min) {
1142 min = a.V[j];
1143 }
1144 }
1145 max += (max-min) / 10;;
1146 T interval = (max - min) / 10;
1147 for (j = 0;j < a.x;j++) {
1148 this->V[int((a.V[j] - min) / interval)] += 1;
1149 }
1150 for (j = 0;j < this->x;j++) {
1151 this->V[j] = this->V[j] / a.x;
1152 }
1153 }
1154 Vector<T> operator +(const Vector<T> &b) {
1155 Vector<T> result;
1156 if (this->x != b.x) return result;
1157 result.Set(this->x);
1158 for (int i = 0;i < this->x;i++) {
1159 result.V[i] = this->V[i] + b.V[i];
1160 }
1161 return result;
1162 }
1163
1164 void operator +=(const Vector<T> &b) {
1165 if (this->x != b.x) return;
1166 for (int i = 0;i < this->x;i++) {
1167 this->V[i] += b.V[i];
1168 }
1169 }
1170 Vector<T> operator -(const Vector<T> &b) {
1171 Vector<T> result;
1172 if (this->x != b.x) return result;
1173 result.Set(this->x);
1174 for (int i = 0;i < this->x;i++) {
1175 result.V[i] = this->V[i] - b.V[i];
1176 }
1177 return result;
1178 }
1179
1180 void operator -=(const Vector<T> &b) {
1181 if (this->x != b.x) return;
1182 for (int i = 0;i < this->x;i++) {
1183 this->V[i] -= b.V[i];
1184 }
1185 }
1186 T operator *(const Vector<T> &a) {
1187 T result = 0;
1188 if (this->x != a.x) return result;
1189 for (int i = 0;i < this->x;i++) {
1190 result += this->V[i] * a.V[i];
1191 }
1192 return result;
1193 }
1194 void Mult(const Matrix<T> &a,T landa) {
1195 this->Set(a.x);
1196 for (int i = 0;i < this->x;i++) {
1197 this->V[i] = a.V[i] * landa;
1198 }
1199 }
1200 void operator *=(T landa) {
1201 for (int i = 0;i < this->x;i++) {
1202 this->V[i] = this->V[i] * landa;
1203 }
1204 }
1205 void operator +=(T landa) {
1206 for (int i = 0;i < this->x;i++) {
1207 this->V[i] += landa;
1208 }
1209 }
1210 void TimesPlus(const Matrix<T> &a, T landa) {
1211 if (this->x != a.x) return;
1212 for (int i = 0;i < this->x;i++) {
1213 this->V[i] += a.V[i] * landa;
1214 }
1215 }
1216 void Zero() {
1217 for (int i = 0;i < this->x;i++) {
1218 this->V[i] = 0;
1219 }
1220 }
1221 void DelError(T error) {
1222 for (int i = 0;i < this->x;i++) {
1223 if (abs(this->V[i]) < error) {
1224 this->V[i] = 0;
1225 }
1226 }
1227 }
1228 void Add(T landa) {
1229 this->V.push_back(landa);
1230 this->x++;
1231 }
1232 void Ins(T landa,int k) {
1233 this->V.insert(this->V.begin()+k,landa);
1234 this->x++;
1235 }
1236 void Del(int k) {
1237 if (k >= this->x) return;
1238 this->V.erase(this->V.begin()+k);
1239 this->x--;
1240 }
1241 int Size() {
1242 return this -> x;
1243 }
1244 void Link(const Vector<T> &a,const Vector<T> &b) {
1245 this->Set(a.x+b.x);
1246 int i;
1247 for (i = 0;i < a.x;i++) {
1248 this->V[i] = a.V[i];
1249 }
1250 for (i = a.x;i < a.x + b.x;i++) {
1251 this->V[i] = b.V[i-a.x];
1252 }
1253 }
1254 void Link(const Vector<T> &a) {
1255 for (int i = 0;i < a.x;i++) {
1256 this->V.push_back(a.V[i]);
1257 }
1258 this->x=a.x + this->x;
1259 }
1260 void GetLeft(const Vector<T> &a,int k) {
1261 if (k > a.x) return;
1262 this->Set(k);
1263 for (int i = 0;i < k;i++) {
1264 this->V[i] = a.V[i];
1265 }
1266 }
1267 void GetMid(const Vector<T> &a,int pos,int k) {
1268 if (pos + k > a.x) return;
1269 this->Set(k);
1270 for (int i = 0;i < k;i++) {
1271 this->V[i] = a.V[pos+i];
1272 }
1273 }
1274 void GetRight(const Vector<T> &a,int k) {
1275 if (k > a.x) return;
1276 this->Set(k);
1277 for (int i = 0;i < k;i++) {
1278 this->V[i] = a.V[a.x-k+i];
1279 }
1280 }
1281 void FillIn(const Vector<T> &a,int k) {
1282 if (k >= this->x) return;
1283 for (int i = 0;i < a.x&&k + i < this->x;i++) {
1284 this->V[k+i] = a.V[i];
1285 }
1286 }
1287 void operator >> (vector<T> &a) {
1288 for (int i = 0;i < this->x;i++) {
1289 a.push_back(this->V[i]);
1290 }
1291 }
1292 void operator <<(const vector<T> &a) {
1293 this->Set(a.size());
1294 for (int i = 0;i < this->x;i++) {
1295 this->V[i] = a[i];
1296 }
1297 }
1298 T Mean() {
1299 T result = 0;
1300 for (int i = 0;i < this->x;i++) {
1301 result += this->V[i];
1302 }
1303 return result/this->x;
1304 }
1305 T Var() {
1306 T avg = this->Mean();
1307 T result = 0;
1308 for (int i = 0;i < this->x;i++) {
1309 result += (this->V[i]-avg)*(this->V[i] - avg);
1310 }
1311 return result / this->x;
1312 }
1313 T AbsMax() {
1314 T result = 0;
1315 for (int i = 0;i < this->x;i++) {
1316 T tmp = abs(this->V[i]);
1317 if (result < tmp) {
1318 result = tmp;
1319 }
1320 }
1321 return result;
1322 }
1323 T AbsMax(int &k) {
1324 T result = 0;
1325 for (int i = 0;i < this->x;i++) {
1326 T tmp = abs(this->V[i]);
1327 if (result < tmp) {
1328 result = tmp;
1329 k = i;
1330 }
1331 }
1332 return result;
1333 }
1334 T AbsMin() {
1335 T result = 0;
1336 for (int i = 0;i < this->x;i++) {
1337 T tmp = abs(this->V[i]);
1338 if (result > tmp) {
1339 result = tmp;
1340 }
1341 }
1342 return result;
1343 }
1344 T AbsMin(int &k) {
1345 T result = 0;
1346 for (int i = 0;i < this->x;i++) {
1347 T tmp = abs(this->V[i]);
1348 if (result > tmp) {
1349 result = tmp;
1350 k = i;
1351 }
1352 }
1353 return result;
1354 }
1355 T Min() {
1356 T result = 0;
1357 for (int i = 0;i < this->x;i++) {
1358 T tmp = this->V[i];
1359 if (result > tmp) {
1360 result = tmp;
1361 }
1362 }
1363 return result;
1364 }
1365 T Max() {
1366 T result = 0;
1367 for (int i = 0;i < this->x;i++) {
1368 T tmp = this->V[i];
1369 if (result < tmp) {
1370 result = tmp;
1371 }
1372 }
1373 return result;
1374 }
1375 T Min(int &k) {
1376 T result = 0;
1377 for (int i = 0;i < this->x;i++) {
1378 T tmp = this->V[i];
1379 if (result > tmp) {
1380 result = tmp;
1381 k = i;
1382 }
1383 }
1384 return result;
1385 }
1386 T Max(int &k) {
1387 T result = 0;
1388 for (int i = 0;i < this->x;i++) {
1389 T tmp = this->V[i];
1390 if (result < tmp) {
1391 result = tmp;
1392 k = i;
1393 }
1394 }
1395 return result;
1396 }
1397 void SmallToBig() {
1398 for (int i = 0;i < this->x;i++) {
1399 for(int j = 0;j < this->x - i - 1;j++){
1400 if (this->V[j + 1] < this->V[j]) {
1401 this->V[j + 1] = this->V[j + 1] + this->V[j];
1402 this->V[j] = this->V[j + 1] - this->V[j];
1403 this->V[j + 1] = this->V[j + 1] - this->V[j];
1404 }
1405 }
1406 }
1407 }
1408 void BigToSmall() {
1409 for (int i = 0;i < this->x;i++) {
1410 for (int j = 0;j < this->x - i - 1;j++) {
1411 if (this->V[j + 1] > this->V[j]) {
1412 this->V[j + 1] = this->V[j + 1] + this->V[j];
1413 this->V[j] = this->V[j + 1] - this->V[j];
1414 this->V[j + 1] = this->V[j + 1] - this->V[j];
1415 }
1416 }
1417 }
1418 }
1419 int Local(T landa) {
1420 int result = 0;
1421 while (landa > this->V[result]) {
1422 result++;
1423 }
1424 return result;
1425 }
1426 T Factorial(T n) {
1427 T result = 1;
1428 while (n > 0) {
1429 result *= n;
1430 n--;
1431 }
1432 return result;
1433 }
1434 void Com(T n) {
1435 this->Set(n+1);
1436 for (int i = 0;i <= n;i++) {
1437 this->V[i] = (Factorial(n))/((Factorial(n-i))*(Factorial(i)));
1438 }
1439 }
1440 void Invers() {
1441 for (int i = 0;i < this->x / 2;i++) {
1442 this->V[i] = this->V[i] + this->V[this->x-1-i];
1443 this->V[this->x - 1 - i] = this->V[i] - this->V[this->x - 1 - i];
1444 this->V[i] = this->V[i] - this->V[this->x - 1 - i];
1445 }
1446 }
1447 void operator <<(int k) {
1448 if (k >= this->x) this->Zero();
1449 if (k <= 0) return;
1450 int i;
1451 for (i = 0;i < this->x - k;i++) {
1452 this->V[i] = this->V[i+k];
1453 }
1454 for (;i < this->x;i++) {
1455 this->V[i] = 0;
1456 }
1457 }
1458 void operator >>(int k) {
1459 if (k >= this->x) this->Zero();
1460 if (k <= 0) return;
1461 int i;
1462 for (i=this->x-k-1;i>=0;i--) {
1463 this->V[i+k] = this->V[i];
1464 }
1465 for (i=0;i < k;i++) {
1466 this->V[i] = 0;
1467 }
1468 }
1469 void GetLine(const Matrix<T> &a,int k) {
1470 if (k >= a.x) {
1471 return;
1472 }
1473 this->V = a.M[k];
1474 this->x = a.M[k].size();
1475 }
1476 void GetColum(const Matrix<T> &a,int k) {
1477 if (k >= a.y) {
1478 return;
1479 }
1480 this->Set(a.x);
1481 for (int i = 0;i < this->x;i++) {
1482 this->V[i] = a.M[k][i];
1483 }
1484 }
1485 void Mult(const Matrix<T> &a,const Vector<T> &b) {
1486 if (a.y != b.x) return;
1487 this->Set(a.x);
1488 for (int i = 0;i < this->x;i++) {
1489 T tmp = 0;
1490 for (int j = 0;j < a.y;j++) {
1491 tmp += a.M[i][j]*b.V[j];
1492 }
1493 this->V[i] = tmp;
1494 }
1495 }
1496 void Mult(const Vector<T> &b, const Matrix<T> &a) {
1497 if (a.x != b.x) return;
1498 this->Set(a.y);
1499 for (int i = 0;i < this->x;i++) {
1500 T tmp = 0;
1501 for (int j = 0;j < a.x;j++) {
1502 tmp += a.M[i][j] * b.V[j];
1503 }
1504 this->V[i] = tmp;
1505 }
1506 }
1507 void Mult(const Vector<T> &a,const Matrix<T> &b, const Vector<T> &c) {
1508 if (a.x != b.x || b.y != c.x) return;
1509 this->Mult(a,b);
1510 T result = 0;
1511 for (int i = 0;i < c.x;i++) {
1512 result += this->V[i] * c.V[i];
1513 }
1514 return result;
1515 }
1516 void GetDiag(const Matrix<T> &a) {
1517 if (a.x != a.y) return;
1518 this->Set(a.x);
1519 for (int i = 0;i < this->x;i++) {
1520 this->V[i] = a.M[i][i];
1521 }
1522 }
1523 void SetDiag(Matrix<T> &a) {
1524 if (a.x != a.y || this->x != a.x) return;
1525 for (int i = 0;i < this->x;i++) {
1526 a.M[i][i] = this->V[i];
1527 }
1528 }
1529
1530 };
1531
1532 template <typename T>
1533 double ED(const Vector<T> &a,const Vector<T> &b) {
1534 double result = 0;
1535 if (a.x != b.x) return;
1536 for (int i = 0;i < a.x;i++) {
1537 result += (a.V[i]-b.V[i])*(a.V[i] - b.V[i]);
1538 }
1539 return sqrt(result);
1540 }
1541
1542 #endif
