home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d3xx / d386 / xlispstat.lha / XLispStat / src2.lzh / XLisp-Stat / matrices1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-10-05  |  16.0 KB  |  658 lines
  1. /* matrices1 - Elementary matrix operations                            */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlproto.h"
  11. #include "xlsproto.h"
  12. #else
  13. #include "xlfun.h"
  14. #include "xlsfun.h"
  15. #endif ANSI
  16. #include "xlsvar.h"
  17.  
  18. /* forward declarations */
  19. #ifdef ANSI
  20. LVAL matmult(void),transpose_list(void);
  21. #else
  22. LVAL matmult(),transpose_list();
  23. #endif ANSI
  24.  
  25. /*************************************************************************/
  26. /**                                                                     **/
  27. /**                       Matrix Data Type                              **/
  28. /**                                                                     **/
  29. /*************************************************************************/
  30.  
  31. /* Many routines here assume that displacedarraydim returns a vector.    */
  32.  
  33. /* is a a matrix? */
  34. int matrixp(a)
  35.      LVAL a;
  36. {
  37.   return(displacedarrayp(a) && (arrayrank(a) == 2));
  38. }
  39.  
  40. /* get a matrix from the argument stack */
  41. LVAL xsgetmatrix()
  42. {
  43.   LVAL m;
  44.  
  45.   m = xlgetarg();
  46.   if (! matrixp(m)) xlerror("not a matrix", m);
  47.   return(m);
  48. }
  49.  
  50. /* number of rows in a matrix */
  51. int numrows(a)
  52.      LVAL a;
  53. {
  54.   return((int) getfixnum(getelement(displacedarraydim(a), 0)));
  55. }
  56.  
  57. /* number of columns in a matrix */
  58. int numcols(a)
  59.      LVAL a;
  60. {
  61.   return((int) getfixnum(getelement(displacedarraydim(a), 1)));
  62. }
  63.  
  64. /*************************************************************************/
  65. /**                                                                     **/
  66. /**                   Matrix Multiplication Functions                   **/
  67. /**                                                                     **/
  68. /*************************************************************************/
  69.  
  70. /* numerical inner product. Result is always of type FLOAT*/
  71. LVAL innerproduct(x, y)
  72.      LVAL x, y;
  73. {
  74.   int n, i;
  75.   double val, xi, yi;
  76.   int complex = FALSE;
  77.   Complex cval, cxi, cyi;
  78.   
  79.   if (! vectorp(x)) xlerror("not a vector", x);
  80.   if (! vectorp(y)) xlerror("not a vector", y);
  81.   if (getsize(x) != getsize(y)) xlfail("vector lengths do not match");
  82.   if (getsize(x) == 0) xlfail("vectors are too short");
  83.   n = getsize(x);
  84.  
  85.   /* check for a complex argument */
  86.   for (i = 0; i < n && ! complex; i++) {
  87.     if (complexp(getelement(x, i))) complex = TRUE;
  88.     if (complexp(getelement(y, i))) complex = TRUE;
  89.   }
  90.   
  91.   if (complex) {
  92.     cval = cart2complex(0.0, 0.0);
  93.     for (i = 0; i < n; i++) {
  94.       cxi = makecomplex(getelement(x, i));
  95.       cyi = makecomplex(getelement(y, i));
  96.       cval = cadd(cval, cmul(cxi, cyi));
  97.     }
  98.     return(cvcomplex(cval));
  99.   }
  100.   else {
  101.     for (val = 0.0, i = 0; i < n; i++) {
  102.       xi = makedouble(getelement(x, i));
  103.       yi = makedouble(getelement(y, i));
  104.       val = f_plus(val, f_times(xi, yi));
  105.     }
  106.     return(cvflonum((FLOTYPE) val));
  107.   }
  108. }
  109.  
  110. /* copy row r of matrix a into the vector x */
  111. void copy_row(a, r, x)
  112.      LVAL a, x;
  113.      int r;
  114. {
  115.   int cols, i;
  116.   LVAL data;
  117.  
  118.   if (! matrixp(a)) xlerror("not a matrix", a);
  119.   if (! vectorp(x)) xlerror("not a vector", x);
  120.   if (numcols(a) != getsize(x)) xlfail("dimensions do not match");
  121.     
  122.   cols = numcols(a);
  123.   data = arraydata(a);
  124.     
  125.   for (i = 0; i < cols; i++)
  126.     setelement(x, i, getelement(data, cols * r + i));
  127. }
  128.  
  129. /* copy column c of matrix a into vector x */
  130. void copy_column(a, c, x)
  131.      LVAL a, x;
  132.      int c;
  133. {
  134.   int rows, cols, i;
  135.   LVAL data;
  136.  
  137.   if (! matrixp(a)) xlerror("not a matrix", a);
  138.   if (! vectorp(x)) xlerror("not a vector", x);
  139.   if (numrows(a) != getsize(x)) xlfail("dimensions do not match");
  140.     
  141.   rows = numrows(a);
  142.   cols = numcols(a);
  143.   data = arraydata(a);
  144.     
  145.   for (i = 0; i < rows; i++)
  146.     setelement(x, i, getelement(data, cols * i + c));
  147. }
  148.  
  149. /* Built in MATMULT function. Result is always of type FLOAT */
  150. static LVAL matmult()
  151. {
  152.   LVAL x, y, val, dim, nn, row, col, mm, result, result_data;
  153.   int i, j, n, m;
  154.   int list_arg = FALSE;
  155.     
  156.   /* protect some pointers */
  157.   xlstkcheck(9);
  158.   xlsave(x);
  159.   xlsave(y);
  160.   xlsave(dim);
  161.   xlsave(nn);
  162.   xlsave(mm);
  163.   xlsave(val);
  164.   xlsave(row);
  165.   xlsave(col);
  166.   xlsave(result);
  167.  
  168.   /* scalar multiplication */
  169.   if (numberp(peekarg(0)) || numberp(peekarg(1))) result = xsrmul();
  170.   else {
  171.     x = xlgetarg();
  172.     y = xlgetarg();
  173.     xllastarg();
  174.     
  175.     /* coerce lists to vectors and check */
  176.     if (! arrayp(x) && consp(x)) { list_arg = TRUE; x = coerce_to_vector(x); }
  177.     if (! arrayp(y) && consp(y)) { list_arg = TRUE; y = coerce_to_vector(y); }
  178.     if ((! vectorp(x)) && (! matrixp(x))) xlbadtype(x);
  179.     if ((! vectorp(y)) && (! matrixp(y))) xlbadtype(y);
  180.  
  181.     /* simple inner product */
  182.     if (simplevectorp(x) && simplevectorp(y)) {
  183.       if (getsize(x) != getsize(y))
  184.     xlfail("dimensions do not match");
  185.       result = innerproduct(x, y);
  186.     }
  187.     else {
  188.  
  189.       /* check dimensions */
  190.       if (simplevectorp(x) && matrixp(y) && (getsize(x) != numrows(y)))
  191.     xlfail("dimensions do not match");
  192.       if (matrixp(x) && simplevectorp(y) && (numcols(x) != getsize(y)))
  193.     xlfail("dimensions do not match");
  194.       if (matrixp(x) && matrixp(y) && (numcols(x) != numrows(y)))
  195.     xlfail("dimensions do not match");
  196.  
  197.       /* compute result dimensions */
  198.       nn = (simplevectorp(x)) ? cvfixnum((FIXTYPE) 1)
  199.                           : getelement(displacedarraydim(x), 0);
  200.       mm = (simplevectorp(y)) ? cvfixnum((FIXTYPE) 1)
  201.                               : getelement(displacedarraydim(y), 1);
  202.       n = getfixnum(nn);
  203.       m = getfixnum(mm);
  204.       dim = list2(nn, mm);
  205.       result = newarray(dim, NIL, NIL);
  206.     
  207.       /* set up the vectors for holding rows and columns for innerproduct */   
  208.       row = (simplevectorp(x)) ? x : newvector(numcols(x));
  209.       col = (simplevectorp(y)) ? y : newvector(numrows(y));
  210.     
  211.       /* compute the result matrix */
  212.       result_data = arraydata(result);
  213.       for (i = 0; i < n; i++)
  214.     for (j = 0; j < m; j++) {
  215.       if (! simplevectorp(x)) copy_row(x, i, row);
  216.       if (! simplevectorp(y)) copy_column(y, j, col);
  217.       val = innerproduct(row, col);
  218.       setelement(result_data, m * i + j, val);
  219.     }
  220.     }
  221.   }
  222.   
  223.   /* reformat result if one of inputs was a sequence, one a matrix */
  224.   if (matrixp(result) && (simplevectorp(x) || simplevectorp(y))) {
  225.     result = (list_arg) ? coerce_to_list(arraydata(result))
  226.                         : arraydata(result);
  227.   }
  228.   
  229.   /* restore the stack frame */
  230.   xlpopn(9);
  231.  
  232.   return(result);
  233. }
  234.  
  235. LVAL xsmatmult()
  236. {
  237.   LVAL args, fcn, result;
  238.  
  239.   if (xlargc <= 2) result = matmult();
  240.   else {
  241.     xlstkcheck(2);
  242.     xlsave1(args);
  243.     xlsave1(fcn);
  244.     fcn = cvsubr(matmult, SUBR, 0);
  245.     args = makearglist(xlargc, xlargv);
  246.     result = reduce(fcn, args, FALSE, NIL);
  247.     xlpopn(2);
  248.   }
  249.   return(result);
  250. }
  251.     
  252. /* Built in CROSS-PRODUCT function */
  253. LVAL xscrossproduct()
  254. {
  255.   LVAL x, val, dim, nn, col_i, col_j, result, result_data;
  256.   int i, j, n;
  257.     
  258.   x = xlgetarg();
  259.   xllastarg();
  260.  
  261.   if (simplevectorp(x)) result = innerproduct(x, x);
  262.   else if (matrixp(x)) {
  263.  
  264.     /* save some pointers */
  265.     xlstkcheck(6);
  266.     xlsave(dim);
  267.     xlsave(nn);
  268.     xlsave(val);
  269.     xlsave(col_i);
  270.     xlsave(col_j);
  271.     xlsave(result);
  272.  
  273.     /* determine dimensions and set up result */
  274.     nn = getelement(displacedarraydim(x), 1);
  275.     n = getfixnum(nn);
  276.     dim = list2(nn, nn);
  277.     result = newarray(dim, NIL, NIL);
  278.     
  279.     /* set up vectors for holding columns for inner product */
  280.     col_i = newvector(numrows(x));
  281.     col_j = newvector(numrows(x));
  282.  
  283.     /* compute the result */
  284.     result_data = arraydata(result);
  285.     for (i = 0; i < n; i++)
  286.       for (j = 0; j < n; j++) {
  287.     copy_column(x, i, col_i);
  288.     copy_column(x, j, col_j);
  289.     val = innerproduct(col_i, col_j);
  290.     setelement(result_data, n * i + j, val);
  291.     }
  292.         
  293.     /* restore the stack frame */
  294.     xlpopn(6);
  295.   }
  296.   else xlbadtype(x);
  297.  
  298.   return(result);
  299. }
  300.  
  301. /* Built in OUTER-PRODUCT finction */
  302. LVAL xsouterproduct()
  303. {
  304.   LVAL x, y, f, dim, val, result, result_data, xe, ye;
  305.   int i, j, n, m;
  306.     
  307.   /* protect some pointers */
  308.   xlstkcheck(5);
  309.   xlsave(x);
  310.   xlsave(y);
  311.   xlsave(dim);
  312.   xlsave(val);
  313.   xlsave(result);
  314.  
  315.   x = xlgetarg();
  316.   x = coerce_to_vector(x);
  317.   y = xlgetarg();
  318.   y = coerce_to_vector(y);
  319.   f = (moreargs()) ? xlgetarg() : NIL;
  320.   if (! simplevectorp(x)) xlerror("not a simple vector", x);
  321.   if (! simplevectorp(y)) xlerror("not a simple vector", y);
  322.   xllastarg();
  323.     
  324.   n = getsize(x);
  325.   m = getsize(y);
  326.   dim = integer_list_2(n, m);
  327.   result = newarray(dim, NIL, NIL);
  328.     
  329.   result_data = arraydata(result);
  330.   for (i = 0; i < n; i++)
  331.     for (j = 0; j < m; j++) {
  332.       xe = getelement(x, i);
  333.       ye = getelement(y, j);
  334.       val = (f == NIL) ? xscallsubr2(xsrmul, xe, ye) : xsfuncall2(f, xe, ye);
  335.       setelement(result_data, m * i + j, val);
  336.     }
  337.         
  338.   /* restore the stack frame */
  339.   xlpopn(5);
  340.  
  341.   return(result);
  342. }
  343.  
  344. /*************************************************************************/
  345. /**                                                                     **/
  346. /**           Matrix Construction and Decomposition Functions           **/
  347. /**                                                                     **/
  348. /*************************************************************************/
  349.  
  350. /* Internal version of DIAGONAL function. Given a matrix returns the */
  351. /* diagonal; given a sequence returns a diagonal matrix.             */
  352. LVAL diagonal(arg)
  353.      LVAL arg;
  354. {
  355.   LVAL next, dim, val, data, result, result_data;
  356.   int n, m, i;
  357.  
  358.   /* protect some pointers */
  359.   xlstkcheck(3);
  360.   xlsave(dim);
  361.   xlsave(val);
  362.   xlsave(result);
  363.     
  364.   if (matrixp(arg)) {
  365.  
  366.     /* extract diagonal from a matrix */
  367.     n = (numrows(arg) < numcols(arg)) ? numrows(arg) : numcols(arg);
  368.     m = numcols(arg);
  369.     result = mklist(n, NIL);
  370.     data = arraydata(arg);
  371.     for (i = 0, next = result; i < n; i++, next = cdr(next))
  372.       rplaca(next, getelement(data, m * i + i));
  373.   }
  374.   else if (sequencep(arg)) {
  375.  
  376.     /* construct a diagonal matrix */
  377.     n = (vectorp(arg)) ? getsize(arg) : llength(arg);
  378.     dim = cvfixnum((FIXTYPE) n);
  379.     dim = list2(dim, dim);
  380.     val = cvfixnum((FIXTYPE) 0);
  381.     result = newarray(dim, s_ielement, val);
  382.     result_data = arraydata(result);
  383.     for (i = 0; i < n; i++)
  384.       setelement(result_data, n * i + i, getnextelement(&arg, i));
  385.   }
  386.   else xlbadtype(arg);
  387.     
  388.   /* restore the stack frame */
  389.   xlpopn(3);
  390.  
  391.   return(result);
  392. }
  393.  
  394. /* Built in DIAGONAL function */
  395. LVAL xsdiagonal()
  396. {
  397.   LVAL arg;
  398.  
  399.   arg = xlgetarg();
  400.   xllastarg();
  401.  
  402.   return(diagonal(arg));
  403. }
  404.  
  405. /* Internal IDENTITY-MATRIX function */
  406. LVAL identitymatrix(n)
  407.      int n;
  408. {
  409.   LVAL result, val;
  410.  
  411.   /* protect some pointers */
  412.   xlstkcheck(2);
  413.   xlsave(val);
  414.   xlsave(result);
  415.  
  416.   val = cvfixnum((FIXTYPE) 1);
  417.   result = mklist(n, val);
  418.   result = diagonal(result);
  419.   
  420.   /* restore the stack frame */
  421.   xlpopn(2);
  422.  
  423.   return(result);
  424. }
  425.  
  426. /* Built in IDENTITY-MATRIX function */
  427. LVAL xsidentitymatrix()
  428. {
  429.   int n;
  430.     
  431.   n = getfixnum(checknonnegint(xlgetarg()));
  432.   xllastarg();
  433.     
  434.   return(identitymatrix(n));
  435. }
  436.  
  437. /* Return a list of rows or columns of a matrix read from the stack */
  438. LVAL facelist(face)
  439.      int face;
  440. {
  441.   LVAL a, result, next, vect, data;
  442.   int rows, cols, i, j;
  443.     
  444.   a = xsgetmatrix();
  445.   xllastarg();
  446.  
  447.   rows = numrows(a);
  448.   cols = numcols(a);
  449.     
  450.   /* protect some pointers */
  451.   xlsave1(result);
  452.  
  453.   data = arraydata(a);
  454.   switch(face) {
  455.   case 0: /* rows */
  456.     result = mklist(rows, NIL);
  457.     for (next = result, i = 0; i < rows; i++, next = cdr(next)) {
  458.       vect = newvector(cols);
  459.       rplaca(next, vect);
  460.       for (j = 0; j < cols; j++) 
  461.     setelement(vect, j, getelement(data, cols * i + j));
  462.     }
  463.     break;
  464.   case 1: /* columns */
  465.     result = mklist(cols, NIL);
  466.     for (next = result, j = 0; j < cols; j++, next = cdr(next)) {
  467.       vect = newvector(rows);
  468.       rplaca(next, vect);
  469.       for (i = 0; i < rows; i++) 
  470.     setelement(vect, i, getelement(data, cols * i + j));
  471.     }
  472.     break;
  473.   default:
  474.     xlfail(" bad face selector");
  475.   }
  476.     
  477.   /* restore the stack frame */
  478.   xlpop();
  479.  
  480.   return(result);
  481. }
  482.  
  483. /* Built in ROW-LIST and COLUMN-LIST functions */
  484. LVAL xsrowlist() { return(facelist(0)); }
  485. LVAL xscolumnlist() { return(facelist(1)); }
  486.  
  487. /* Bind list of sequences or matrices along rows or columns */
  488. LVAL xsbindfaces(face)
  489.      int face;
  490. {
  491.   LVAL next, data, dim, result, result_data;
  492.   int totalsize, rows, cols, c, r, n, i, j;
  493.   
  494.   /* protect some pointers */
  495.   xlstkcheck(3);
  496.   xlsave(data);
  497.   xlsave(dim);
  498.   xlsave(result);
  499.   
  500.   /* Check the first argument and establish size of the binding face */
  501.   next = peekarg(0);
  502.   switch (face) {
  503.   case 0:
  504.     if (matrixp(next)) cols = numcols(next);
  505.     else if (sequencep(next)) cols = seqlen(next);
  506.     else if (! compoundp(next)) cols = 1;
  507.     else xlbadtype(next);
  508.     break;
  509.   case 1:
  510.     if (matrixp(next)) rows = numrows(next);
  511.     else if (sequencep(next)) rows = seqlen(next);
  512.     else if (! compoundp(next)) rows = 1;
  513.     else xlbadtype(next);
  514.     break;
  515.   }
  516.  
  517.   /* Pass through the arguments on the stack to determine the result size */
  518.   n = xlargc;
  519.   for (i = 0, totalsize = 0; i < n; i++) {
  520.     next = peekarg(i);
  521.     if (matrixp(next)) {
  522.       c = numcols(next);
  523.       r = numrows(next); 
  524.     }
  525.     else if (sequencep(next))
  526.       switch (face) {
  527.       case 0:  c = seqlen(next); r = 1; break;
  528.       case 1:  c = 1; r = seqlen(next); break;
  529.       }
  530.     else if (! compoundp(next)) {
  531.       c = 1;
  532.       r = 1;
  533.     }
  534.     else xlbadtype(next);
  535.  
  536.     switch (face) {
  537.     case 0:
  538.       if (c != cols) xlfail("dimensions do not match");
  539.       else totalsize += r;
  540.       break;
  541.     case 1:
  542.       if (r != rows) xlfail("dimensions do not match");
  543.       else totalsize += c;
  544.     }
  545.   }
  546.   
  547.   /* set up the result matrix */
  548.   dim = newvector(2);
  549.     switch (face) {
  550.     case 0:
  551.       setelement(dim, 0, cvfixnum((FIXTYPE) totalsize));
  552.       setelement(dim, 1, cvfixnum((FIXTYPE) cols));
  553.       break;
  554.     case 1:
  555.       setelement(dim, 0, cvfixnum((FIXTYPE) rows));
  556.       setelement(dim, 1, cvfixnum((FIXTYPE) totalsize));
  557.       break;
  558.     }
  559.   result = newarray(dim, NIL, NIL);
  560.   result_data = arraydata(result);
  561.  
  562.   /* compute the result */
  563.   for (r = 0, c = 0; moreargs();) {
  564.     next = xlgetarg();
  565.     if (matrixp(next)) {
  566.       rows = numrows(next);
  567.       cols = numcols(next);
  568.       data = arraydata(next);
  569.     }
  570.     else {
  571.       switch (face) {
  572.       case 0: rows = 1; break;
  573.       case 1: cols = 1; break;
  574.       }
  575.       data = coerce_to_vector(next);
  576.     }
  577.     switch (face) {
  578.     case 0:
  579.       for (i = 0; i < rows; i++, r++) 
  580.     for (j = 0; j < cols; j++)
  581.       setelement(result_data, cols * r + j,
  582.              getelement(data, cols * i + j));
  583.       break;
  584.     case 1:
  585.       for (j = 0; j < cols; j++, c++)
  586.     for (i = 0; i < rows; i++) 
  587.       setelement(result_data, totalsize * i + c,
  588.              getelement(data, cols * i + j));
  589.       break;
  590.     }
  591.   }
  592.   
  593.   /* restore the stack frame */
  594.   xlpopn(3);
  595.   
  596.   return(result);
  597. }
  598.  
  599. /* Built in BIND-ROWS and BIND-COLUMNS functions */
  600. LVAL xsbindrows() { return(xsbindfaces(0)); }
  601. LVAL xsbindcols() { return(xsbindfaces(1)); }
  602.  
  603. /* Built in TRANSPOSE function */
  604. /* FORWARD LVAL transpose_list(); above JKL *?
  605. LVAL xstranspose()
  606. {
  607.   LVAL m, perm, result;
  608.  
  609.   if (consp(peekarg(0))) return(transpose_list());
  610.  
  611.   m = xsgetmatrix();
  612.   xllastarg();
  613.  
  614.   xlsave1(perm);
  615.   perm = newvector(2);
  616.   setelement(perm, 0, cvfixnum((FIXTYPE) 1));
  617.   setelement(perm, 1, cvfixnum((FIXTYPE) 0));
  618.   result = xscallsubr2(xspermutearray, m, perm);
  619.   xlpop();
  620.  
  621.   return(result);
  622. }
  623. /* above JKL
  624. extern LVAL copylist();
  625. */
  626. LOCAL LVAL transpose_list()
  627. {
  628.   LVAL list, result, nextr, row, nextl;
  629.   int m, n;
  630.   
  631.   list = xlgalist();
  632.   xllastarg();
  633.   
  634.   xlstkcheck(2);
  635.   xlsave(result);
  636.   xlprotect(list);
  637.   
  638.   list = copylist(list);
  639.   m = llength(list);
  640.   if (! consp(car(list))) xlerror("not a list", car(list));
  641.   n = llength(car(list));
  642.   
  643.   result = mklist(n, NIL);
  644.   for (nextr = result; consp(nextr); nextr = cdr(nextr)) {
  645.     row = mklist(m, NIL);
  646.     rplaca(nextr, row);
  647.     for (nextl = list; consp(nextl); nextl = cdr(nextl)) {
  648.       if (!consp(car(nextl))) xlerror("not a list", car(nextl));
  649.       rplaca(row, car(car(nextl)));
  650.       row = cdr(row);
  651.       rplaca(nextl, cdr(car(nextl)));
  652.     }
  653.   }
  654.   
  655.   xlpopn(2);
  656.   return(result);
  657. }
  658.