BSD 3 development

[unix-history] / usr / src / cmd / apl / a7.c

#include "apl.h"

ex_iprod()
{
       register i, j;
       struct item *p, *q, *r;
       int param[10], ipr1();

       param[0] = exop[*pcp++];
       param[1] = exop[*pcp++];
       p = fetch2();
       q = sp[-2];
       if(p->type != DA || q->type != DA)
               error("iprod T");
       bidx(p);
       idx.rank--;
       param[2] = idx.dim[idx.rank];
       if(param[2] != q->dim[0])
               error("inner prod C");
       param[3] = q->size/param[2];
       for(i=1; i<q->rank; i++)
               idx.dim[idx.rank++] = q->dim[i];
       r = newdat(DA, idx.rank, size());
       copy(IN, idx.dim, r->dim, idx.rank);
       param[4] = 0;
       param[5] = 0;
       param[6] = p->datap;
       param[7] = q->datap;
       param[8] = r->datap;
       param[9] = p->size;
       forloop(ipr1, param);
       pop();
       pop();
       push(r);
}

ipr1(param)
int param[];
{
       register i, dk;
       int lk, a, b;
       data *dp1, *dp2, *dp3;
       data (*f1)(), (*f2)(), d;

       f1 = param[0];
       f2 = param[1];
       dk = param[2];
       lk = param[3];
       a = param[4];
       b = param[5];
       dp1 = param[6];
       dp2 = param[7];
       dp3 = param[8];
       a =+ dk;
       b =+ (dk * lk);
       for(i=0; i<dk; i++) {
               a--;
               b =- lk;
               d = (*f2)(dp1[a], dp2[b]);
               if(i == 0)
                       datum = d; else
                       datum = (*f1)(d, datum);
       }
       *dp3++ = datum;
       param[8] = dp3;
       param[5]++;
       if(param[5] >= lk) {
               param[5] = 0;
               param[4] =+ dk;
               if(param[4] >= param[9])
                       param[4] = 0;
       }
}

ex_oprod()
{
       register i, j;
       register data *dp;
       struct item *p, *q, *r;
       data *dp1, *dp2;
       data (*f)();

       f = exop[*pcp++];
       p = fetch2();
       q = sp[-2];
       if(p->type != DA || q->type != DA)
               error("oprod T");
       bidx(p);
       for(i=0; i<q->rank; i++)
               idx.dim[idx.rank++] = q->dim[i];
       r = newdat(DA, idx.rank, size());
       copy(IN, idx.dim, r->dim, idx.rank);
       dp = r->datap;
       dp1 = p->datap;
       for(i=0; i<p->size; i++) {
               datum = *dp1++;
               dp2 = q->datap;
               for(j=0; j<q->size; j++)
                       *dp++ = (*f)(datum, *dp2++);
       }
       pop();
       pop();
       push(r);
}