+#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);
+}