home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / yamp / virtop.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  40.5 KB  |  1,647 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21. /////////////////// matrix functions ---- non-member functions
  22.  
  23. int Setwid(int w)
  24. {
  25.     static int wid = 6;
  26.     wid = (w < 0) ? wid : w;
  27.     return wid;
  28. }
  29. int Setdec(int d)
  30. {
  31.     static int dec = 3;
  32.     dec = (d < 0) ? dec : d;
  33.     if (d >= Getwid()) dec = Getwid() - 1;
  34.     dec = (dec <= 0) ? 0 : dec;
  35.     return dec;
  36. }
  37. int Getwid(void)
  38. {
  39.     static int wid = 6;
  40.     wid = Setwid(- 1);
  41.     return wid;
  42. }
  43. int Getdec(void)
  44. {
  45.     static int dec = 3;
  46.     dec = Setdec(- 1);
  47.     return dec;
  48. }
  49.  
  50.  
  51. VMatrix &Reada(char *fid)     /* this reads an ascii matrix */
  52. {
  53.     unsigned int chh[MAXCHARS], *chc;
  54.     int i = 0, j = 0, x1 = 0, x2 = 0;
  55.     FILE *file;
  56.     char mname[MAXCHARS];
  57.     float x = 0;
  58.     
  59.     
  60.     file = fopen(fid, "r");
  61.     if (!file) Dispatch->Nrerror("READA: error while opening file");
  62.     
  63.     for (i = 0; i < 80; i++) chh[i] = 0;
  64.     
  65.     i = 0;
  66.     
  67.     do {     /* read until first carriage return */
  68.         j = fgetc(file);
  69.         chh[i] = j;
  70.         i++;
  71.     } while (j != ((int) 0x0A) && !ferror(file) && (i < 80));
  72.     
  73.     if (ferror(file)) Dispatch->Nrerror(
  74.             "READA: error while reading matrix name" );
  75.  
  76.  
  77.     chh[i] = (int) '\0';     /* null terminate the string */
  78.  
  79.     if (chh[0] != (int) 0x0A) {
  80.         strcpy(mname, ((const char *) chh));
  81.         chc = chh;
  82.         for (j = 1; j < i - 1; j++) strcat(mname, ((const char *) (++ chc)));
  83.     }
  84.     else {
  85.         strcpy(mname, " ");
  86.     }
  87.  
  88.     fscanf(file, "%d%d", &x1, &x2);
  89.  
  90.     VMatrix *mat = new VMatrix(mname, x1, x2);
  91.     for (i = 1; i <= x1; i++) {
  92.         for (j = 1; j <= x2; j++) {
  93.             fscanf(file, "%f", &x);
  94.             mat->M(i, j) = (double) x;
  95.         }
  96.     }
  97.     fclose(file);
  98.     Dispatch->Push(*mat);
  99.     delete mat;
  100.     return Dispatch->ReturnMat();
  101. }     /* reada */
  102.  
  103. VMatrix &Readb(char *fid)     /* this reads an binary matrix */
  104. {
  105.     int i = 0, j = 0, x1 = 0, x2 = 0;
  106.     FILE *file;
  107.     char name[MAXCHARS];
  108.     double x = 0;
  109.  
  110.     file = fopen(fid, "rb");
  111.     if (!file) Dispatch->Nrerror("READB: error while opening file");
  112.  
  113.     fread(&i, sizeof (int), 1, file);
  114.     fgets(name, (i + 1), file);
  115.     fread(&x1, sizeof (int), 1, file);
  116.     fread(&x2, sizeof (int), 1, file);
  117.  
  118.     VMatrix *mat = new VMatrix(name, x1, x2);
  119.     for (i = 1; i <= x1; i++) {
  120.         for (j = 1; j <= x2; j++) {
  121.             fread(&x, sizeof (double), 1, file);
  122.             mat->M(i, j) = x;
  123.         }
  124.     }
  125.     fclose(file);
  126.     Dispatch->Push(*mat);
  127.     delete mat;
  128.     return Dispatch->ReturnMat();
  129. }     /* readb */
  130.  
  131. void Writea(char *fid, VMatrix &mat)     /* write a matrix in an ascii file */
  132. {
  133.     int i, j;
  134.     FILE *file;
  135.  
  136.     file = fopen(fid, "w");
  137.     if (!file) Dispatch->Nrerror("WRITEA: unable to open file");
  138.  
  139.     mat.Garbage("Writea");
  140.  
  141.     fprintf(file, "%s\n", mat.StringAddress());
  142.     fprintf(file, "%d %d \n", mat.r, mat.c);
  143.     for (i = 1; i <= mat.r; i++) {
  144.         for (j = 1; j <= mat.c; j++)
  145.             fprintf(file, "%lf ", mat.m(i, j));
  146.         fprintf(file, "\n");
  147.     }
  148.     fclose(file);
  149. }     /* writea */
  150.  
  151. void heapsort(VMatrix &a)
  152. {
  153.     int n = a.r;
  154.     int s0, s1, s2, s3, s4, s5, s4m1, s5m1;
  155.     double t1, ts;
  156.     s0 = s1 = n;     /* Shell-Metzger sort, PC-World, May 1986 */
  157.     s1 /= 2;
  158.     while (s1 > 0) {
  159.         s2 = s0 - s1;
  160.         for (s3 = 1; s3 <= s2; s3++) {
  161.             s4 = s3;
  162.             while (s4 > 0) {
  163.                 s5 = s4 + s1;
  164.                 s4m1 = s4;
  165.                 s5m1 = s5;
  166.                 if (a.m(s4m1, 1) > a.m(s5m1, 1)) {
  167.                     t1 = a.m(s4m1, 1);
  168.                     a.M(s4m1, 1) = a.m(s5m1, 1);
  169.                     a.M(s5m1, 1) = t1;
  170.                     ts = a.m(s4m1, 2);
  171.                     a.M(s4m1, 2) = a.m(s5m1, 2);
  172.                     a.M(s5m1, 2) = ts;
  173.                     s4 = s4 - s1;
  174.                 }
  175.                 else {
  176.                     s4 = 0;
  177.                 }
  178.             }
  179.         }
  180.         s1 /= 2;
  181.     }
  182. }
  183.  
  184.  
  185.  
  186. VMatrix &MSort(VMatrix &a, int col, int order)
  187. {
  188.     a.Garbage();
  189.     VMatrix *sorted = new VMatrix("sorted(", a.r, a.c);
  190.     if (!order) {
  191.         // sort column col, then reorder other rows
  192.         VMatrix *temp = new VMatrix("temp", a.r, 2);
  193.         for (int i = 1; i <= a.r; i++) {
  194.             temp->M(i, 1) = a.m(i, col);
  195.             temp->M(i, 2) = (double) i;
  196.         }
  197.         heapsort(*temp);
  198.         for (i = 1; i <= a.r; i++) {
  199.             for (int j = 1; j <= a.c; j++)
  200.                 sorted->M(i, j) = a.m(((int) temp->m(i, 2)), j);
  201.         }
  202.         delete temp;
  203.     }
  204.     else {
  205.         // sort row col, then reorder other columns
  206.         VMatrix *temp = new VMatrix("temp", a.c, 2);
  207.         for (int i = 1; i <= a.c; i++) {
  208.             temp->M(i, 1) = a.m(col, i);
  209.             temp->M(i, 2) = (double) i;
  210.         }
  211.         heapsort(*temp);
  212.         for (i = 1; i <= a.c; i++) {
  213.             for (int j = 1; j <= a.r; j++)
  214.                 sorted->M(j, i) = a.m(j, ((int) temp->m(i, 2)));
  215.         }
  216.         delete temp;
  217.     }
  218.     strtype name = sorted->Getname(*sorted);
  219.     name = name + a.Getname(a) + ")";
  220.     sorted->Nameit(name);
  221.     Dispatch->Push(*sorted);
  222.     delete sorted;
  223.     return Dispatch->ReturnMat();
  224. }
  225.  
  226.  
  227.  
  228.  
  229.  
  230. VMatrix &Submat(VMatrix &a, int r, int c, int lr, int lc)
  231. /* extract a submatrix of a into b*/
  232. {
  233.     a.Garbage("Submat");
  234.  
  235.     if ((r > a.r) || (c > a.c))
  236.         a.Nrerror("SUBMAT: index out of range");
  237.  
  238.     VMatrix *b = new VMatrix("b", r - lr + 1, c - lc + 1);
  239.     strtype temp = a.Getname(a);
  240.     b->Nameit(temp);
  241.     int r2 = r - lr + 1;
  242.     int c2 = c - lc + 1;
  243.     for (int i = 1; i <= r2; i++) {
  244.         for (int j = 1; j <= c2; j++)
  245.             b-> M(i, j) = a.m(lr - 1+ i, lc - 1+ j);
  246.     }
  247.     Dispatch->Push(*b);
  248.     delete b;
  249.     return Dispatch->ReturnMat();
  250. }     /* submat */
  251.  
  252.  
  253. VMatrix &Ch(VMatrix &a, VMatrix &b)     /* horizontal concatination */
  254. {
  255.     a.Garbage("Ch");
  256.     b.Garbage("Ch");
  257.     int adim = a.c;
  258.     int bdim = b.c;
  259.     int k = adim + bdim;
  260.  
  261.     if (a.r != b.r)
  262.         a.Nrerror("CH: matrices have different number of rows");
  263.  
  264.     VMatrix *c = new VMatrix("(", a.r, k);
  265.     if (!c) a.Nrerror("CH: failed to create temp storage");
  266.  
  267.     strtype mname = c->Getname(*c) + a.Getname(a) + "||" + b.Getname(b) +")";
  268.     c->Nameit(mname);
  269.  
  270.     for (int i = 1; i <= a.r; i++) {
  271.         for (int j = 1; j <= a.c; j++)
  272.             c->M(i, j) = a.m(i, j);
  273.     }
  274.     int ii;
  275.     for (i = 1, ii = 1; i <= b.r; i++, ii++) {
  276.         for (int j = 1, jj = 1; j <= b.c; j++, jj++)
  277.             c->M(ii, jj + a.c) = b.m(i,j);
  278.     }
  279.     Dispatch->Push(*c);
  280.     delete c;
  281.     return Dispatch->ReturnMat();
  282. }     /* ch */
  283.  
  284. VMatrix &Cv(VMatrix &a, VMatrix &b)     /* vertical concatination */
  285. {
  286.     a.Garbage("Cv");
  287.     b.Garbage("Cv");
  288.     int adim = a.r;
  289.     int bdim = b.r;
  290.     int k = adim + bdim;
  291.  
  292.     if (a.c != b.c)
  293.         a.Nrerror("CV: matrices have different number of columns");
  294.  
  295.     VMatrix *c = new VMatrix("(", k, a.c);
  296.  
  297.     strtype mname = c->Getname(*c) + a.Getname(a) + "//" + b.Getname(b) +")";
  298.     c->Nameit(mname);
  299.  
  300.     for (int i = 1; i <= a.r; i++) {
  301.         for (int j = 1; j <= a.c; j++)
  302.             c->M(i, j) = a.m(i, j);
  303.     }
  304.     int ii;
  305.     for (i = 1, ii = 1; i <= b.r; i++, ii++) {
  306.         for (int j = 1, jj = 1; j <= b.c; j++, jj++)
  307.             c->M(ii + a.r, jj) = b.m(i,j);
  308.     }
  309.     Dispatch->Push(*c);
  310.     delete c;
  311.     return Dispatch->ReturnMat();
  312. }     /* cv */
  313.  
  314.  
  315. //////////////////  math funtions
  316.  
  317. VMatrix& Tran(VMatrix &ROp)
  318. {
  319.     ROp.Garbage("Tran");
  320.     VMatrix *temp = new VMatrix("t", ROp.c, ROp.r);
  321.     strtype newname = temp->Getname(ROp);
  322.     for (int i = 1; i <= ROp.r; i++) {
  323.         for (int j = 1; j <= ROp.c; j++) temp->M(j, i) = ROp.m(i, j);
  324.     }
  325.     newname = newname + "'";
  326.     temp->Nameit(newname);
  327.     Dispatch->Push(*temp);
  328.     delete temp;
  329.     return Dispatch->ReturnMat();
  330. }
  331. VMatrix& Mexp(VMatrix &ROp)
  332. {
  333.     // take exponent of all elements
  334.     ROp.Garbage("Mexp");
  335.     VMatrix *temp = new VMatrix("exp(", ROp.r, ROp.c);
  336.     for (int i = 1; i <= ROp.r; i++) {
  337.         for (int j = 1; j <= ROp.c; j++) temp->M(i, j) = exp(ROp.m(i, j));
  338.     }
  339.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  340.     temp->Nameit(newname);
  341.     Dispatch->Push(*temp);
  342.     delete temp;
  343.     return Dispatch->ReturnMat();
  344. }
  345. VMatrix& Mabs(VMatrix &ROp)
  346. {
  347.     // take absolute values of all elements
  348.     ROp.Garbage("Mabs");
  349.     VMatrix *temp = new VMatrix("abs(", ROp.r, ROp.c);
  350.     for (int i = 1; i <= ROp.r; i++) {
  351.         for (int j = 1; j <= ROp.c; j++) temp->M(i, j) = fabs(ROp.m(i, j));
  352.     }
  353.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  354.     temp->Nameit(newname);
  355.     Dispatch->Push(*temp);
  356.     delete temp;
  357.     return Dispatch->ReturnMat();
  358. }
  359. VMatrix& Mlog(VMatrix &ROp)
  360. {
  361.     // take logarithm of all elements
  362.     ROp.Garbage("Mlog");
  363.     VMatrix *temp = new VMatrix("log(", ROp.r, ROp.c);
  364.     for (int i = 1; i <= ROp.r; i++) {
  365.         for (int j = 1; j <= ROp.c; j++) {
  366.             double R = ROp.m(i, j);
  367.             temp->M(i, j) = (R > 0.0) ? log(R) :
  368.             Dispatch->Nrerror("Mlog: zero or negative argument to log");
  369.         }
  370.     }
  371.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  372.     temp->Nameit(newname);
  373.     Dispatch->Push(*temp);
  374.     delete temp;
  375.     return Dispatch->ReturnMat();
  376. }
  377. VMatrix& Msqrt(VMatrix &ROp)
  378. {
  379.     // take sqrt of all elements
  380.     ROp.Garbage("Msqrt");
  381.     VMatrix *temp = new VMatrix("sqrt(", ROp.r, ROp.c);
  382.     for (int i = 1; i <= ROp.r; i++) {
  383.         for (int j = 1; j <= ROp.c; j++) {
  384.             double R = ROp.m(i, j);
  385.             temp->M(i, j) = (R >= 0.0) ? sqrt(R) :
  386.             Dispatch->Nrerror("Msqrt: zero or negative argument to sqrt");
  387.         }
  388.     }
  389.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  390.     temp->Nameit(newname);
  391.     Dispatch->Push(*temp);
  392.     delete temp;
  393.     return Dispatch->ReturnMat();
  394. }
  395. double Trace(VMatrix &ROp)
  396. {
  397.     ROp.Garbage("Trace");
  398.     double trace = 0;
  399.     if (ROp.r != ROp.c)
  400.         ROp.Nrerror("Trace: matrix not square in trace");
  401.     for (int i = 1; i <= ROp.r; i++) trace += ROp.m(i, i);
  402.     return trace;
  403. }
  404. VMatrix& Ident(int n)
  405. {
  406.     VMatrix *temp = new VMatrix("I(", n, n);
  407.     for (int i = 1; i <= n; i++) {
  408.         for (int j = 1; j <= n; j++) temp-> M(i, j) = (i == j) ? 1.0 : 0.0;
  409.     }
  410.  
  411.     char buffer[MAXCHARS] = { '\0' };
  412.     gcvt(((double) n), NDECS, buffer);
  413.     strtype string = buffer;
  414.     string = temp->Getname(*temp) + string + ")";
  415.  
  416.     temp->Nameit(string);
  417.     Dispatch->Push(*temp);
  418.     delete temp;
  419.     return Dispatch->ReturnMat();
  420. }
  421. VMatrix& Helm(int n)
  422. {
  423.     VMatrix *temp = new VMatrix("Helm(", n, n);
  424.     for (int i = 1; i <= n; i++) {
  425.         for (int j = 1; j <= n; j++) {
  426.             double d = 1.0 / sqrt((double) n);
  427.             double den = (j > 1) ? 1.0 / sqrt((double) j*(j - 1)) : d;
  428.             if (j > 1 && j < i) d = 0;
  429.             if (j > 1 && j == i) d = -den * ((double) (j - 1));
  430.             if (j > 1 && j > i) d = den;
  431.             temp->M(i, j) = d;
  432.         }
  433.     }
  434.     char buffer[MAXCHARS] = { '\0' };
  435.     gcvt(((double) n), NDECS, buffer);
  436.     strtype string = buffer;
  437.     string = temp->Getname(*temp) + string + ")";
  438.     temp->Nameit(string);
  439.     Dispatch->Push(*temp);
  440.     delete temp;
  441.     return Dispatch->ReturnMat();
  442. }
  443. VMatrix& Fill(int r, int c, double d)
  444. {
  445.     VMatrix *temp = new VMatrix("j", r, c);
  446.     for (int i = 1; i <= r; i++) {
  447.         for (int j = 1; j <= c; j++) temp-> M(i, j) = d;
  448.     }
  449.     Dispatch->Push(*temp);
  450.     delete temp;
  451.     return Dispatch->ReturnMat();
  452. }
  453. VMatrix& Index(int lower, int upper, int step)
  454. {
  455.     if (step == 0) Dispatch-> Nrerror("Index: step is zero");
  456.     int cnter = 0;
  457.     if (lower < upper) {
  458.         for (int i = lower; i <= upper; cnter++, i += step);
  459.     }
  460.     else {
  461.         for (int i = lower; i >= upper; cnter++, i += step);
  462.     }
  463.     if (cnter == 0)
  464.         Dispatch->Nrerror(
  465.                 "Index: trying to allocate a matrix with zero rows");
  466.     VMatrix *temp = new VMatrix("Index", cnter, 1);
  467.  
  468.     cnter = 1;
  469.     if (lower < upper) {
  470.         for (int i = lower; i <= upper; i += step)
  471.             temp->M(cnter++, 1) = (double) i;
  472.     }
  473.     else {
  474.         for (int i = lower; i >= upper; i += step)
  475.             temp->M(cnter++, 1) = (double) i;
  476.     }
  477.  
  478.     Dispatch->Push(*temp);
  479.     delete temp;
  480.     return Dispatch->ReturnMat();
  481. }
  482.  
  483. VMatrix& Sweep(int k1, int k2, VMatrix& ROp)     /* sweep out a row */
  484. {
  485.     long double eps = ZERO, d;
  486.     int i, j, k, n, it;
  487.  
  488.     ROp.Garbage("Sweep");
  489.     if (ROp.r != ROp.c) ROp.Nrerror("SWEEP: matrix a is not square");
  490.  
  491.     n = ROp.r;
  492.     if (k2 < k1) { k = k1; k1 = k2; k2 = k; }
  493.  
  494.     for (k = k1; k <= k2; k++) {
  495.         if (fabs(ROp.m(k, k)) < eps)
  496.             for (it = 1; it <= n; it++) {
  497.                 ROp.M(it, k) = 0;
  498.                 ROp.M(k, it) = 0;
  499.         }
  500.         else {
  501.             d = 1.0 / ROp.m(k, k);
  502.             ROp.M(k, k) = d;
  503.             for (i = 1; i <= n; i++) { if (i != k)
  504.                 ROp.M(i, k) = ROp.m(i, k) *(double) - d; }
  505.             for (j = 1; j <= n; j++) { if (j != k)
  506.                 ROp.M(k, j) = ROp.m(k, j) *(double) d; }
  507.                 for (i = 1; i <= n; i++) {
  508.                     if (i != k) {
  509.                         for (j = 1; j <= n; j++) {
  510.                             if (j != k)
  511.                                 ROp.M(i, j) = ROp.m(i, j) + (ROp.m(i, k)) *(ROp.m(k,j))/d;
  512.  
  513.                         }
  514.                     }
  515.                 }
  516.         }
  517.     }
  518.     strtype newname = "sweep(";
  519.     newname = newname + ROp.Getname(ROp) + ")";
  520.     ROp.Nameit(newname);
  521.     Dispatch->Push(ROp);
  522.     return Dispatch->ReturnMat();
  523. }     /*sweep*/
  524.  
  525. VMatrix& Inv(VMatrix &ROp)     /*matrix inversion using sweep */
  526. {
  527.     struct mat *b;
  528.     ROp.Garbage("Inv");
  529.     Dispatch->Inclevel();
  530.  
  531.     if (ROp.c != ROp.r) ROp.Nrerror("INVERSE: matrix a not square");
  532.     strtype newname = "invs(";
  533.     newname = newname + ROp.Getname(ROp) + ")";
  534.     VMatrix *temp = new VMatrix("t", ROp.c, ROp.r);
  535.  
  536.     *temp = ROp;
  537.     *temp = Sweep(1, temp->r, *temp);
  538.  
  539.     temp->Nameit(newname);
  540.     Dispatch->Push(*temp);
  541.     delete temp;
  542.     return Dispatch->DecReturn();
  543. }     /* inverse */
  544.  
  545. VMatrix &Kron(VMatrix &a, VMatrix &b)
  546. {
  547.     a.Garbage("Kron");
  548.     b.Garbage("Kron");
  549.  
  550.     int cr = a.r*b.r;
  551.     int cc = a.c*b.c;
  552.     VMatrix *c = new VMatrix("(", cr, cc);
  553.     strtype mname = c->Getname(*c) + a.Getname(a) + "@" + b.Getname(b) + ")";
  554.     c->Nameit(mname);
  555.  
  556.     for (int i = 1; i <= a.r; i++) {
  557.         for (int j = 1; j <= a.c; j++) {
  558.             for (int k = 1; k <= b.r; k++) {
  559.                 for (int l = 1; l <= b.c; l++) {
  560.                     c->M((i - 1) *b.r + k, (j - 1) *b.c + l) = a.m(i, j) *b.m(k,l);
  561.                 }
  562.             }
  563.         }
  564.     }
  565.     Dispatch->Push(*c);
  566.     delete c;
  567.     return Dispatch->ReturnMat();
  568. }     /* kron */
  569.  
  570.  
  571. void Pivot(VMatrix &Data, int RefRow,
  572.     double &Determ, enum boolean &DetEqualsZero)
  573. {
  574.     DetEqualsZero = TTRUE;
  575.     int NewRow = RefRow;
  576.     while (DetEqualsZero && (NewRow < Data.r)) {
  577.         NewRow++;
  578.         if (fabs(Data.m(NewRow, RefRow)) > ZERO) {
  579.             for (int i = 1; i <= Data.r; i++) {
  580.                 double temp = Data.m(NewRow, i);
  581.                 Data.M(NewRow, i) = Data.m(RefRow, i);
  582.                 Data.M(RefRow, i) = temp;
  583.             }
  584.             DetEqualsZero = FFALSE;
  585.             Determ = -Determ;     /* row swap adjustment to det */
  586.         }
  587.     }
  588. }
  589. double Det(VMatrix &Datain)
  590. {
  591.     Datain.Garbage("Det");
  592.  
  593.     double Determ = 1;
  594.     if (Datain.r == Datain.c) {
  595.         Dispatch->Inclevel();
  596.         int Dimen = Datain.r;
  597.         VMatrix *Data = new VMatrix("data", Dimen, Dimen);
  598.  
  599.         *Data = Datain;
  600.         long double Coef;
  601.         int Row, RefRow = 0;
  602.         enum boolean DetEqualsZero = FFALSE;
  603.  
  604.         while (!(DetEqualsZero) && (RefRow < Dimen - 1)) {
  605.             RefRow++;
  606.             if (fabs(Data->m(RefRow, RefRow)) < ZERO)
  607.                 Pivot(*Data, RefRow, Determ, DetEqualsZero);
  608.             if (!(DetEqualsZero))
  609.                 for (Row = RefRow + 1; Row <= Dimen; Row++)
  610.                     if (fabs(Data->m(Row, RefRow)) > ZERO) {
  611.                         Coef = -((long double) Data->m(Row, RefRow)) /
  612.                         ((long double) Data->m(RefRow, RefRow));
  613.                         for (int i = 1; i <= Dimen; i++)
  614.                             Data->M(Row, i) = Data->m(Row, i) +
  615.                                 (double) (Coef*((long double) Data->m(RefRow,i)));
  616.                     }
  617.                     Determ *= Data->m(RefRow, RefRow);
  618.         }
  619.         Determ = (DetEqualsZero) ? 0 : Determ * Data-> m(Dimen, Dimen);
  620.         delete Data;
  621.         Dispatch->Declevel();
  622.     }
  623.     else Datain.Nrerror("Det: Matrix is not square\n");
  624.     return Determ;
  625. }
  626.  
  627. //--------------- cholesky decomposition ----------------------
  628. // ROp is symmetric, returns upper triangular S where ROp = S'S
  629. //-------------------------------------------------------------
  630. #define EPS 1.0e-7
  631. VMatrix& Cholesky(VMatrix& ROp)
  632. {
  633.     int nr = ROp.r;
  634.  
  635.     ROp.Garbage("Cholesky");
  636.     for (int i = 1; i <= nr; i++) {
  637.         for (int j = 1; j < i; j++)
  638.             if (fabs(ROp.m(i, j) - ROp.m(j, i)) > EPS)
  639.                 Dispatch->Nrerror("Cholesky: Input matrix is not symmetric");
  640.     }
  641.     VMatrix *temp = new VMatrix();
  642.     *temp = Fill(nr, nr, 0.0);
  643.     for (i = 1; i <= nr; i++) {
  644.         for (int j = i; j <= nr; j++) {
  645.             double sum = 0.0;
  646.             for (int k = 1; k < i; k++) {
  647.                 sum += temp->m(k, i) * temp->m(k, j);
  648.             }
  649.             if (i == j) temp-> M(i, j) = (ROp.m(i, j) < sum) ?
  650.                 Dispatch->Nrerror("Cholesky: negative pivot") :
  651.                 sqrt(ROp.m(i, j) - sum);
  652.             else temp->M(i, j) = (fabs(temp->m(i, i)) < EPS) ?
  653.                 Dispatch->Nrerror("Cholesky: zero or negative pivot") :
  654.                 (ROp.m(i, j) - sum) / temp->m(i, i);
  655.         }
  656.     }
  657.     strtype newname = "half(";
  658.     newname = newname + ROp.Getname(ROp) + ")";
  659.     temp->Nameit(newname);
  660.     Dispatch->Push(*temp);
  661.     delete temp;
  662.     return Dispatch->ReturnMat();
  663. }
  664.  
  665.  
  666. void QR(VMatrix& ROp, VMatrix& Q, VMatrix& R, boolean makeq)
  667. {
  668.     ROp.Garbage("QR");
  669.     Q.Garbage("QR");
  670.     R.Garbage("QR");
  671.  
  672.     Dispatch->Inclevel();
  673.     int r = ROp.r;
  674.     int c = ROp.c;
  675.  
  676.     Q = ROp;
  677.     R = Fill(c, c, 0.0);
  678.     double s;
  679.  
  680.     for (int j = 1; j <= c; j++) {
  681.         double sigma = 0.0;
  682.         for (int i = j; i <= r; i++) sigma += Q.m(i, j) *Q.m(i, j);
  683.         if (fabs(sigma) <= 1.0e-10) Dispatch->Nrerror("QR: singular X");
  684.         R.M(j, j) = s = (Q.m(j, j) < 0) ? sqrt(sigma) : - sqrt(sigma);
  685.         double beta = 1.0 /(s*Q.m(j, j) - sigma);
  686.         Q.M(j, j) = Q.m(j, j) - s;
  687.         for (int k = j + 1; k <= c; k++) {
  688.             double sum = 0.0;
  689.             for (int i = j; i <= r; i++) sum += Q.m(i, j) *Q.m(i, k);
  690.             sum *= beta;
  691.             for (i = j; i <= r; i++) Q.M(i, k) = Q.m(i, k) + Q.m(i, j) *sum;
  692.             R.M(j, k) = Q.m(j, k);
  693.         }
  694.     }
  695.  
  696.     if (makeq) Q = ROp*Ginv(R);
  697.  
  698.     Q.Nameit("Q");
  699.     R.Nameit("R");
  700.     Dispatch->Declevel();
  701. }
  702. #undef EPS
  703.  
  704.  
  705. //------------------- Singular Valued Decomposition ----------
  706. // from EISPACK SVD
  707. //------------------------------------------------------------
  708. long double safehypot(long double a1, long double b1)
  709.     // function to get hypotenuse numbers a1 and b1
  710.     // where a=log(a1) and b=log(b1)
  711. {
  712.             long double del, sum;
  713.             long double a = 2.0*log(fabs(a1));
  714.             long double b = 2.0*log(fabs(b1));
  715.             del = (a > b) ? a : b;
  716.             // add a1^2 + b1^2, but in logarithms
  717.             sum = log(exp(a - del) + exp(b - del)) + del;
  718.             return (exp(0.5*sum));
  719. }
  720.  
  721. void svdtemp(VMatrix &a, VMatrix &w, boolean matu, VMatrix &u,
  722.     boolean matv, VMatrix &v, int &ierr, VMatrix &rv1)
  723. {
  724.     Dispatch->Inclevel();
  725.     a.Garbage("SVD");
  726.     w.Garbage("SVD");
  727.     u.Garbage("SVD");
  728.     v.Garbage("SVD");
  729.     rv1.Garbage("SVD");
  730.  
  731.     int m = a.r;
  732.     int n = a.c;
  733.  
  734.     int i, j, k, l, ii, i1, kk, k1, ll, l1, its, mn;
  735.     //      VMatrix a("a",m,n),w("w",n,1),
  736.     //              u("u",m,m),v("v",n,n),rv1("rv1",n,1);
  737.     long double c, f, g, h, s, x, y, z, eps, scale, machep, fss;
  738.     //      boolean matu,matv;
  739.  
  740.     //       ********** Machep is a machine dependent parameter specifying
  741.     //              the relative precision of floating point aritmetic.
  742.     //
  743.     //       for pc doubles, range is 1.6^-308 to 1.6^308
  744.     //       **************
  745.  
  746.     machep = pow(1.6, -308.0);
  747.  
  748.     ierr = 0;
  749.  
  750.     u = a;
  751.     //     ********householder reduction to bi diagonal form ********
  752.     g = 0.0;
  753.     scale = 0.0;
  754.     x = 0.0;
  755.  
  756.     for (i = 1; i <= n; i++) {
  757.         l = i + 1;
  758.         rv1.M(i, 1) = scale*g;
  759.         g = 0.0;
  760.         s = 0.0;
  761.         scale = 0.0;
  762.         if (i > m) goto l210;
  763.  
  764.         for (k = i; k <= m; k++) scale += fabs(u.m(k, i));
  765.  
  766.         if (scale == 0.0) goto l210;
  767.  
  768.         for (k = i; k <= m; k++) {
  769.             u.M(k, i) = u.m(k, i) / scale;
  770.             s += u.m(k, i) *u.m(k, i);
  771.         }
  772.  
  773.         f = u.m(i, i);
  774.         g = -((f >= 0.0) ? fabs(sqrt(s)) : - fabs(sqrt(s)));
  775.         h = f*g - s;
  776.         u.M(i, i) = f - g;
  777.         if (i == n) goto l190;
  778.  
  779.         for (j = l; j <= n; j++) {
  780.             s = 0.0;
  781.             for (k = i; k <= m; k++) s += u.m(k, i) *u.m(k, j);
  782.             f = s / h;
  783.             for (k = i; k <= m; k++) u.M(k, j) = u.m(k, j) + f*u.m(k, i);
  784.         }
  785.         l190 : for (k = i; k <= m; k++) u.M(k, i) = scale*u.m(k, i);
  786.  
  787.         l210 : w.M(i, 1) = scale*g;
  788.         g = 0.0;
  789.         s = 0.0;
  790.         scale = 0.0;
  791.         if (i > m || i == n) goto l290;
  792.  
  793.         for (k = l; k <= n; k++) scale += fabs(u.m(i, k));
  794.  
  795.         if (scale == 0.0) goto l290;
  796.  
  797.         for (k = l; k <= n; k++) {
  798.             u.M(i, k) = u.m(i, k) / scale;
  799.             s += u.m(i, k) *u.m(i, k);
  800.         }
  801.         f = u.m(i, l);
  802.         g = -((f >= 0.0) ? fabs(sqrt(s)) : - fabs(sqrt(s)));
  803.         h = f*g - s;
  804.         u.M(i, l) = f - g;
  805.  
  806.         for (k = l; k <= n; k++) rv1.M(k, 1) = u.m(i, k) / h;
  807.  
  808.         if (i == m) goto l270;
  809.  
  810.         for (j = l; j <= m; j++) {
  811.             s = 0.0;
  812.             for (k = l; k <= n; k++) s += u.m(j, k) *u.m(i, k);
  813.             for (k = l; k <= n; k++) u.M(j, k) = u.m(j, k) + s*rv1.m(k, 1);
  814.         }
  815.  
  816.         l270 : for (k = l; k <= n; k++) u.M(i, k) = scale*u.m(i, k);
  817.  
  818.         l290 : x = (x > fabs(w.m(i, 1)) + fabs(rv1.m(i, 1))) ? x :
  819.                         fabs(w.m(i, 1)) + fabs(rv1.m(i, 1));
  820.     }
  821.  
  822.     //       ********** accumulation of right-hand transformations ********
  823.  
  824.     if (!matv) goto l410;
  825.  
  826.     for (ii = 1; ii <= n; ii++) {
  827.         i = n + 1- ii;
  828.         if (i == n) goto l390;
  829.         if (g == 0.0) goto l360;
  830.  
  831.         //               double division avoids possible underflow
  832.         for (j = l; j <= n; j++) v.M(j, i) = (u.m(i, j) / u.m(i, l)) / g;
  833.  
  834.         for (j = l; j <= n; j++) {
  835.             s = 0.0;
  836.             for (k = l; k <= n; k++) s += u.m(i, k) *v.m(k, j);
  837.             for (k = l; k <= n; k++) v.M(k, j) = v.m(k, j) + s*v.m(k, i);
  838.         }
  839.  
  840.         l360 : for (j = l; j <= n; j++) {
  841.             v.M(i, j) = 0.0;
  842.             v.M(j, i) = 0.0;
  843.         }
  844.  
  845.         l390 : v.M(i, i) = 1.0;
  846.         g = rv1.m(i, 1);
  847.         l = i;
  848.     }
  849.     //       *************accumulation of left-hand transformations
  850.     l410 : if (!matu) goto l510;
  851.     //       ************* for i= min(m,n) step -1 until 1 do ******
  852.     mn = n;
  853.     if (m < n) mn = m;
  854.  
  855.     for (ii = 1; ii <= mn; ii++) {
  856.         i = mn + 1- ii;
  857.         l = i + 1;
  858.         g = w.m(i, 1);
  859.         if (i == n) goto l430;
  860.  
  861.         for (j = l; j <= n; j++) u.M(i, j) = 0.0;
  862.  
  863.         l430 : if (g == 0.0) goto l475;
  864.         if (i == mn) goto l460;
  865.  
  866.         for (j = l; j <= n; j++) {
  867.             s = 0.0;
  868.  
  869.             for (k = l; k <= m; k++) s += u.m(k, i) *u.m(k, j);
  870.             //double division avoids possible underflow
  871.             f =(s / u.m(i, i)) / g;
  872.             for (k = i; k <= m; k++) u.M(k, j) = u.m(k, j) + f*u.m(k, i);
  873.         }
  874.  
  875.         l460 : for (j = i; j <= m; j++) u.M(j, i) = u.m(j, i) / g;
  876.  
  877.         goto l490;
  878.  
  879.         l475 : for (j = i; j <= m; j++) u.M(j, i) = 0.0;
  880.  
  881.         l490 : u.M(i, i) = u.m(i, i) + 1.0;
  882.     }
  883.  
  884.     //       ********** diagonalization of the bidiagonal form ***********
  885.  
  886.     l510 : eps = machep*x;
  887.     //       ********** for k=n step -1 until 1 do ***********************
  888.     for (kk = 1; kk <= n; kk++) {
  889.         k1 = n - kk;
  890.         k = k1 + 1;
  891.         its = 0;
  892.         //       ********** test for splitting .**********
  893.         //                 for l=k step -1 until 1 do
  894.         l520 : for (ll = 1; ll <= k; ll++) {
  895.             l1 = k - ll;
  896.             l = l1 + 1;
  897.             if (fabs(rv1.m(l, 1)) <= eps) goto l565;
  898.             //  rv1(1) is always zero, so there is no exit
  899.             //  through the bottom of the loop *********
  900.             if (fabs(w.m(l1, 1)) <= eps) goto l540;
  901.         }
  902.         //  ************* cancellation of rv1(l) if l greater than 1******
  903.         l540 : c = 0.0;
  904.         s = 1.0;
  905.  
  906.         for (i = l; i <= k; i++) {
  907.             f = s*rv1.m(i, 1);
  908.             rv1.M(i, 1) = c*rv1.m(i, 1);
  909.             if (fabs(f) <= eps) goto l565;
  910.             g = w.m(i, 1);
  911.             h = safehypot(f,g);
  912.             w.M(i, 1) = h;
  913.             c = g / h;
  914.             s = -f / h;
  915.             if (matu) {     //go to 560
  916.                 for (j = 1; j <= m; j++) {
  917.                     y = u.m(j, l1);
  918.                     z = u.m(j, i);
  919.                     u.M(j, l1) = y*c + z*s;
  920.                     u.M(j, i) = -y*s + z*c;
  921.                 }
  922.             }
  923.         }
  924.         //       ************** test for convergence ********************
  925.         l565 : z = w.m(k, 1);
  926.         if (l == k) goto l650;
  927.         //       *************** if shift from bottom 2 by 2 minor *******
  928.         if (its == 50) goto l1000;
  929.         its = its + 1;
  930.         x = w.m(l, 1);
  931.         y = w.m(k1, 1);
  932.         g = rv1.m(k1, 1);
  933.         h = rv1.m(k, 1);
  934.         f =((y - z) *(y + z) + (g - h) *(g + h)) / (2.0*h*y);
  935.         g = safehypot(f,1.0);
  936.         fss =(f >= 0.0) ? fabs(g) : - fabs(g);
  937.         f =((x - z) *(x + z) + h*(y / (f + fss) - h)) / x;
  938.         //       **************** next qr transformation ***************
  939.         c = 1.0;
  940.         s = 1.0;
  941.         //
  942.         for (i1 = l; i1 <= k1; i1++) {
  943.             i = i1 + 1;
  944.             g = rv1.m(i, 1);
  945.             y = w.m(i, 1);
  946.             h = s*g;
  947.             g = c*g;
  948.             z = safehypot(f,h);
  949.             rv1.M(i1, 1) = z;
  950.             c = f / z;
  951.             s = h / z;
  952.             f = x*c + g*s;
  953.             g = -x*s + g*c;
  954.             h = y*s;
  955.             y = y*c;
  956.             if (matv) {     // goto l575;
  957.                 for (j = 1; j <= n; j++) {
  958.                     x = v.m(j, i1);
  959.                     z = v.m(j, i);
  960.                     v.M(j, i1) = x*c + z*s;
  961.                     v.M(j, i) = -x*s + z*c;
  962.                 }
  963.             }
  964.             z = safehypot(f,h);
  965.             w.M(i1, 1) = z;
  966.             //  ************ rotation can be arbitrary if z is zero *******
  967.             if (z != 0.0) {     // goto l580;
  968.                 c = f / z;
  969.                 s = h / z;
  970.             }
  971.             // l580:
  972.             f = c*g + s*y;
  973.             x = -s*g + c*y;
  974.             if (matu) {     //goto l600;
  975.                 for (j = 1; j <= m; j++) {
  976.                     y = u.m(j, i1);
  977.                     z = u.m(j, i);
  978.                     u.M(j, i1) = y*c + z*s;
  979.                     u.M(j, i) = -y*s + z*c;
  980.                 }
  981.             }
  982.  
  983.             l600 : x = x; }     //continue
  984.  
  985.         rv1.M(l, 1) = 0.0;
  986.         rv1.M(k, 1) = f;
  987.         w.M(k, 1) = x;
  988.         goto l520;
  989.         //       ************** convergence **************
  990.         l650 : if (z >= 0.0) goto l700;
  991.         //       *************  w(k) is made non-negative *********
  992.         w.M(k, 1) = -z;
  993.         if (!matv) goto l700;
  994.  
  995.         for (j = 1; j <= n; j++) v.M(j, k) = -v.m(j, k);
  996.  
  997.         l700 : x = x; }
  998.  
  999.  
  1000.     ierr = 0;
  1001.     Dispatch->Declevel();
  1002.     return;
  1003.     //        ***************** set error -- no convergence to a
  1004.     //              singular value after 30 iterations
  1005.     l1000 : ierr = k;
  1006.     Dispatch->Declevel();
  1007.     return;
  1008. }
  1009.  
  1010. int Svd(VMatrix &A, VMatrix &U, VMatrix &S, VMatrix &V,
  1011.     boolean makeu, boolean makev)
  1012. {
  1013.  
  1014.     Dispatch->Inclevel();
  1015.  
  1016.     A.Garbage("SVD");
  1017.     VMatrix aa, uu, ss, vv, rv1;
  1018.     int ierr = 0;
  1019.  
  1020.     if (A.r < A.c) aa = Tran(A);
  1021.     else aa = A;
  1022.  
  1023.     uu = Fill(aa.r, aa.r, 0.0);
  1024.     vv = Fill(aa.c, aa.c, 0.0);
  1025.     ss = Fill(aa.c, 1, 0.0);
  1026.     rv1 = Fill(aa.c, 1, 0.0);
  1027.  
  1028.  
  1029.     svdtemp(aa, ss, makeu, uu, makev, vv, ierr, rv1);
  1030.  
  1031.     if (A.r < A.c) {
  1032.         if (makeu) U = vv;
  1033.         if (makev) V = uu;
  1034.     } else {
  1035.         if (makeu) U = uu;
  1036.         if (makev) V = vv;
  1037.     }
  1038.     S = ss;
  1039.     Dispatch->Declevel();
  1040.     return ierr;
  1041. }
  1042. //---------------- end svd ------------------------------------
  1043.  
  1044. VMatrix& Ginv(VMatrix& ROp)
  1045. {
  1046.     ROp.Garbage("Ginv");
  1047.     Dispatch->Inclevel();
  1048.  
  1049.     VMatrix *u = new VMatrix;
  1050.     VMatrix *s = new VMatrix;
  1051.     VMatrix *v = new VMatrix;
  1052.     VMatrix *g = new VMatrix;
  1053.  
  1054.     Svd(ROp, *u, *s, *v);
  1055.     for (int i = 1; i <= s-> r; i++) {
  1056.         double t = s->m(i, 1);
  1057.         s->M(i, 1) = (fabs(t) <= 0.0) ? 0.0 : 1.0 / t;
  1058.     }
  1059.     *g = *v *Diag(*s) *Tran(*u);
  1060.     strtype newname = "ginv(";
  1061.     newname = newname + ROp.Getname(ROp) + ")";
  1062.     g->Nameit(newname);
  1063.     Dispatch->Push(*g);
  1064.     delete u;
  1065.     delete s;
  1066.     delete v;
  1067.     delete g;
  1068.     return Dispatch->DecReturn();
  1069. }
  1070.  
  1071. //-------------------------- fft ------------------------------
  1072. // de Boor (1980) SIAM J SCI. STAT. COMPUT. V1 no 1, pp 173-178
  1073. // and NewMat by R. G. Davies a C++ matrix Package
  1074. //-------------------------------------------------------------
  1075.  
  1076. static void cossin(int n, int d, double& c, double& s)
  1077.     // calculate cos(twopi*n/d) and sin(twopi*n/d)
  1078.         // minimise roundoff error
  1079.         {
  1080.             long n4 = n * 4;
  1081.             int sector =(int) floor((double) n4 / (double) d + 0.5);
  1082.             n4 -= sector * d;
  1083.             if (sector < 0) sector = 3 - (3 - sector) % 4; else sector %= 4;
  1084.             double ratio = 1.5707963267948966192 * (double) n4 / (double) d;
  1085.  
  1086.             switch (sector) {
  1087.                 case 0 : c = cos(ratio);
  1088.                 s = sin(ratio);
  1089.                 break;
  1090.                 case 1 : c = -sin(ratio);
  1091.                 s = cos(ratio);
  1092.                 break;
  1093.                 case 2 : c = -cos(ratio);
  1094.                 s = -sin(ratio);
  1095.                 break;
  1096.                 case 3 : c = sin(ratio);
  1097.                 s = -cos(ratio);
  1098.                 break;
  1099.             }
  1100. }
  1101.  
  1102.  
  1103. static void fftstep(VMatrix &AB, VMatrix &XY,
  1104.     int after, int now, int before)
  1105. {
  1106.  
  1107.     int gamma = after * before;
  1108.     int delta = now * after;
  1109.     double r_value, i_value, ridx_value, iidx_value, temp;
  1110.     int j, x1, ia, a, a1, x2, ib, a2, idx, in;
  1111.  
  1112.     double r_arg = 1.0, i_arg = 0.0;
  1113.  
  1114.     int x = 1;
  1115.     int m = AB.r - gamma;
  1116.  
  1117.     for (j = 0; j < now; j++) {
  1118.         a = 1;
  1119.         x1 = x;
  1120.         x += after;
  1121.         for ( ia = 0; ia < after; ia++) {
  1122.             // generate sins & cosines explicitly rather than iteratively
  1123.             // for more accuracy; but slower
  1124.             cossin(- (j*after + ia), delta, r_arg, i_arg);
  1125.             a1 = a++;
  1126.             x2 = x1++;
  1127.             if (now == 2) {
  1128.                 ib = before;
  1129.                 while (ib--) {
  1130.                     a2 = m + a1;
  1131.                     a1 += after;
  1132.                     idx = a2 - gamma;
  1133.                     r_value = AB.m(a2, 1);
  1134.                     i_value = AB.m(a2, 2);
  1135.                     ridx_value = AB.m(idx,1);
  1136.                     iidx_value = AB.m(idx,2);
  1137.                     XY.M(x2, 1) = r_value*r_arg - i_value*i_arg + ridx_value;
  1138.                     XY.M(x2, 2) = r_value*i_arg + i_value*r_arg + iidx_value;
  1139.                     x2 += delta;
  1140.                 }
  1141.             }
  1142.             else {
  1143.                 ib = before;
  1144.                 while (ib--) {
  1145.                     a2 = m + a1;
  1146.                     a1 += after;
  1147.                     r_value = AB.m(a2, 1);
  1148.                     i_value = AB.m(a2, 2);
  1149.                     in = now - 1;
  1150.                     while (in--) {
  1151.                         a2 -= gamma;
  1152.                         ridx_value = AB.m(a2, 1);
  1153.                         iidx_value = AB.m(a2, 2);
  1154.                         temp = r_value;
  1155.                         r_value = r_value*r_arg - i_value*i_arg + ridx_value;
  1156.                         i_value = temp*i_arg + i_value*r_arg + iidx_value;
  1157.                     }
  1158.                     XY.M(x2, 1) = r_value;
  1159.                     XY.M(x2, 2) = i_value;
  1160.                     x2 += delta;
  1161.                 }
  1162.             }
  1163.         }
  1164.     }
  1165. }
  1166.  
  1167.  
  1168. VMatrix& Fft(VMatrix &ROp, boolean INZEE)
  1169. // if INZEE == TTRUE  then calculate fft
  1170. // if INZEE == FFALSE then calculate inverse fft
  1171. {
  1172.     ROp.Garbage("FFT");
  1173.     if (ROp.c > 2) Dispatch->Nrerror("FFT: Input has more than two columns");
  1174.  
  1175.     int n = ROp.r;     // length of arrays
  1176.     const int nextmx = 26;
  1177.     int prime[26] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
  1178.         41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
  1179.         89, 97, 101 };
  1180.     int after = 1;
  1181.     int before = n;
  1182.     int next = 0;
  1183.     boolean inzee = TTRUE;
  1184.  
  1185.     Dispatch->Inclevel();
  1186.  
  1187.     VMatrix *AB = new VMatrix("AB", n, 2);
  1188.     VMatrix *XY = new VMatrix("XY", n, 2);
  1189.  
  1190.     if (ROp.c == 1) {
  1191.         for (int i = 1; i <= n; i++) {
  1192.             AB->M(i, 1) = ROp.m(i, 1);
  1193.             AB->M(i, 2) = 0.0;
  1194.         }
  1195.     } else *AB = ROp;
  1196.     *XY = Fill(n, 2, 0.0);
  1197.  
  1198.     if (INZEE == FFALSE) {
  1199.         // take complex conjugate for ifft
  1200.         for (int i = 1; i <= n; i++) AB-> M(i, 2) = -AB->m(i, 2);
  1201.     }
  1202.  
  1203.     do {
  1204.         int now, b1;
  1205.         for (;;) {
  1206.             if (next < nextmx) now = prime[next];
  1207.             b1 = before / now;
  1208.             if (b1 * now == before) break;
  1209.             next++;
  1210.             now += 2;
  1211.         }
  1212.         before = b1;
  1213.  
  1214.         if (inzee) fftstep(*AB, *XY, after, now, before);
  1215.         else fftstep(*XY, *AB, after, now, before);
  1216.  
  1217.         inzee =(inzee == TTRUE) ? FFALSE : TTRUE;
  1218.         after *= now;
  1219.     }
  1220.     while (before != 1);
  1221.  
  1222.     if (inzee == TTRUE) *XY = *AB;
  1223.     // divide by n for ifft
  1224.     if (INZEE == FFALSE) *XY = (*XY) / ((double) n);
  1225.  
  1226.     strtype newname = (INZEE == TTRUE) ? "fft( " : "ifft(";
  1227.     newname = newname + ROp.Getname(ROp) + ")";
  1228.     XY->Nameit(newname);
  1229.  
  1230.     Dispatch->Push(*XY);
  1231.     delete XY;
  1232.     delete AB;
  1233.     return Dispatch->DecReturn();
  1234. }
  1235.  
  1236.  
  1237. ////////////////////  reshaping functions
  1238.  
  1239. VMatrix& Vec(VMatrix& ROp)
  1240. {     // send columns of ROp to a vector
  1241.     ROp.Garbage("Vec");
  1242.  
  1243.     long lrc = ((long) ROp.r) *((long) ROp.c);
  1244.     if (lrc >= 32768 || lrc < 1)
  1245.         Dispatch->Nrerror("Vec: vec produces invalid indices");
  1246.     int r = ROp.r*ROp.c;
  1247.     int c = 1;
  1248.     VMatrix *temp = new VMatrix("t", r, c);
  1249.     int cnter = 1;
  1250.     for (int j = 1; j <= ROp.c; j++)
  1251.         for (int i = 1; i <= ROp.r; i++) temp->M(cnter++, 1) = ROp.m(i, j);
  1252.  
  1253.     strtype newname = "vec(";
  1254.     newname = newname + ROp.Getname(ROp) + ")";
  1255.     temp->Nameit(newname);
  1256.     Dispatch->Push(*temp);
  1257.     delete temp;
  1258.     return Dispatch->ReturnMat();
  1259. }
  1260.  
  1261. VMatrix& Vecdiag(VMatrix& ROp)
  1262. {     // Extract the diagonal from a square matrix and put in a vector
  1263.     ROp.Garbage("Vecdiag");
  1264.     if (ROp.r != ROp.c)
  1265.         Dispatch->Nrerror("Vecdiag: ROp is not square");
  1266.  
  1267.     VMatrix *temp = new VMatrix("vdiag(", ROp.r, 1);
  1268.     for (int i = 1; i <= ROp.r; i++) temp->M(i, 1) = ROp.m(i, i);
  1269.  
  1270.     strtype newname = "vdiag(";
  1271.     newname = newname + ROp.Getname(ROp) + ")";
  1272.     temp->Nameit(newname);
  1273.     Dispatch->Push(*temp);
  1274.     delete temp;
  1275.     return Dispatch->ReturnMat();
  1276. }
  1277.  
  1278. VMatrix& Diag(VMatrix& ROp)
  1279. {     // make a diagonal matrix from a vector or the principal diagonal of
  1280.  
  1281.     // a square matrix, off diagonal elements are zero
  1282.     ROp.Garbage("Diag");
  1283.  
  1284.     if (ROp.r != ROp.c && ROp.c != 1)
  1285.         Dispatch->Nrerror(
  1286.                 "Diag: ROp is not square or is not a column vector"    );
  1287.     Dispatch->Inclevel();
  1288.  
  1289.     VMatrix *temp = new VMatrix;
  1290.     *temp = Fill(ROp.r, ROp.r, 0.0);
  1291.  
  1292.     int rr = ROp.r;
  1293.     int cc = ROp.c;
  1294.     for (int i = 1; i <= rr; i++)
  1295.         temp->M(i, i) = (rr == cc) ? ROp.m(i, i) : ROp.m(i, 1);
  1296.         strtype newname = "diag(";
  1297.         newname = newname + ROp.Getname(ROp) + ")";
  1298.         temp->Nameit(newname);
  1299.         Dispatch->Push(*temp);
  1300.         delete temp;
  1301.         return Dispatch->DecReturn();
  1302. }
  1303.  
  1304. VMatrix& Shape(VMatrix& ROp, int nrows)
  1305. {     // reshape a matrix into n rows, nrows must divide r*c without a
  1306.     // remainder
  1307.     ROp.Garbage("Shape");
  1308.  
  1309.     long nr = (long) nrows;
  1310.     long lr = (long) ROp.r;
  1311.     long lc = (long) ROp.c;
  1312.  
  1313.     if (lr*lc % nr)
  1314.         Dispatch->Nrerror("Shape: nrows divides r*c with a remainder");
  1315.  
  1316.     Dispatch->Inclevel();
  1317.     VMatrix *temp = new VMatrix;
  1318.     *temp = ROp;
  1319.  
  1320.     long lrc = (lr*lc) / nr;
  1321.     if (nr >= 32768 || nr < 1 || lrc >= 32768 || lrc < 1)
  1322.         Dispatch->Nrerror("Shape: reshape produces invalid indices");
  1323.     temp->r = nrows;
  1324.     temp->c = (int) lrc;
  1325.  
  1326.     strtype newname = "shape(";
  1327.     newname = newname + ROp.Getname(ROp) + ")";
  1328.     temp->Nameit(newname);
  1329.     Dispatch->Push(*temp);
  1330.     delete temp;
  1331.     return Dispatch->DecReturn();
  1332. }
  1333.  
  1334.  
  1335. ///////////////////////////// summation functions // maxs and mins
  1336.  
  1337. //// typedef enum margins { ALL, ROWS, COLUMNS } margins;
  1338.  
  1339. VMatrix& Sum(VMatrix& ROp, margins marg)
  1340. {     // sum(ROp, ALL) returns 1x1
  1341.     // sum(ROp, ROWS) returns 1xc
  1342.     // sum(ROp, COLUMNS) returns cx1
  1343.     ROp.Garbage("Sum");
  1344.  
  1345.     int r, c;
  1346.     double sum;
  1347.     switch (marg) {
  1348.         case ALL : r = 1; c = 1; break;
  1349.         case ROWS : r = 1; c = ROp.c; break;
  1350.         case COLUMNS : r = ROp.r; c = 1; break;
  1351.         default : Dispatch->Nrerror("Sum: invalid margin specification");
  1352.     }
  1353.     VMatrix *temp = new VMatrix("t", r, c);
  1354.  
  1355.     int i, j;
  1356.     switch (marg) {
  1357.         case ALL : sum = 0.0;
  1358.         for (i = 1; i <= ROp.r; i++)
  1359.             for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j);
  1360.         temp->M(1, 1) = sum;
  1361.         break;
  1362.         case ROWS : for (j = 1; j <= c; j++) {
  1363.             sum = 0.0;
  1364.             for (i = 1; i <= ROp.r; i++) sum += ROp.m(i, j);
  1365.             temp->M(1, j) = sum;
  1366.         }
  1367.         break;
  1368.         case COLUMNS : for (i = 1; i <= r; i++) {
  1369.             sum = 0.0;
  1370.             for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j);
  1371.             temp->M(i, 1) = sum;
  1372.         }
  1373.         break;
  1374.     }
  1375.  
  1376.     strtype newname = "sum(";
  1377.     newname = newname + ROp.Getname(ROp) + ")";
  1378.     temp->Nameit(newname);
  1379.     Dispatch->Push(*temp);
  1380.     delete temp;
  1381.     return Dispatch->ReturnMat();
  1382. }
  1383.  
  1384. VMatrix& Sumsq(VMatrix& ROp, margins marg)
  1385. {     // sumsq(ROp, ALL) returns 1x1
  1386.     // sumsq(ROp, ROWS) returns 1xc
  1387.     // sumsq(ROp, COLUMNS) returns cx1
  1388.     ROp.Garbage("Sum");
  1389.  
  1390.     int r, c;
  1391.     double sum;
  1392.     switch (marg) {
  1393.         case ALL : r = 1; c = 1; break;
  1394.         case ROWS : r = 1; c = ROp.c; break;
  1395.         case COLUMNS : r = ROp.r; c = 1; break;
  1396.         default : Dispatch->Nrerror("Sumsq: invalid margin specification");
  1397.     }
  1398.     VMatrix *temp = new VMatrix("t", r, c);
  1399.  
  1400.     int i, j;
  1401.     switch (marg) {
  1402.         case ALL : sum = 0.0;
  1403.         for (i = 1; i <= ROp.r; i++)
  1404.             for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j) *ROp.m(i, j);
  1405.         temp->M(1, 1) = sum;
  1406.         break;
  1407.         case ROWS : for (j = 1; j <= c; j++) {
  1408.             sum = 0.0;
  1409.             for (i = 1; i <= ROp.r; i++) sum += ROp.m(i, j) *ROp.m(i, j);
  1410.             temp->M(1, j) = sum;
  1411.         }
  1412.         break;
  1413.         case COLUMNS : for (i = 1; i <= r; i++) {
  1414.             sum = 0.0;
  1415.             for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j) *ROp.m(i, j);
  1416.             temp->M(i, 1) = sum;
  1417.         }
  1418.         break;
  1419.     }
  1420.  
  1421.     strtype newname = "sumsq(";
  1422.     newname = newname + ROp.Getname(ROp) + ")";
  1423.     temp->Nameit(newname);
  1424.     Dispatch->Push(*temp);
  1425.     delete temp;
  1426.     return Dispatch->ReturnMat();
  1427. }
  1428.  
  1429. VMatrix& Cusum(VMatrix& ROp)
  1430. {     // cumulative sum along rows
  1431.     ROp.Garbage("Cusum");
  1432.     VMatrix *temp = new VMatrix("t", ROp.r, ROp.c);
  1433.  
  1434.     double sum = 0.0;
  1435.     for (int i = 1; i <= ROp.r; i++) {
  1436.         for (int j = 1; j <= ROp.c; j++) {
  1437.             sum += ROp.m(i, j);
  1438.             temp->M(i, j) = sum;
  1439.         }
  1440.     }
  1441.     strtype newname = "cusum(";
  1442.     newname = newname + ROp.Getname(ROp) + ")";
  1443.     temp->Nameit(newname);
  1444.     Dispatch->Push(*temp);
  1445.     delete temp;
  1446.     return Dispatch->ReturnMat();
  1447. }
  1448.  
  1449. VMatrix& Mmax(VMatrix& ROp, margins marg)
  1450. {     // Mmax(ROp, ALL) returns 3x1
  1451.     // Mmax(ROp, ROWS) returns 3xc
  1452.     // Mmax(ROp, COLUMNS) returns cx3
  1453.     ROp.Garbage("Sum");
  1454.  
  1455.     int r, c;
  1456.     switch (marg) {
  1457.         case ALL : r = 3; c = 1; break;
  1458.         case ROWS : r = 3; c = ROp.c; break;
  1459.         case COLUMNS : r = ROp.r; c = 3; break;
  1460.         default : Dispatch->Nrerror("Mmax: invalid margin specification");
  1461.     }
  1462.     VMatrix *temp = new VMatrix("t", r, c);
  1463.  
  1464.     int maxr, maxc, i, j;
  1465.     double mmax;
  1466.     switch (marg) {
  1467.         case ALL : mmax = ROp.m(1, 1);
  1468.         maxr = 1; maxc = 1;
  1469.         for (i = 1; i <= ROp.r; i++)
  1470.             for (j = 1; j <= ROp.c; j++)
  1471.                 mmax = (mmax < ROp.m(i, j)) ?
  1472.                     maxr = i, maxc = j, ROp.m(i, j) : mmax;
  1473.                     temp->M(1, 1) = (double) maxr;
  1474.                     temp->M(2, 1) = (double) maxc;
  1475.                     temp->M(3, 1) = mmax;
  1476.                     break;
  1477.                     case ROWS : for (j = 1; j <= c; j++) {
  1478.                         mmax = ROp.m(1, j);
  1479.                         maxr = 1; maxc = j;
  1480.                         for (i = 1; i <= ROp.r; i++)
  1481.                             mmax = (mmax < ROp.m(i, j)) ? maxr = i, ROp.m(i,j): mmax;
  1482.                             temp->M(1, j) = (double) maxr;
  1483.                             temp->M(2, j) = (double) maxc;
  1484.                             temp->M(3, j) = mmax;
  1485.                     }
  1486.                     break;
  1487.                     case COLUMNS : for (i = 1; i <= r; i++) {
  1488.                         mmax = ROp.m(i, 1);
  1489.                         maxr = i; maxc = 1;
  1490.                         for (j = 1; j <= ROp.c; j++)
  1491.                             mmax = (mmax < ROp.m(i, j)) ? maxc = j, ROp.m(i,j):    mmax;
  1492.                             temp->M(i, 1) = (double) maxr;
  1493.                             temp->M(i, 2) = (double) maxc;
  1494.                             temp->M(i, 3) = mmax;
  1495.                     }
  1496.                     break;
  1497.     }
  1498.  
  1499.     strtype newname = "max(";
  1500.     newname = newname + ROp.Getname(ROp) + ")";
  1501.     temp->Nameit(newname);
  1502.     Dispatch->Push(*temp);
  1503.     delete temp;
  1504.     return Dispatch->ReturnMat();
  1505. }
  1506.  
  1507.  
  1508.  
  1509. VMatrix& Mmin(VMatrix& ROp, margins marg)
  1510. {     // Mmin(ROp, ALL) returns 3x1
  1511.     // Mmin(ROp, ROWS) returns 3xc
  1512.     // Mmin(ROp, COLUMNS) returns cx3
  1513.     ROp.Garbage("Sum");
  1514.  
  1515.     int r, c;
  1516.     switch (marg) {
  1517.         case ALL : r = 3; c = 1; break;
  1518.         case ROWS : r = 3; c = ROp.c; break;
  1519.         case COLUMNS : r = ROp.r; c = 3; break;
  1520.         default : Dispatch->Nrerror("Mmin: invalid margin specification");
  1521.     }
  1522.     VMatrix *temp = new VMatrix("t", r, c);
  1523.  
  1524.     int maxr, maxc, i, j;
  1525.     double mmin;
  1526.     switch (marg) {
  1527.         case ALL : mmin = ROp.m(1, 1);
  1528.         maxr = 1; maxc = 1;
  1529.         for (i = 1; i <= ROp.r; i++)
  1530.             for (j = 1; j <= ROp.c; j++)
  1531.                 mmin = (mmin > ROp.m(i, j)) ?
  1532.                     maxr = i, maxc = j, ROp.m(i, j) : mmin;
  1533.                     temp->M(1, 1) = (double) maxr;
  1534.                     temp->M(2, 1) = (double) maxc;
  1535.                     temp->M(3, 1) = mmin;
  1536.                     break;
  1537.                     case ROWS : for (j = 1; j <= c; j++) {
  1538.                         mmin = ROp.m(1, j);
  1539.                         maxr = 1; maxc = j;
  1540.                         for (i = 1; i <= ROp.r; i++)
  1541.                             mmin = (mmin > ROp.m(i, j)) ? maxr = i, ROp.m(i,j): mmin;
  1542.                             temp->M(1, j) = (double) maxr;
  1543.                             temp->M(2, j) = (double) maxc;
  1544.                             temp->M(3, j) = mmin;
  1545.                     }
  1546.                     break;
  1547.                     case COLUMNS : for (i = 1; i <= r; i++) {
  1548.                         mmin = ROp.m(i, 1);
  1549.                         maxr = i; maxc = 1;
  1550.                         for (j = 1; j <= ROp.c; j++)
  1551.                             mmin = (mmin > ROp.m(i, j)) ? maxc = j, ROp.m(i,j):    mmin;
  1552.                             temp->M(i, 1) = (double) maxr;
  1553.                             temp->M(i, 2) = (double) maxc;
  1554.                             temp->M(i, 3) = mmin;
  1555.                     }
  1556.                     break;
  1557.     }
  1558.  
  1559.     strtype newname = "min(";
  1560.     newname = newname + ROp.Getname(ROp) + ")";
  1561.     temp->Nameit(newname);
  1562.     Dispatch->Push(*temp);
  1563.     delete temp;
  1564.     return Dispatch->ReturnMat();
  1565. }
  1566.  
  1567. //////////////////////////////// elemenatary row and column operations
  1568. ////////////////// note these functions operate directly on ROp
  1569.  
  1570. void Crow(VMatrix& ROp, int row, double c)
  1571. {     // multiply a row by a constant
  1572.     ROp.Garbage("Crow");
  1573.  
  1574.     if (row < 1 || row > ROp.r)
  1575.         Dispatch->Nrerror("Crow: row index out of range");
  1576.  
  1577.     for (int i = 1; i <= ROp.c; i++) ROp.M(row, i) = c*ROp.m(row, i);
  1578.     return;
  1579. }
  1580.  
  1581. void Srow(VMatrix &ROp, int row1, int row2)
  1582. {     // swap rows
  1583.     ROp.Garbage("Srow");
  1584.     if (row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
  1585.         Dispatch->Nrerror("Srow: row index out of range");
  1586.  
  1587.     double tmp = 0;
  1588.     for (int i = 1; i <= ROp.c; i++) {
  1589.         tmp = ROp.m(row1, i);
  1590.         ROp.M(row1, i) = ROp.m(row2, i);
  1591.         ROp.M(row2, i) = tmp;
  1592.     }
  1593.     return;
  1594. }
  1595.  
  1596. void Lrow(VMatrix &ROp, int row1, int row2, double c)
  1597. {     // add c*r1 to r2
  1598.  
  1599.     ROp.Garbage("Srow");
  1600.     if (row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
  1601.         Dispatch->Nrerror("Lrow: row index out of range");
  1602.  
  1603.     for (int i = 1; i <= ROp.c; i++)
  1604.         ROp.M(row2, i) = ROp.m(row2, i) + c*ROp.m(row1, i);
  1605.  
  1606.     return;
  1607. }
  1608.  
  1609. void Ccol(VMatrix& ROp, int col, double c)
  1610. {     // multiply a col by a constant
  1611.     ROp.Garbage("Ccol");
  1612.  
  1613.     if (col < 1 || col > ROp.c)
  1614.         Dispatch->Nrerror("Ccol: col index out of range");
  1615.  
  1616.     for (int i = 1; i <= ROp.r; i++) ROp.M(i, col) = c*ROp.m(i, col);
  1617.     return;
  1618. }
  1619.  
  1620. void Scol(VMatrix &ROp, int col1, int col2)
  1621. {     // swap cols
  1622.     ROp.Garbage("Scol");
  1623.     if (col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
  1624.         Dispatch->Nrerror("Scol: col index out of range");
  1625.  
  1626.     double tmp = 0;
  1627.     for (int i = 1; i <= ROp.r; i++) {
  1628.         tmp = ROp.m(i, col1);
  1629.         ROp.M(i, col1) = ROp.m(i, col2);
  1630.         ROp.M(i, col2) = tmp;
  1631.     }
  1632.     return;
  1633. }
  1634.  
  1635. void Lcol(VMatrix &ROp, int col1, int col2, double c)
  1636. {     // add c*c1 to c2
  1637.  
  1638.     ROp.Garbage("Scol");
  1639.     if (col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
  1640.         Dispatch->Nrerror("Lcol: col index out of range");
  1641.  
  1642.     for (int i = 1; i <= ROp.r; i++)
  1643.         ROp.M(i, col2) = ROp.m(i, col2) + c*ROp.m(i, col1);
  1644.  
  1645.     return;
  1646. }
  1647.