Facebook
From Lousy Lemur, 6 Years ago, written in C++.
Embed
Download Paste or View Raw
Hits: 244
  1. int n = 100;
  2. int minRange = 0;
  3. int maxRange = 1000000000;
  4.  
  5. LIB "poly.lib";
  6. ring r = 0, x, dp;
  7. int i, j;
  8. poly zero = 0;
  9. poly mlt, p, ratio;
  10. matrix m[n][n];
  11.  
  12. proc generateMatrix() {
  13.         for(i = 1; i <= n; i++) {
  14.                 for(j = 1; j <= n;  j++) {
  15.                         m[i, j] = random(minRange, maxRange);
  16.                 }
  17.         }
  18. }
  19.  
  20. proc getDiagonal() {
  21.         poly res = 1;
  22.         for(i = 1; i <= n; i++)  {
  23.                 res = res * m[i, i];
  24.         }
  25.         return(res[1][1] / ratio);
  26. }
  27.  
  28. proc getNotZero(int cl) {
  29.         for(i = cl; i <= n; i++) {
  30.                 if(m[i, cl] != zero) {
  31.                         return(i);
  32.                 }
  33.         }
  34.         return(-1);
  35. }
  36.  
  37. proc changeRow(int first, int second) {
  38.         poly tmp;
  39.         for(i = first; i <= n; i++) {
  40.                 tmp = m[first, i];
  41.                 m[first, i] = m[second, i];
  42.                 m[second, i] = tmp;
  43.         }
  44. }
  45.  
  46. proc getLCM(int col) {
  47.         int res = 1;
  48.         for(i = col; i <= n; i++) {
  49.                 res = lcm(res, int(m[i, col]));
  50.         }
  51.         return(res);
  52. }
  53.  
  54. proc gaussMatrix() {
  55.         int i, j, col, row;
  56.         ratio = 1;
  57.         mlt = 1;
  58.         p = 1;
  59.         for(col = 1; col <= n; col++) {
  60.                 row = getNotZero(col);
  61.                 if(row == -1) {
  62.                         ratio = 0;
  63.                         return;
  64.                 } else {
  65.                         if(col != row) {
  66.                                 changeRow(col, row);
  67.                                 ratio = ratio * (-1);
  68.                         }
  69.                         mlt = getLCM(col);
  70.                         p = mlt / m[col, col];
  71.                         ratio = ratio * p;
  72.                         m = multrow(m, col, p);
  73.                         for(j = col + 1; j <= n; j++) {
  74.                                 if(m[j, col] != zero) {
  75.                                         p = mlt / m[j, col];
  76.                                         ratio = ratio * p;
  77.                                                 for(i = col; i <= n; i++) {
  78.                                                 m[j, i] = m[j, i] * p - m[col, i];
  79.                                         }
  80.                                 }
  81.                         }
  82.                 }
  83.         }
  84. }
  85.  
  86. proc solve() {
  87.         gaussMatrix();
  88.         if(ratio == 0) {
  89.                 return(0);
  90.         } else {
  91.                 return(getDiagonal());
  92.         }
  93. }
  94.  
  95.  
  96. generateMatrix();
  97. det(m);
  98. solve();
  99. exit;
  100.