|
|
1.1 ! root 1: static char Sccsid[] = "ae.c @(#)ae.c 1.1 10/1/82 Berkeley "; ! 2: #include "apl.h" ! 3: ! 4: ! 5: char base_com[] = {ADD, MUL}; ! 6: ! 7: ex_base() ! 8: { ! 9: struct item *extend(); ! 10: register struct item *p, *q; ! 11: register i; ! 12: char *savpcp; ! 13: data d1, d2; ! 14: ! 15: p = fetch2(); ! 16: q = sp[-2]; ! 17: if(scalar(p)){ ! 18: if(q->rank > 0) ! 19: i = q->dim[0]; ! 20: else ! 21: i = q->size; ! 22: q = extend(DA, i, p->datap[0]); ! 23: pop(); ! 24: *sp++ = p = q; ! 25: q = sp[-2]; ! 26: } ! 27: d1 = p->datap[p->size-1]; ! 28: p->datap[p->size-1] = 1.0; ! 29: for(i=p->size-2; i>= 0; i--){ ! 30: d2 = p->datap[i]; ! 31: p->datap[i] = d1; ! 32: d1 *= d2; ! 33: } ! 34: savpcp = pcp; ! 35: pcp = base_com; ! 36: ex_iprod(); ! 37: pcp = savpcp; ! 38: } ! 39: ! 40: ex_rep() ! 41: { ! 42: register struct item *p, *q; ! 43: struct item *r; ! 44: double d1, d2, d3; ! 45: data *p1, *p2, *p3; ! 46: ! 47: p = fetch2(); ! 48: q = sp[-2]; ! 49: /* ! 50: * first map 1 element vectors to scalars: ! 51: * ! 52: if(scalar(p)) ! 53: p->rank = 0; ! 54: if(scalar(q)) ! 55: q->rank = 0; ! 56: */ ! 57: r = newdat(DA, p->rank+q->rank, p->size*q->size); ! 58: copy(IN, p->dim, r->dim, p->rank); ! 59: copy(IN, q->dim, r->dim+p->rank, q->rank); ! 60: p3 = &r->datap[r->size]; ! 61: for(p1 = &p->datap[p->size]; p1 > p->datap; ){ ! 62: d1 = *--p1; ! 63: if(d1 == 0.0) ! 64: d1 = 1.0e38; /* all else goes here */ ! 65: for(p2 = &q->datap[q->size]; p2 > q->datap; ){ ! 66: d2 = *--p2; ! 67: d3 = d2 /= d1; ! 68: *p2 = d2 = floor(d2); ! 69: *--p3 = (d3 - d2)*d1; ! 70: } ! 71: } ! 72: pop(); ! 73: pop(); ! 74: *sp++ = r; ! 75: } ! 76: ! 77: /* ! 78: * scalar -- return true if arg is a scalar ! 79: */ ! 80: scalar(aip) ! 81: struct item *aip; ! 82: { ! 83: return(aip->size == 1); ! 84: } ! 85: ! 86: ! 87: struct item * ! 88: extend(ty, n, d) ! 89: data d; ! 90: { ! 91: register i; ! 92: register struct item *q; ! 93: ! 94: if(ty != DA) ! 95: error("extend T"); ! 96: q = newdat(ty, 1, n); ! 97: for(i=0; i<n; i++) ! 98: q->datap[i] = d; ! 99: return(q); ! 100: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.