added [mtx_cumprod] and [mtx_polyval] to evaluate polynomials

svn path=/trunk/externals/iem/iemmatrix/; revision=13685

6 files changed, 403 insertions, 66 deletions

diff --git a/abs/mtx_polyval-help.pd b/abs/mtx_polyval-help.pd
new file mode 100644
index 0000000..9a6ead3
--- /dev/null
+++ b/abs/mtx_polyval-help.pd
@@ -0,0 +1,14 @@
+#N canvas 0 0 514 300 10;
+#X obj 47 142 mtx_polyval;
+#X msg 122 100 matrix 1 5 1 -2 3 2 1;
+#X msg 43 73 matrix 1 3 0 1 2;
+#X obj 47 174 mtx_print;
+#X text 185 152 evaluate polynomial;
+#X text 181 178 1*x^4 -2*x^3 +3*x^2 +2*x^1 +1*x^0;
+#X text 45 13 [mtx_polyval] evaluates polynomial with given coefficients
+;
+#X text 94 33 at given values for x;
+#X text 252 241 Franz Zotter \, 2010;
+#X connect 0 0 3 0;
+#X connect 1 0 0 1;
+#X connect 2 0 0 0;
diff --git a/abs/mtx_polyval.pd b/abs/mtx_polyval.pd
index d63045d..69c719b 100644
--- a/abs/mtx_polyval.pd
+++ b/abs/mtx_polyval.pd
@@ -1,74 +1,64 @@
-#N canvas 96 236 810 668 10;
+#N canvas 0 326 645 532 10;
 #X obj 50 32 inlet;
 #X obj 310 34 inlet;
 #X text 354 32 coeffs;
 #X text 94 32 values;
 #X obj 50 60 mtx_:;
-#X obj 310 64 t a a;
-#X obj 342 86 mtx_size;
-#X msg 123 52 matrix 1 2 1 2;
-#X obj 273 547 mtx;
-#X obj 82 112 mtx_size;
-#X obj 272 245 pack f 4;
-#X obj 53 358 mtx;
-#X obj 53 379 mtx_.*;
-#X obj 103 379 t a;
-#X msg 92 316 1 \$1;
-#X obj 92 337 mtx_ones;
-#X obj 53 231 until;
-#X obj 53 273 f;
-#X obj 82 273 + 1;
-#X obj 115 437 list split 3;
-#X obj 155 459 list prepend;
-#X msg 236 437 row \$1;
-#X obj 155 480 list trim;
-#X msg 272 268 size \$1 \$2;
-#X obj 92 357 t a a b;
-#X msg 148 357 1;
-#X obj 53 293 t b f;
-#X msg 71 252 2;
-#X obj 53 210 t f b;
-#X obj 53 188 - 1;
-#X obj 53 169 f 2;
-#X obj 50 86 t b a;
-#X obj 286 588 mtx_print;
+#X obj 310 56 t a a;
+#X obj 342 78 mtx_size;
+#X obj 82 160 mtx_repmat 1 1;
+#X obj 82 135 spigot;
+#X obj 121 113 > 1;
+#X obj 50 115 t a a;
+#X obj 177 117 - 1;
+#X obj 50 184 mtx_size;
+#X obj 50 227 mtx_ones;
+#X obj 106 261 spigot;
+#X obj 50 261 spigot;
+#X obj 116 226 == 1;
+#X obj 50 389 mtx_*;
+#X msg 50 205 \$1 1;
+#X msg 177 137 1 \$1;
+#X obj 50 87 mtx_transpose;
+#X obj 106 280 mtx_concat row;
+#X obj 310 156 mtx_transpose;
+#X obj 50 427 outlet;
+#X text 181 319 powers of x:;
+#X text 192 335 [x1^0 x1^1 x1^2 ...;
+#X text 197 350 x2^0 x2^1 x2^2 ...;
+#X text 199 365 x3^0 x3^1 x3^2 ...;
+#X text 320 198 coefficients of polynomials;
+#X text 315 240 a1 b1 c1 d1 ...;
+#X text 309 223 [a0 b0 c0 d0 ...;
+#X text 315 255 a2 b3 c3 d3 ...;
+#X obj 50 319 mtx_cumprod row;
+#X obj 310 128 mtx_reverse col;
 #X connect 0 0 4 0;
 #X connect 1 0 5 0;
-#X connect 4 0 31 0;
+#X connect 4 0 20 0;
+#X connect 5 0 33 0;
 #X connect 5 1 6 0;
