1 |
|
2 |
|
3 |
|
4 |
|
5 |
|
6 |
|
7 |
|
8 |
|
9 |
|
10 |
|
11 |
|
12 |
|
13 |
|
14 |
|
15 |
|
16 |
|
17 |
|
18 |
|
19 |
|
20 |
|
21 | Eval_contract = ->
|
22 | push(cadr(p1))
|
23 | Eval()
|
24 | if (cddr(p1) == symbol(NIL))
|
25 | push_integer(1)
|
26 | push_integer(2)
|
27 | else
|
28 | push(caddr(p1))
|
29 | Eval()
|
30 | push(cadddr(p1))
|
31 | Eval()
|
32 | contract()
|
33 |
|
34 | contract = ->
|
35 | save()
|
36 | yycontract()
|
37 | restore()
|
38 |
|
39 | yycontract = ->
|
40 | h = 0
|
41 | i = 0
|
42 | j = 0
|
43 | k = 0
|
44 | l = 0
|
45 | m = 0
|
46 | n = 0
|
47 | ndim = 0
|
48 | nelem = 0
|
49 | ai = []
|
50 | an = []
|
51 |
|
52 | p3 = pop()
|
53 | p2 = pop()
|
54 | p1 = pop()
|
55 |
|
56 | if (!istensor(p1))
|
57 | if (!iszero(p1))
|
58 | stop("contract: tensor expected, 1st arg is not a tensor")
|
59 | push(zero)
|
60 | return
|
61 |
|
62 | push(p2)
|
63 | l = pop_integer()
|
64 |
|
65 | push(p3)
|
66 | m = pop_integer()
|
67 |
|
68 | ndim = p1.tensor.ndim
|
69 |
|
70 | if (l < 1 || l > ndim || m < 1 || m > ndim || l == m \
|
71 | || p1.tensor.dim[l - 1] != p1.tensor.dim[m - 1])
|
72 | stop("contract: index out of range")
|
73 |
|
74 | l--
|
75 | m--
|
76 |
|
77 | n = p1.tensor.dim[l]
|
78 |
|
79 |
|
80 |
|
81 | nelem = 1
|
82 | for i in [0...ndim]
|
83 | if (i != l && i != m)
|
84 | nelem *= p1.tensor.dim[i]
|
85 |
|
86 | p2 = alloc_tensor(nelem)
|
87 |
|
88 | p2.tensor.ndim = ndim - 2
|
89 |
|
90 | j = 0
|
91 | for i in [0...ndim]
|
92 | if (i != l && i != m)
|
93 | p2.tensor.dim[j++] = p1.tensor.dim[i]
|
94 |
|
95 | a = p1.tensor.elem
|
96 | b = p2.tensor.elem
|
97 |
|
98 | for i in [0...ndim]
|
99 | ai[i] = 0
|
100 | an[i] = p1.tensor.dim[i]
|
101 |
|
102 | for i in [0...nelem]
|
103 | push(zero)
|
104 | for j in [0...n]
|
105 | ai[l] = j
|
106 | ai[m] = j
|
107 | h = 0
|
108 | for k in [0...ndim]
|
109 | h = (h * an[k]) + ai[k]
|
110 | push(a[h])
|
111 | add()
|
112 | b[i] = pop()
|
113 | for j in [(ndim - 1)..0]
|
114 | if (j == l || j == m)
|
115 | continue
|
116 | if (++ai[j] < an[j])
|
117 | break
|
118 | ai[j] = 0
|
119 |
|
120 | if (nelem == 1)
|
121 | push(b[0])
|
122 | else
|
123 | push(p2)
|
124 |
|