1 /*
2   Teem: Tools to process and visualize scientific data and images             .
3   Copyright (C) 2012, 2011, 2010, 2009  University of Chicago
4   Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
5   Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah
6 
7   This library is free software; you can redistribute it and/or
8   modify it under the terms of the GNU Lesser General Public License
9   (LGPL) as published by the Free Software Foundation; either
10   version 2.1 of the License, or (at your option) any later version.
11   The terms of redistributing and/or modifying this software also
12   include exceptions to the LGPL that facilitate static linking.
13 
14   This library is distributed in the hope that it will be useful,
15   but WITHOUT ANY WARRANTY; without even the implied warranty of
16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17   Lesser General Public License for more details.
18 
19   You should have received a copy of the GNU Lesser General Public License
20   along with this library; if not, write to Free Software Foundation, Inc.,
21   51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22 */
23 
24 #ifndef ELLMACROS_HAS_BEEN_INCLUDED
25 #define ELLMACROS_HAS_BEEN_INCLUDED
26 
27 #ifdef __cplusplus
28 extern "C" {
29 #endif
30 
31 /*
32 ******** ELL_SWAP2, ELL_SWAP3
33 **
34 ** used to interchange 2 or 3 values, using the given temp variable
35 ** SWAP2: a <=> b
36 ** SWAP3: a <= b, b <= c, c <= a
37 */
38 #define ELL_SWAP2(a, b, t)    ((t)=(a),(a)=(b),(b)=(t))
39 #define ELL_SWAP3(a, b, c, t) ((t)=(a),(a)=(b),(b)=(c),(c)=(t))
40 
41 /*
42 ******** ELL_SORT3
43 **
44 ** sorts v0, v1, v2 in descending order, using given temp variable t,
45 */
46 #define ELL_SORT3(v0, v1, v2, t)             \
47   if (v0 > v1) {                             \
48     if (v1 < v2) {                           \
49       if (v0 > v2) { ELL_SWAP2(v1, v2, t); } \
50       else { ELL_SWAP3(v0, v2, v1, t); }     \
51     }                                        \
52   }                                          \
53   else {                                     \
54     if (v1 > v2) {                           \
55       if (v0 > v2) { ELL_SWAP2(v0, v1, t); } \
56       else { ELL_SWAP3(v0, v1, v2, t); }     \
57     }                                        \
58     else {                                   \
59       ELL_SWAP2(v0, v2, t);                  \
60     }                                        \
61   }
62 
63 /*
64 ******** ELL_MAX3_IDX
65 **
66 ** returns 0, 1, 2, to indicate which of the three arguments is largest
67 */
68 #define ELL_MAX3_IDX(v0, v1, v2) \
69   (v0 > v1                       \
70    ? (v1 > v2                    \
71       ? 0                        \
72       : (v0 > v2                 \
73          ? 0                     \
74          : 2))                   \
75    : (v2 > v1                    \
76       ? 2                        \
77       : 1))
78 
79 /*
80 ******** ELL_MIN3_IDX
81 **
82 ** returns 0, 1, 2, to indicate which of the three arguments is smallest
83 */
84 #define ELL_MIN3_IDX(v0, v1, v2) \
85   (v0 < v1                       \
86    ? (v1 < v2                    \
87       ? 0                        \
88       : (v0 < v2                 \
89          ? 0                     \
90          : 2))                   \
91    : (v2 < v1                    \
92       ? 2                        \
93       : 1))
94 
95 #define ELL_2V_SET(v, a, b) \
96   ((v)[0]=(a), (v)[1]=(b))
97 
98 #define ELL_2V_SET_TT(v, TT, a, b) \
99   ((v)[0] = AIR_CAST(TT, (a)), \
100    (v)[1] = AIR_CAST(TT, (b)))
101 
102 #define ELL_2V_COPY(v2, v1) \
103   ((v2)[0] = (v1)[0], (v2)[1] = (v1)[1])
104 
105 #define ELL_2V_INCR(v2, v1) \
106   ((v2)[0] += (v1)[0],      \
107    (v2)[1] += (v1)[1])
108 
109 #define ELL_2V_LERP(v3, w, v1, v2)            \
110   ((v3)[0] = AIR_LERP((w), (v1)[0], (v2)[0]), \
111    (v3)[1] = AIR_LERP((w), (v1)[1], (v2)[1]))
112 
113 #define ELL_2V_LERP_TT(v3, TT, w, v1, v2)                      \
114   ((v3)[0] = AIR_CAST(TT, AIR_LERP((w), (v1)[0], (v2)[0])),    \
115    (v3)[1] = AIR_CAST(TT, AIR_LERP((w), (v1)[1], (v2)[1])))
116 
117 #define ELL_2V_ADD2(v3, v1, v2) \
118   ((v3)[0] = (v1)[0] + (v2)[0], \
119    (v3)[1] = (v1)[1] + (v2)[1])
120 
121 #define ELL_2V_ADD3(v4, v1, v2, v3)       \
122   ((v4)[0] = (v1)[0] + (v2)[0] + (v3)[0], \
123    (v4)[1] = (v1)[1] + (v2)[1] + (v3)[1])
124 
125 #define ELL_2V_ADD4(v5, v1, v2, v3, v4)       \
126   ((v5)[0] = (v1)[0] + (v2)[0] + (v3)[0] + (v4)[0], \
127    (v5)[1] = (v1)[1] + (v2)[1] + (v3)[1] + (v4)[1])
128 
129 #define ELL_2V_SUB(v3, v1, v2)  \
130   ((v3)[0] = (v1)[0] - (v2)[0], \
131    (v3)[1] = (v1)[1] - (v2)[1])
132 
133 #define ELL_2V_COPY(v2, v1) \
134   ((v2)[0] = (v1)[0], (v2)[1] = (v1)[1])
135 
136 #define ELL_2V_COPY_TT(v2, TYPE, v1)  \
137   ((v2)[0] = AIR_CAST(TYPE, (v1)[0]), \
138    (v2)[1] = AIR_CAST(TYPE, (v1)[1]))
139 
140 #define ELL_2V_DOT(v1, v2) ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1])
141 
142 #define ELL_2V_LEN(v) (sqrt(ELL_2V_DOT((v),(v))))
143 
144 #define ELL_2V_SCALE(v2, a, v1) \
145   ((v2)[0] = (a)*(v1)[0],       \
146    (v2)[1] = (a)*(v1)[1])
147 
148 #define ELL_2V_SCALE_ADD2(v2, s0, v0, s1, v1) \
149   ((v2)[0] = (s0)*(v0)[0] + (s1)*(v1)[0],     \
150    (v2)[1] = (s0)*(v0)[1] + (s1)*(v1)[1])
151 
152 #define ELL_2V_NORM(v2, v1, length) \
153   (length = ELL_2V_LEN(v1), ELL_2V_SCALE(v2, 1.0/length, v1))
154 
155 /*
156 ** the 2x2 matrix-related macros assume that the matrix indexing is:
157 ** 0  1
158 ** 2  3
159 */
160 
161 #define _ELL_2M_DET(a,b,c,d) ((a)*(d) - (b)*(c))
162 
163 #define ELL_2M_DET(m) _ELL_2M_DET((m)[0],(m)[1],(m)[2],(m)[3])
164 
165 #define ELL_2M_TRANSPOSE(m2, m1)                \
166   ((m2)[0] = (m1)[0],                           \
167    (m2)[1] = (m1)[2],                           \
168    (m2)[2] = (m1)[1],                           \
169    (m2)[3] = (m1)[3])
170 
171 #define ELL_2M_MUL(m3, m1, m2)                  \
172   ((m3)[0] = (m1)[0]*(m2)[0] + (m1)[1]*(m2)[2], \
173    (m3)[1] = (m1)[0]*(m2)[1] + (m1)[1]*(m2)[3], \
174                                                 \
175    (m3)[2] = (m1)[2]*(m2)[0] + (m1)[3]*(m2)[2], \
176    (m3)[3] = (m1)[2]*(m2)[1] + (m1)[3]*(m2)[3])
177 
178 #define ELL_2M_ROTATE_SET(m, th)                \
179   (ELL_2V_SET((m)+ 0,  cos(th) , -sin(th)),     \
180    ELL_2V_SET((m)+ 2, +sin(th) ,  cos(th)))
181 
182 #define ELL_2MV_MUL(v2, m, v1)                  \
183   ((v2)[0] = (m)[0]*(v1)[0] + (m)[1]*(v1)[1],   \
184    (v2)[1] = (m)[2]*(v1)[0] + (m)[3]*(v1)[1])
185 
186 /*
187 ** the 3x3 matrix-related macros assume that the matrix indexing is:
188 ** 0  1  2
189 ** 3  4  5
190 ** 6  7  8
191 */
192 
193 #define ELL_3V_SET(v, a, b, c) \
194   ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
195 
196 #define ELL_3V_ZERO_SET(v)  ELL_3V_SET(v, 0, 0, 0)
197 
198 #define ELL_3V_SET_TT(v, TT, a, b, c) \
199   ((v)[0] = AIR_CAST(TT, (a)), \
200    (v)[1] = AIR_CAST(TT, (b)), \
201    (v)[2] = AIR_CAST(TT, (c)))
202 
203 #define ELL_3V_GET(a, b, c, v) \
204   ((a) = (v)[0], (b) = (v)[1], (c) = (v)[2])
205 
206 #define ELL_3V_EQUAL(a, b) \
207   ((a)[0] == (b)[0] && (a)[1] == (b)[1] && (a)[2] == (b)[2])
208 
209 #define ELL_3V_COPY(v2, v1) \
210   ((v2)[0] = (v1)[0], (v2)[1] = (v1)[1], (v2)[2] = (v1)[2])
211 
212 #define ELL_3V_COPY_TT(v2, TYPE, v1) \
213   ((v2)[0] = AIR_CAST(TYPE, (v1)[0]), \
214    (v2)[1] = AIR_CAST(TYPE, (v1)[1]), \
215    (v2)[2] = AIR_CAST(TYPE, (v1)[2]))
216 
217 #define ELL_3V_INCR(v2, v1) \
218   ((v2)[0] += (v1)[0],      \
219    (v2)[1] += (v1)[1],      \
220    (v2)[2] += (v1)[2])
221 
222 #define ELL_3V_LERP(v3, w, v1, v2)            \
223   ((v3)[0] = AIR_LERP((w), (v1)[0], (v2)[0]), \
224    (v3)[1] = AIR_LERP((w), (v1)[1], (v2)[1]), \
225    (v3)[2] = AIR_LERP((w), (v1)[2], (v2)[2]))
226 
227 #define ELL_3V_LERP_TT(v3, TT, w, v1, v2)                      \
228   ((v3)[0] = AIR_CAST(TT, AIR_LERP((w), (v1)[0], (v2)[0])),    \
229    (v3)[1] = AIR_CAST(TT, AIR_LERP((w), (v1)[1], (v2)[1])),    \
230    (v3)[2] = AIR_CAST(TT, AIR_LERP((w), (v1)[2], (v2)[2])))
231 
232 #define ELL_3V_ADD2(v3, v1, v2) \
233   ((v3)[0] = (v1)[0] + (v2)[0], \
234    (v3)[1] = (v1)[1] + (v2)[1], \
235    (v3)[2] = (v1)[2] + (v2)[2])
236 
237 #define ELL_3V_ADD3(v4, v1, v2, v3)       \
238   ((v4)[0] = (v1)[0] + (v2)[0] + (v3)[0], \
239    (v4)[1] = (v1)[1] + (v2)[1] + (v3)[1], \
240    (v4)[2] = (v1)[2] + (v2)[2] + (v3)[2])
241 
242 #define ELL_3V_ADD4(v5, v1, v2, v3, v4)       \
243   ((v5)[0] = (v1)[0] + (v2)[0] + (v3)[0] + (v4)[0], \
244    (v5)[1] = (v1)[1] + (v2)[1] + (v3)[1] + (v4)[1], \
245    (v5)[2] = (v1)[2] + (v2)[2] + (v3)[2] + (v4)[2])
246 
247 #define ELL_3V_SUB(v3, v1, v2)  \
248   ((v3)[0] = (v1)[0] - (v2)[0], \
249    (v3)[1] = (v1)[1] - (v2)[1], \
250    (v3)[2] = (v1)[2] - (v2)[2])
251 
252 #define ELL_3V_DOT(v1, v2) \
253   ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
254 
255 #define ELL_3V_SCALE(v2, a, v1) \
256   ((v2)[0] = (a)*(v1)[0],       \
257    (v2)[1] = (a)*(v1)[1],       \
258    (v2)[2] = (a)*(v1)[2])
259 
260 #define ELL_3V_SCALE_TT(v2, TT, a, v1)   \
261   ((v2)[0] = AIR_CAST(TT, (a)*(v1)[0]), \
262    (v2)[1] = AIR_CAST(TT, (a)*(v1)[1]), \
263    (v2)[2] = AIR_CAST(TT, (a)*(v1)[2]))
264 
265 #define ELL_3V_SCALE_INCR(v2, s0, v0) \
266   ((v2)[0] += (s0)*(v0)[0], \
267    (v2)[1] += (s0)*(v0)[1], \
268    (v2)[2] += (s0)*(v0)[2])
269 
270 #define ELL_3V_SCALE_INCR_TT(v2, TT, s0, v0) \
271   ((v2)[0] += AIR_CAST(TT, (s0)*(v0)[0]), \
272    (v2)[1] += AIR_CAST(TT, (s0)*(v0)[1]), \
273    (v2)[2] += AIR_CAST(TT, (s0)*(v0)[2]))
274 
275 #define ELL_3V_SCALE_ADD2(v2, s0, v0, s1, v1) \
276   ((v2)[0] = (s0)*(v0)[0] + (s1)*(v1)[0],     \
277    (v2)[1] = (s0)*(v0)[1] + (s1)*(v1)[1],     \
278    (v2)[2] = (s0)*(v0)[2] + (s1)*(v1)[2])
279 
280 #define ELL_3V_SCALE_ADD2_TT(v2, TT, s0, v0, s1, v1) \
281   ((v2)[0] = AIR_CAST(TT, (s0)*(v0)[0] + (s1)*(v1)[0]),     \
282    (v2)[1] = AIR_CAST(TT, (s0)*(v0)[1] + (s1)*(v1)[1]),     \
283    (v2)[2] = AIR_CAST(TT, (s0)*(v0)[2] + (s1)*(v1)[2]))
284 
285 #define ELL_3V_SCALE_INCR2(v2, s0, v0, s1, v1) \
286   ((v2)[0] += (s0)*(v0)[0] + (s1)*(v1)[0],     \
287    (v2)[1] += (s0)*(v0)[1] + (s1)*(v1)[1],     \
288    (v2)[2] += (s0)*(v0)[2] + (s1)*(v1)[2])
289 
290 #define ELL_3V_SCALE_ADD3(v3, s0, v0, s1, v1, s2, v2)     \
291   ((v3)[0] = (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0],  \
292    (v3)[1] = (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1],  \
293    (v3)[2] = (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2])
294 
295 #define ELL_3V_SCALE_ADD3_TT(v3, TT, s0, v0, s1, v1, s2, v2)     \
296   ((v3)[0] = AIR_CAST(TT, (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0]),  \
297    (v3)[1] = AIR_CAST(TT, (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1]),  \
298    (v3)[2] = AIR_CAST(TT, (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2]))
299 
300 #define ELL_3V_SCALE_ADD4(v4, s0, v0, s1, v1, s2, v2, s3, v3)           \
301   ((v4)[0] = (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0] + (s3)*(v3)[0], \
302    (v4)[1] = (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1] + (s3)*(v3)[1], \
303    (v4)[2] = (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2] + (s3)*(v3)[2])
304 
305 #define ELL_3V_SCALE_ADD6(v6, s0, v0, s1, v1, s2, v2,     \
306                               s3, v3, s4, v4, s5, v5)     \
307   ((v6)[0] = (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0]   \
308            + (s3)*(v3)[0] + (s4)*(v4)[0] + (s5)*(v5)[0],  \
309    (v6)[1] = (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1]   \
310            + (s3)*(v3)[1] + (s4)*(v4)[1] + (s5)*(v5)[1],  \
311    (v6)[2] = (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2]   \
312            + (s3)*(v3)[2] + (s4)*(v4)[2] + (s5)*(v5)[2])
313 
314 #define ELL_3V_SCALE_INCR3(v3, s0, v0, s1, v1, s2, v2)     \
315   ((v3)[0] += (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0],  \
316    (v3)[1] += (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1],  \
317    (v3)[2] += (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2])
318 
319 #define ELL_3V_LEN(v) (sqrt(ELL_3V_DOT((v),(v))))
320 
321 #define ELL_3V_DIST(a, b)                    \
322   sqrt(((a)[0] - (b)[0])*((a)[0] - (b)[0]) + \
323        ((a)[1] - (b)[1])*((a)[1] - (b)[1]) + \
324        ((a)[2] - (b)[2])*((a)[2] - (b)[2]))
325 
326 #define ELL_3V_NORM(v2, v1, length) \
327   (length = ELL_3V_LEN(v1), ELL_3V_SCALE(v2, 1.0/length, v1))
328 
329 #define ELL_3V_NORM_TT(v2, TT, v1, length) \
330   (length = AIR_CAST(TT, ELL_3V_LEN(v1)), \
331    ELL_3V_SCALE_TT(v2, TT, 1.0/length, v1))
332 
333 #define ELL_3V_CROSS(v3, v1, v2) \
334   ((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
335    (v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
336    (v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
337 
338 #define ELL_3V_MIN(v3,v1,v2) (         \
339   (v3)[0] = AIR_MIN((v1)[0], (v2)[0]), \
340   (v3)[1] = AIR_MIN((v1)[1], (v2)[1]), \
341   (v3)[2] = AIR_MIN((v1)[2], (v2)[2]))
342 
343 #define ELL_3V_MAX(v3,v1,v2) (         \
344   (v3)[0] = AIR_MAX((v1)[0], (v2)[0]), \
345   (v3)[1] = AIR_MAX((v1)[1], (v2)[1]), \
346   (v3)[2] = AIR_MAX((v1)[2], (v2)[2]))
347 
348 #define ELL_3V_EXISTS(v) \
349    (AIR_EXISTS((v)[0]) && AIR_EXISTS((v)[1]) && AIR_EXISTS((v)[2]))
350 
351 #define ELL_3V_AFFINE(v,i,x,I,o,O) ( \
352   (v)[0] = AIR_AFFINE((i), (x), (I), (o)[0], (O)[0]), \
353   (v)[1] = AIR_AFFINE((i), (x), (I), (o)[1], (O)[1]), \
354   (v)[2] = AIR_AFFINE((i), (x), (I), (o)[2], (O)[2]))
355 
356 #define ELL_3V_ABS(v2,v1) ( \
357   (v2)[0] = AIR_ABS((v1)[0]), \
358   (v2)[1] = AIR_ABS((v1)[1]), \
359   (v2)[2] = AIR_ABS((v1)[2]))
360 
361 #define ELL_3V_NAN_SET(v) ( \
362   (v)[0] = AIR_NAN, \
363   (v)[1] = AIR_NAN, \
364   (v)[2] = AIR_NAN)
365 
366 #define ELL_3M_EQUAL(m1, m2) \
367   ((m1)[0] == (m2)[0] &&     \
368    (m1)[1] == (m2)[1] &&     \
369    (m1)[2] == (m2)[2] &&     \
370    (m1)[3] == (m2)[3] &&     \
371    (m1)[4] == (m2)[4] &&     \
372    (m1)[5] == (m2)[5] &&     \
373    (m1)[6] == (m2)[6] &&     \
374    (m1)[7] == (m2)[7] &&     \
375    (m1)[8] == (m2)[8])
376 
377 #define ELL_3M_SET(m, a, b, c, d, e, f, g, h, i) \
378   (ELL_3V_SET(m + 0*3, a, b, c),                 \
379    ELL_3V_SET(m + 1*3, d, e, f),                 \
380    ELL_3V_SET(m + 2*3, g, h, i))
381 
382 #define ELL_3M_SCALE(m2, s, m1) \
383   (ELL_3V_SCALE((m2)+0, (s), (m1)+0), \
384    ELL_3V_SCALE((m2)+3, (s), (m1)+3), \
385    ELL_3V_SCALE((m2)+6, (s), (m1)+6))
386 
387 #define ELL_3M_SCALE_INCR(m2, s, m1) \
388   (ELL_3V_SCALE_INCR((m2)+0, (s), (m1)+0), \
389    ELL_3V_SCALE_INCR((m2)+3, (s), (m1)+3), \
390    ELL_3V_SCALE_INCR((m2)+6, (s), (m1)+6))
391 
392 #define ELL_3M_SCALE_ADD2(m2, s0, m0, s1, m1) \
393   (ELL_3V_SCALE_ADD2((m2)+0, (s0), (m0)+0, (s1), (m1)+0), \
394    ELL_3V_SCALE_ADD2((m2)+3, (s0), (m0)+3, (s1), (m1)+3), \
395    ELL_3V_SCALE_ADD2((m2)+6, (s0), (m0)+6, (s1), (m1)+6))
396 
397 #define ELL_3M_LERP(m3, w, m1, m2)           \
398   (ELL_3V_LERP((m3)+0, (w), (m1)+0, (m2)+0), \
399    ELL_3V_LERP((m3)+3, (w), (m1)+3, (m2)+3), \
400    ELL_3V_LERP((m3)+6, (w), (m1)+6, (m2)+6))
401 
402 #define ELL_3M_ADD2(m3, m1, m2) \
403   ((m3)[0] = (m1)[0] + (m2)[0],  \
404    (m3)[1] = (m1)[1] + (m2)[1],  \
405    (m3)[2] = (m1)[2] + (m2)[2],  \
406    (m3)[3] = (m1)[3] + (m2)[3],  \
407    (m3)[4] = (m1)[4] + (m2)[4],  \
408    (m3)[5] = (m1)[5] + (m2)[5],  \
409    (m3)[6] = (m1)[6] + (m2)[6],  \
410    (m3)[7] = (m1)[7] + (m2)[7],  \
411    (m3)[8] = (m1)[8] + (m2)[8])
412 
413 #define ELL_3M_SUB(m3, m1, m2) \
414   ((m3)[0] = (m1)[0] - (m2)[0],  \
415    (m3)[1] = (m1)[1] - (m2)[1],  \
416    (m3)[2] = (m1)[2] - (m2)[2],  \
417    (m3)[3] = (m1)[3] - (m2)[3],  \
418    (m3)[4] = (m1)[4] - (m2)[4],  \
419    (m3)[5] = (m1)[5] - (m2)[5],  \
420    (m3)[6] = (m1)[6] - (m2)[6],  \
421    (m3)[7] = (m1)[7] - (m2)[7],  \
422    (m3)[8] = (m1)[8] - (m2)[8])
423 
424 #define ELL_3M_SCALE_ADD3(m3, s0, m0, s1, m1, s2, m2) \
425   (ELL_3V_SCALE_ADD3((m3)+0, (s0), (m0)+0, (s1), (m1)+0, (s2), (m2)+0), \
426    ELL_3V_SCALE_ADD3((m3)+3, (s0), (m0)+3, (s1), (m1)+3, (s2), (m2)+3), \
427    ELL_3V_SCALE_ADD3((m3)+6, (s0), (m0)+6, (s1), (m1)+6, (s2), (m2)+6))
428 
429 #define ELL_3M_COPY(m2, m1) \
430   (ELL_3V_COPY((m2)+0, (m1)+0), \
431    ELL_3V_COPY((m2)+3, (m1)+3), \
432    ELL_3V_COPY((m2)+6, (m1)+6))
433 
434 #define ELL_3M_COPY_TT(m2, TYPE, m1) \
435   (ELL_3V_COPY_TT((m2)+0, TYPE, (m1)+0), \
436    ELL_3V_COPY_TT((m2)+3, TYPE, (m1)+3), \
437    ELL_3V_COPY_TT((m2)+6, TYPE, (m1)+6))
438 
439 #define ELL_3M_IDENTITY_SET(m) \
440   (ELL_3V_SET((m)+0,  1 ,  0 ,  0), \
441    ELL_3V_SET((m)+3,  0 ,  1 ,  0), \
442    ELL_3V_SET((m)+6,  0 ,  0 ,  1))
443 
444 #define ELL_3M_EXISTS(m) \
445   (ELL_3V_EXISTS((m) + 0) \
446    && ELL_3V_EXISTS((m) + 3) \
447    && ELL_3V_EXISTS((m) + 6))
448 
449 #define ELL_3M_ZERO_SET(m) \
450   (ELL_3V_SET((m)+0,  0 ,  0 ,  0), \
451    ELL_3V_SET((m)+3,  0 ,  0 ,  0), \
452    ELL_3V_SET((m)+6,  0 ,  0 ,  0))
453 
454 #define ELL_3M_NAN_SET(m) \
455   (ELL_3V_NAN_SET((m)+0), \
456    ELL_3V_NAN_SET((m)+3), \
457    ELL_3V_NAN_SET((m)+6))
458 
459 #define ELL_3M_DIAG_SET(m, a, b, c) \
460   ((m)[0] = (a), (m)[4] = (b), (m)[8] = (c))
461 
462 #define ELL_3M_TRANSPOSE(m2, m1) \
463   ((m2)[0] = (m1)[0],            \
464    (m2)[1] = (m1)[3],            \
465    (m2)[2] = (m1)[6],            \
466    (m2)[3] = (m1)[1],            \
467    (m2)[4] = (m1)[4],            \
468    (m2)[5] = (m1)[7],            \
469    (m2)[6] = (m1)[2],            \
470    (m2)[7] = (m1)[5],            \
471    (m2)[8] = (m1)[8])
472 
473 #define ELL_3M_TRANSPOSE_IP(m, t) \
474   (ELL_SWAP2((m)[1],(m)[3],(t)),  \
475    ELL_SWAP2((m)[2],(m)[6],(t)),  \
476    ELL_SWAP2((m)[5],(m)[7],(t)))
477 
478 #define ELL_3M_TRACE(m) ((m)[0] + (m)[4] + (m)[8])
479 
480 #define ELL_3M_FROB(m) \
481   (sqrt(ELL_3V_DOT((m)+0, (m)+0) + \
482         ELL_3V_DOT((m)+3, (m)+3) + \
483         ELL_3V_DOT((m)+6, (m)+6)))
484 
485 #define _ELL_3M_DET(a,b,c,d,e,f,g,h,i) \
486   (  (a)*(e)*(i) \
487    + (d)*(h)*(c) \
488    + (g)*(b)*(f) \
489    - (g)*(e)*(c) \
490    - (d)*(b)*(i) \
491    - (a)*(h)*(f))
492 
493 #define ELL_3M_DET(m) _ELL_3M_DET((m)[0],(m)[1],(m)[2],\
494                                   (m)[3],(m)[4],(m)[5],\
495                                   (m)[6],(m)[7],(m)[8])
496 
497 #define ELL_3MV_COL0_GET(v, m) \
498   (ELL_3V_SET((v), (m)[0], (m)[3], (m)[6]))
499 
500 #define ELL_3MV_COL1_GET(v, m) \
501   (ELL_3V_SET((v), (m)[1], (m)[4], (m)[7]))
502 
503 #define ELL_3MV_COL2_GET(v, m) \
504   (ELL_3V_SET((v), (m)[2], (m)[5], (m)[8]))
505 
506 #define ELL_3MV_ROW0_GET(v, m) \
507   (ELL_3V_SET((v), (m)[0], (m)[1], (m)[2]))
508 
509 #define ELL_3MV_ROW1_GET(v, m) \
510   (ELL_3V_SET((v), (m)[3], (m)[4], (m)[5]))
511 
512 #define ELL_3MV_ROW2_GET(v, m) \
513   (ELL_3V_SET((v), (m)[6], (m)[7], (m)[8]))
514 
515 #define ELL_3MV_COL0_SET(m, v) \
516   (ELL_3V_GET((m)[0], (m)[3], (m)[6], (v)))
517 
518 #define ELL_3MV_COL1_SET(m, v) \
519   (ELL_3V_GET((m)[1], (m)[4], (m)[7], (v)))
520 
521 #define ELL_3MV_COL2_SET(m, v) \
522   (ELL_3V_GET((m)[2], (m)[5], (m)[8], (v)))
523 
524 #define ELL_3MV_ROW0_SET(m, v) \
525   (ELL_3V_GET((m)[0], (m)[1], (m)[2], (v)))
526 
527 #define ELL_3MV_ROW1_SET(m, v) \
528   (ELL_3V_GET((m)[3], (m)[4], (m)[5], (v)))
529 
530 #define ELL_3MV_ROW2_SET(m, v) \
531   (ELL_3V_GET((m)[6], (m)[7], (m)[8], (v)))
532 
533 #define ELL_3MV_OUTER(m, v1, v2) \
534   (ELL_3V_SCALE((m)+0, (v1)[0], (v2)), \
535    ELL_3V_SCALE((m)+3, (v1)[1], (v2)), \
536    ELL_3V_SCALE((m)+6, (v1)[2], (v2)))
537 
538 #define ELL_3MV_OUTER_TT(m, T, v1, v2)  \
539   (ELL_3V_SCALE_TT((m)+0, T, (v1)[0], (v2)),        \
540    ELL_3V_SCALE_TT((m)+3, T, (v1)[1], (v2)),         \
541    ELL_3V_SCALE_TT((m)+6, T, (v1)[2], (v2)))
542 
543 #define ELL_3MV_OUTER_INCR(m, v1, v2) \
544   (ELL_3V_SCALE_INCR((m)+0, (v1)[0], (v2)), \
545    ELL_3V_SCALE_INCR((m)+3, (v1)[1], (v2)), \
546    ELL_3V_SCALE_INCR((m)+6, (v1)[2], (v2)))
547 
548 #define ELL_3MV_SCALE_OUTER_INCR(m, s, v1, v2) \
549   (ELL_3V_SCALE_INCR((m)+0, (s)*(v1)[0], (v2)), \
550    ELL_3V_SCALE_INCR((m)+3, (s)*(v1)[1], (v2)), \
551    ELL_3V_SCALE_INCR((m)+6, (s)*(v1)[2], (v2)))
552 
553 #define ELL_3MV_MUL(v2, m, v1) \
554   ((v2)[0] = (m)[0]*(v1)[0] + (m)[1]*(v1)[1] + (m)[2]*(v1)[2], \
555    (v2)[1] = (m)[3]*(v1)[0] + (m)[4]*(v1)[1] + (m)[5]*(v1)[2], \
556    (v2)[2] = (m)[6]*(v1)[0] + (m)[7]*(v1)[1] + (m)[8]*(v1)[2])
557 
558 #define ELL_3MV_CONTR(m, v) \
559   ((m)[0]*(v)[0]*(v)[0] + (m)[1]*(v)[1]*(v)[0] + (m)[2]*(v)[2]*(v)[0] + \
560    (m)[3]*(v)[0]*(v)[1] + (m)[4]*(v)[1]*(v)[1] + (m)[5]*(v)[2]*(v)[1] + \
561    (m)[6]*(v)[0]*(v)[2] + (m)[7]*(v)[1]*(v)[2] + (m)[8]*(v)[2]*(v)[2])
562 
563 #define ELL_3MV_CONTR2(u, m, v)                                         \
564   ((m)[0]*(v)[0]*(u)[0] + (m)[1]*(v)[1]*(u)[0] + (m)[2]*(v)[2]*(u)[0] + \
565    (m)[3]*(v)[0]*(u)[1] + (m)[4]*(v)[1]*(u)[1] + (m)[5]*(v)[2]*(u)[1] + \
566    (m)[6]*(v)[0]*(u)[2] + (m)[7]*(v)[1]*(u)[2] + (m)[8]*(v)[2]*(u)[2])
567 
568 #define ELL_3MV_MUL_TT(v2, TT, m, v1) \
569   ((v2)[0] = AIR_CAST(TT, (m)[0]*(v1)[0] + (m)[1]*(v1)[1] + (m)[2]*(v1)[2]), \
570    (v2)[1] = AIR_CAST(TT, (m)[3]*(v1)[0] + (m)[4]*(v1)[1] + (m)[5]*(v1)[2]), \
571    (v2)[2] = AIR_CAST(TT, (m)[6]*(v1)[0] + (m)[7]*(v1)[1] + (m)[8]*(v1)[2]))
572 
573 #define ELL_3MV_TMUL(v2, m, v1) \
574   ((v2)[0] = (m)[0]*(v1)[0] + (m)[3]*(v1)[1] + (m)[6]*(v1)[2], \
575    (v2)[1] = (m)[1]*(v1)[0] + (m)[4]*(v1)[1] + (m)[7]*(v1)[2], \
576    (v2)[2] = (m)[2]*(v1)[0] + (m)[5]*(v1)[1] + (m)[8]*(v1)[2])
577 
578 #define ELL_3M_MUL(m3, m1, m2)                                    \
579   ((m3)[0] = (m1)[0]*(m2)[0] + (m1)[1]*(m2)[3] + (m1)[2]*(m2)[6], \
580    (m3)[1] = (m1)[0]*(m2)[1] + (m1)[1]*(m2)[4] + (m1)[2]*(m2)[7], \
581    (m3)[2] = (m1)[0]*(m2)[2] + (m1)[1]*(m2)[5] + (m1)[2]*(m2)[8], \
582                                                                   \
583    (m3)[3] = (m1)[3]*(m2)[0] + (m1)[4]*(m2)[3] + (m1)[5]*(m2)[6], \
584    (m3)[4] = (m1)[3]*(m2)[1] + (m1)[4]*(m2)[4] + (m1)[5]*(m2)[7], \
585    (m3)[5] = (m1)[3]*(m2)[2] + (m1)[4]*(m2)[5] + (m1)[5]*(m2)[8], \
586                                                                   \
587    (m3)[6] = (m1)[6]*(m2)[0] + (m1)[7]*(m2)[3] + (m1)[8]*(m2)[6], \
588    (m3)[7] = (m1)[6]*(m2)[1] + (m1)[7]*(m2)[4] + (m1)[8]*(m2)[7], \
589    (m3)[8] = (m1)[6]*(m2)[2] + (m1)[7]*(m2)[5] + (m1)[8]*(m2)[8])
590 
591 #define ELL_3M_MUL_TT(m3, TT, m1, m2)                                    \
592   ((m3)[0] = AIR_CAST(TT, (m1)[0]*(m2)[0]+(m1)[1]*(m2)[3]+(m1)[2]*(m2)[6]), \
593    (m3)[1] = AIR_CAST(TT, (m1)[0]*(m2)[1]+(m1)[1]*(m2)[4]+(m1)[2]*(m2)[7]), \
594    (m3)[2] = AIR_CAST(TT, (m1)[0]*(m2)[2]+(m1)[1]*(m2)[5]+(m1)[2]*(m2)[8]), \
595                                                                   \
596    (m3)[3] = AIR_CAST(TT, (m1)[3]*(m2)[0]+(m1)[4]*(m2)[3]+(m1)[5]*(m2)[6]), \
597    (m3)[4] = AIR_CAST(TT, (m1)[3]*(m2)[1]+(m1)[4]*(m2)[4]+(m1)[5]*(m2)[7]), \
598    (m3)[5] = AIR_CAST(TT, (m1)[3]*(m2)[2]+(m1)[4]*(m2)[5]+(m1)[5]*(m2)[8]), \
599                                                                   \
600    (m3)[6] = AIR_CAST(TT, (m1)[6]*(m2)[0]+(m1)[7]*(m2)[3]+(m1)[8]*(m2)[6]), \
601    (m3)[7] = AIR_CAST(TT, (m1)[6]*(m2)[1]+(m1)[7]*(m2)[4]+(m1)[8]*(m2)[7]), \
602    (m3)[8] = AIR_CAST(TT, (m1)[6]*(m2)[2]+(m1)[7]*(m2)[5]+(m1)[8]*(m2)[8]))
603 
604 #define ELL_3M_INV(m2, m1, det)                                   \
605   ((det) = ELL_3M_DET(m1),                                        \
606    (m2)[0] =  _ELL_2M_DET((m1)[4],(m1)[5],(m1)[7],(m1)[8])/(det), \
607    (m2)[1] = -_ELL_2M_DET((m1)[1],(m1)[2],(m1)[7],(m1)[8])/(det), \
608    (m2)[2] =  _ELL_2M_DET((m1)[1],(m1)[2],(m1)[4],(m1)[5])/(det), \
609    (m2)[3] = -_ELL_2M_DET((m1)[3],(m1)[5],(m1)[6],(m1)[8])/(det), \
610    (m2)[4] =  _ELL_2M_DET((m1)[0],(m1)[2],(m1)[6],(m1)[8])/(det), \
611    (m2)[5] = -_ELL_2M_DET((m1)[0],(m1)[2],(m1)[3],(m1)[5])/(det), \
612    (m2)[6] =  _ELL_2M_DET((m1)[3],(m1)[4],(m1)[6],(m1)[7])/(det), \
613    (m2)[7] = -_ELL_2M_DET((m1)[0],(m1)[1],(m1)[6],(m1)[7])/(det), \
614    (m2)[8] =  _ELL_2M_DET((m1)[0],(m1)[1],(m1)[3],(m1)[4])/(det))
615 
616 #define ELL_3M_SCALE_SET(m, x, y, z)  \
617   (ELL_3V_SET((m)+ 0, (x),  0 ,  0 ), \
618    ELL_3V_SET((m)+ 3,  0 , (y),  0 ), \
619    ELL_3V_SET((m)+ 6,  0 ,  0 , (z)))
620 
621 #define ELL_3M_ROTATE_X_SET(m, th)               \
622   (ELL_3V_SET((m)+ 0,  1 ,     0    ,     0   ), \
623    ELL_3V_SET((m)+ 3,  0 ,  cos(th) , -sin(th)), \
624    ELL_3V_SET((m)+ 6,  0 , +sin(th) ,  cos(th)))
625 
626 #define ELL_3M_ROTATE_Y_SET(m, th)               \
627   (ELL_3V_SET((m)+ 0,  cos(th) ,  0 , +sin(th)), \
628    ELL_3V_SET((m)+ 3,     0    ,  1 ,     0   ), \
629    ELL_3V_SET((m)+ 6, -sin(th) ,  0 ,  cos(th)))
630 
631 #define ELL_3M_ROTATE_Z_SET(m, th)               \
632   (ELL_3V_SET((m)+ 0,  cos(th) , -sin(th) ,  0), \
633    ELL_3V_SET((m)+ 3, +sin(th) ,  cos(th) ,  0), \
634    ELL_3V_SET((m)+ 6,     0    ,     0    ,  1))
635 
636 /*
637 ** the 4x4 matrix-related macros assume that the matrix indexing is:
638 **
639 **  0   1   2   3
640 **  4   5   6   7
641 **  8   9  10  11
642 ** 12  13  14  15
643 */
644 
645 #define ELL_4V_SET(v, a, b, c, d) \
646   ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c), (v)[3] = (d))
647 
648 #define ELL_4V_ZERO_SET(v)  ELL_4V_SET(v, 0, 0, 0, 0)
649 
650 #define ELL_4V_NAN_SET(v) ( \
651   (v)[0] = AIR_NAN, \
652   (v)[1] = AIR_NAN, \
653   (v)[2] = AIR_NAN, \
654   (v)[3] = AIR_NAN)
655 
656 #define ELL_4V_SET_TT(v, TT, a, b, c, d) \
657   ((v)[0] = AIR_CAST(TT, (a)), \
658    (v)[1] = AIR_CAST(TT, (b)), \
659    (v)[2] = AIR_CAST(TT, (c)), \
660    (v)[3] = AIR_CAST(TT, (d)))
661 
662 #define ELL_4V_GET(a, b, c, d, v) \
663   ((a) = (v)[0], (b) = (v)[1], (c) = (v)[2], (d) = (v)[3])
664 
665 #define ELL_4V_EQUAL(a, b) \
666   ((a)[0]==(b)[0] && (a)[1]==(b)[1] && (a)[2]==(b)[2] && (a)[3]==(b)[3])
667 
668 #define ELL_4V_COPY(v2, v1) \
669   ((v2)[0] = (v1)[0],       \
670    (v2)[1] = (v1)[1],       \
671    (v2)[2] = (v1)[2],       \
672    (v2)[3] = (v1)[3])
673 
674 #define ELL_4V_COPY_TT(v2, TT, v1)  \
675   ((v2)[0] = AIR_CAST(TT, (v1)[0]), \
676    (v2)[1] = AIR_CAST(TT, (v1)[1]), \
677    (v2)[2] = AIR_CAST(TT, (v1)[2]), \
678    (v2)[3] = AIR_CAST(TT, (v1)[3]))
679 
680 #define ELL_4V_INCR(v2, v1) \
681   ((v2)[0] += (v1)[0],      \
682    (v2)[1] += (v1)[1],      \
683    (v2)[2] += (v1)[2],      \
684    (v2)[3] += (v1)[3])
685 
686 #define ELL_4V_ADD2(v3, v1, v2)  \
687   ((v3)[0] = (v1)[0] + (v2)[0], \
688    (v3)[1] = (v1)[1] + (v2)[1], \
689    (v3)[2] = (v1)[2] + (v2)[2], \
690    (v3)[3] = (v1)[3] + (v2)[3])
691 
692 #define ELL_4V_SUB(v3, v1, v2)  \
693   ((v3)[0] = (v1)[0] - (v2)[0], \
694    (v3)[1] = (v1)[1] - (v2)[1], \
695    (v3)[2] = (v1)[2] - (v2)[2], \
696    (v3)[3] = (v1)[3] - (v2)[3])
697 
698 #define ELL_4V_DOT(v1, v2) \
699   ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2] + (v1)[3]*(v2)[3])
700 
701 #define ELL_4V_SCALE(v2, a, v1) \
702   ((v2)[0] = (v1)[0]*a, (v2)[1] = (v1)[1]*a, \
703    (v2)[2] = (v1)[2]*a, (v2)[3] = (v1)[3]*a)
704 
705 #define ELL_4V_SCALE_TT(v2, TT, a, v1)   \
706   ((v2)[0] = AIR_CAST(TT, (v1)[0]*(a)),  \
707    (v2)[1] = AIR_CAST(TT, (v1)[1]*(a)),  \
708    (v2)[2] = AIR_CAST(TT, (v1)[2]*(a)),  \
709    (v2)[3] = AIR_CAST(TT, (v1)[3]*(a)))
710 
711 #define ELL_4V_SCALE_ADD2(v2, s0, v0, s1, v1) \
712   ((v2)[0] = (s0)*(v0)[0] + (s1)*(v1)[0],    \
713    (v2)[1] = (s0)*(v0)[1] + (s1)*(v1)[1],    \
714    (v2)[2] = (s0)*(v0)[2] + (s1)*(v1)[2],    \
715    (v2)[3] = (s0)*(v0)[3] + (s1)*(v1)[3])
716 
717 #define ELL_4V_SCALE_ADD3(v, s0, v0, s1, v1, s2, v2)    \
718   ((v)[0] = (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0], \
719    (v)[1] = (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1], \
720    (v)[2] = (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2], \
721    (v)[3] = (s0)*(v0)[3] + (s1)*(v1)[3] + (s2)*(v2)[3])
722 
723 #define ELL_4V_SCALE_ADD4(v, s0, v0, s1, v1, s2, v2, s3, v3)           \
724   ((v)[0] = (s0)*(v0)[0] + (s1)*(v1)[0] + (s2)*(v2)[0] + (s3)*(v3)[0], \
725    (v)[1] = (s0)*(v0)[1] + (s1)*(v1)[1] + (s2)*(v2)[1] + (s3)*(v3)[1], \
726    (v)[2] = (s0)*(v0)[2] + (s1)*(v1)[2] + (s2)*(v2)[2] + (s3)*(v3)[2], \
727    (v)[3] = (s0)*(v0)[3] + (s1)*(v1)[3] + (s2)*(v2)[3] + (s3)*(v3)[3])
728 
729 #define ELL_4V_SCALE_INCR(v2, s0, v0) \
730   ((v2)[0] += (s0)*(v0)[0], \
731    (v2)[1] += (s0)*(v0)[1], \
732    (v2)[2] += (s0)*(v0)[2], \
733    (v2)[3] += (s0)*(v0)[3])
734 
735 #define ELL_4V_LEN(v) (sqrt(ELL_4V_DOT((v),(v))))
736 
737 #define ELL_4V_NORM(v2, v1, length) \
738   (length = ELL_4V_LEN(v1), ELL_4V_SCALE(v2, 1.0/length, v1))
739 
740 #define ELL_4V_LERP(v3, w, v1, v2)            \
741   ((v3)[0] = AIR_LERP((w), (v1)[0], (v2)[0]), \
742    (v3)[1] = AIR_LERP((w), (v1)[1], (v2)[1]), \
743    (v3)[2] = AIR_LERP((w), (v1)[2], (v2)[2]), \
744    (v3)[3] = AIR_LERP((w), (v1)[3], (v2)[3]))
745 
746 #define ELL_4V_LERP_TT(v3, TT, w, v1, v2)                      \
747   ((v3)[0] = AIR_CAST(TT, AIR_LERP((w), (v1)[0], (v2)[0])),    \
748    (v3)[1] = AIR_CAST(TT, AIR_LERP((w), (v1)[1], (v2)[1])),    \
749    (v3)[2] = AIR_CAST(TT, AIR_LERP((w), (v1)[2], (v2)[2])),    \
750    (v3)[3] = AIR_CAST(TT, AIR_LERP((w), (v1)[3], (v2)[3])))
751 
752 #define ELL_4V_EXISTS(v) \
753    (AIR_EXISTS((v)[0]) && AIR_EXISTS((v)[1]) \
754     && AIR_EXISTS((v)[2]) && AIR_EXISTS((v)[3]))
755 
756 #define ELL_4M_EQUAL(m1, m2) \
757   ((m1)[ 0] == (m2)[ 0] &&   \
758    (m1)[ 1] == (m2)[ 1] &&   \
759    (m1)[ 2] == (m2)[ 2] &&   \
760    (m1)[ 3] == (m2)[ 3] &&   \
761    (m1)[ 4] == (m2)[ 4] &&   \
762    (m1)[ 5] == (m2)[ 5] &&   \
763    (m1)[ 6] == (m2)[ 6] &&   \
764    (m1)[ 7] == (m2)[ 7] &&   \
765    (m1)[ 8] == (m2)[ 8] &&   \
766    (m1)[ 9] == (m2)[ 9] &&   \
767    (m1)[10] == (m2)[10] &&   \
768    (m1)[11] == (m2)[11] &&   \
769    (m1)[12] == (m2)[12] &&   \
770    (m1)[13] == (m2)[13] &&   \
771    (m1)[14] == (m2)[14] &&   \
772    (m1)[15] == (m2)[15])
773 
774 #define ELL_4M_ADD2(m3, m1, m2)            \
775   (ELL_4V_ADD2((m3)+ 0, (m1)+ 0, (m2)+ 0), \
776    ELL_4V_ADD2((m3)+ 4, (m1)+ 4, (m2)+ 4), \
777    ELL_4V_ADD2((m3)+ 8, (m1)+ 8, (m2)+ 8), \
778    ELL_4V_ADD2((m3)+12, (m1)+12, (m2)+12))
779 
780 #define ELL_4M_SUB(m3, m1, m2)            \
781   (ELL_4V_SUB((m3)+ 0, (m1)+ 0, (m2)+ 0), \
782    ELL_4V_SUB((m3)+ 4, (m1)+ 4, (m2)+ 4), \
783    ELL_4V_SUB((m3)+ 8, (m1)+ 8, (m2)+ 8), \
784    ELL_4V_SUB((m3)+12, (m1)+12, (m2)+12))
785 
786 #define ELL_4M_SCALE(m2, a, m1)         \
787   (ELL_4V_SCALE((m2)+ 0, (a), (m1)+ 0), \
788    ELL_4V_SCALE((m2)+ 4, (a), (m1)+ 4), \
789    ELL_4V_SCALE((m2)+ 8, (a), (m1)+ 8), \
790    ELL_4V_SCALE((m2)+12, (a), (m1)+12))
791 
792 #define ELL_4M_COPY(m2, m1)     \
793   (ELL_4V_COPY((m2)+ 0, (m1)+ 0), \
794    ELL_4V_COPY((m2)+ 4, (m1)+ 4), \
795    ELL_4V_COPY((m2)+ 8, (m1)+ 8), \
796    ELL_4V_COPY((m2)+12, (m1)+12))
797 
798 #define ELL_4M_COPY_TT(m2, TT, m1)     \
799   (ELL_4V_COPY_TT((m2)+ 0, TT, (m1)+ 0), \
800    ELL_4V_COPY_TT((m2)+ 4, TT, (m1)+ 4), \
801    ELL_4V_COPY_TT((m2)+ 8, TT, (m1)+ 8), \
802    ELL_4V_COPY_TT((m2)+12, TT, (m1)+12))
803 
804 #define ELL_4M_TRANSPOSE(m2, m1) \
805   ((m2)[ 0] = (m1)[ 0],          \
806    (m2)[ 1] = (m1)[ 4],          \
807    (m2)[ 2] = (m1)[ 8],          \
808    (m2)[ 3] = (m1)[12],          \
809    (m2)[ 4] = (m1)[ 1],          \
810    (m2)[ 5] = (m1)[ 5],          \
811    (m2)[ 6] = (m1)[ 9],          \
812    (m2)[ 7] = (m1)[13],          \
813    (m2)[ 8] = (m1)[ 2],          \
814    (m2)[ 9] = (m1)[ 6],          \
815    (m2)[10] = (m1)[10],          \
816    (m2)[11] = (m1)[14],          \
817    (m2)[12] = (m1)[ 3],          \
818    (m2)[13] = (m1)[ 7],          \
819    (m2)[14] = (m1)[11],          \
820    (m2)[15] = (m1)[15])
821 
822 #define ELL_4M_TRANSPOSE_TT(m2, TT, m1) \
823   ((m2)[ 0] = AIR_CAST(TT, (m1)[ 0]),   \
824    (m2)[ 1] = AIR_CAST(TT, (m1)[ 4]),   \
825    (m2)[ 2] = AIR_CAST(TT, (m1)[ 8]),   \
826    (m2)[ 3] = AIR_CAST(TT, (m1)[12]),   \
827    (m2)[ 4] = AIR_CAST(TT, (m1)[ 1]),   \
828    (m2)[ 5] = AIR_CAST(TT, (m1)[ 5]),   \
829    (m2)[ 6] = AIR_CAST(TT, (m1)[ 9]),   \
830    (m2)[ 7] = AIR_CAST(TT, (m1)[13]),   \
831    (m2)[ 8] = AIR_CAST(TT, (m1)[ 2]),   \
832    (m2)[ 9] = AIR_CAST(TT, (m1)[ 6]),   \
833    (m2)[10] = AIR_CAST(TT, (m1)[10]),   \
834    (m2)[11] = AIR_CAST(TT, (m1)[14]),   \
835    (m2)[12] = AIR_CAST(TT, (m1)[ 3]),   \
836    (m2)[13] = AIR_CAST(TT, (m1)[ 7]),   \
837    (m2)[14] = AIR_CAST(TT, (m1)[11]),   \
838    (m2)[15] = AIR_CAST(TT, (m1)[15]))
839 
840 #define ELL_4M_TRANSPOSE_IP(m, t)   \
841   (ELL_SWAP2((m)[ 1],(m)[ 4],(t)),  \
842    ELL_SWAP2((m)[ 2],(m)[ 8],(t)),  \
843    ELL_SWAP2((m)[ 3],(m)[12],(t)),  \
844    ELL_SWAP2((m)[ 6],(m)[ 9],(t)),  \
845    ELL_SWAP2((m)[ 7],(m)[13],(t)),  \
846    ELL_SWAP2((m)[11],(m)[14],(t)))
847 
848 #define ELL_4MV_ROW0_GET(v, m) \
849   (ELL_4V_SET((v), (m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3]))
850 
851 #define ELL_4MV_ROW1_GET(v, m) \
852   (ELL_4V_SET((v), (m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7]))
853 
854 #define ELL_4MV_ROW2_GET(v, m) \
855   (ELL_4V_SET((v), (m)[ 8], (m)[ 9], (m)[10], (m)[11]))
856 
857 #define ELL_4MV_ROW3_GET(v, m) \
858   (ELL_4V_SET((v), (m)[12], (m)[13], (m)[14], (m)[15]))
859 
860 #define ELL_4MV_COL0_GET(v, m) \
861   (ELL_4V_SET((v), (m)[ 0], (m)[ 4], (m)[ 8], (m)[12]))
862 
863 #define ELL_4MV_COL1_GET(v, m) \
864   (ELL_4V_SET((v), (m)[ 1], (m)[ 5], (m)[ 9], (m)[13]))
865 
866 #define ELL_4MV_COL2_GET(v, m) \
867   (ELL_4V_SET((v), (m)[ 2], (m)[ 6], (m)[10], (m)[14]))
868 
869 #define ELL_4MV_COL3_GET(v, m) \
870   (ELL_4V_SET((v), (m)[ 3], (m)[ 7], (m)[11], (m)[15]))
871 
872 #define ELL_4MV_ROW0_SET(m, v) \
873   (ELL_4V_GET((m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3], (v)))
874 
875 #define ELL_4MV_ROW1_SET(m, v) \
876   (ELL_4V_GET((m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7], (v)))
877 
878 #define ELL_4MV_ROW2_SET(m, v) \
879   (ELL_4V_GET((m)[ 8], (m)[ 9], (m)[10], (m)[11], (v)))
880 
881 #define ELL_4MV_ROW3_SET(m, v) \
882   (ELL_4V_GET((m)[12], (m)[13], (m)[14], (m)[15], (v)))
883 
884 #define ELL_4MV_COL0_SET(m, v) \
885   (ELL_4V_GET((m)[ 0], (m)[ 4], (m)[ 8], (m)[12], (v)))
886 
887 #define ELL_4MV_COL1_SET(m, v) \
888   (ELL_4V_GET((m)[ 1], (m)[ 5], (m)[ 9], (m)[13], (v)))
889 
890 #define ELL_4MV_COL2_SET(m, v) \
891   (ELL_4V_GET((m)[ 2], (m)[ 6], (m)[10], (m)[14], (v)))
892 
893 #define ELL_4MV_COL3_SET(m, v) \
894   (ELL_4V_GET((m)[ 3], (m)[ 7], (m)[11], (m)[15], (v)))
895 
896 #define ELL_4MV_MUL(v2, m, v1)                                              \
897   ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 1]*(v1)[1]+(m)[ 2]*(v1)[2]+(m)[ 3]*(v1)[3], \
898    (v2)[1]=(m)[ 4]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 6]*(v1)[2]+(m)[ 7]*(v1)[3], \
899    (v2)[2]=(m)[ 8]*(v1)[0]+(m)[ 9]*(v1)[1]+(m)[10]*(v1)[2]+(m)[11]*(v1)[3], \
900    (v2)[3]=(m)[12]*(v1)[0]+(m)[13]*(v1)[1]+(m)[14]*(v1)[2]+(m)[15]*(v1)[3])
901 
902 #define ELL_4MV_MUL_TT(v2, TT, m, v1)                       \
903   ((v2)[0]=AIR_CAST(TT, ((m)[ 0]*(v1)[0]+(m)[ 1]*(v1)[1]    \
904                         +(m)[ 2]*(v1)[2]+(m)[ 3]*(v1)[3])), \
905    (v2)[1]=AIR_CAST(TT, ((m)[ 4]*(v1)[0]+(m)[ 5]*(v1)[1]    \
906                         +(m)[ 6]*(v1)[2]+(m)[ 7]*(v1)[3])), \
907    (v2)[2]=AIR_CAST(TT, ((m)[ 8]*(v1)[0]+(m)[ 9]*(v1)[1]    \
908                         +(m)[10]*(v1)[2]+(m)[11]*(v1)[3])), \
909    (v2)[3]=AIR_CAST(TT, ((m)[12]*(v1)[0]+(m)[13]*(v1)[1]    \
910                         +(m)[14]*(v1)[2]+(m)[15]*(v1)[3])))
911 
912 #define ELL_4MV_TMUL(v2, m, v1)                                             \
913   ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 4]*(v1)[1]+(m)[ 8]*(v1)[2]+(m)[12]*(v1)[3], \
914    (v2)[1]=(m)[ 1]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 9]*(v1)[2]+(m)[13]*(v1)[3], \
915    (v2)[2]=(m)[ 2]*(v1)[0]+(m)[ 6]*(v1)[1]+(m)[10]*(v1)[2]+(m)[14]*(v1)[3], \
916    (v2)[3]=(m)[ 3]*(v1)[0]+(m)[ 7]*(v1)[1]+(m)[11]*(v1)[2]+(m)[15]*(v1)[3])
917 
918 #define ELL_34V_HOMOG(v2, v1) \
919   ((v2)[0] = (v1)[0]/(v1)[3], \
920    (v2)[1] = (v1)[1]/(v1)[3], \
921    (v2)[2] = (v1)[2]/(v1)[3])
922 
923 /* 3-vector v2 = 4x4 homog coord matrix m times 3-vector v1,
924    using 4-vector tv as tmp */
925 #define ELL_4M3V_HOMOG_MUL(v2, m, v1, tv)                               \
926   ((tv)[0]=(m)[ 0]*(v1)[0]+(m)[ 1]*(v1)[1]+(m)[ 2]*(v1)[2]+(m)[ 3],     \
927    (tv)[1]=(m)[ 4]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 6]*(v1)[2]+(m)[ 7],     \
928    (tv)[2]=(m)[ 8]*(v1)[0]+(m)[ 9]*(v1)[1]+(m)[10]*(v1)[2]+(m)[11],     \
929    (tv)[3]=(m)[12]*(v1)[0]+(m)[13]*(v1)[1]+(m)[14]*(v1)[2]+(m)[15]);    \
930   ELL_34V_HOMOG(v2, tv)
931 
932 #define ELL_34V_HOMOG_TT(v2, TT, v1) \
933   ((v2)[0] = AIR_CAST(TT, (v1)[0]/(v1)[3]), \
934    (v2)[1] = AIR_CAST(TT, (v1)[1]/(v1)[3]), \
935    (v2)[2] = AIR_CAST(TT, (v1)[2]/(v1)[3]))
936 
937 #define ELL_4V_HOMOG(v2, v1)  \
938   ((v2)[0] = (v1)[0]/(v1)[3], \
939    (v2)[1] = (v1)[1]/(v1)[3], \
940    (v2)[2] = (v1)[2]/(v1)[3], \
941    (v2)[3] = 1.0)
942 
943 /*
944 ** These macros are intended to be used as aids with homogeneous transforms
945 */
946 
947 #define ELL_4M_COLS_SET(m, a, b, c, d)                 \
948   (ELL_4V_SET((m)+ 0, (a)[0], (b)[0], (c)[0], (d)[0]), \
949    ELL_4V_SET((m)+ 4, (a)[1], (b)[1], (c)[1], (d)[1]), \
950    ELL_4V_SET((m)+ 8, (a)[2], (b)[2], (c)[2], (d)[2]), \
951    ELL_4V_SET((m)+12, (a)[3], (b)[3], (c)[3], (d)[3]))
952 
953 #define ELL_4M_ROWS_SET(m, a, b, c, d)  \
954   (ELL_4V_COPY((m)+ 0, a),              \
955    ELL_4V_COPY((m)+ 4, b),              \
956    ELL_4V_COPY((m)+ 8, c),              \
957    ELL_4V_COPY((m)+12, d))
958 
959 #define ELL_4M_IDENTITY_SET(m) \
960   (ELL_4V_SET((m)+ 0,  1 ,  0 ,  0 , 0), \
961    ELL_4V_SET((m)+ 4,  0 ,  1 ,  0 , 0), \
962    ELL_4V_SET((m)+ 8,  0 ,  0 ,  1 , 0), \
963    ELL_4V_SET((m)+12,  0 ,  0 ,  0 , 1))
964 
965 #define ELL_4M_EXISTS(m) \
966   (ELL_4V_EXISTS((m) + 0) \
967    && ELL_4V_EXISTS((m) + 4) \
968    && ELL_4V_EXISTS((m) + 8) \
969    && ELL_4V_EXISTS((m) + 12))
970 
971 #define ELL_4M_ZERO_SET(m) \
972   (ELL_4V_SET((m)+ 0,  0 ,  0 ,  0 , 0), \
973    ELL_4V_SET((m)+ 4,  0 ,  0 ,  0 , 0), \
974    ELL_4V_SET((m)+ 8,  0 ,  0 ,  0 , 0), \
975    ELL_4V_SET((m)+12,  0 ,  0 ,  0 , 0))
976 
977 #define ELL_4M_SCALE_SET(m, x, y, z)     \
978   (ELL_4V_SET((m)+ 0, (x),  0 ,  0 , 0), \
979    ELL_4V_SET((m)+ 4,  0 , (y),  0 , 0), \
980    ELL_4V_SET((m)+ 8,  0 ,  0 , (z), 0), \
981    ELL_4V_SET((m)+12,  0 ,  0 ,  0 , 1))
982 
983 #define ELL_4M_TRANSLATE_SET(m, x, y, z) \
984   (ELL_4V_SET((m)+ 0,  1 ,  0 ,  0 , (x)), \
985    ELL_4V_SET((m)+ 4,  0 ,  1 ,  0 , (y)), \
986    ELL_4V_SET((m)+ 8,  0 ,  0 ,  1 , (z)), \
987    ELL_4V_SET((m)+12,  0 ,  0 ,  0 ,  1))
988 
989 #define ELL_4M_ROTATE_X_SET(m, th)                   \
990   (ELL_4V_SET((m)+ 0,  1 ,     0    ,     0    , 0), \
991    ELL_4V_SET((m)+ 4,  0 ,  cos(th) , -sin(th) , 0), \
992    ELL_4V_SET((m)+ 8,  0 , +sin(th) ,  cos(th) , 0), \
993    ELL_4V_SET((m)+12,  0 ,     0    ,     0    , 1))
994 
995 #define ELL_4M_ROTATE_Y_SET(m, th)                   \
996   (ELL_4V_SET((m)+ 0,  cos(th) ,  0 , +sin(th) , 0), \
997    ELL_4V_SET((m)+ 4,     0    ,  1 ,     0    , 0), \
998    ELL_4V_SET((m)+ 8, -sin(th) ,  0 ,  cos(th) , 0), \
999    ELL_4V_SET((m)+12,     0    ,  0 ,     0    , 1))
1000 
1001 #define ELL_4M_ROTATE_Z_SET(m, th)                   \
1002   (ELL_4V_SET((m)+ 0,  cos(th) , -sin(th) ,  0 , 0), \
1003    ELL_4V_SET((m)+ 4, +sin(th) ,  cos(th) ,  0 , 0), \
1004    ELL_4V_SET((m)+ 8,     0    ,     0    ,  1 , 0), \
1005    ELL_4V_SET((m)+12,     0    ,     0    ,  0 , 1))
1006 
1007 #define ELL_4M_NAN_SET(m)   \
1008   (ELL_4V_NAN_SET((m)+ 0),  \
1009    ELL_4V_NAN_SET((m)+ 4),  \
1010    ELL_4V_NAN_SET((m)+ 8),  \
1011    ELL_4V_NAN_SET((m)+ 12))
1012 
1013 #define ELL_4M_MUL(n, l, m)                                                 \
1014   ((n)[ 0]=(l)[ 0]*(m)[ 0]+(l)[ 1]*(m)[ 4]+(l)[ 2]*(m)[ 8]+(l)[ 3]*(m)[12], \
1015    (n)[ 1]=(l)[ 0]*(m)[ 1]+(l)[ 1]*(m)[ 5]+(l)[ 2]*(m)[ 9]+(l)[ 3]*(m)[13], \
1016    (n)[ 2]=(l)[ 0]*(m)[ 2]+(l)[ 1]*(m)[ 6]+(l)[ 2]*(m)[10]+(l)[ 3]*(m)[14], \
1017    (n)[ 3]=(l)[ 0]*(m)[ 3]+(l)[ 1]*(m)[ 7]+(l)[ 2]*(m)[11]+(l)[ 3]*(m)[15], \
1018                                                                             \
1019    (n)[ 4]=(l)[ 4]*(m)[ 0]+(l)[ 5]*(m)[ 4]+(l)[ 6]*(m)[ 8]+(l)[ 7]*(m)[12], \
1020    (n)[ 5]=(l)[ 4]*(m)[ 1]+(l)[ 5]*(m)[ 5]+(l)[ 6]*(m)[ 9]+(l)[ 7]*(m)[13], \
1021    (n)[ 6]=(l)[ 4]*(m)[ 2]+(l)[ 5]*(m)[ 6]+(l)[ 6]*(m)[10]+(l)[ 7]*(m)[14], \
1022    (n)[ 7]=(l)[ 4]*(m)[ 3]+(l)[ 5]*(m)[ 7]+(l)[ 6]*(m)[11]+(l)[ 7]*(m)[15], \
1023                                                                             \
1024    (n)[ 8]=(l)[ 8]*(m)[ 0]+(l)[ 9]*(m)[ 4]+(l)[10]*(m)[ 8]+(l)[11]*(m)[12], \
1025    (n)[ 9]=(l)[ 8]*(m)[ 1]+(l)[ 9]*(m)[ 5]+(l)[10]*(m)[ 9]+(l)[11]*(m)[13], \
1026    (n)[10]=(l)[ 8]*(m)[ 2]+(l)[ 9]*(m)[ 6]+(l)[10]*(m)[10]+(l)[11]*(m)[14], \
1027    (n)[11]=(l)[ 8]*(m)[ 3]+(l)[ 9]*(m)[ 7]+(l)[10]*(m)[11]+(l)[11]*(m)[15], \
1028                                                                             \
1029    (n)[12]=(l)[12]*(m)[ 0]+(l)[13]*(m)[ 4]+(l)[14]*(m)[ 8]+(l)[15]*(m)[12], \
1030    (n)[13]=(l)[12]*(m)[ 1]+(l)[13]*(m)[ 5]+(l)[14]*(m)[ 9]+(l)[15]*(m)[13], \
1031    (n)[14]=(l)[12]*(m)[ 2]+(l)[13]*(m)[ 6]+(l)[14]*(m)[10]+(l)[15]*(m)[14], \
1032    (n)[15]=(l)[12]*(m)[ 3]+(l)[13]*(m)[ 7]+(l)[14]*(m)[11]+(l)[15]*(m)[15])
1033 
1034 #define ELL_34M_EXTRACT(m, l) \
1035   ((m)[0] = (l)[ 0], (m)[1] = (l)[ 1], (m)[2] = (l)[ 2], \
1036    (m)[3] = (l)[ 4], (m)[4] = (l)[ 5], (m)[5] = (l)[ 6], \
1037    (m)[6] = (l)[ 8], (m)[7] = (l)[ 9], (m)[8] = (l)[10])
1038 
1039 #define ELL_43M_INSET(l, m) \
1040   ((l)[ 0] = (m)[0], (l)[ 1] = (m)[1], (l)[ 2] = (m)[2], (l)[ 3] = 0, \
1041    (l)[ 4] = (m)[3], (l)[ 5] = (m)[4], (l)[ 6] = (m)[5], (l)[ 7] = 0, \
1042    (l)[ 8] = (m)[6], (l)[ 9] = (m)[7], (l)[10] = (m)[8], (l)[11] = 0, \
1043    (l)[12] =   0   , (l)[13] =   0   , (l)[14] =   0   , (l)[15] = 1)
1044 
1045 #define ELL_4M_FROB(m) \
1046   (sqrt(ELL_4V_DOT((m)+ 0, (m)+ 0) + \
1047         ELL_4V_DOT((m)+ 4, (m)+ 4) + \
1048         ELL_4V_DOT((m)+ 8, (m)+ 8) + \
1049         ELL_4V_DOT((m)+12, (m)+12)))
1050 
1051 #define ELL_4M_DET(m) \
1052   (  (m)[ 0] * _ELL_3M_DET((m)[ 5], (m)[ 6], (m)[ 7], \
1053                            (m)[ 9], (m)[10], (m)[11], \
1054                            (m)[13], (m)[14], (m)[15]) \
1055    - (m)[ 1] * _ELL_3M_DET((m)[ 4], (m)[ 6], (m)[ 7], \
1056                            (m)[ 8], (m)[10], (m)[11], \
1057                            (m)[12], (m)[14], (m)[15]) \
1058    + (m)[ 2] * _ELL_3M_DET((m)[ 4], (m)[ 5], (m)[ 7], \
1059                            (m)[ 8], (m)[ 9], (m)[11], \
1060                            (m)[12], (m)[13], (m)[15]) \
1061    - (m)[ 3] * _ELL_3M_DET((m)[ 4], (m)[ 5], (m)[ 6], \
1062                            (m)[ 8], (m)[ 9], (m)[10], \
1063                            (m)[12], (m)[13], (m)[14]))
1064 
1065 #define ELL_4MV_OUTER(m, v1, v2)        \
1066   (ELL_4V_SCALE((m)+ 0, (v1)[0], (v2)), \
1067    ELL_4V_SCALE((m)+ 4, (v1)[1], (v2)), \
1068    ELL_4V_SCALE((m)+ 8, (v1)[2], (v2)), \
1069    ELL_4V_SCALE((m)+12, (v1)[3], (v2)))
1070 
1071 #define ELL_4MV_OUTER_TT(m, TT, v1, v2) \
1072   (ELL_4V_SCALE_TT((m)+ 0, TT, (v1)[0], (v2)),   \
1073    ELL_4V_SCALE_TT((m)+ 4, TT, (v1)[1], (v2)),    \
1074    ELL_4V_SCALE_TT((m)+ 8, TT, (v1)[2], (v2)),     \
1075    ELL_4V_SCALE_TT((m)+12, TT, (v1)[3], (v2)))
1076 
1077 #define ELL_4MV_OUTER_INCR(m, v1, v2)        \
1078   (ELL_4V_SCALE_INCR((m)+ 0, (v1)[0], (v2)), \
1079    ELL_4V_SCALE_INCR((m)+ 4, (v1)[1], (v2)), \
1080    ELL_4V_SCALE_INCR((m)+ 8, (v1)[2], (v2)), \
1081    ELL_4V_SCALE_INCR((m)+12, (v1)[3], (v2)))
1082 
1083 #define ELL_Q_MUL(q3, q1, q2)                                            \
1084   ELL_4V_SET((q3),                                                       \
1085   (q1)[0]*(q2)[0] - (q1)[1]*(q2)[1] - (q1)[2]*(q2)[2] - (q1)[3]*(q2)[3], \
1086   (q1)[0]*(q2)[1] + (q1)[1]*(q2)[0] + (q1)[2]*(q2)[3] - (q1)[3]*(q2)[2], \
1087   (q1)[0]*(q2)[2] - (q1)[1]*(q2)[3] + (q1)[2]*(q2)[0] + (q1)[3]*(q2)[1], \
1088   (q1)[0]*(q2)[3] + (q1)[1]*(q2)[2] - (q1)[2]*(q2)[1] + (q1)[3]*(q2)[0])
1089 
1090 #define ELL_Q_CONJ(q2, q1) \
1091   ELL_4V_SET((q2), (q1)[0], -(q1)[1], -(q1)[2], -(q1)[3])
1092 
1093 #define ELL_Q_INV(i, q, n)                                            \
1094  (n = ELL_4V_DOT(q, q),                                               \
1095   ELL_4V_SET((i), (q)[0]/(n), -(q)[1]/(n), -(q)[2]/(n), -(q)[3]/(n)))
1096 
1097 #define ELL_Q_TO_3M(m, q)                                                   \
1098  (ELL_3V_SET((m)+0,                                                         \
1099              (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
1100              2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
1101              2*((q)[1]*(q)[3] + (q)[0]*(q)[2])),                            \
1102   ELL_3V_SET((m)+3,                                                         \
1103              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
1104              (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
1105              2*((q)[2]*(q)[3] - (q)[0]*(q)[1])),                            \
1106   ELL_3V_SET((m)+6,                                                         \
1107              2*((q)[1]*(q)[3] - (q)[0]*(q)[2]),                             \
1108              2*((q)[2]*(q)[3] + (q)[0]*(q)[1]),                             \
1109              (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3]))
1110 
1111 #define ELL_Q_TO_4M(m, q)                                                   \
1112  (ELL_4V_SET((m)+0,                                                         \
1113              (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
1114              2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
1115              2*((q)[1]*(q)[3] + (q)[0]*(q)[2]),                             \
1116              0),                                                            \
1117   ELL_4V_SET((m)+4,                                                         \
1118              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
1119              (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
1120              2*((q)[2]*(q)[3] - (q)[0]*(q)[1]),                             \
1121              0),                                                            \
1122   ELL_4V_SET((m)+8,                                                         \
1123              2*((q)[1]*(q)[3] - (q)[0]*(q)[2]),                             \
1124              2*((q)[2]*(q)[3] + (q)[0]*(q)[1]),                             \
1125              (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3], \
1126              0),                                                            \
1127   ELL_4V_SET((m)+12, 0, 0, 0, 1))
1128 
1129 #define ELL_5V_SET(v, a, b, c, d, e) \
1130   ((v)[0]=(a), (v)[1]=(b), (v)[2]=(c), (v)[3]=(d), (v)[4]=(e))
1131 
1132 #define ELL_5V_COPY(v, w) \
1133   ((v)[0]=(w)[0], (v)[1]=(w)[1], (v)[2]=(w)[2], (v)[3]=(w)[3], (v)[4]=(w)[4])
1134 
1135 #define ELL_6V_SET(v, a, b, c, d, e, f) \
1136   ((v)[0]=(a), (v)[1]=(b), (v)[2]=(c), (v)[3]=(d), (v)[4]=(e), (v)[5]=(f))
1137 
1138 #define ELL_6V_ZERO_SET(v) ELL_6V_SET(v, 0, 0, 0, 0, 0, 0)
1139 
1140 #define ELL_6V_COPY(v2, v1) \
1141   ((v2)[0] = (v1)[0],       \
1142    (v2)[1] = (v1)[1],       \
1143    (v2)[2] = (v1)[2],       \
1144    (v2)[3] = (v1)[3],       \
1145    (v2)[4] = (v1)[4],       \
1146    (v2)[5] = (v1)[5])       \
1147 
1148 #define ELL_6V_INCR(v2, v1) \
1149   ((v2)[0] += (v1)[0],       \
1150    (v2)[1] += (v1)[1],       \
1151    (v2)[2] += (v1)[2],       \
1152    (v2)[3] += (v1)[3],       \
1153    (v2)[4] += (v1)[4],       \
1154    (v2)[5] += (v1)[5])
1155 
1156 #define ELL_6V_SCALE_INCR2(v2, s0, v0, s1, v1) \
1157   ((v2)[0] += (s0)*(v0)[0] + (s1)*(v1)[0],     \
1158    (v2)[1] += (s0)*(v0)[1] + (s1)*(v1)[1],     \
1159    (v2)[2] += (s0)*(v0)[2] + (s1)*(v1)[2],     \
1160    (v2)[3] += (s0)*(v0)[3] + (s1)*(v1)[3],     \
1161    (v2)[4] += (s0)*(v0)[4] + (s1)*(v1)[4],     \
1162    (v2)[5] += (s0)*(v0)[5] + (s1)*(v1)[5])
1163 
1164 #define ELL_6V_SCALE_INCR(v2, a, v1) \
1165   ((v2)[0] += (a)*(v1)[0],       \
1166    (v2)[1] += (a)*(v1)[1],       \
1167    (v2)[2] += (a)*(v1)[2],       \
1168    (v2)[3] += (a)*(v1)[3],       \
1169    (v2)[4] += (a)*(v1)[4],       \
1170    (v2)[5] += (a)*(v1)[5])
1171 
1172 #define ELL_6V_SCALE(v2, a, v1) \
1173   ((v2)[0] = (a)*(v1)[0],       \
1174    (v2)[1] = (a)*(v1)[1],       \
1175    (v2)[2] = (a)*(v1)[2],       \
1176    (v2)[3] = (a)*(v1)[3],       \
1177    (v2)[4] = (a)*(v1)[4],       \
1178    (v2)[5] = (a)*(v1)[5])
1179 
1180 #define ELL_6V_DOT(v1, v2) \
1181   ((v1)[0]*(v2)[0] +       \
1182    (v1)[1]*(v2)[1] +       \
1183    (v1)[2]*(v2)[2] +       \
1184    (v1)[3]*(v2)[3] +       \
1185    (v1)[4]*(v2)[4] +       \
1186    (v1)[5]*(v2)[5])
1187 
1188 #define ELL_6V_SUB(m3, m1, m2) \
1189   ((m3)[0] = (m1)[0] - (m2)[0],  \
1190    (m3)[1] = (m1)[1] - (m2)[1],  \
1191    (m3)[2] = (m1)[2] - (m2)[2],  \
1192    (m3)[3] = (m1)[3] - (m2)[3],  \
1193    (m3)[4] = (m1)[4] - (m2)[4],  \
1194    (m3)[5] = (m1)[5] - (m2)[5])
1195 
1196 #define ELL_6V_ADD2(m3, m1, m2) \
1197   ((m3)[0] = (m1)[0] + (m2)[0],  \
1198    (m3)[1] = (m1)[1] + (m2)[1],  \
1199    (m3)[2] = (m1)[2] + (m2)[2],  \
1200    (m3)[3] = (m1)[3] + (m2)[3],  \
1201    (m3)[4] = (m1)[4] + (m2)[4],  \
1202    (m3)[5] = (m1)[5] + (m2)[5])
1203 
1204 #define ELL_9V_SET(v, a, b, c, d, e, f, g, h, i) \
1205   ((v)[0]=(a), (v)[1]=(b), (v)[2]=(c), \
1206    (v)[3]=(d), (v)[4]=(e), (v)[5]=(f), \
1207    (v)[6]=(g), (v)[7]=(h), (v)[8]=(i))
1208 
1209 #define ELL_9V_COPY(v2, v1) \
1210   ((v2)[0] = (v1)[0],       \
1211    (v2)[1] = (v1)[1],       \
1212    (v2)[2] = (v1)[2],       \
1213    (v2)[3] = (v1)[3],       \
1214    (v2)[4] = (v1)[4],       \
1215    (v2)[5] = (v1)[5],       \
1216    (v2)[6] = (v1)[6],       \
1217    (v2)[7] = (v1)[7],       \
1218    (v2)[8] = (v1)[8])
1219 
1220 #define ELL_9V_SUB(v3, v1, v2)   \
1221   ((v3)[0] = (v1)[0] - (v2)[0],  \
1222    (v3)[1] = (v1)[1] - (v2)[1],  \
1223    (v3)[2] = (v1)[2] - (v2)[2],  \
1224    (v3)[3] = (v1)[3] - (v2)[3],  \
1225    (v3)[4] = (v1)[4] - (v2)[4],  \
1226    (v3)[5] = (v1)[5] - (v2)[5],  \
1227    (v3)[6] = (v1)[6] - (v2)[6],  \
1228    (v3)[7] = (v1)[7] - (v2)[7],  \
1229    (v3)[8] = (v1)[8] - (v2)[8])
1230 
1231 #define ELL_9V_DOT(v1, v2) \
1232   ((v1)[0]*(v2)[0] +       \
1233    (v1)[1]*(v2)[1] +       \
1234    (v1)[2]*(v2)[2] +       \
1235    (v1)[3]*(v2)[3] +       \
1236    (v1)[4]*(v2)[4] +       \
1237    (v1)[5]*(v2)[5] +       \
1238    (v1)[6]*(v2)[6] +       \
1239    (v1)[7]*(v2)[7] +       \
1240    (v1)[8]*(v2)[8])
1241 
1242 #define ELL_9V_LEN(v) (sqrt(ELL_9V_DOT((v),(v))))
1243 
1244 #define ELL_10V_ZERO_SET(v)                         \
1245   ((v)[0]=0, (v)[1]=0, (v)[2]=0,                    \
1246    (v)[3]=0, (v)[4]=0, (v)[5]=0,                    \
1247    (v)[6]=0, (v)[7]=0, (v)[8]=0, (v)[9]=0)
1248 
1249 #define ELL_10V_SET(v, a, b, c, d, e, f, g, h, i, j)      \
1250   ((v)[0]=(a), (v)[1]=(b), (v)[2]=(c),                    \
1251    (v)[3]=(d), (v)[4]=(e), (v)[5]=(f),                    \
1252    (v)[6]=(g), (v)[7]=(h), (v)[8]=(i), (v)[9]=(j))
1253 
1254 #define ELL_10V_COPY(v2, v1) \
1255   ((v2)[0] = (v1)[0],       \
1256    (v2)[1] = (v1)[1],       \
1257    (v2)[2] = (v1)[2],       \
1258    (v2)[3] = (v1)[3],       \
1259    (v2)[4] = (v1)[4],       \
1260    (v2)[5] = (v1)[5],       \
1261    (v2)[6] = (v1)[6],       \
1262    (v2)[7] = (v1)[7],       \
1263    (v2)[8] = (v1)[8],       \
1264    (v2)[9] = (v1)[9])
1265 
1266 #define ELL_10V_SCALE(v2, a, v1) \
1267   ((v2)[0] = (a)*(v1)[0],        \
1268    (v2)[1] = (a)*(v1)[1],        \
1269    (v2)[2] = (a)*(v1)[2],        \
1270    (v2)[3] = (a)*(v1)[3],        \
1271    (v2)[4] = (a)*(v1)[4],        \
1272    (v2)[5] = (a)*(v1)[5],        \
1273    (v2)[6] = (a)*(v1)[6],        \
1274    (v2)[7] = (a)*(v1)[7],        \
1275    (v2)[8] = (a)*(v1)[8],        \
1276    (v2)[9] = (a)*(v1)[9])
1277 
1278 #ifdef __cplusplus
1279    }
1280 #endif
1281 
1282 #endif /* ELLMACROS_HAS_BEEN_INCLUDED */
1283