matrix
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