Commit | Line | Data |
---|---|---|
0f4556f1 C |
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 | } |