source: trunk/src/lol/math/matrix.h @ 1138

Last change on this file since 1138 was 1138, checked in by sam, 8 years ago

math: some simplifications in the magic vector templates, and some macros
removed because they were used only once.

  • Property svn:keywords set to Id
File size: 35.5 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11//
12// The Matrix classes
13// ------------------
14//
15
16#if !defined __LOL_MATH_MATRIX_H__
17#define __LOL_MATH_MATRIX_H__
18
19#include <stdint.h>
20#include <cmath>
21#if !defined __ANDROID__
22#   include <iostream>
23#endif
24
25namespace lol
26{
27
28class half;
29class real;
30
31#define VECTOR_TYPES(tname, suffix) \
32    template <typename T> struct tname; \
33    typedef tname<half> f16##suffix; \
34    typedef tname<float> suffix; \
35    typedef tname<double> f64##suffix; \
36    typedef tname<int8_t> i8##suffix; \
37    typedef tname<uint8_t> u8##suffix; \
38    typedef tname<int16_t> i16##suffix; \
39    typedef tname<uint16_t> u16##suffix; \
40    typedef tname<int32_t> i##suffix; \
41    typedef tname<uint32_t> u##suffix; \
42    typedef tname<int64_t> i64##suffix; \
43    typedef tname<uint64_t> u64##suffix;
44
45VECTOR_TYPES(Vec2, vec2)
46VECTOR_TYPES(Cmplx, cmplx)
47VECTOR_TYPES(Vec3, vec3)
48VECTOR_TYPES(Vec4, vec4)
49VECTOR_TYPES(Quat, quat)
50VECTOR_TYPES(Mat4, mat4)
51
52/*
53 * Magic vector swizzling (part 1/2)
54 */
55
56template<typename T, int N> struct XVec2
57{
58    inline Vec2<T> operator =(Vec2<T> that);
59
60    static int const I = (N >> 2) & 3;
61    static int const J = (N >> 0) & 3;
62    T ptr[1 + (I > J ? I : J)];
63};
64
65template<typename T, int N> struct XVec3
66{
67    inline Vec3<T> operator =(Vec3<T> that);
68
69    static int const I = (N >> 4) & 3;
70    static int const J = (N >> 2) & 3;
71    static int const K = (N >> 0) & 3;
72    T ptr[1 + (I > J ? I > K ? I : K
73                     : J > K ? J : K)];
74};
75
76template<typename T, int N> struct XVec4
77{
78    inline Vec4<T> operator =(Vec4<T> that);
79
80    static int const I = (N >> 6) & 3;
81    static int const J = (N >> 4) & 3;
82    static int const K = (N >> 2) & 3;
83    static int const L = (N >> 0) & 3;
84    T ptr[1 + (I > J ? I > K ? I > L ? I : L : K > L ? K : L
85                     : J > K ? J > L ? J : L : K > L ? K : L)];
86};
87
88/*
89 * Helper macros for vector type member functions
90 */
91
92#define MEMBER_OPS() \
93    inline T& operator[](int n) { return *(&x + n); } \
94    inline T const& operator[](int n) const { return *(&x + n); } \
95    \
96    void printf() const;
97
98#define OTHER_MEMBER_OPS(tname) \
99    template<typename U> \
100    inline operator tname<U>() const \
101    { \
102        tname<U> ret; \
103        for (size_t n = 0; n < sizeof(*this) / sizeof(T); n++) \
104            ret[n] = static_cast<U>((*this)[n]); \
105        return ret; \
106    }
107
108/*
109 * 2-element vectors
110 */
111
112template <typename T> struct Vec2
113{
114    inline Vec2() {}
115    inline Vec2(T X, T Y) : x(X), y(Y) {}
116
117    explicit inline Vec2(T X) : x(X), y(X) {}
118
119    template<int N>
120    inline Vec2(XVec2<T, N> const &v)
121      : x(v.ptr[v.I]), y(v.ptr[v.J]) {}
122
123    template<typename U, int N>
124    explicit inline Vec2(XVec2<U, N> const &v)
125      : x(v.ptr[v.I]), y(v.ptr[v.J]) {}
126
127    MEMBER_OPS()
128    OTHER_MEMBER_OPS(Vec2)
129
130#if !defined __ANDROID__
131    template<typename U>
132    friend std::ostream &operator<<(std::ostream &stream, Vec2<U> const &v);
133#endif
134
135    union
136    {
137        struct { T x, y; };
138        struct { T r, g; };
139        struct { T s, t; };
140
141        XVec2<T,0> xx, rr, ss;
142        XVec2<T,1> xy, rg, st;
143        XVec2<T,4> yx, gr, ts;
144        XVec2<T,5> yy, gg, tt;
145        XVec3<T,0> xxx, rrr, sss;
146        XVec3<T,1> xxy, rrg, sst;
147        XVec3<T,4> xyx, rgr, sts;
148        XVec3<T,5> xyy, rgg, stt;
149        XVec3<T,16> yxx, grr, tss;
150        XVec3<T,17> yxy, grg, tst;
151        XVec3<T,20> yyx, ggr, tts;
152        XVec3<T,21> yyy, ggg, ttt;
153        XVec4<T,0> xxxx, rrrr, ssss;
154        XVec4<T,1> xxxy, rrrg, ssst;
155        XVec4<T,4> xxyx, rrgr, ssts;
156        XVec4<T,5> xxyy, rrgg, sstt;
157        XVec4<T,16> xyxx, rgrr, stss;
158        XVec4<T,17> xyxy, rgrg, stst;
159        XVec4<T,20> xyyx, rggr, stts;
160        XVec4<T,21> xyyy, rggg, sttt;
161        XVec4<T,64> yxxx, grrr, tsss;
162        XVec4<T,65> yxxy, grrg, tsst;
163        XVec4<T,68> yxyx, grgr, tsts;
164        XVec4<T,69> yxyy, grgg, tstt;
165        XVec4<T,80> yyxx, ggrr, ttss;
166        XVec4<T,81> yyxy, ggrg, ttst;
167        XVec4<T,84> yyyx, gggr, ttts;
168        XVec4<T,85> yyyy, gggg, tttt;
169    };
170};
171
172/*
173 * 2-element complexes
174 */
175
176template <typename T> struct Cmplx
177{
178    inline Cmplx() {}
179    inline Cmplx(T X) : x(X), y(0) {}
180    inline Cmplx(T X, T Y) : x(X), y(Y) {}
181
182    MEMBER_OPS()
183
184    inline Cmplx<T> operator *(Cmplx<T> const &val) const
185    {
186        return Cmplx<T>(x * val.x - y * val.y, x * val.y + y * val.x);
187    }
188
189    inline Cmplx<T> operator *=(Cmplx<T> const &val)
190    {
191        return *this = (*this) * val;
192    }
193
194    inline Cmplx<T> operator ~() const
195    {
196        return Cmplx<T>(x, -y);
197    }
198
199    inline T norm() const { return len(*this); }
200#if !defined __ANDROID__
201    template<typename U>
202    friend std::ostream &operator<<(std::ostream &stream, Cmplx<U> const &v);
203#endif
204
205    T x, y;
206};
207
208template<typename T>
209static inline Cmplx<T> re(Cmplx<T> const &val)
210{
211    return ~val / sqlen(val);
212}
213
214template<typename T>
215static inline Cmplx<T> operator /(T a, Cmplx<T> const &b)
216{
217    return a * re(b);
218}
219
220template<typename T>
221static inline Cmplx<T> operator /(Cmplx<T> a, Cmplx<T> const &b)
222{
223    return a * re(b);
224}
225
226template<typename T>
227static inline bool operator ==(Cmplx<T> const &a, T b)
228{
229    return (a.x == b) && !a.y;
230}
231
232template<typename T>
233static inline bool operator !=(Cmplx<T> const &a, T b)
234{
235    return (a.x != b) || a.y;
236}
237
238template<typename T>
239static inline bool operator ==(T a, Cmplx<T> const &b) { return b == a; }
240
241template<typename T>
242static inline bool operator !=(T a, Cmplx<T> const &b) { return b != a; }
243
244/*
245 * 3-element vectors
246 */
247
248template <typename T> struct Vec3
249{
250    inline Vec3() {}
251    inline Vec3(T X, T Y, T Z) : x(X), y(Y), z(Z) {}
252    inline Vec3(Vec2<T> XY, T Z) : x(XY.x), y(XY.y), z(Z) {}
253    inline Vec3(T X, Vec2<T> YZ) : x(X), y(YZ.x), z(YZ.y) {}
254
255    explicit inline Vec3(T X) : x(X), y(X), z(X) {}
256
257    template<int N>
258    inline Vec3(XVec3<T, N> const &v)
259      : x(v.ptr[v.I]), y(v.ptr[v.J]), z(v.ptr[v.K]) {}
260
261    template<typename U, int N>
262    explicit inline Vec3(XVec3<U, N> const &v)
263      : x(v.ptr[v.I]), y(v.ptr[v.J]), z(v.ptr[v.K]) {}
264
265    MEMBER_OPS()
266    OTHER_MEMBER_OPS(Vec3)
267
268    template<typename U>
269    friend Vec3<U> cross(Vec3<U>, Vec3<U>);
270
271#if !defined __ANDROID__
272    template<typename U>
273    friend std::ostream &operator<<(std::ostream &stream, Vec3<U> const &v);
274#endif
275
276    union
277    {
278        struct { T x, y, z; };
279        struct { T r, g, b; };
280        struct { T s, t, p; };
281
282        XVec2<T,0> xx, rr, ss;
283        XVec2<T,1> xy, rg, st;
284        XVec2<T,2> xz, rb, sp;
285        XVec2<T,4> yx, gr, ts;
286        XVec2<T,5> yy, gg, tt;
287        XVec2<T,6> yz, gb, tp;
288        XVec2<T,8> zx, br, ps;
289        XVec2<T,9> zy, bg, pt;
290        XVec2<T,10> zz, bb, pp;
291        XVec3<T,0> xxx, rrr, sss;
292        XVec3<T,1> xxy, rrg, sst;
293        XVec3<T,2> xxz, rrb, ssp;
294        XVec3<T,4> xyx, rgr, sts;
295        XVec3<T,5> xyy, rgg, stt;
296        XVec3<T,6> xyz, rgb, stp;
297        XVec3<T,8> xzx, rbr, sps;
298        XVec3<T,9> xzy, rbg, spt;
299        XVec3<T,10> xzz, rbb, spp;
300        XVec3<T,16> yxx, grr, tss;
301        XVec3<T,17> yxy, grg, tst;
302        XVec3<T,18> yxz, grb, tsp;
303        XVec3<T,20> yyx, ggr, tts;
304        XVec3<T,21> yyy, ggg, ttt;
305        XVec3<T,22> yyz, ggb, ttp;
306        XVec3<T,24> yzx, gbr, tps;
307        XVec3<T,25> yzy, gbg, tpt;
308        XVec3<T,26> yzz, gbb, tpp;
309        XVec3<T,32> zxx, brr, pss;
310        XVec3<T,33> zxy, brg, pst;
311        XVec3<T,34> zxz, brb, psp;
312        XVec3<T,36> zyx, bgr, pts;
313        XVec3<T,37> zyy, bgg, ptt;
314        XVec3<T,38> zyz, bgb, ptp;
315        XVec3<T,40> zzx, bbr, pps;
316        XVec3<T,41> zzy, bbg, ppt;
317        XVec3<T,42> zzz, bbb, ppp;
318        XVec4<T,0> xxxx, rrrr, ssss;
319        XVec4<T,1> xxxy, rrrg, ssst;
320        XVec4<T,2> xxxz, rrrb, sssp;
321        XVec4<T,4> xxyx, rrgr, ssts;
322        XVec4<T,5> xxyy, rrgg, sstt;
323        XVec4<T,6> xxyz, rrgb, sstp;
324        XVec4<T,8> xxzx, rrbr, ssps;
325        XVec4<T,9> xxzy, rrbg, sspt;
326        XVec4<T,10> xxzz, rrbb, sspp;
327        XVec4<T,16> xyxx, rgrr, stss;
328        XVec4<T,17> xyxy, rgrg, stst;
329        XVec4<T,18> xyxz, rgrb, stsp;
330        XVec4<T,20> xyyx, rggr, stts;
331        XVec4<T,21> xyyy, rggg, sttt;
332        XVec4<T,22> xyyz, rggb, sttp;
333        XVec4<T,24> xyzx, rgbr, stps;
334        XVec4<T,25> xyzy, rgbg, stpt;
335        XVec4<T,26> xyzz, rgbb, stpp;
336        XVec4<T,32> xzxx, rbrr, spss;
337        XVec4<T,33> xzxy, rbrg, spst;
338        XVec4<T,34> xzxz, rbrb, spsp;
339        XVec4<T,36> xzyx, rbgr, spts;
340        XVec4<T,37> xzyy, rbgg, sptt;
341        XVec4<T,38> xzyz, rbgb, sptp;
342        XVec4<T,40> xzzx, rbbr, spps;
343        XVec4<T,41> xzzy, rbbg, sppt;
344        XVec4<T,42> xzzz, rbbb, sppp;
345        XVec4<T,64> yxxx, grrr, tsss;
346        XVec4<T,65> yxxy, grrg, tsst;
347        XVec4<T,66> yxxz, grrb, tssp;
348        XVec4<T,68> yxyx, grgr, tsts;
349        XVec4<T,69> yxyy, grgg, tstt;
350        XVec4<T,70> yxyz, grgb, tstp;
351        XVec4<T,72> yxzx, grbr, tsps;
352        XVec4<T,73> yxzy, grbg, tspt;
353        XVec4<T,74> yxzz, grbb, tspp;
354        XVec4<T,80> yyxx, ggrr, ttss;
355        XVec4<T,81> yyxy, ggrg, ttst;
356        XVec4<T,82> yyxz, ggrb, ttsp;
357        XVec4<T,84> yyyx, gggr, ttts;
358        XVec4<T,85> yyyy, gggg, tttt;
359        XVec4<T,86> yyyz, gggb, tttp;
360        XVec4<T,88> yyzx, ggbr, ttps;
361        XVec4<T,89> yyzy, ggbg, ttpt;
362        XVec4<T,90> yyzz, ggbb, ttpp;
363        XVec4<T,96> yzxx, gbrr, tpss;
364        XVec4<T,97> yzxy, gbrg, tpst;
365        XVec4<T,98> yzxz, gbrb, tpsp;
366        XVec4<T,100> yzyx, gbgr, tpts;
367        XVec4<T,101> yzyy, gbgg, tptt;
368        XVec4<T,102> yzyz, gbgb, tptp;
369        XVec4<T,104> yzzx, gbbr, tpps;
370        XVec4<T,105> yzzy, gbbg, tppt;
371        XVec4<T,106> yzzz, gbbb, tppp;
372        XVec4<T,128> zxxx, brrr, psss;
373        XVec4<T,129> zxxy, brrg, psst;
374        XVec4<T,130> zxxz, brrb, pssp;
375        XVec4<T,132> zxyx, brgr, psts;
376        XVec4<T,133> zxyy, brgg, pstt;
377        XVec4<T,134> zxyz, brgb, pstp;
378        XVec4<T,136> zxzx, brbr, psps;
379        XVec4<T,137> zxzy, brbg, pspt;
380        XVec4<T,138> zxzz, brbb, pspp;
381        XVec4<T,144> zyxx, bgrr, ptss;
382        XVec4<T,145> zyxy, bgrg, ptst;
383        XVec4<T,146> zyxz, bgrb, ptsp;
384        XVec4<T,148> zyyx, bggr, ptts;
385        XVec4<T,149> zyyy, bggg, pttt;
386        XVec4<T,150> zyyz, bggb, pttp;
387        XVec4<T,152> zyzx, bgbr, ptps;
388        XVec4<T,153> zyzy, bgbg, ptpt;
389        XVec4<T,154> zyzz, bgbb, ptpp;
390        XVec4<T,160> zzxx, bbrr, ppss;
391        XVec4<T,161> zzxy, bbrg, ppst;
392        XVec4<T,162> zzxz, bbrb, ppsp;
393        XVec4<T,164> zzyx, bbgr, ppts;
394        XVec4<T,165> zzyy, bbgg, pptt;
395        XVec4<T,166> zzyz, bbgb, pptp;
396        XVec4<T,168> zzzx, bbbr, ppps;
397        XVec4<T,169> zzzy, bbbg, pppt;
398        XVec4<T,170> zzzz, bbbb, pppp;
399    };
400};
401
402/*
403 * 4-element vectors
404 */
405
406template <typename T> struct Vec4
407{
408    inline Vec4() {}
409    inline Vec4(T X, T Y, T Z, T W) : x(X), y(Y), z(Z), w(W) {}
410    inline Vec4(Vec2<T> XY, T Z, T W) : x(XY.x), y(XY.y), z(Z), w(W) {}
411    inline Vec4(T X, Vec2<T> YZ, T W) : x(X), y(YZ.x), z(YZ.y), w(W) {}
412    inline Vec4(T X, T Y, Vec2<T> ZW) : x(X), y(Y), z(ZW.x), w(ZW.y) {}
413    inline Vec4(Vec2<T> XY, Vec2<T> ZW) : x(XY.x), y(XY.y), z(ZW.x), w(ZW.y) {}
414    inline Vec4(Vec3<T> XYZ, T W) : x(XYZ.x), y(XYZ.y), z(XYZ.z), w(W) {}
415    inline Vec4(T X, Vec3<T> YZW) : x(X), y(YZW.x), z(YZW.y), w(YZW.z) {}
416
417    explicit inline Vec4(T X) : x(X), y(X), z(X), w(X) {}
418
419    template<int N>
420    inline Vec4(XVec4<T, N> const &v)
421      : x(v.ptr[v.I]), y(v.ptr[v.J]), z(v.ptr[v.K]), w(v.ptr[v.L]) {}
422
423    template<typename U, int N>
424    explicit inline Vec4(XVec4<U, N> const &v)
425      : x(v.ptr[v.I]), y(v.ptr[v.J]), z(v.ptr[v.K]), w(v.ptr[v.L]) {}
426
427    MEMBER_OPS()
428    OTHER_MEMBER_OPS(Vec4)
429
430#if !defined __ANDROID__
431    template<typename U>
432    friend std::ostream &operator<<(std::ostream &stream, Vec4<U> const &v);
433#endif
434
435    union
436    {
437        struct { T x, y, z, w; };
438        struct { T r, g, b, a; };
439        struct { T s, t, p, q; };
440
441        XVec2<T,0> xx, rr, ss;
442        XVec2<T,1> xy, rg, st;
443        XVec2<T,2> xz, rb, sp;
444        XVec2<T,3> xw, ra, sq;
445        XVec2<T,4> yx, gr, ts;
446        XVec2<T,5> yy, gg, tt;
447        XVec2<T,6> yz, gb, tp;
448        XVec2<T,7> yw, ga, tq;
449        XVec2<T,8> zx, br, ps;
450        XVec2<T,9> zy, bg, pt;
451        XVec2<T,10> zz, bb, pp;
452        XVec2<T,11> zw, ba, pq;
453        XVec2<T,12> wx, ar, qs;
454        XVec2<T,13> wy, ag, qt;
455        XVec2<T,14> wz, ab, qp;
456        XVec2<T,15> ww, aa, qq;
457        XVec3<T,0> xxx, rrr, sss;
458        XVec3<T,1> xxy, rrg, sst;
459        XVec3<T,2> xxz, rrb, ssp;
460        XVec3<T,3> xxw, rra, ssq;
461        XVec3<T,4> xyx, rgr, sts;
462        XVec3<T,5> xyy, rgg, stt;
463        XVec3<T,6> xyz, rgb, stp;
464        XVec3<T,7> xyw, rga, stq;
465        XVec3<T,8> xzx, rbr, sps;
466        XVec3<T,9> xzy, rbg, spt;
467        XVec3<T,10> xzz, rbb, spp;
468        XVec3<T,11> xzw, rba, spq;
469        XVec3<T,12> xwx, rar, sqs;
470        XVec3<T,13> xwy, rag, sqt;
471        XVec3<T,14> xwz, rab, sqp;
472        XVec3<T,15> xww, raa, sqq;
473        XVec3<T,16> yxx, grr, tss;
474        XVec3<T,17> yxy, grg, tst;
475        XVec3<T,18> yxz, grb, tsp;
476        XVec3<T,19> yxw, gra, tsq;
477        XVec3<T,20> yyx, ggr, tts;
478        XVec3<T,21> yyy, ggg, ttt;
479        XVec3<T,22> yyz, ggb, ttp;
480        XVec3<T,23> yyw, gga, ttq;
481        XVec3<T,24> yzx, gbr, tps;
482        XVec3<T,25> yzy, gbg, tpt;
483        XVec3<T,26> yzz, gbb, tpp;
484        XVec3<T,27> yzw, gba, tpq;
485        XVec3<T,28> ywx, gar, tqs;
486        XVec3<T,29> ywy, gag, tqt;
487        XVec3<T,30> ywz, gab, tqp;
488        XVec3<T,31> yww, gaa, tqq;
489        XVec3<T,32> zxx, brr, pss;
490        XVec3<T,33> zxy, brg, pst;
491        XVec3<T,34> zxz, brb, psp;
492        XVec3<T,35> zxw, bra, psq;
493        XVec3<T,36> zyx, bgr, pts;
494        XVec3<T,37> zyy, bgg, ptt;
495        XVec3<T,38> zyz, bgb, ptp;
496        XVec3<T,39> zyw, bga, ptq;
497        XVec3<T,40> zzx, bbr, pps;
498        XVec3<T,41> zzy, bbg, ppt;
499        XVec3<T,42> zzz, bbb, ppp;
500        XVec3<T,43> zzw, bba, ppq;
501        XVec3<T,44> zwx, bar, pqs;
502        XVec3<T,45> zwy, bag, pqt;
503        XVec3<T,46> zwz, bab, pqp;
504        XVec3<T,47> zww, baa, pqq;
505        XVec3<T,48> wxx, arr, qss;
506        XVec3<T,49> wxy, arg, qst;
507        XVec3<T,50> wxz, arb, qsp;
508        XVec3<T,51> wxw, ara, qsq;
509        XVec3<T,52> wyx, agr, qts;
510        XVec3<T,53> wyy, agg, qtt;
511        XVec3<T,54> wyz, agb, qtp;
512        XVec3<T,55> wyw, aga, qtq;
513        XVec3<T,56> wzx, abr, qps;
514        XVec3<T,57> wzy, abg, qpt;
515        XVec3<T,58> wzz, abb, qpp;
516        XVec3<T,59> wzw, aba, qpq;
517        XVec3<T,60> wwx, aar, qqs;
518        XVec3<T,61> wwy, aag, qqt;
519        XVec3<T,62> wwz, aab, qqp;
520        XVec3<T,63> www, aaa, qqq;
521        XVec4<T,0> xxxx, rrrr, ssss;
522        XVec4<T,1> xxxy, rrrg, ssst;
523        XVec4<T,2> xxxz, rrrb, sssp;
524        XVec4<T,3> xxxw, rrra, sssq;
525        XVec4<T,4> xxyx, rrgr, ssts;
526        XVec4<T,5> xxyy, rrgg, sstt;
527        XVec4<T,6> xxyz, rrgb, sstp;
528        XVec4<T,7> xxyw, rrga, sstq;
529        XVec4<T,8> xxzx, rrbr, ssps;
530        XVec4<T,9> xxzy, rrbg, sspt;
531        XVec4<T,10> xxzz, rrbb, sspp;
532        XVec4<T,11> xxzw, rrba, sspq;
533        XVec4<T,12> xxwx, rrar, ssqs;
534        XVec4<T,13> xxwy, rrag, ssqt;
535        XVec4<T,14> xxwz, rrab, ssqp;
536        XVec4<T,15> xxww, rraa, ssqq;
537        XVec4<T,16> xyxx, rgrr, stss;
538        XVec4<T,17> xyxy, rgrg, stst;
539        XVec4<T,18> xyxz, rgrb, stsp;
540        XVec4<T,19> xyxw, rgra, stsq;
541        XVec4<T,20> xyyx, rggr, stts;
542        XVec4<T,21> xyyy, rggg, sttt;
543        XVec4<T,22> xyyz, rggb, sttp;
544        XVec4<T,23> xyyw, rgga, sttq;
545        XVec4<T,24> xyzx, rgbr, stps;
546        XVec4<T,25> xyzy, rgbg, stpt;
547        XVec4<T,26> xyzz, rgbb, stpp;
548        XVec4<T,27> xyzw, rgba, stpq;
549        XVec4<T,28> xywx, rgar, stqs;
550        XVec4<T,29> xywy, rgag, stqt;
551        XVec4<T,30> xywz, rgab, stqp;
552        XVec4<T,31> xyww, rgaa, stqq;
553        XVec4<T,32> xzxx, rbrr, spss;
554        XVec4<T,33> xzxy, rbrg, spst;
555        XVec4<T,34> xzxz, rbrb, spsp;
556        XVec4<T,35> xzxw, rbra, spsq;
557        XVec4<T,36> xzyx, rbgr, spts;
558        XVec4<T,37> xzyy, rbgg, sptt;
559        XVec4<T,38> xzyz, rbgb, sptp;
560        XVec4<T,39> xzyw, rbga, sptq;
561        XVec4<T,40> xzzx, rbbr, spps;
562        XVec4<T,41> xzzy, rbbg, sppt;
563        XVec4<T,42> xzzz, rbbb, sppp;
564        XVec4<T,43> xzzw, rbba, sppq;
565        XVec4<T,44> xzwx, rbar, spqs;
566        XVec4<T,45> xzwy, rbag, spqt;
567        XVec4<T,46> xzwz, rbab, spqp;
568        XVec4<T,47> xzww, rbaa, spqq;
569        XVec4<T,48> xwxx, rarr, sqss;
570        XVec4<T,49> xwxy, rarg, sqst;
571        XVec4<T,50> xwxz, rarb, sqsp;
572        XVec4<T,51> xwxw, rara, sqsq;
573        XVec4<T,52> xwyx, ragr, sqts;
574        XVec4<T,53> xwyy, ragg, sqtt;
575        XVec4<T,54> xwyz, ragb, sqtp;
576        XVec4<T,55> xwyw, raga, sqtq;
577        XVec4<T,56> xwzx, rabr, sqps;
578        XVec4<T,57> xwzy, rabg, sqpt;
579        XVec4<T,58> xwzz, rabb, sqpp;
580        XVec4<T,59> xwzw, raba, sqpq;
581        XVec4<T,60> xwwx, raar, sqqs;
582        XVec4<T,61> xwwy, raag, sqqt;
583        XVec4<T,62> xwwz, raab, sqqp;
584        XVec4<T,63> xwww, raaa, sqqq;
585        XVec4<T,64> yxxx, grrr, tsss;
586        XVec4<T,65> yxxy, grrg, tsst;
587        XVec4<T,66> yxxz, grrb, tssp;
588        XVec4<T,67> yxxw, grra, tssq;
589        XVec4<T,68> yxyx, grgr, tsts;
590        XVec4<T,69> yxyy, grgg, tstt;
591        XVec4<T,70> yxyz, grgb, tstp;
592        XVec4<T,71> yxyw, grga, tstq;
593        XVec4<T,72> yxzx, grbr, tsps;
594        XVec4<T,73> yxzy, grbg, tspt;
595        XVec4<T,74> yxzz, grbb, tspp;
596        XVec4<T,75> yxzw, grba, tspq;
597        XVec4<T,76> yxwx, grar, tsqs;
598        XVec4<T,77> yxwy, grag, tsqt;
599        XVec4<T,78> yxwz, grab, tsqp;
600        XVec4<T,79> yxww, graa, tsqq;
601        XVec4<T,80> yyxx, ggrr, ttss;
602        XVec4<T,81> yyxy, ggrg, ttst;
603        XVec4<T,82> yyxz, ggrb, ttsp;
604        XVec4<T,83> yyxw, ggra, ttsq;
605        XVec4<T,84> yyyx, gggr, ttts;
606        XVec4<T,85> yyyy, gggg, tttt;
607        XVec4<T,86> yyyz, gggb, tttp;
608        XVec4<T,87> yyyw, ggga, tttq;
609        XVec4<T,88> yyzx, ggbr, ttps;
610        XVec4<T,89> yyzy, ggbg, ttpt;
611        XVec4<T,90> yyzz, ggbb, ttpp;
612        XVec4<T,91> yyzw, ggba, ttpq;
613        XVec4<T,92> yywx, ggar, ttqs;
614        XVec4<T,93> yywy, ggag, ttqt;
615        XVec4<T,94> yywz, ggab, ttqp;
616        XVec4<T,95> yyww, ggaa, ttqq;
617        XVec4<T,96> yzxx, gbrr, tpss;
618        XVec4<T,97> yzxy, gbrg, tpst;
619        XVec4<T,98> yzxz, gbrb, tpsp;
620        XVec4<T,99> yzxw, gbra, tpsq;
621        XVec4<T,100> yzyx, gbgr, tpts;
622        XVec4<T,101> yzyy, gbgg, tptt;
623        XVec4<T,102> yzyz, gbgb, tptp;
624        XVec4<T,103> yzyw, gbga, tptq;
625        XVec4<T,104> yzzx, gbbr, tpps;
626        XVec4<T,105> yzzy, gbbg, tppt;
627        XVec4<T,106> yzzz, gbbb, tppp;
628        XVec4<T,107> yzzw, gbba, tppq;
629        XVec4<T,108> yzwx, gbar, tpqs;
630        XVec4<T,109> yzwy, gbag, tpqt;
631        XVec4<T,110> yzwz, gbab, tpqp;
632        XVec4<T,111> yzww, gbaa, tpqq;
633        XVec4<T,112> ywxx, garr, tqss;
634        XVec4<T,113> ywxy, garg, tqst;
635        XVec4<T,114> ywxz, garb, tqsp;
636        XVec4<T,115> ywxw, gara, tqsq;
637        XVec4<T,116> ywyx, gagr, tqts;
638        XVec4<T,117> ywyy, gagg, tqtt;
639        XVec4<T,118> ywyz, gagb, tqtp;
640        XVec4<T,119> ywyw, gaga, tqtq;
641        XVec4<T,120> ywzx, gabr, tqps;
642        XVec4<T,121> ywzy, gabg, tqpt;
643        XVec4<T,122> ywzz, gabb, tqpp;
644        XVec4<T,123> ywzw, gaba, tqpq;
645        XVec4<T,124> ywwx, gaar, tqqs;
646        XVec4<T,125> ywwy, gaag, tqqt;
647        XVec4<T,126> ywwz, gaab, tqqp;
648        XVec4<T,127> ywww, gaaa, tqqq;
649        XVec4<T,128> zxxx, brrr, psss;
650        XVec4<T,129> zxxy, brrg, psst;
651        XVec4<T,130> zxxz, brrb, pssp;
652        XVec4<T,131> zxxw, brra, pssq;
653        XVec4<T,132> zxyx, brgr, psts;
654        XVec4<T,133> zxyy, brgg, pstt;
655        XVec4<T,134> zxyz, brgb, pstp;
656        XVec4<T,135> zxyw, brga, pstq;
657        XVec4<T,136> zxzx, brbr, psps;
658        XVec4<T,137> zxzy, brbg, pspt;
659        XVec4<T,138> zxzz, brbb, pspp;
660        XVec4<T,139> zxzw, brba, pspq;
661        XVec4<T,140> zxwx, brar, psqs;
662        XVec4<T,141> zxwy, brag, psqt;
663        XVec4<T,142> zxwz, brab, psqp;
664        XVec4<T,143> zxww, braa, psqq;
665        XVec4<T,144> zyxx, bgrr, ptss;
666        XVec4<T,145> zyxy, bgrg, ptst;
667        XVec4<T,146> zyxz, bgrb, ptsp;
668        XVec4<T,147> zyxw, bgra, ptsq;
669        XVec4<T,148> zyyx, bggr, ptts;
670        XVec4<T,149> zyyy, bggg, pttt;
671        XVec4<T,150> zyyz, bggb, pttp;
672        XVec4<T,151> zyyw, bgga, pttq;
673        XVec4<T,152> zyzx, bgbr, ptps;
674        XVec4<T,153> zyzy, bgbg, ptpt;
675        XVec4<T,154> zyzz, bgbb, ptpp;
676        XVec4<T,155> zyzw, bgba, ptpq;
677        XVec4<T,156> zywx, bgar, ptqs;
678        XVec4<T,157> zywy, bgag, ptqt;
679        XVec4<T,158> zywz, bgab, ptqp;
680        XVec4<T,159> zyww, bgaa, ptqq;
681        XVec4<T,160> zzxx, bbrr, ppss;
682        XVec4<T,161> zzxy, bbrg, ppst;
683        XVec4<T,162> zzxz, bbrb, ppsp;
684        XVec4<T,163> zzxw, bbra, ppsq;
685        XVec4<T,164> zzyx, bbgr, ppts;
686        XVec4<T,165> zzyy, bbgg, pptt;
687        XVec4<T,166> zzyz, bbgb, pptp;
688        XVec4<T,167> zzyw, bbga, pptq;
689        XVec4<T,168> zzzx, bbbr, ppps;
690        XVec4<T,169> zzzy, bbbg, pppt;
691        XVec4<T,170> zzzz, bbbb, pppp;
692        XVec4<T,171> zzzw, bbba, pppq;
693        XVec4<T,172> zzwx, bbar, ppqs;
694        XVec4<T,173> zzwy, bbag, ppqt;
695        XVec4<T,174> zzwz, bbab, ppqp;
696        XVec4<T,175> zzww, bbaa, ppqq;
697        XVec4<T,176> zwxx, barr, pqss;
698        XVec4<T,177> zwxy, barg, pqst;
699        XVec4<T,178> zwxz, barb, pqsp;
700        XVec4<T,179> zwxw, bara, pqsq;
701        XVec4<T,180> zwyx, bagr, pqts;
702        XVec4<T,181> zwyy, bagg, pqtt;
703        XVec4<T,182> zwyz, bagb, pqtp;
704        XVec4<T,183> zwyw, baga, pqtq;
705        XVec4<T,184> zwzx, babr, pqps;
706        XVec4<T,185> zwzy, babg, pqpt;
707        XVec4<T,186> zwzz, babb, pqpp;
708        XVec4<T,187> zwzw, baba, pqpq;
709        XVec4<T,188> zwwx, baar, pqqs;
710        XVec4<T,189> zwwy, baag, pqqt;
711        XVec4<T,190> zwwz, baab, pqqp;
712        XVec4<T,191> zwww, baaa, pqqq;
713        XVec4<T,192> wxxx, arrr, qsss;
714        XVec4<T,193> wxxy, arrg, qsst;
715        XVec4<T,194> wxxz, arrb, qssp;
716        XVec4<T,195> wxxw, arra, qssq;
717        XVec4<T,196> wxyx, argr, qsts;
718        XVec4<T,197> wxyy, argg, qstt;
719        XVec4<T,198> wxyz, argb, qstp;
720        XVec4<T,199> wxyw, arga, qstq;
721        XVec4<T,200> wxzx, arbr, qsps;
722        XVec4<T,201> wxzy, arbg, qspt;
723        XVec4<T,202> wxzz, arbb, qspp;
724        XVec4<T,203> wxzw, arba, qspq;
725        XVec4<T,204> wxwx, arar, qsqs;
726        XVec4<T,205> wxwy, arag, qsqt;
727        XVec4<T,206> wxwz, arab, qsqp;
728        XVec4<T,207> wxww, araa, qsqq;
729        XVec4<T,208> wyxx, agrr, qtss;
730        XVec4<T,209> wyxy, agrg, qtst;
731        XVec4<T,210> wyxz, agrb, qtsp;
732        XVec4<T,211> wyxw, agra, qtsq;
733        XVec4<T,212> wyyx, aggr, qtts;
734        XVec4<T,213> wyyy, aggg, qttt;
735        XVec4<T,214> wyyz, aggb, qttp;
736        XVec4<T,215> wyyw, agga, qttq;
737        XVec4<T,216> wyzx, agbr, qtps;
738        XVec4<T,217> wyzy, agbg, qtpt;
739        XVec4<T,218> wyzz, agbb, qtpp;
740        XVec4<T,219> wyzw, agba, qtpq;
741        XVec4<T,220> wywx, agar, qtqs;
742        XVec4<T,221> wywy, agag, qtqt;
743        XVec4<T,222> wywz, agab, qtqp;
744        XVec4<T,223> wyww, agaa, qtqq;
745        XVec4<T,224> wzxx, abrr, qpss;
746        XVec4<T,225> wzxy, abrg, qpst;
747        XVec4<T,226> wzxz, abrb, qpsp;
748        XVec4<T,227> wzxw, abra, qpsq;
749        XVec4<T,228> wzyx, abgr, qpts;
750        XVec4<T,229> wzyy, abgg, qptt;
751        XVec4<T,230> wzyz, abgb, qptp;
752        XVec4<T,231> wzyw, abga, qptq;
753        XVec4<T,232> wzzx, abbr, qpps;
754        XVec4<T,233> wzzy, abbg, qppt;
755        XVec4<T,234> wzzz, abbb, qppp;
756        XVec4<T,235> wzzw, abba, qppq;
757        XVec4<T,236> wzwx, abar, qpqs;
758        XVec4<T,237> wzwy, abag, qpqt;
759        XVec4<T,238> wzwz, abab, qpqp;
760        XVec4<T,239> wzww, abaa, qpqq;
761        XVec4<T,240> wwxx, aarr, qqss;
762        XVec4<T,241> wwxy, aarg, qqst;
763        XVec4<T,242> wwxz, aarb, qqsp;
764        XVec4<T,243> wwxw, aara, qqsq;
765        XVec4<T,244> wwyx, aagr, qqts;
766        XVec4<T,245> wwyy, aagg, qqtt;
767        XVec4<T,246> wwyz, aagb, qqtp;
768        XVec4<T,247> wwyw, aaga, qqtq;
769        XVec4<T,248> wwzx, aabr, qqps;
770        XVec4<T,249> wwzy, aabg, qqpt;
771        XVec4<T,250> wwzz, aabb, qqpp;
772        XVec4<T,251> wwzw, aaba, qqpq;
773        XVec4<T,252> wwwx, aaar, qqqs;
774        XVec4<T,253> wwwy, aaag, qqqt;
775        XVec4<T,254> wwwz, aaab, qqqp;
776        XVec4<T,255> wwww, aaaa, qqqq;
777    };
778};
779
780/*
781 * 4-element quaternions
782 */
783
784template <typename T> struct Quat
785{
786    inline Quat() {}
787    inline Quat(T X) : x(0), y(0), z(0), w(X) {}
788    inline Quat(T X, T Y, T Z, T W) : x(X), y(Y), z(Z), w(W) {}
789
790    Quat(Mat4<T> const &m);
791
792    MEMBER_OPS()
793
794    inline Quat<T> operator *(Quat<T> const &val) const
795    {
796        Quat<T> ret;
797        Vec3<T> v1(x, y, z);
798        Vec3<T> v2(val.x, val.y, val.z);
799        Vec3<T> v3 = cross(v1, v2) + w * v2 + val.w * v1;
800        ret.x = v3.x;
801        ret.y = v3.y;
802        ret.z = v3.z;
803        ret.w = w * val.w - dot(v1, v2);
804        return ret;
805    }
806
807    inline Quat<T> operator *=(Quat<T> const &val)
808    {
809        return *this = (*this) * val;
810    }
811
812    inline Quat<T> operator ~() const
813    {
814        Quat<T> ret;
815        for (int n = 0; n < 3; n++)
816            ret[n] = -(*this)[n];
817        ret[3] = (*this)[3];
818        return ret;
819    }
820
821#if !defined __ANDROID__
822    template<typename U>
823    friend std::ostream &operator<<(std::ostream &stream, Quat<U> const &v);
824#endif
825
826    T x, y, z, w;
827};
828
829template<typename T>
830inline T norm(Quat<T> const &val)
831{
832    return sqlen(val);
833}
834
835template<typename T>
836static inline Quat<T> re(Quat<T> const &val)
837{
838    return ~val / norm(val);
839}
840
841template<typename T>
842static inline Quat<T> operator /(T x, Quat<T> const &y)
843{
844    return x * re(y);
845}
846
847template<typename T>
848static inline Quat<T> operator /(Quat<T> x, Quat<T> const &y)
849{
850    return x * re(y);
851}
852
853/*
854 * Common operators for all vector types, including quaternions
855 */
856
857#define VECTOR_OP(tname, op, tprefix, T) \
858    tprefix \
859    static inline tname<T> operator op(tname<T> const &a, tname<T> const &b) \
860    { \
861        tname<T> ret; \
862        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
863            ret[n] = a[n] op b[n]; \
864        return ret; \
865    } \
866    \
867    tprefix \
868    static inline tname<T> operator op##=(tname<T> &a, tname<T> const &b) \
869    { \
870        return a = a op b; \
871    }
872
873#define BOOL_OP(tname, op, op2, ret, tprefix, T) \
874    tprefix \
875    static inline bool operator op(tname<T> const &a, tname<T> const &b) \
876    { \
877        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
878            if (!(a[n] op2 b[n])) \
879                return !ret; \
880        return ret; \
881    }
882
883#define SCALAR_OP(tname, op, tprefix, T) \
884    tprefix \
885    static inline tname<T> operator op(tname<T> const &a, T const &val) \
886    { \
887        tname<T> ret; \
888        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
889            ret[n] = a[n] op val; \
890        return ret; \
891    } \
892    \
893    tprefix \
894    static inline tname<T> operator op(T const &val, tname<T> const &a) \
895    { \
896        tname<T> ret; \
897        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
898            ret[n] = a[n] op val; \
899        return ret; \
900    } \
901    \
902    tprefix \
903    static inline tname<T> operator op##=(tname<T> &a, T const &val) \
904    { \
905        return a = a op val; \
906    }
907
908#define SCALAR_PROMOTE_OP(tname, op, U) \
909    template<typename T> \
910    static inline tname<U> operator op(U const &val, tname<T> const &a) \
911    { \
912        tname<U> ret; \
913        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
914            ret[n] = val op a[n]; \
915        return ret; \
916    }
917
918#define GLOBAL_OPS(tname, tprefix, T) \
919    SCALAR_OP(tname, *, tprefix, T) \
920    SCALAR_OP(tname, /, tprefix, T) \
921    \
922    VECTOR_OP(tname, -, tprefix, T) \
923    VECTOR_OP(tname, +, tprefix, T) \
924    \
925    BOOL_OP(tname, ==, ==, true, tprefix, T) \
926    BOOL_OP(tname, !=, ==, false, tprefix, T) \
927    \
928    tprefix \
929    static inline tname<T> operator -(tname<T> const &a) \
930    { \
931        tname<T> ret; \
932        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
933            ret[n] = -a[n]; \
934        return ret; \
935    } \
936    \
937    tprefix \
938    static inline T sqlen(tname<T> const &a) \
939    { \
940        T acc = 0; \
941        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
942            acc += a[n] * a[n]; \
943        return acc; \
944    } \
945    \
946    tprefix \
947    static inline double len(tname<T> const &a) \
948    { \
949        using namespace std; \
950        return sqrt((double)sqlen(a)); \
951    } \
952    \
953    tprefix \
954    static inline T dot(tname<T> const &a, tname<T> const &b) \
955    { \
956        T ret = 0; \
957        for (size_t n = 0; n < sizeof(a) / sizeof(T); n++) \
958            ret += a[n] * b[n]; \
959        return ret; \
960    } \
961    \
962    tprefix \
963    static inline tname<T> normalize(tname<T> const &val) \
964    { \
965        T norm = len(val); \
966        return norm ? val / norm : val * (T)0; \
967    }
968
969#define GLOBAL_TEMPLATE_OPS(tname) \
970    GLOBAL_OPS(tname, template<typename T>, T)
971
972#define ALL_GLOBAL_OPS(tname, tprefix, T) \
973    VECTOR_OP(tname, *, tprefix, T) \
974    VECTOR_OP(tname, /, tprefix, T) \
975    \
976    GLOBAL_OPS(tname, tprefix, T) \
977    \
978    BOOL_OP(tname, <=, <=, true, tprefix, T) \
979    BOOL_OP(tname, >=, >=, true, tprefix, T) \
980    BOOL_OP(tname, <, <, true, tprefix, T) \
981    BOOL_OP(tname, >, >, true, tprefix, T)
982
983/* FIXME: a few problems need to be fixed before we can use "half" here. It
984 * will probably never work until we switch to C++11 because it's not really
985 * a POD class. */
986#define GLOBAL_TYPED_OPS(tname) \
987    ALL_GLOBAL_OPS(tname, , float) \
988    ALL_GLOBAL_OPS(tname, , double) \
989    ALL_GLOBAL_OPS(tname, , int8_t) \
990    ALL_GLOBAL_OPS(tname, , uint8_t) \
991    ALL_GLOBAL_OPS(tname, , int16_t) \
992    ALL_GLOBAL_OPS(tname, , uint16_t) \
993    ALL_GLOBAL_OPS(tname, , int32_t) \
994    ALL_GLOBAL_OPS(tname, , uint32_t) \
995    ALL_GLOBAL_OPS(tname, , int64_t) \
996    ALL_GLOBAL_OPS(tname, , uint64_t)
997
998GLOBAL_TEMPLATE_OPS(Cmplx)
999GLOBAL_TEMPLATE_OPS(Quat)
1000
1001GLOBAL_TYPED_OPS(Vec2)
1002GLOBAL_TYPED_OPS(Vec3)
1003GLOBAL_TYPED_OPS(Vec4)
1004
1005/*
1006 * Magic vector swizzling (part 2/2)
1007 */
1008
1009template<typename T, int N>
1010inline Vec2<T> XVec2<T, N>::operator =(Vec2<T> that)
1011{
1012    ptr[I] = that.x; ptr[J] = that.y;
1013    return *this;
1014}
1015
1016template<typename T, int N>
1017inline Vec3<T> XVec3<T, N>::operator =(Vec3<T> that)
1018{
1019    ptr[I] = that.x; ptr[J] = that.y; ptr[K] = that.z;
1020    return *this;
1021}
1022
1023template<typename T, int N>
1024inline Vec4<T> XVec4<T, N>::operator =(Vec4<T> that)
1025{
1026    ptr[I] = that.x; ptr[J] = that.y; ptr[K] = that.z; ptr[L] = that.w;
1027    return *this;
1028}
1029
1030/*
1031 * 4×4-element matrices
1032 */
1033
1034template <typename T> struct Mat4
1035{
1036    inline Mat4() {}
1037    explicit inline Mat4(T val)
1038    {
1039        for (int j = 0; j < 4; j++)
1040            for (int i = 0; i < 4; i++)
1041                v[i][j] = (i == j) ? val : 0;
1042    }
1043    inline Mat4(Vec4<T> v0, Vec4<T> v1, Vec4<T> v2, Vec4<T> v3)
1044    {
1045        v[0] = v0; v[1] = v1; v[2] = v2; v[3] = v3;
1046    }
1047
1048    inline Vec4<T>& operator[](int n) { return v[n]; }
1049    inline Vec4<T> const& operator[](int n) const { return v[n]; }
1050
1051    T det() const;
1052    Mat4<T> invert() const;
1053
1054    /* Helpers for transformation matrices */
1055    static Mat4<T> translate(T x, T y, T z);
1056    static Mat4<T> translate(Vec3<T> v);
1057    static Mat4<T> rotate(T angle, T x, T y, T z);
1058    static Mat4<T> rotate(T angle, Vec3<T> v);
1059    static Mat4<T> rotate(Quat<T> q);
1060
1061    static inline Mat4<T> translate(Mat4<T> mat, Vec3<T> v)
1062    {
1063        return translate(v) * mat;
1064    }
1065
1066    static inline Mat4<T> rotate(Mat4<T> mat, T angle, Vec3<T> v)
1067    {
1068        return rotate(angle, v) * mat;
1069    }
1070
1071    /* Helpers for view matrices */
1072    static Mat4<T> lookat(Vec3<T> eye, Vec3<T> center, Vec3<T> up);
1073
1074    /* Helpers for projection matrices */
1075    static Mat4<T> ortho(T left, T right, T bottom, T top, T near, T far);
1076    static Mat4<T> frustum(T left, T right, T bottom, T top, T near, T far);
1077    static Mat4<T> perspective(T fov_y, T width, T height, T near, T far);
1078
1079    void printf() const;
1080
1081#if !defined __ANDROID__
1082    template<class U>
1083    friend std::ostream &operator<<(std::ostream &stream, Mat4<U> const &m);
1084#endif
1085
1086    inline Mat4<T> operator +(Mat4<T> const val) const
1087    {
1088        Mat4<T> ret;
1089        for (int j = 0; j < 4; j++)
1090            for (int i = 0; i < 4; i++)
1091                ret[i][j] = v[i][j] + val[i][j];
1092        return ret;
1093    }
1094
1095    inline Mat4<T> operator +=(Mat4<T> const val)
1096    {
1097        return *this = *this + val;
1098    }
1099
1100    inline Mat4<T> operator -(Mat4<T> const val) const
1101    {
1102        Mat4<T> ret;
1103        for (int j = 0; j < 4; j++)
1104            for (int i = 0; i < 4; i++)
1105                ret[i][j] = v[i][j] - val[i][j];
1106        return ret;
1107    }
1108
1109    inline Mat4<T> operator -=(Mat4<T> const val)
1110    {
1111        return *this = *this - val;
1112    }
1113
1114    inline Mat4<T> operator *(Mat4<T> const val) const
1115    {
1116        Mat4<T> ret;
1117        for (int j = 0; j < 4; j++)
1118            for (int i = 0; i < 4; i++)
1119            {
1120                T tmp = 0;
1121                for (int k = 0; k < 4; k++)
1122                    tmp += v[k][j] * val[i][k];
1123                ret[i][j] = tmp;
1124            }
1125        return ret;
1126    }
1127
1128    inline Mat4<T> operator *=(Mat4<T> const val)
1129    {
1130        return *this = *this * val;
1131    }
1132
1133    inline Vec4<T> operator *(Vec4<T> const val) const
1134    {
1135        Vec4<T> ret;
1136        for (int j = 0; j < 4; j++)
1137        {
1138            T tmp = 0;
1139            for (int i = 0; i < 4; i++)
1140                tmp += v[i][j] * val[i];
1141            ret[j] = tmp;
1142        }
1143        return ret;
1144    }
1145
1146    Vec4<T> v[4];
1147};
1148
1149/*
1150 * Arbitrarily-sized square matrices; for now this only supports
1151 * naive inversion and is used for the Remez inversion method.
1152 */
1153
1154template<int N, typename T> struct Mat
1155{
1156    inline Mat<N, T>() {}
1157
1158    Mat(T x)
1159    {
1160        for (int j = 0; j < N; j++)
1161            for (int i = 0; i < N; i++)
1162                if (i == j)
1163                    m[i][j] = x;
1164                else
1165                    m[i][j] = 0;
1166    }
1167
1168    /* Naive matrix inversion */
1169    Mat<N, T> inv() const
1170    {
1171        Mat a = *this, b((T)1);
1172
1173        /* Inversion method: iterate through all columns and make sure
1174         * all the terms are 1 on the diagonal and 0 everywhere else */
1175        for (int i = 0; i < N; i++)
1176        {
1177            /* If the expected coefficient is zero, add one of
1178             * the other lines. The first we meet will do. */
1179            if (!a.m[i][i])
1180            {
1181                for (int j = i + 1; j < N; j++)
1182                {
1183                    if (!a.m[i][j])
1184                        continue;
1185                    /* Add row j to row i */
1186                    for (int n = 0; n < N; n++)
1187                    {
1188                        a.m[n][i] += a.m[n][j];
1189                        b.m[n][i] += b.m[n][j];
1190                    }
1191                    break;
1192                }
1193            }
1194
1195            /* Now we know the diagonal term is non-zero. Get its inverse
1196             * and use that to nullify all other terms in the column */
1197            T x = (T)1 / a.m[i][i];
1198            for (int j = 0; j < N; j++)
1199            {
1200                if (j == i)
1201                    continue;
1202                T mul = x * a.m[i][j];
1203                for (int n = 0; n < N; n++)
1204                {
1205                    a.m[n][j] -= mul * a.m[n][i];
1206                    b.m[n][j] -= mul * b.m[n][i];
1207                }
1208            }
1209
1210            /* Finally, ensure the diagonal term is 1 */
1211            for (int n = 0; n < N; n++)
1212            {
1213                a.m[n][i] *= x;
1214                b.m[n][i] *= x;
1215            }
1216        }
1217
1218        return b;
1219    }
1220
1221    T m[N][N];
1222};
1223
1224} /* namespace lol */
1225
1226#endif // __LOL_MATH_MATRIX_H__
1227
Note: See TracBrowser for help on using the repository browser.