-#X connect 6 1 10 1;
-#X connect 6 1 30 1;
-#X connect 7 0 4 0;
-#X connect 8 0 32 0;
-#X connect 9 1 10 0;
-#X connect 9 1 14 0;
-#X connect 10 0 23 0;
-#X connect 11 0 12 0;
-#X connect 12 0 13 0;
-#X connect 13 0 12 1;
-#X connect 13 0 19 0;
-#X connect 14 0 15 0;
-#X connect 15 0 24 0;
-#X connect 16 0 17 0;
-#X connect 17 0 18 0;
-#X connect 17 0 26 0;
-#X connect 18 0 17 1;
-#X connect 19 1 20 0;
-#X connect 20 0 22 0;
-#X connect 21 0 20 1;
-#X connect 22 0 8 0;
-#X connect 23 0 8 0;
-#X connect 24 0 12 1;
-#X connect 24 1 19 0;
-#X connect 24 2 25 0;
-#X connect 25 0 21 0;
-#X connect 26 0 11 0;
-#X connect 26 1 21 0;
-#X connect 27 0 17 1;
-#X connect 28 0 16 0;
-#X connect 28 1 27 0;
-#X connect 29 0 28 0;
-#X connect 30 0 29 0;
-#X connect 31 0 29 0;
-#X connect 31 1 9 0;
-#X connect 31 1 11 1;
+#X connect 6 1 9 0;
+#X connect 6 1 11 0;
+#X connect 6 1 16 0;
+#X connect 7 0 21 1;
+#X connect 8 0 7 0;
+#X connect 9 0 8 1;
+#X connect 9 0 14 1;
+#X connect 10 0 12 0;
+#X connect 10 1 8 0;
+#X connect 11 0 19 0;
+#X connect 12 0 18 0;
+#X connect 13 0 14 0;
+#X connect 13 0 15 0;
+#X connect 14 0 21 0;
+#X connect 15 0 32 0;
+#X connect 16 0 15 1;
+#X connect 17 0 23 0;
+#X connect 18 0 13 0;
+#X connect 19 0 7 1;
+#X connect 20 0 10 0;
+#X connect 21 0 32 0;
+#X connect 22 0 17 1;
+#X connect 32 0 17 0;
+#X connect 33 0 22 0;
diff --git a/doc/mtx_cumprod-help.pd b/doc/mtx_cumprod-help.pd
new file mode 100644
index 0000000..d4f2e02
--- /dev/null
+++ b/doc/mtx_cumprod-help.pd
@@ -0,0 +1,60 @@
+#N canvas 27 49 670 511 10;
+#X text 458 17 updated for;
+#X obj 546 17 iemmatrix 0.2;
+#X obj 595 43 matrix;
+#X text 465 42 see also help for;
+#X obj 12 189 t a a;
+#X msg 12 124 bang;
+#X msg 265 390 mode row;
+#X msg 265 364 mode column;
+#X msg 352 363 mode col;
+#X msg 418 363 mode :;
+#X obj 263 416 t a;
+#X obj 12 147 mtx_ones 3 4;
+#X obj 17 379 t a a;
+#X msg 17 305 bang;
+#X obj 17 331 mtx_ones 2 2;
+#X obj 20 211 mtx_print original;
+#X obj 46 400 mtx_print original;
+#X text 410 444 see also:;
+#X msg 298 336 direction -1;
+#X obj 12 235 t a a a;
+#X obj 409 466 mtx_cumsum;
+#X text 90 15 [mtx_cumprod];
+#X text 47 34 cumulative product in row/col direction;
+#X text 18 63 you can calculate factorials or geometric series \; you
+can also do this into the reverse direction;
+#X obj 12 168 mtx_* 2;
+#X obj 17 357 mtx_* 2;
+#X obj 12 257 mtx_cumprod row;
+#X obj 138 258 mtx_cumprod col;
+#X obj 269 258 mtx_cumprod row -1;
+#X obj 17 442 mtx_cumprod row;
+#X obj 17 471 mtx_print cumprod;
+#X obj 12 285 mtx_print cumprod-row;
+#X obj 169 285 mtx_print cumprod-col;
+#X obj 330 285 mtx_print cumprod-col-reverse;
+#X text 377 484 Franz Zotter \, 2010;
+#X connect 4 0 19 0;
+#X connect 4 1 15 0;
+#X connect 5 0 11 0;
+#X connect 6 0 10 0;
+#X connect 7 0 10 0;
+#X connect 8 0 10 0;
+#X connect 9 0 10 0;
+#X connect 10 0 29 0;
+#X connect 11 0 24 0;
+#X connect 12 0 29 0;
+#X connect 12 1 16 0;
+#X connect 13 0 14 0;
+#X connect 14 0 25 0;
+#X connect 18 0 10 0;
+#X connect 19 0 26 0;
+#X connect 19 1 27 0;
+#X connect 19 2 28 0;
+#X connect 24 0 4 0;
+#X connect 25 0 12 0;
+#X connect 26 0 31 0;
+#X connect 27 0 32 0;
+#X connect 28 0 33 0;
+#X connect 29 0 30 0;
diff --git a/src/iemmatrix_sources.c b/src/iemmatrix_sources.c
index 5e05345..2e6b4f4 100644
--- a/src/iemmatrix_sources.c
+++ b/src/iemmatrix_sources.c
@@ -24,6 +24,7 @@ void iemmatrix_sources_setup(void)
 	iemtx_concat_setup(); /* mtx_concat.c */
 	iemtx_conv_setup(); /* mtx_conv.c */
 	iemtx_cos_setup(); /* mtx_cos.c */
