home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 13 / MA_Cover_13.bin / source / c / apl / a7.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-11-23  |  2.3 KB  |  125 lines
  1. #include "apl.h"
  2.  
  3. ex_iprod()
  4. {
  5.     int i, j, param[10], ipr1();
  6.     struct item *p, *q, *r;
  7.     data (*fn)();
  8.  
  9.     param[0] = exop[*pcp++];
  10.     param[1] = exop[*pcp++];
  11.     p = fetch2();
  12.     q = sp[-2];
  13.     if(p->type != DA || q->type != DA) error("iprod T");
  14.     /*
  15.      * extend scalars to match corresponding arg
  16.      */
  17.     if(scalar(p)) {
  18.         if(scalar(q)){
  19.             r = newdat(DA, 0, 1);
  20.             fn = param[1];
  21.             r->datap[0] = (*fn)(p->datap[0], q->datap[0]);
  22.             goto out;
  23.         }
  24.         r = extend(DA, q->dim[0], p->datap[0]);
  25.         pop();
  26.         *sp++ = p = r;
  27.     }
  28.     if(scalar(q)){
  29.         r = extend(DA,p->dim[p->rank - 1], q->datap[0]);
  30.         aplfree(sp[-2]);
  31.         sp[-2] = q = r;
  32.     }
  33.     bidx(p);
  34.     idx.rank--;
  35.     param[2] = idx.dim[idx.rank];
  36.     if((param[2] != q->dim[0])) error("inner prod C");
  37.     param[3] = q->size/param[2];
  38.     for(i=1; i<q->rank; i++) idx.dim[idx.rank++] = q->dim[i];
  39.     r = newdat(DA, idx.rank, size());
  40.     copy(IN, idx.dim, r->dim, idx.rank);
  41.     param[4] = 0;
  42.     param[5] = 0;
  43.     param[6] = p->datap;
  44.     param[7] = q->datap;
  45.     param[8] = r->datap;
  46.     param[9] = p->size;
  47.     forloop(ipr1, param);
  48. out:
  49.     pop();
  50.     pop();
  51.     /*
  52.      * KLUDGE (we need the dim[0]'s for above stuff to work)
  53.      */
  54.     if(r->rank == 1 && r->size == 1) r->rank = 0;
  55.     *sp++ = r;
  56. }
  57.  
  58. ipr1(param)
  59. int param[];
  60. {
  61.     int i, dk, lk, a, b;
  62.     data *dp1, *dp2, *dp3;
  63.     data (*f1)(), (*f2)(), d;
  64.  
  65.     f1 = param[0];
  66.     f2 = param[1];
  67.     dk = param[2];
  68.     lk = param[3];
  69.     a = param[4];
  70.     b = param[5];
  71.     dp1 = param[6];
  72.     dp2 = param[7];
  73.     dp3 = param[8];
  74.     a += dk;
  75.     b += (dk * lk);
  76.     for(i=0; i<dk; i++) {
  77.         a--;
  78.         b -= lk;
  79.         d = (*f2)(dp1[a], dp2[b]);
  80.         if(i == 0) datum = d;
  81.         else datum = (*f1)(d, datum);
  82.     }
  83.     *dp3++ = datum;
  84.     param[8] = dp3;
  85.     param[5]++;
  86.     if(param[5] >= lk) {
  87.         param[5] = 0;
  88.         param[4] += dk;
  89.         if(param[4] >= param[9]) param[4] = 0;
  90.     }
  91. }
  92.  
  93. ex_oprod()
  94. {
  95.     int i, j;
  96.     data *dp, *dp1, *dp2;
  97.     struct item *p, *q, *r;
  98.     data (*f)();
  99.  
  100.     f = (data *)exop[*pcp++];
  101.     p = fetch2();
  102.     q = sp[-2];
  103.     if(p->type != DA || q->type != DA) error("oprod T");
  104.     /*
  105.      * collapse 1 element vectors to scalars
  106.      *
  107.     if(scalar(p)) p->rank = 0;
  108.     if(scalar(q)) q->rank = 0;
  109.     */
  110.     bidx(p);
  111.     for(i=0; i<q->rank; i++) idx.dim[idx.rank++] = q->dim[i];
  112.     r = newdat(DA, idx.rank, size());
  113.     copy(IN, idx.dim, r->dim, idx.rank);
  114.     dp = r->datap;
  115.     dp1 = p->datap;
  116.     for(i=0; i<p->size; i++) {
  117.         datum = *dp1++;
  118.         dp2 = q->datap;
  119.         for(j=0; j<q->size; j++) *dp++ = (*f)(datum, *dp2++);
  120.     }
  121.     pop();
  122.     pop();
  123.     *sp++ = r;
  124. }
  125.