+	iemtx_cumprod_setup(); /* mtx_cumprod.c */
 	iemtx_cumsum_setup(); /* mtx_cumsum.c */
 	iemtx_dbtopow_setup(); /* mtx_dbtopow.c */
 	iemtx_dbtorms_setup(); /* mtx_dbtorms.c */
diff --git a/src/iemmatrix_sources.h b/src/iemmatrix_sources.h
index 3e56eba..df6dc7b 100644
--- a/src/iemmatrix_sources.h
+++ b/src/iemmatrix_sources.h
@@ -22,6 +22,7 @@ void iemtx_colon_setup(void); /* mtx_colon.c */
 void iemtx_concat_setup(void); /* mtx_concat.c */
 void iemtx_conv_setup(void); /* mtx_conv.c */
 void iemtx_cos_setup(void); /* mtx_cos.c */
+void iemtx_cumprod_setup(void); /* mtx_cumprod.c */
 void iemtx_cumsum_setup(void); /* mtx_cumsum.c */
 void iemtx_dbtopow_setup(void); /* mtx_dbtopow.c */
 void iemtx_dbtorms_setup(void); /* mtx_dbtorms.c */
diff --git a/src/mtx_cumprod.c b/src/mtx_cumprod.c
new file mode 100644
index 0000000..03f0c9f
--- /dev/null
+++ b/src/mtx_cumprod.c
@@ -0,0 +1,271 @@
+/*
+ *  iemmatrix
+ *
+ *  objects for manipulating simple matrices
+ *  mostly refering to matlab/octave matrix functions
+ *
+ * Copyright (c) 2005, Franz Zotter
+ * IEM, Graz, Austria
+ *
+ * For information on usage and redistribution, and for a DISCLAIMER OF ALL
+ * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
+ *
+ */
+
+
+#include "iemmatrix.h"
+
+static t_class *mtx_cumprod_class;
+static t_symbol *row_sym;
+static t_symbol *col_sym;
+static t_symbol *col_sym2;
+
+typedef struct _MTXCumprod_ MTXCumprod;
+struct _MTXCumprod_
+{
+   t_object x_obj;
+   int rows;
+   int columns;
+   int size;
+   int cumprod_direction;
+   t_symbol *cumprod_mode;
+
+   t_outlet *list_outlet;
+
+   t_atom *list_out;
+   t_atom *list_in;
+   t_float *x;
+   t_float *y;
+};
+
+static void deleteMTXCumprod (MTXCumprod *mtx_cumprod_obj) 
+{
+   if (mtx_cumprod_obj->list_out)
+      freebytes (mtx_cumprod_obj->list_out, sizeof(t_atom)*(mtx_cumprod_obj->size+2));
+   if (mtx_cumprod_obj->x)
+      freebytes (mtx_cumprod_obj->x, sizeof(t_float)*(mtx_cumprod_obj->size));
+   if (mtx_cumprod_obj->y)
+      freebytes (mtx_cumprod_obj->y, sizeof(t_float)*(mtx_cumprod_obj->size));
+}
+
+static void mTXSetCumprodDirection (MTXCumprod *mtx_cumprod_obj, t_float c_dir)
+{
+   int direction = (int) c_dir;
+   mtx_cumprod_obj->cumprod_direction = (direction==-1)?direction:1;
+}
+
+static void mTXSetCumprodMode (MTXCumprod *mtx_cumprod_obj, t_symbol *m_sym)
+{
+   mtx_cumprod_obj->cumprod_mode = m_sym;
+}
+
+static void *newMTXCumprod (t_symbol *s, int argc, t_atom *argv)
+{
+   MTXCumprod *mtx_cumprod_obj = (MTXCumprod *) pd_new (mtx_cumprod_class);
+   mTXSetCumprodMode (mtx_cumprod_obj, gensym(":"));
+   mTXSetCumprodDirection (mtx_cumprod_obj, 1.0f);
+   if (argc>=1) {
+      if (argv[0].a_type == A_SYMBOL) {
+	 mTXSetCumprodMode (mtx_cumprod_obj, atom_getsymbol (argv));
+	 if (argc>=2) {
+	    if (argv[1].a_type != A_SYMBOL)
+	       mTXSetCumprodDirection (mtx_cumprod_obj, atom_getfloat (argv+1));
+	    else
+	       post("mtx_cumprod: 2nd arg ignored. supposed to be float");
+	 }
+      }
+      else {
+	 mTXSetCumprodDirection (mtx_cumprod_obj, atom_getfloat (argv));
+	 if (argc>=2) {
+	    if (argv[1].a_type == A_SYMBOL)
+	       mTXSetCumprodMode (mtx_cumprod_obj, atom_getsymbol (argv+1));
+	    else
+	       post("mtx_cumprod: 2nd arg ignored. supposed to be symbolic, e.g. \"row\", \"col\", \":\"");
+	 }
+      }
+   }
+
+   mtx_cumprod_obj->list_outlet = outlet_new (&mtx_cumprod_obj->x_obj, gensym("matrix"));
+   return ((void *) mtx_cumprod_obj);
+} 
+
+static void mTXCumprodBang (MTXCumprod *mtx_cumprod_obj)
+{
+   if (mtx_cumprod_obj->list_out) 
+      outlet_anything(mtx_cumprod_obj->list_outlet, gensym("matrix"), 
+	    mtx_cumprod_obj->size+2, mtx_cumprod_obj->list_out);
+}
+
+static void writeFloatIntoList (int n, t_atom *l, t_float *f) 
+{
+   for (;n--;f++, l++) 
+      SETFLOAT (l, *f);
+}
+static void readFloatFromList (int n, t_atom *l, t_float *f) 
+{
+   while (n--) 
+      *f++ = atom_getfloat (l++);
+}
+static void readFloatFromListModulo (int n, int m, t_atom *l, t_float *f) 
+{
+   t_atom *ptr = l;
+   int count1, count2;
+   n /= m;
+   count1 = m;
+   while (count1--) 
+      for (count2 = n, ptr = l++; count2--; ptr += m, f++) 
+	 *f = atom_getfloat (ptr);
+}
+static void writeFloatIntoListModulo (int n, int m, t_atom *l, t_float *f) 
+{
+   t_atom *ptr = l;
+   int count1, count2;
+   n /= m;
+   count1 = m;
+   while (count1--) 
+      for (count2 = n, ptr = l++; count2--; ptr += m, f++) 
+	 SETFLOAT(ptr,*f);
+}
+
+static void cumProd (int n, t_float *x, t_float *y)
+{
+   t_float accu = 1.0f;
+   for (;n--; x++, y++) {
+      accu *= *x;
+      *y = accu;
+   }
+}
+static void cumProdReverse (int n, t_float *x, t_float *y)
+{
+   t_float accu = 1.0f;
+   for (;n--; x--, y--) {
+      accu *= *x;
+      *y = accu;
+   }
+}
+
+static void mTXCumprodMatrix (MTXCumprod *mtx_cumprod_obj, t_symbol *s, 
+      int argc, t_atom *argv)
+{
+   int rows = atom_getint (argv++);
+   int columns = atom_getint (argv++);
+   int size = rows * columns;
+   int list_size = argc - 2;
+   t_atom *list_ptr = argv;
+   t_atom *list_out = mtx_cumprod_obj->list_out;
+   t_float *x = mtx_cumprod_obj->x;
+   t_float *y = mtx_cumprod_obj->y;
+   int count;
+
+   /* size check */
+   if (!size) {
+      post("mtx_cumprod: invalid dimensions");
+      return;
+   }
+   else if (list_size<size) {
+      post("mtx_cumprod: sparse matrix not yet supported: use \"mtx_check\"");
+      return;
+   }
+   else if ((!x)||(!list_out)||(!y)) {
+      if (!x)
+	 x = (t_float *) getbytes (sizeof (t_float) * (size));
+      if (!y)
+	 y = (t_float *) getbytes (sizeof (t_float) * (size));
+      if (!list_out)
+	 list_out = (t_atom *) getbytes (sizeof (t_atom) * (size+2));
+   }
+   else if (size != mtx_cumprod_obj->size) {
+      x = (t_float *) resizebytes (x,
+	    sizeof (t_float) * (mtx_cumprod_obj->size),
+	    sizeof (t_float) * (size));
+      y = (t_float *) resizebytes (y,
+	    sizeof (t_float) * (mtx_cumprod_obj->size),
+	    sizeof (t_float) * (size));
+      list_out = (t_atom *) resizebytes (list_out,
+	    sizeof (t_atom) * (mtx_cumprod_obj->size+2),
+	    sizeof (t_atom) * (size + 2));
+   }
+   mtx_cumprod_obj->size = size;
+   mtx_cumprod_obj->rows = rows;
+   mtx_cumprod_obj->columns = columns;
+   mtx_cumprod_obj->list_out = list_out;
+   mtx_cumprod_obj->x = x;
+   mtx_cumprod_obj->y = y;
+
+   /* main part */
+   /* reading matrix from inlet */
+   if ((mtx_cumprod_obj->cumprod_mode == col_sym) ||
+	(mtx_cumprod_obj->cumprod_mode == col_sym2)) {
+      readFloatFromListModulo (size, columns, list_ptr, x);
+      columns = mtx_cumprod_obj->rows;
+      rows = mtx_cumprod_obj->columns;
+   }
+   else
+      readFloatFromList (size, list_ptr, x);
+   
+   /* calculating cumprod */
+   if (mtx_cumprod_obj->cumprod_direction == -1) {
+      if ((mtx_cumprod_obj->cumprod_mode == row_sym) ||
+	    (mtx_cumprod_obj->cumprod_mode == col_sym) ||
+	    (mtx_cumprod_obj->cumprod_mode == col_sym2)) {
+	 x += columns-1;
+	 y += columns-1;
+
+	 for (count = rows; count--; x += columns, y += columns)
+	    cumProdReverse (columns,x,y);
+      }
+      else {
+	 x += size-1;
+	 y += size-1;
+	 cumProdReverse (size, x, y);
+      }
+   }
+   else if ((mtx_cumprod_obj->cumprod_mode == row_sym) ||
+	 (mtx_cumprod_obj->cumprod_mode == col_sym) ||
+	 (mtx_cumprod_obj->cumprod_mode == col_sym2))
+      for (count = rows; count--; x += columns, y += columns)
+	 cumProd (columns,x,y);
+   else
+      cumProd (size, x, y);
+
+   x = mtx_cumprod_obj->x;
+   y = mtx_cumprod_obj->y;
+
+   /* writing matrix to outlet */
+   if ((mtx_cumprod_obj->cumprod_mode == col_sym) ||
+	(mtx_cumprod_obj->cumprod_mode == col_sym2)) {
+      columns = mtx_cumprod_obj->columns;
+      rows = mtx_cumprod_obj->rows;
+      writeFloatIntoListModulo (size, columns, list_out+2, y);
+   }
+   else
+      writeFloatIntoList (size, list_out+2, y);
+
+   SETSYMBOL(list_out, gensym("matrix"));
+   SETFLOAT(list_out, rows);
+   SETFLOAT(&list_out[1], columns);
+   outlet_anything(mtx_cumprod_obj->list_outlet, gensym("matrix"), 
+	 mtx_cumprod_obj->size+2, list_out);
+}
+
+void mtx_cumprod_setup (void)
+{
+   mtx_cumprod_class = class_new 
+      (gensym("mtx_cumprod"),
+       (t_newmethod) newMTXCumprod,
+       (t_method) deleteMTXCumprod,
+       sizeof (MTXCumprod),
+       CLASS_DEFAULT, A_GIMME, 0);
+   class_addbang (mtx_cumprod_class, (t_method) mTXCumprodBang);
+   class_addmethod (mtx_cumprod_class, (t_method) mTXCumprodMatrix, gensym("matrix"), A_GIMME,0);
+   class_addmethod (mtx_cumprod_class, (t_method) mTXSetCumprodMode, gensym("mode"), A_DEFSYMBOL,0);
+   class_addmethod (mtx_cumprod_class, (t_method) mTXSetCumprodDirection, gensym("direction"), A_DEFFLOAT,0);
+
+   row_sym = gensym("row");
+   col_sym = gensym("col");
+   col_sym2 = gensym("column");
+}
+
+void iemtx_cumprod_setup(void){
+  mtx_cumprod_setup();
+}