source: trunk/src/lol/math/vector.h @ 2056

Last change on this file since 2056 was 2056, checked in by sam, 10 years ago

math: reimplement min(), max(), abs() and fmod() in the lol:: namespace to
avoid conflicts with the C++ stdlib.

  • Property svn:keywords set to Id
File size: 68.3 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2012 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 vector, complex, quaternion and matrix classes
13// --------------------------------------------------
14//
15
16#if !defined __LOL_MATH_VECTOR_H__
17#define __LOL_MATH_VECTOR_H__
18
19#include <stdint.h>
20#include <ostream>
21
22#include "lol/math/math.h"
23#include "lol/math/half.h"
24#include "lol/math/real.h"
25
26namespace lol
27{
28
29/* Some compilers do not support const members in anonymous unions. So
30 * far, GCC (>= 4.6), CLang (3.0) and Visual Studio (>= 2010) appear to
31 * work properly. */
32#undef LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
33#if defined __GNUC__ && (__GNUC__ < 4 || (__GNUC__ == 4 && __GNUC_MINOR__ < 6))
34#   define LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS 1
35#endif
36
37#define DECLARE_VECTOR_TYPEDEFS(tname, suffix) \
38    template <typename T> struct tname; \
39    typedef tname<half> f16##suffix; \
40    typedef tname<float> suffix; \
41    typedef tname<double> d##suffix; \
42    typedef tname<ldouble> f128##suffix; \
43    typedef tname<int8_t> i8##suffix; \
44    typedef tname<uint8_t> u8##suffix; \
45    typedef tname<int16_t> i16##suffix; \
46    typedef tname<uint16_t> u16##suffix; \
47    typedef tname<int32_t> i##suffix; \
48    typedef tname<uint32_t> u##suffix; \
49    typedef tname<int64_t> i64##suffix; \
50    typedef tname<uint64_t> u64##suffix; \
51    typedef tname<real> r##suffix; \
52
53DECLARE_VECTOR_TYPEDEFS(Vec2, vec2)
54DECLARE_VECTOR_TYPEDEFS(Cmplx, cmplx)
55DECLARE_VECTOR_TYPEDEFS(Vec3, vec3)
56DECLARE_VECTOR_TYPEDEFS(Vec4, vec4)
57DECLARE_VECTOR_TYPEDEFS(Quat, quat)
58DECLARE_VECTOR_TYPEDEFS(Mat2, mat2)
59DECLARE_VECTOR_TYPEDEFS(Mat3, mat3)
60DECLARE_VECTOR_TYPEDEFS(Mat4, mat4)
61
62/*
63 * HLSL/Cg-compliant type names.
64 */
65
66typedef vec2 float2;
67typedef vec3 float3;
68typedef vec4 float4;
69typedef mat2 float2x2;
70typedef mat3 float3x3;
71typedef mat4 float4x4;
72
73typedef ivec2 int2;
74typedef ivec3 int3;
75typedef ivec4 int4;
76typedef imat2 int2x2;
77typedef imat3 int3x3;
78typedef imat4 int4x4;
79
80/*
81 * Magic vector swizzling (part 1/2)
82 * These vectors are empty, but thanks to static_cast we can take their
83 * address and access the vector of T's that they are union'ed with. We
84 * use static_cast instead of reinterpret_cast because there is a stronger
85 * guarantee (by the standard) that the address will stay the same across
86 * casts.
87 */
88
89template<typename T, int N> struct XVec2
90{
91    inline Vec2<T> operator =(Vec2<T> const &that);
92
93    inline T& operator[](size_t n)
94    {
95        int i = (N >> (4 * (1 - n))) & 3;
96        return static_cast<T*>(static_cast<void*>(this))[i];
97    }
98    inline T const& operator[](size_t n) const
99    {
100        int i = (N >> (4 * (1 - n))) & 3;
101        return static_cast<T const*>(static_cast<void const *>(this))[i];
102    }
103};
104
105template<typename T, int N> struct XVec3
106{
107    inline Vec3<T> operator =(Vec3<T> const &that);
108
109    inline T& operator[](size_t n)
110    {
111        int i = (N >> (4 * (2 - n))) & 3;
112        return static_cast<T*>(static_cast<void*>(this))[i];
113    }
114    inline T const& operator[](size_t n) const
115    {
116        int i = (N >> (4 * (2 - n))) & 3;
117        return static_cast<T const*>(static_cast<void const *>(this))[i];
118    }
119};
120
121template<typename T, int N> struct XVec4
122{
123    inline Vec4<T> operator =(Vec4<T> const &that);
124
125    inline T& operator[](size_t n)
126    {
127        int i = (N >> (4 * (3 - n))) & 3;
128        return static_cast<T*>(static_cast<void*>(this))[i];
129    }
130    inline T const& operator[](size_t n) const
131    {
132        int i = (N >> (4 * (3 - n))) & 3;
133        return static_cast<T const*>(static_cast<void const *>(this))[i];
134    }
135};
136
137/*
138 * Helper macro for vector type member functions
139 */
140
141#define DECLARE_MEMBER_OPS(tname, first) \
142    inline T& operator[](size_t n) { return *(&this->first + n); } \
143    inline T const& operator[](size_t n) const { return *(&this->first + n); } \
144    \
145    /* Visual Studio insists on having an assignment operator. */ \
146    inline tname<T> const & operator =(tname<T> const &that) \
147    { \
148        for (size_t n = 0; n < sizeof(*this) / sizeof(T); n++) \
149            (*this)[n] = that[n]; \
150        return *this; \
151    } \
152    \
153    template<typename U> \
154    inline operator tname<U>() const \
155    { \
156        tname<U> ret; \
157        for (size_t n = 0; n < sizeof(*this) / sizeof(T); n++) \
158            ret[n] = static_cast<U>((*this)[n]); \
159        return ret; \
160    } \
161    \
162    void printf() const;
163
164/*
165 * 2-element vectors
166 */
167
168template <typename T> struct BVec2
169{
170    explicit inline BVec2() {}
171    explicit inline BVec2(T X, T Y) : x(X), y(Y) {}
172
173    union
174    {
175        struct { T x, y; };
176        struct { T r, g; };
177        struct { T s, t; };
178
179#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
180#   define const /* disabled */
181#endif
182        XVec2<T,0x00> const xx, rr, ss;
183        XVec2<T,0x01> const xy, rg, st; /* lvalue */
184        XVec2<T,0x10> const yx, gr, ts; /* lvalue */
185        XVec2<T,0x11> const yy, gg, tt;
186
187        XVec3<T,0x000> const xxx, rrr, sss;
188        XVec3<T,0x001> const xxy, rrg, sst;
189        XVec3<T,0x010> const xyx, rgr, sts;
190        XVec3<T,0x011> const xyy, rgg, stt;
191        XVec3<T,0x100> const yxx, grr, tss;
192        XVec3<T,0x101> const yxy, grg, tst;
193        XVec3<T,0x110> const yyx, ggr, tts;
194        XVec3<T,0x111> const yyy, ggg, ttt;
195
196        XVec4<T,0x0000> const xxxx, rrrr, ssss;
197        XVec4<T,0x0001> const xxxy, rrrg, ssst;
198        XVec4<T,0x0010> const xxyx, rrgr, ssts;
199        XVec4<T,0x0011> const xxyy, rrgg, sstt;
200        XVec4<T,0x0100> const xyxx, rgrr, stss;
201        XVec4<T,0x0101> const xyxy, rgrg, stst;
202        XVec4<T,0x0110> const xyyx, rggr, stts;
203        XVec4<T,0x0111> const xyyy, rggg, sttt;
204        XVec4<T,0x1000> const yxxx, grrr, tsss;
205        XVec4<T,0x1001> const yxxy, grrg, tsst;
206        XVec4<T,0x1010> const yxyx, grgr, tsts;
207        XVec4<T,0x1011> const yxyy, grgg, tstt;
208        XVec4<T,0x1100> const yyxx, ggrr, ttss;
209        XVec4<T,0x1101> const yyxy, ggrg, ttst;
210        XVec4<T,0x1110> const yyyx, gggr, ttts;
211        XVec4<T,0x1111> const yyyy, gggg, tttt;
212#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
213#   undef const
214#endif
215    };
216};
217
218template <> struct BVec2<half>
219{
220    explicit inline BVec2() {}
221    explicit inline BVec2(half X, half Y) : x(X), y(Y) {}
222
223    half x, y;
224};
225
226template <> struct BVec2<real>
227{
228    explicit inline BVec2() {}
229    explicit inline BVec2(real X, real Y) : x(X), y(Y) {}
230
231    real x, y;
232};
233
234template <typename T> struct Vec2 : BVec2<T>
235{
236    inline Vec2() {}
237    inline Vec2(T X, T Y) : BVec2<T>(X, Y) {}
238
239    explicit inline Vec2(T X) : BVec2<T>(X, X) {}
240
241    template<int N>
242    inline Vec2(XVec2<T, N> const &v)
243      : BVec2<T>(v[0], v[1]) {}
244
245    template<typename U, int N>
246    explicit inline Vec2(XVec2<U, N> const &v)
247      : BVec2<T>(v[0], v[1]) {}
248
249    DECLARE_MEMBER_OPS(Vec2, x)
250
251    template<typename U>
252    friend std::ostream &operator<<(std::ostream &stream, Vec2<U> const &v);
253};
254
255/*
256 * 2-element complexes
257 */
258
259template <typename T> struct Cmplx
260{
261    inline Cmplx() {}
262    inline Cmplx(T X) : x(X), y(0) {}
263    inline Cmplx(T X, T Y) : x(X), y(Y) {}
264
265    DECLARE_MEMBER_OPS(Cmplx, x)
266
267    inline Cmplx<T> operator *(Cmplx<T> const &val) const
268    {
269        return Cmplx<T>(x * val.x - y * val.y, x * val.y + y * val.x);
270    }
271
272    inline Cmplx<T> operator *=(Cmplx<T> const &val)
273    {
274        return *this = (*this) * val;
275    }
276
277    inline Cmplx<T> operator ~() const
278    {
279        return Cmplx<T>(x, -y);
280    }
281
282    inline T norm() const { return length(*this); }
283    template<typename U>
284    friend std::ostream &operator<<(std::ostream &stream, Cmplx<U> const &v);
285
286    T x, y;
287};
288
289template<typename T>
290static inline Cmplx<T> re(Cmplx<T> const &val)
291{
292    return ~val / sqlength(val);
293}
294
295template<typename T>
296static inline Cmplx<T> operator /(T a, Cmplx<T> const &b)
297{
298    return a * re(b);
299}
300
301template<typename T>
302static inline Cmplx<T> operator /(Cmplx<T> a, Cmplx<T> const &b)
303{
304    return a * re(b);
305}
306
307template<typename T>
308static inline bool operator ==(Cmplx<T> const &a, T b)
309{
310    return (a.x == b) && !a.y;
311}
312
313template<typename T>
314static inline bool operator !=(Cmplx<T> const &a, T b)
315{
316    return (a.x != b) || a.y;
317}
318
319template<typename T>
320static inline bool operator ==(T a, Cmplx<T> const &b) { return b == a; }
321
322template<typename T>
323static inline bool operator !=(T a, Cmplx<T> const &b) { return b != a; }
324
325/*
326 * 3-element vectors
327 */
328
329template <typename T> struct BVec3
330{
331    explicit inline BVec3() {}
332    explicit inline BVec3(T X, T Y, T Z) : x(X), y(Y), z(Z) {}
333
334    union
335    {
336        struct { T x, y, z; };
337        struct { T r, g, b; };
338        struct { T s, t, p; };
339
340#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
341#   define const /* disabled */
342#endif
343        XVec2<T,0x00> const xx, rr, ss;
344        XVec2<T,0x01> const xy, rg, st; /* lvalue */
345        XVec2<T,0x02> const xz, rb, sp; /* lvalue */
346        XVec2<T,0x10> const yx, gr, ts; /* lvalue */
347        XVec2<T,0x11> const yy, gg, tt;
348        XVec2<T,0x12> const yz, gb, tp; /* lvalue */
349        XVec2<T,0x20> const zx, br, ps; /* lvalue */
350        XVec2<T,0x21> const zy, bg, pt; /* lvalue */
351        XVec2<T,0x22> const zz, bb, pp;
352
353        XVec3<T,0x000> const xxx, rrr, sss;
354        XVec3<T,0x001> const xxy, rrg, sst;
355        XVec3<T,0x002> const xxz, rrb, ssp;
356        XVec3<T,0x010> const xyx, rgr, sts;
357        XVec3<T,0x011> const xyy, rgg, stt;
358        XVec3<T,0x012> const xyz, rgb, stp; /* lvalue */
359        XVec3<T,0x020> const xzx, rbr, sps;
360        XVec3<T,0x021> const xzy, rbg, spt; /* lvalue */
361        XVec3<T,0x022> const xzz, rbb, spp;
362        XVec3<T,0x100> const yxx, grr, tss;
363        XVec3<T,0x101> const yxy, grg, tst;
364        XVec3<T,0x102> const yxz, grb, tsp; /* lvalue */
365        XVec3<T,0x110> const yyx, ggr, tts;
366        XVec3<T,0x111> const yyy, ggg, ttt;
367        XVec3<T,0x112> const yyz, ggb, ttp;
368        XVec3<T,0x120> const yzx, gbr, tps; /* lvalue */
369        XVec3<T,0x121> const yzy, gbg, tpt;
370        XVec3<T,0x122> const yzz, gbb, tpp;
371        XVec3<T,0x200> const zxx, brr, pss;
372        XVec3<T,0x201> const zxy, brg, pst; /* lvalue */
373        XVec3<T,0x202> const zxz, brb, psp;
374        XVec3<T,0x210> const zyx, bgr, pts; /* lvalue */
375        XVec3<T,0x211> const zyy, bgg, ptt;
376        XVec3<T,0x212> const zyz, bgb, ptp;
377        XVec3<T,0x220> const zzx, bbr, pps;
378        XVec3<T,0x221> const zzy, bbg, ppt;
379        XVec3<T,0x222> const zzz, bbb, ppp;
380
381        XVec4<T,0x0000> const xxxx, rrrr, ssss;
382        XVec4<T,0x0001> const xxxy, rrrg, ssst;
383        XVec4<T,0x0002> const xxxz, rrrb, sssp;
384        XVec4<T,0x0010> const xxyx, rrgr, ssts;
385        XVec4<T,0x0011> const xxyy, rrgg, sstt;
386        XVec4<T,0x0012> const xxyz, rrgb, sstp;
387        XVec4<T,0x0020> const xxzx, rrbr, ssps;
388        XVec4<T,0x0021> const xxzy, rrbg, sspt;
389        XVec4<T,0x0022> const xxzz, rrbb, sspp;
390        XVec4<T,0x0100> const xyxx, rgrr, stss;
391        XVec4<T,0x0101> const xyxy, rgrg, stst;
392        XVec4<T,0x0102> const xyxz, rgrb, stsp;
393        XVec4<T,0x0110> const xyyx, rggr, stts;
394        XVec4<T,0x0111> const xyyy, rggg, sttt;
395        XVec4<T,0x0112> const xyyz, rggb, sttp;
396        XVec4<T,0x0120> const xyzx, rgbr, stps;
397        XVec4<T,0x0121> const xyzy, rgbg, stpt;
398        XVec4<T,0x0122> const xyzz, rgbb, stpp;
399        XVec4<T,0x0200> const xzxx, rbrr, spss;
400        XVec4<T,0x0201> const xzxy, rbrg, spst;
401        XVec4<T,0x0202> const xzxz, rbrb, spsp;
402        XVec4<T,0x0210> const xzyx, rbgr, spts;
403        XVec4<T,0x0211> const xzyy, rbgg, sptt;
404        XVec4<T,0x0212> const xzyz, rbgb, sptp;
405        XVec4<T,0x0220> const xzzx, rbbr, spps;
406        XVec4<T,0x0221> const xzzy, rbbg, sppt;
407        XVec4<T,0x0222> const xzzz, rbbb, sppp;
408        XVec4<T,0x1000> const yxxx, grrr, tsss;
409        XVec4<T,0x1001> const yxxy, grrg, tsst;
410        XVec4<T,0x1002> const yxxz, grrb, tssp;
411        XVec4<T,0x1010> const yxyx, grgr, tsts;
412        XVec4<T,0x1011> const yxyy, grgg, tstt;
413        XVec4<T,0x1012> const yxyz, grgb, tstp;
414        XVec4<T,0x1020> const yxzx, grbr, tsps;
415        XVec4<T,0x1021> const yxzy, grbg, tspt;
416        XVec4<T,0x1022> const yxzz, grbb, tspp;
417        XVec4<T,0x1100> const yyxx, ggrr, ttss;
418        XVec4<T,0x1101> const yyxy, ggrg, ttst;
419        XVec4<T,0x1102> const yyxz, ggrb, ttsp;
420        XVec4<T,0x1110> const yyyx, gggr, ttts;
421        XVec4<T,0x1111> const yyyy, gggg, tttt;
422        XVec4<T,0x1112> const yyyz, gggb, tttp;
423        XVec4<T,0x1120> const yyzx, ggbr, ttps;
424        XVec4<T,0x1121> const yyzy, ggbg, ttpt;
425        XVec4<T,0x1122> const yyzz, ggbb, ttpp;
426        XVec4<T,0x1200> const yzxx, gbrr, tpss;
427        XVec4<T,0x1201> const yzxy, gbrg, tpst;
428        XVec4<T,0x1202> const yzxz, gbrb, tpsp;
429        XVec4<T,0x1210> const yzyx, gbgr, tpts;
430        XVec4<T,0x1211> const yzyy, gbgg, tptt;
431        XVec4<T,0x1212> const yzyz, gbgb, tptp;
432        XVec4<T,0x1220> const yzzx, gbbr, tpps;
433        XVec4<T,0x1221> const yzzy, gbbg, tppt;
434        XVec4<T,0x1222> const yzzz, gbbb, tppp;
435        XVec4<T,0x2000> const zxxx, brrr, psss;
436        XVec4<T,0x2001> const zxxy, brrg, psst;
437        XVec4<T,0x2002> const zxxz, brrb, pssp;
438        XVec4<T,0x2010> const zxyx, brgr, psts;
439        XVec4<T,0x2011> const zxyy, brgg, pstt;
440        XVec4<T,0x2012> const zxyz, brgb, pstp;
441        XVec4<T,0x2020> const zxzx, brbr, psps;
442        XVec4<T,0x2021> const zxzy, brbg, pspt;
443        XVec4<T,0x2022> const zxzz, brbb, pspp;
444        XVec4<T,0x2100> const zyxx, bgrr, ptss;
445        XVec4<T,0x2101> const zyxy, bgrg, ptst;
446        XVec4<T,0x2102> const zyxz, bgrb, ptsp;
447        XVec4<T,0x2110> const zyyx, bggr, ptts;
448        XVec4<T,0x2111> const zyyy, bggg, pttt;
449        XVec4<T,0x2112> const zyyz, bggb, pttp;
450        XVec4<T,0x2120> const zyzx, bgbr, ptps;
451        XVec4<T,0x2121> const zyzy, bgbg, ptpt;
452        XVec4<T,0x2122> const zyzz, bgbb, ptpp;
453        XVec4<T,0x2200> const zzxx, bbrr, ppss;
454        XVec4<T,0x2201> const zzxy, bbrg, ppst;
455        XVec4<T,0x2202> const zzxz, bbrb, ppsp;
456        XVec4<T,0x2210> const zzyx, bbgr, ppts;
457        XVec4<T,0x2211> const zzyy, bbgg, pptt;
458        XVec4<T,0x2212> const zzyz, bbgb, pptp;
459        XVec4<T,0x2220> const zzzx, bbbr, ppps;
460        XVec4<T,0x2221> const zzzy, bbbg, pppt;
461        XVec4<T,0x2222> const zzzz, bbbb, pppp;
462#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
463#   undef const
464#endif
465    };
466};
467
468template <> struct BVec3<half>
469{
470    explicit inline BVec3() {}
471    explicit inline BVec3(half X, half Y, half Z) : x(X), y(Y), z(Z) {}
472
473    half x, y, z;
474};
475
476template <> struct BVec3<real>
477{
478    explicit inline BVec3() {}
479    explicit inline BVec3(real X, real Y, real Z) : x(X), y(Y), z(Z) {}
480
481    real x, y, z;
482};
483
484template <typename T> struct Vec3 : BVec3<T>
485{
486    inline Vec3() {}
487    inline Vec3(T X, T Y, T Z) : BVec3<T>(X, Y, Z) {}
488    inline Vec3(Vec2<T> XY, T Z) : BVec3<T>(XY.x, XY.y, Z) {}
489    inline Vec3(T X, Vec2<T> YZ) : BVec3<T>(X, YZ.x, YZ.y) {}
490
491    explicit inline Vec3(T X) : BVec3<T>(X, X, X) {}
492
493    template<int N>
494    inline Vec3(XVec3<T, N> const &v)
495      : BVec3<T>(v[0], v[1], v[2]) {}
496
497    template<typename U, int N>
498    explicit inline Vec3(XVec3<U, N> const &v)
499      : BVec3<T>(v[0], v[1], v[2]) {}
500
501    static Vec3<T> toeuler(Quat<T> const &q);
502
503    DECLARE_MEMBER_OPS(Vec3, x)
504
505    template<typename U>
506    friend std::ostream &operator<<(std::ostream &stream, Vec3<U> const &v);
507};
508
509/*
510 * 4-element vectors
511 */
512
513template <typename T> struct BVec4
514{
515    explicit inline BVec4() {}
516    explicit inline BVec4(T X, T Y, T Z, T W) : x(X), y(Y), z(Z), w(W) {}
517
518    union
519    {
520        struct { T x, y, z, w; };
521        struct { T r, g, b, a; };
522        struct { T s, t, p, q; };
523
524#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
525#   define const /* disabled */
526#endif
527        XVec2<T,0x00> const xx, rr, ss;
528        XVec2<T,0x01> const xy, rg, st; /* lvalue */
529        XVec2<T,0x02> const xz, rb, sp; /* lvalue */
530        XVec2<T,0x03> const xw, ra, sq; /* lvalue */
531        XVec2<T,0x10> const yx, gr, ts; /* lvalue */
532        XVec2<T,0x11> const yy, gg, tt;
533        XVec2<T,0x12> const yz, gb, tp; /* lvalue */
534        XVec2<T,0x13> const yw, ga, tq; /* lvalue */
535        XVec2<T,0x20> const zx, br, ps; /* lvalue */
536        XVec2<T,0x21> const zy, bg, pt; /* lvalue */
537        XVec2<T,0x22> const zz, bb, pp;
538        XVec2<T,0x23> const zw, ba, pq; /* lvalue */
539        XVec2<T,0x30> const wx, ar, qs; /* lvalue */
540        XVec2<T,0x31> const wy, ag, qt; /* lvalue */
541        XVec2<T,0x32> const wz, ab, qp; /* lvalue */
542        XVec2<T,0x33> const ww, aa, qq;
543
544        XVec3<T,0x000> const xxx, rrr, sss;
545        XVec3<T,0x001> const xxy, rrg, sst;
546        XVec3<T,0x002> const xxz, rrb, ssp;
547        XVec3<T,0x003> const xxw, rra, ssq;
548        XVec3<T,0x010> const xyx, rgr, sts;
549        XVec3<T,0x011> const xyy, rgg, stt;
550        XVec3<T,0x012> const xyz, rgb, stp; /* lvalue */
551        XVec3<T,0x013> const xyw, rga, stq; /* lvalue */
552        XVec3<T,0x020> const xzx, rbr, sps;
553        XVec3<T,0x021> const xzy, rbg, spt; /* lvalue */
554        XVec3<T,0x022> const xzz, rbb, spp;
555        XVec3<T,0x023> const xzw, rba, spq; /* lvalue */
556        XVec3<T,0x030> const xwx, rar, sqs;
557        XVec3<T,0x031> const xwy, rag, sqt; /* lvalue */
558        XVec3<T,0x032> const xwz, rab, sqp; /* lvalue */
559        XVec3<T,0x033> const xww, raa, sqq;
560        XVec3<T,0x100> const yxx, grr, tss;
561        XVec3<T,0x101> const yxy, grg, tst;
562        XVec3<T,0x102> const yxz, grb, tsp; /* lvalue */
563        XVec3<T,0x103> const yxw, gra, tsq; /* lvalue */
564        XVec3<T,0x110> const yyx, ggr, tts;
565        XVec3<T,0x111> const yyy, ggg, ttt;
566        XVec3<T,0x112> const yyz, ggb, ttp;
567        XVec3<T,0x113> const yyw, gga, ttq;
568        XVec3<T,0x120> const yzx, gbr, tps; /* lvalue */
569        XVec3<T,0x121> const yzy, gbg, tpt;
570        XVec3<T,0x122> const yzz, gbb, tpp;
571        XVec3<T,0x123> const yzw, gba, tpq; /* lvalue */
572        XVec3<T,0x130> const ywx, gar, tqs; /* lvalue */
573        XVec3<T,0x131> const ywy, gag, tqt;
574        XVec3<T,0x132> const ywz, gab, tqp; /* lvalue */
575        XVec3<T,0x133> const yww, gaa, tqq;
576        XVec3<T,0x200> const zxx, brr, pss;
577        XVec3<T,0x201> const zxy, brg, pst; /* lvalue */
578        XVec3<T,0x202> const zxz, brb, psp;
579        XVec3<T,0x203> const zxw, bra, psq; /* lvalue */
580        XVec3<T,0x210> const zyx, bgr, pts; /* lvalue */
581        XVec3<T,0x211> const zyy, bgg, ptt;
582        XVec3<T,0x212> const zyz, bgb, ptp;
583        XVec3<T,0x213> const zyw, bga, ptq; /* lvalue */
584        XVec3<T,0x220> const zzx, bbr, pps;
585        XVec3<T,0x221> const zzy, bbg, ppt;
586        XVec3<T,0x222> const zzz, bbb, ppp;
587        XVec3<T,0x223> const zzw, bba, ppq;
588        XVec3<T,0x230> const zwx, bar, pqs; /* lvalue */
589        XVec3<T,0x231> const zwy, bag, pqt; /* lvalue */
590        XVec3<T,0x232> const zwz, bab, pqp;
591        XVec3<T,0x233> const zww, baa, pqq;
592        XVec3<T,0x300> const wxx, arr, qss;
593        XVec3<T,0x301> const wxy, arg, qst; /* lvalue */
594        XVec3<T,0x302> const wxz, arb, qsp; /* lvalue */
595        XVec3<T,0x303> const wxw, ara, qsq;
596        XVec3<T,0x310> const wyx, agr, qts; /* lvalue */
597        XVec3<T,0x311> const wyy, agg, qtt;
598        XVec3<T,0x312> const wyz, agb, qtp; /* lvalue */
599        XVec3<T,0x313> const wyw, aga, qtq;
600        XVec3<T,0x320> const wzx, abr, qps; /* lvalue */
601        XVec3<T,0x321> const wzy, abg, qpt; /* lvalue */
602        XVec3<T,0x322> const wzz, abb, qpp;
603        XVec3<T,0x323> const wzw, aba, qpq;
604        XVec3<T,0x330> const wwx, aar, qqs;
605        XVec3<T,0x331> const wwy, aag, qqt;
606        XVec3<T,0x332> const wwz, aab, qqp;
607        XVec3<T,0x333> const www, aaa, qqq;
608
609        XVec4<T,0x0000> const xxxx, rrrr, ssss;
610        XVec4<T,0x0001> const xxxy, rrrg, ssst;
611        XVec4<T,0x0002> const xxxz, rrrb, sssp;
612        XVec4<T,0x0003> const xxxw, rrra, sssq;
613        XVec4<T,0x0010> const xxyx, rrgr, ssts;
614        XVec4<T,0x0011> const xxyy, rrgg, sstt;
615        XVec4<T,0x0012> const xxyz, rrgb, sstp;
616        XVec4<T,0x0013> const xxyw, rrga, sstq;
617        XVec4<T,0x0020> const xxzx, rrbr, ssps;
618        XVec4<T,0x0021> const xxzy, rrbg, sspt;
619        XVec4<T,0x0022> const xxzz, rrbb, sspp;
620        XVec4<T,0x0023> const xxzw, rrba, sspq;
621        XVec4<T,0x0030> const xxwx, rrar, ssqs;
622        XVec4<T,0x0031> const xxwy, rrag, ssqt;
623        XVec4<T,0x0032> const xxwz, rrab, ssqp;
624        XVec4<T,0x0033> const xxww, rraa, ssqq;
625        XVec4<T,0x0100> const xyxx, rgrr, stss;
626        XVec4<T,0x0101> const xyxy, rgrg, stst;
627        XVec4<T,0x0102> const xyxz, rgrb, stsp;
628        XVec4<T,0x0103> const xyxw, rgra, stsq;
629        XVec4<T,0x0110> const xyyx, rggr, stts;
630        XVec4<T,0x0111> const xyyy, rggg, sttt;
631        XVec4<T,0x0112> const xyyz, rggb, sttp;
632        XVec4<T,0x0113> const xyyw, rgga, sttq;
633        XVec4<T,0x0120> const xyzx, rgbr, stps;
634        XVec4<T,0x0121> const xyzy, rgbg, stpt;
635        XVec4<T,0x0122> const xyzz, rgbb, stpp;
636        XVec4<T,0x0123> const xyzw, rgba, stpq; /* lvalue */
637        XVec4<T,0x0130> const xywx, rgar, stqs;
638        XVec4<T,0x0131> const xywy, rgag, stqt;
639        XVec4<T,0x0132> const xywz, rgab, stqp; /* lvalue */
640        XVec4<T,0x0133> const xyww, rgaa, stqq;
641        XVec4<T,0x0200> const xzxx, rbrr, spss;
642        XVec4<T,0x0201> const xzxy, rbrg, spst;
643        XVec4<T,0x0202> const xzxz, rbrb, spsp;
644        XVec4<T,0x0203> const xzxw, rbra, spsq;
645        XVec4<T,0x0210> const xzyx, rbgr, spts;
646        XVec4<T,0x0211> const xzyy, rbgg, sptt;
647        XVec4<T,0x0212> const xzyz, rbgb, sptp;
648        XVec4<T,0x0213> const xzyw, rbga, sptq; /* lvalue */
649        XVec4<T,0x0220> const xzzx, rbbr, spps;
650        XVec4<T,0x0221> const xzzy, rbbg, sppt;
651        XVec4<T,0x0222> const xzzz, rbbb, sppp;
652        XVec4<T,0x0223> const xzzw, rbba, sppq;
653        XVec4<T,0x0230> const xzwx, rbar, spqs;
654        XVec4<T,0x0231> const xzwy, rbag, spqt; /* lvalue */
655        XVec4<T,0x0232> const xzwz, rbab, spqp;
656        XVec4<T,0x0233> const xzww, rbaa, spqq;
657        XVec4<T,0x0300> const xwxx, rarr, sqss;
658        XVec4<T,0x0301> const xwxy, rarg, sqst;
659        XVec4<T,0x0302> const xwxz, rarb, sqsp;
660        XVec4<T,0x0303> const xwxw, rara, sqsq;
661        XVec4<T,0x0310> const xwyx, ragr, sqts;
662        XVec4<T,0x0311> const xwyy, ragg, sqtt;
663        XVec4<T,0x0312> const xwyz, ragb, sqtp; /* lvalue */
664        XVec4<T,0x0313> const xwyw, raga, sqtq;
665        XVec4<T,0x0320> const xwzx, rabr, sqps;
666        XVec4<T,0x0321> const xwzy, rabg, sqpt; /* lvalue */
667        XVec4<T,0x0322> const xwzz, rabb, sqpp;
668        XVec4<T,0x0323> const xwzw, raba, sqpq;
669        XVec4<T,0x0330> const xwwx, raar, sqqs;
670        XVec4<T,0x0331> const xwwy, raag, sqqt;
671        XVec4<T,0x0332> const xwwz, raab, sqqp;
672        XVec4<T,0x0333> const xwww, raaa, sqqq;
673        XVec4<T,0x1000> const yxxx, grrr, tsss;
674        XVec4<T,0x1001> const yxxy, grrg, tsst;
675        XVec4<T,0x1002> const yxxz, grrb, tssp;
676        XVec4<T,0x1003> const yxxw, grra, tssq;
677        XVec4<T,0x1010> const yxyx, grgr, tsts;
678        XVec4<T,0x1011> const yxyy, grgg, tstt;
679        XVec4<T,0x1012> const yxyz, grgb, tstp;
680        XVec4<T,0x1013> const yxyw, grga, tstq;
681        XVec4<T,0x1020> const yxzx, grbr, tsps;
682        XVec4<T,0x1021> const yxzy, grbg, tspt;
683        XVec4<T,0x1022> const yxzz, grbb, tspp;
684        XVec4<T,0x1023> const yxzw, grba, tspq; /* lvalue */
685        XVec4<T,0x1030> const yxwx, grar, tsqs;
686        XVec4<T,0x1031> const yxwy, grag, tsqt;
687        XVec4<T,0x1032> const yxwz, grab, tsqp; /* lvalue */
688        XVec4<T,0x1033> const yxww, graa, tsqq;
689        XVec4<T,0x1100> const yyxx, ggrr, ttss;
690        XVec4<T,0x1101> const yyxy, ggrg, ttst;
691        XVec4<T,0x1102> const yyxz, ggrb, ttsp;
692        XVec4<T,0x1103> const yyxw, ggra, ttsq;
693        XVec4<T,0x1110> const yyyx, gggr, ttts;
694        XVec4<T,0x1111> const yyyy, gggg, tttt;
695        XVec4<T,0x1112> const yyyz, gggb, tttp;
696        XVec4<T,0x1113> const yyyw, ggga, tttq;
697        XVec4<T,0x1120> const yyzx, ggbr, ttps;
698        XVec4<T,0x1121> const yyzy, ggbg, ttpt;
699        XVec4<T,0x1122> const yyzz, ggbb, ttpp;
700        XVec4<T,0x1123> const yyzw, ggba, ttpq;
701        XVec4<T,0x1130> const yywx, ggar, ttqs;
702        XVec4<T,0x1131> const yywy, ggag, ttqt;
703        XVec4<T,0x1132> const yywz, ggab, ttqp;
704        XVec4<T,0x1133> const yyww, ggaa, ttqq;
705        XVec4<T,0x1200> const yzxx, gbrr, tpss;
706        XVec4<T,0x1201> const yzxy, gbrg, tpst;
707        XVec4<T,0x1202> const yzxz, gbrb, tpsp;
708        XVec4<T,0x1203> const yzxw, gbra, tpsq; /* lvalue */
709        XVec4<T,0x1210> const yzyx, gbgr, tpts;
710        XVec4<T,0x1211> const yzyy, gbgg, tptt;
711        XVec4<T,0x1212> const yzyz, gbgb, tptp;
712        XVec4<T,0x1213> const yzyw, gbga, tptq;
713        XVec4<T,0x1220> const yzzx, gbbr, tpps;
714        XVec4<T,0x1221> const yzzy, gbbg, tppt;
715        XVec4<T,0x1222> const yzzz, gbbb, tppp;
716        XVec4<T,0x1223> const yzzw, gbba, tppq;
717        XVec4<T,0x1230> const yzwx, gbar, tpqs; /* lvalue */
718        XVec4<T,0x1231> const yzwy, gbag, tpqt;
719        XVec4<T,0x1232> const yzwz, gbab, tpqp;
720        XVec4<T,0x1233> const yzww, gbaa, tpqq;
721        XVec4<T,0x1300> const ywxx, garr, tqss;
722        XVec4<T,0x1301> const ywxy, garg, tqst;
723        XVec4<T,0x1302> const ywxz, garb, tqsp; /* lvalue */
724        XVec4<T,0x1303> const ywxw, gara, tqsq;
725        XVec4<T,0x1310> const ywyx, gagr, tqts;
726        XVec4<T,0x1311> const ywyy, gagg, tqtt;
727        XVec4<T,0x1312> const ywyz, gagb, tqtp;
728        XVec4<T,0x1313> const ywyw, gaga, tqtq;
729        XVec4<T,0x1320> const ywzx, gabr, tqps; /* lvalue */
730        XVec4<T,0x1321> const ywzy, gabg, tqpt;
731        XVec4<T,0x1322> const ywzz, gabb, tqpp;
732        XVec4<T,0x1323> const ywzw, gaba, tqpq;
733        XVec4<T,0x1330> const ywwx, gaar, tqqs;
734        XVec4<T,0x1331> const ywwy, gaag, tqqt;
735        XVec4<T,0x1332> const ywwz, gaab, tqqp;
736        XVec4<T,0x1333> const ywww, gaaa, tqqq;
737        XVec4<T,0x2000> const zxxx, brrr, psss;
738        XVec4<T,0x2001> const zxxy, brrg, psst;
739        XVec4<T,0x2002> const zxxz, brrb, pssp;
740        XVec4<T,0x2003> const zxxw, brra, pssq;
741        XVec4<T,0x2010> const zxyx, brgr, psts;
742        XVec4<T,0x2011> const zxyy, brgg, pstt;
743        XVec4<T,0x2012> const zxyz, brgb, pstp;
744        XVec4<T,0x2013> const zxyw, brga, pstq; /* lvalue */
745        XVec4<T,0x2020> const zxzx, brbr, psps;
746        XVec4<T,0x2021> const zxzy, brbg, pspt;
747        XVec4<T,0x2022> const zxzz, brbb, pspp;
748        XVec4<T,0x2023> const zxzw, brba, pspq;
749        XVec4<T,0x2030> const zxwx, brar, psqs;
750        XVec4<T,0x2031> const zxwy, brag, psqt; /* lvalue */
751        XVec4<T,0x2032> const zxwz, brab, psqp;
752        XVec4<T,0x2033> const zxww, braa, psqq;
753        XVec4<T,0x2100> const zyxx, bgrr, ptss;
754        XVec4<T,0x2101> const zyxy, bgrg, ptst;
755        XVec4<T,0x2102> const zyxz, bgrb, ptsp;
756        XVec4<T,0x2103> const zyxw, bgra, ptsq; /* lvalue */
757        XVec4<T,0x2110> const zyyx, bggr, ptts;
758        XVec4<T,0x2111> const zyyy, bggg, pttt;
759        XVec4<T,0x2112> const zyyz, bggb, pttp;
760        XVec4<T,0x2113> const zyyw, bgga, pttq;
761        XVec4<T,0x2120> const zyzx, bgbr, ptps;
762        XVec4<T,0x2121> const zyzy, bgbg, ptpt;
763        XVec4<T,0x2122> const zyzz, bgbb, ptpp;
764        XVec4<T,0x2123> const zyzw, bgba, ptpq;
765        XVec4<T,0x2130> const zywx, bgar, ptqs; /* lvalue */
766        XVec4<T,0x2131> const zywy, bgag, ptqt;
767        XVec4<T,0x2132> const zywz, bgab, ptqp;
768        XVec4<T,0x2133> const zyww, bgaa, ptqq;
769        XVec4<T,0x2200> const zzxx, bbrr, ppss;
770        XVec4<T,0x2201> const zzxy, bbrg, ppst;
771        XVec4<T,0x2202> const zzxz, bbrb, ppsp;
772        XVec4<T,0x2203> const zzxw, bbra, ppsq;
773        XVec4<T,0x2210> const zzyx, bbgr, ppts;
774        XVec4<T,0x2211> const zzyy, bbgg, pptt;
775        XVec4<T,0x2212> const zzyz, bbgb, pptp;
776        XVec4<T,0x2213> const zzyw, bbga, pptq;
777        XVec4<T,0x2220> const zzzx, bbbr, ppps;
778        XVec4<T,0x2221> const zzzy, bbbg, pppt;
779        XVec4<T,0x2222> const zzzz, bbbb, pppp;
780        XVec4<T,0x2223> const zzzw, bbba, pppq;
781        XVec4<T,0x2230> const zzwx, bbar, ppqs;
782        XVec4<T,0x2231> const zzwy, bbag, ppqt;
783        XVec4<T,0x2232> const zzwz, bbab, ppqp;
784        XVec4<T,0x2233> const zzww, bbaa, ppqq;
785        XVec4<T,0x2300> const zwxx, barr, pqss;
786        XVec4<T,0x2301> const zwxy, barg, pqst; /* lvalue */
787        XVec4<T,0x2302> const zwxz, barb, pqsp;
788        XVec4<T,0x2303> const zwxw, bara, pqsq;
789        XVec4<T,0x2310> const zwyx, bagr, pqts; /* lvalue */
790        XVec4<T,0x2311> const zwyy, bagg, pqtt;
791        XVec4<T,0x2312> const zwyz, bagb, pqtp;
792        XVec4<T,0x2313> const zwyw, baga, pqtq;
793        XVec4<T,0x2320> const zwzx, babr, pqps;
794        XVec4<T,0x2321> const zwzy, babg, pqpt;
795        XVec4<T,0x2322> const zwzz, babb, pqpp;
796        XVec4<T,0x2323> const zwzw, baba, pqpq;
797        XVec4<T,0x2330> const zwwx, baar, pqqs;
798        XVec4<T,0x2331> const zwwy, baag, pqqt;
799        XVec4<T,0x2332> const zwwz, baab, pqqp;
800        XVec4<T,0x2333> const zwww, baaa, pqqq;
801        XVec4<T,0x3000> const wxxx, arrr, qsss;
802        XVec4<T,0x3001> const wxxy, arrg, qsst;
803        XVec4<T,0x3002> const wxxz, arrb, qssp;
804        XVec4<T,0x3003> const wxxw, arra, qssq;
805        XVec4<T,0x3010> const wxyx, argr, qsts;
806        XVec4<T,0x3011> const wxyy, argg, qstt;
807        XVec4<T,0x3012> const wxyz, argb, qstp; /* lvalue */
808        XVec4<T,0x3013> const wxyw, arga, qstq;
809        XVec4<T,0x3020> const wxzx, arbr, qsps;
810        XVec4<T,0x3021> const wxzy, arbg, qspt; /* lvalue */
811        XVec4<T,0x3022> const wxzz, arbb, qspp;
812        XVec4<T,0x3023> const wxzw, arba, qspq;
813        XVec4<T,0x3030> const wxwx, arar, qsqs;
814        XVec4<T,0x3031> const wxwy, arag, qsqt;
815        XVec4<T,0x3032> const wxwz, arab, qsqp;
816        XVec4<T,0x3033> const wxww, araa, qsqq;
817        XVec4<T,0x3100> const wyxx, agrr, qtss;
818        XVec4<T,0x3101> const wyxy, agrg, qtst;
819        XVec4<T,0x3102> const wyxz, agrb, qtsp; /* lvalue */
820        XVec4<T,0x3103> const wyxw, agra, qtsq;
821        XVec4<T,0x3110> const wyyx, aggr, qtts;
822        XVec4<T,0x3111> const wyyy, aggg, qttt;
823        XVec4<T,0x3112> const wyyz, aggb, qttp;
824        XVec4<T,0x3113> const wyyw, agga, qttq;
825        XVec4<T,0x3120> const wyzx, agbr, qtps; /* lvalue */
826        XVec4<T,0x3121> const wyzy, agbg, qtpt;
827        XVec4<T,0x3122> const wyzz, agbb, qtpp;
828        XVec4<T,0x3123> const wyzw, agba, qtpq;
829        XVec4<T,0x3130> const wywx, agar, qtqs;
830        XVec4<T,0x3131> const wywy, agag, qtqt;
831        XVec4<T,0x3132> const wywz, agab, qtqp;
832        XVec4<T,0x3133> const wyww, agaa, qtqq;
833        XVec4<T,0x3200> const wzxx, abrr, qpss;
834        XVec4<T,0x3201> const wzxy, abrg, qpst; /* lvalue */
835        XVec4<T,0x3202> const wzxz, abrb, qpsp;
836        XVec4<T,0x3203> const wzxw, abra, qpsq;
837        XVec4<T,0x3210> const wzyx, abgr, qpts; /* lvalue */
838        XVec4<T,0x3211> const wzyy, abgg, qptt;
839        XVec4<T,0x3212> const wzyz, abgb, qptp;
840        XVec4<T,0x3213> const wzyw, abga, qptq;
841        XVec4<T,0x3220> const wzzx, abbr, qpps;
842        XVec4<T,0x3221> const wzzy, abbg, qppt;
843        XVec4<T,0x3222> const wzzz, abbb, qppp;
844        XVec4<T,0x3223> const wzzw, abba, qppq;
845        XVec4<T,0x3230> const wzwx, abar, qpqs;
846        XVec4<T,0x3231> const wzwy, abag, qpqt;
847        XVec4<T,0x3232> const wzwz, abab, qpqp;
848        XVec4<T,0x3233> const wzww, abaa, qpqq;
849        XVec4<T,0x3300> const wwxx, aarr, qqss;
850        XVec4<T,0x3301> const wwxy, aarg, qqst;
851        XVec4<T,0x3302> const wwxz, aarb, qqsp;
852        XVec4<T,0x3303> const wwxw, aara, qqsq;
853        XVec4<T,0x3310> const wwyx, aagr, qqts;
854        XVec4<T,0x3311> const wwyy, aagg, qqtt;
855        XVec4<T,0x3312> const wwyz, aagb, qqtp;
856        XVec4<T,0x3313> const wwyw, aaga, qqtq;
857        XVec4<T,0x3320> const wwzx, aabr, qqps;
858        XVec4<T,0x3321> const wwzy, aabg, qqpt;
859        XVec4<T,0x3322> const wwzz, aabb, qqpp;
860        XVec4<T,0x3323> const wwzw, aaba, qqpq;
861        XVec4<T,0x3330> const wwwx, aaar, qqqs;
862        XVec4<T,0x3331> const wwwy, aaag, qqqt;
863        XVec4<T,0x3332> const wwwz, aaab, qqqp;
864        XVec4<T,0x3333> const wwww, aaaa, qqqq;
865#if LOL_NO_CONST_MEMBERS_IN_ANONYMOUS_UNIONS
866#   undef const
867#endif
868    };
869};
870
871template <> struct BVec4<half>
872{
873    explicit inline BVec4() {}
874    explicit inline BVec4(half X, half Y, half Z, half W)
875     : x(X), y(Y), z(Z), w(W) {}
876
877    half x, y, z, w;
878};
879
880template <> struct BVec4<real>
881{
882    explicit inline BVec4() {}
883    explicit inline BVec4(real X, real Y, real Z, real W)
884     : x(X), y(Y), z(Z), w(W) {}
885
886    real x, y, z, w;
887};
888
889template <typename T> struct Vec4 : BVec4<T>
890{
891    inline Vec4() {}
892    inline Vec4(T X, T Y, T Z, T W) : BVec4<T>(X, Y, Z, W) {}
893    inline Vec4(Vec2<T> XY, T Z, T W) : BVec4<T>(XY.x, XY.y, Z, W) {}
894    inline Vec4(T X, Vec2<T> YZ, T W) : BVec4<T>(X, YZ.x, YZ.y, W) {}
895    inline Vec4(T X, T Y, Vec2<T> ZW) : BVec4<T>(X, Y, ZW.x, ZW.y) {}
896    inline Vec4(Vec2<T> XY, Vec2<T> ZW) : BVec4<T>(XY.x, XY.y, ZW.x, ZW.y) {}
897    inline Vec4(Vec3<T> XYZ, T W) : BVec4<T>(XYZ.x, XYZ.y, XYZ.z, W) {}
898    inline Vec4(T X, Vec3<T> YZW) : BVec4<T>(X, YZW.x, YZW.y, YZW.z) {}
899
900    explicit inline Vec4(T X) : BVec4<T>(X, X, X, X) {}
901
902    template<int N>
903    inline Vec4(XVec4<T, N> const &v)
904      : BVec4<T>(v[0], v[1], v[2], v[3]) {}
905
906    template<typename U, int N>
907    explicit inline Vec4(XVec4<U, N> const &v)
908      : BVec4<T>(v[0], v[1], v[2], v[3]) {}
909
910    DECLARE_MEMBER_OPS(Vec4, x)
911
912    template<typename U>
913    friend std::ostream &operator<<(std::ostream &stream, Vec4<U> const &v);
914};
915
916/*
917 * 4-element quaternions
918 */
919
920template <typename T> struct Quat
921{
922    inline Quat() {}
923    inline Quat(T W) : w(W),  x(0), y(0), z(0) {}
924    inline Quat(T W, T X, T Y, T Z) : w(W), x(X), y(Y), z(Z) {}
925
926    Quat(Mat3<T> const &m);
927    Quat(Mat4<T> const &m);
928
929    DECLARE_MEMBER_OPS(Quat, w)
930
931    static Quat<T> rotate(T angle, T x, T y, T z);
932    static Quat<T> rotate(T angle, Vec3<T> const &v);
933
934    /* Convert from Euler angles. The axes in fromeuler_xyx are
935     * x, then y', then x", ie. the axes are attached to the model.
936     * If you want to rotate around static axes, just reverse the order
937     * of the arguments. */
938    static Quat<T> fromeuler_xyx(Vec3<T> const &v);
939    static Quat<T> fromeuler_xzx(Vec3<T> const &v);
940    static Quat<T> fromeuler_yxy(Vec3<T> const &v);
941    static Quat<T> fromeuler_yzy(Vec3<T> const &v);
942    static Quat<T> fromeuler_zxz(Vec3<T> const &v);
943    static Quat<T> fromeuler_zyz(Vec3<T> const &v);
944    static Quat<T> fromeuler_xyx(T phi, T theta, T psi);
945    static Quat<T> fromeuler_xzx(T phi, T theta, T psi);
946    static Quat<T> fromeuler_yxy(T phi, T theta, T psi);
947    static Quat<T> fromeuler_yzy(T phi, T theta, T psi);
948    static Quat<T> fromeuler_zxz(T phi, T theta, T psi);
949    static Quat<T> fromeuler_zyz(T phi, T theta, T psi);
950
951    /* Convert from Tait-Bryan angles (incorrectly called Euler angles,
952     * but since everyone does it…). The axes in fromeuler_xyz are
953     * x, then y', then z", ie. the axes are attached to the model.
954     * If you want to apply yaw around x, pitch around y, and roll
955     * around z, use fromeuler_xyz.
956     * If you want to rotate around static axes, reverse the order in
957     * the function name (_zyx instead of _xyz) AND reverse the order
958     * of the arguments. */
959    static Quat<T> fromeuler_xyz(Vec3<T> const &v);
960    static Quat<T> fromeuler_xzy(Vec3<T> const &v);
961    static Quat<T> fromeuler_yxz(Vec3<T> const &v);
962    static Quat<T> fromeuler_yzx(Vec3<T> const &v);
963    static Quat<T> fromeuler_zxy(Vec3<T> const &v);
964    static Quat<T> fromeuler_zyx(Vec3<T> const &v);
965    static Quat<T> fromeuler_xyz(T phi, T theta, T psi);
966    static Quat<T> fromeuler_xzy(T phi, T theta, T psi);
967    static Quat<T> fromeuler_yxz(T phi, T theta, T psi);
968    static Quat<T> fromeuler_yzx(T phi, T theta, T psi);
969    static Quat<T> fromeuler_zxy(T phi, T theta, T psi);
970    static Quat<T> fromeuler_zyx(T phi, T theta, T psi);
971
972    inline Quat<T> operator *(Quat<T> const &val) const;
973
974    inline Quat<T> operator *=(Quat<T> const &val)
975    {
976        return *this = (*this) * val;
977    }
978
979    inline Quat<T> operator ~() const
980    {
981        return Quat<T>(w, -x, -y, -z);
982    }
983
984    inline Vec3<T> transform(Vec3<T> const &v)
985    {
986        Quat<T> p = Quat<T>(0, v.x, v.y, v.z);
987        Quat<T> q = *this * p / *this;
988        return Vec3<T>(q.x, q.y, q.z);
989    }
990
991    template<typename U>
992    friend std::ostream &operator<<(std::ostream &stream, Quat<U> const &v);
993
994    /* XXX: storage order is wxyz, unlike vectors! */
995    T w, x, y, z;
996};
997
998template<typename T>
999inline T norm(Quat<T> const &val)
1000{
1001    return sqlength(val);
1002}
1003
1004template<typename T>
1005static inline Quat<T> re(Quat<T> const &val)
1006{
1007    return ~val / norm(val);
1008}
1009
1010template<typename T>
1011static inline Quat<T> operator /(T x, Quat<T> const &y)
1012{
1013    return x * re(y);
1014}
1015
1016template<typename T>
1017static inline Quat<T> operator /(Quat<T> const &x, Quat<T> const &y)
1018{
1019    return x * re(y);
1020}
1021
1022template<typename T>
1023extern Quat<T> slerp(Quat<T> const &qa, Quat<T> const &qb, T f);
1024
1025/*
1026 * Common operators for all vector types, including quaternions
1027 */
1028
1029/*
1030 * vec +(vec, vec)   (also complex & quaternion)
1031 * vec -(vec, vec)   (also complex & quaternion)
1032 * vec *(vec, vec)
1033 * vec /(vec, vec)
1034 */
1035#define DECLARE_VECTOR_VECTOR_COERCE_OP(tname, op, tprefix, t1, t2, tf) \
1036    tprefix \
1037    inline tname<tf> operator op(tname<t1> const &a, tname<t2> const &b) \
1038    { \
1039        tname<tf> ret; \
1040        for (size_t n = 0; n < sizeof(a) / sizeof(t1); n++) \
1041            ret[n] = a[n] op b[n]; \
1042        return ret; \
1043    }
1044
1045/*
1046 * vec +=(vec, vec)   (also complex & quaternion)
1047 * vec -=(vec, vec)   (also complex & quaternion)
1048 * vec *=(vec, vec)
1049 * vec /=(vec, vec)
1050 */
1051#define DECLARE_VECTOR_VECTOR_OP(tname, op, tprefix, type) \
1052    tprefix \
1053    inline tname<type> operator op##=(tname<type> &a, tname<type> const &b) \
1054    { \
1055        return a = a op b; \
1056    }
1057
1058/*
1059 * vec min(vec, vec)     (also max, fmod)
1060 * vec min(vec, scalar)  (also max, fmod)
1061 * vec min(scalar, vec)  (also max, fmod)
1062 */
1063#define DECLARE_VECTOR_MINMAX_OP(tname, op, tprefix, type) \
1064    tprefix \
1065    inline tname<type> op(tname<type> const &a, tname<type> const &b) \
1066    { \
1067        using lol::op; \
1068        tname<type> ret; \
1069        for (size_t n = 0; n < sizeof(a) / sizeof(type); n++) \
1070            ret[n] = op(a[n], b[n]); \
1071        return ret; \
1072    } \
1073    \
1074    tprefix \
1075    inline tname<type> op(tname<type> const &a, type const &b) \
1076    { \
1077        using lol::op; \
1078        tname<type> ret; \
1079        for (size_t n = 0; n < sizeof(a) / sizeof(type); n++) \
1080            ret[n] = op(a[n], b); \
1081        return ret; \
1082    } \
1083    \
1084    tprefix \
1085    inline tname<type> op(type const &a, tname<type> const &b) \
1086    { \
1087        using lol::op; \
1088        tname<type> ret; \
1089        for (size_t n = 0; n < sizeof(b) / sizeof(type); n++) \
1090            ret[n] = op(a, b[n]); \
1091        return ret; \
1092    }
1093
1094/*
1095 * vec clamp(vec, vec, vec)
1096 * vec clamp(vec, vec, scalar)
1097 * vec clamp(vec, scalar, vec)
1098 * vec clamp(vec, scalar, scalar)
1099 */
1100#define DECLARE_VECTOR_CLAMP_OP(tname, tprefix, type) \
1101    tprefix \
1102    inline tname<type> clamp(tname<type> const &x, \
1103                             tname<type> const &a, tname<type> const &b) \
1104    { \
1105        return max(min(x, b), a); \
1106    } \
1107    \
1108    tprefix \
1109    inline tname<type> clamp(tname<type> const &x, \
1110                             type const &a, tname<type> const &b) \
1111    { \
1112        return max(min(x, b), a); \
1113    } \
1114    \
1115    tprefix \
1116    inline tname<type> clamp(tname<type> const &x, \
1117                             tname<type> const &a, type const &b) \
1118    { \
1119        return max(min(x, b), a); \
1120    } \
1121    \
1122    tprefix \
1123    inline tname<type> clamp(tname<type> const &x, \
1124                             type const &a, type const &b) \
1125    { \
1126        return max(min(x, b), a); \
1127    }
1128
1129/*
1130 * bool ==(vec, vec)   (also complex & quaternion)
1131 * bool !=(vec, vec)   (also complex & quaternion)
1132 * bool >=(vec, vec)
1133 * bool <=(vec, vec)
1134 * bool >(vec, vec)
1135 * bool <(vec, vec)
1136 */
1137#define DECLARE_VECTOR_VECTOR_BOOLOP(tname, op, op2, ret, tprefix, t1, t2) \
1138    tprefix \
1139    inline bool operator op(tname<t1> const &a, tname<t2> const &b) \
1140    { \
1141        for (size_t n = 0; n < sizeof(a) / sizeof(t1); n++) \
1142            if (!(a[n] op2 b[n])) \
1143                return !ret; \
1144        return ret; \
1145    }
1146
1147/*
1148 * vec *(vec, scalar)   (also complex & quaternion)
1149 * vec /(vec, scalar)   (also complex & quaternion)
1150 */
1151#define DECLARE_VECTOR_SCALAR_COERCE_OP(tname, op, tprefix, t1, t2, tf) \
1152    tprefix \
1153    inline tname<tf> operator op(tname<t1> const &a, t2 const &val) \
1154    { \
1155        tname<tf> ret; \
1156        for (size_t n = 0; n < sizeof(a) / sizeof(t1); n++) \
1157            ret[n] = a[n] op val; \
1158        return ret; \
1159    }
1160
1161/*
1162 * vec *(scalar, vec)   (also complex & quaternion)
1163 * vec /(scalar, vec)   (NOT for complex & quaternion!)
1164 */
1165#define DECLARE_SCALAR_VECTOR_COERCE_OP(tname, op, tprefix, t1, t2, tf) \
1166    tprefix \
1167    inline tname<tf> operator op(t1 const &val, tname<t2> const &a) \
1168    { \
1169        tname<tf> ret; \
1170        for (size_t n = 0; n < sizeof(a) / sizeof(t2); n++) \
1171            ret[n] = a[n] op val; \
1172        return ret; \
1173    }
1174
1175/*
1176 * vec *=(vec, scalar)   (also complex & quaternion)
1177 * vec /=(vec, scalar)   (also complex & quaternion)
1178 */
1179#define DECLARE_VECTOR_SCALAR_OP(tname, op, tprefix, type) \
1180    tprefix \
1181    inline tname<type> operator op##=(tname<type> &a, type const &val) \
1182    { \
1183        return a = a op val; \
1184    }
1185
1186#define DECLARE_UNARY_OPS(tname, tprefix, type) \
1187    tprefix \
1188    inline tname<type> operator -(tname<type> const &a) \
1189    { \
1190        tname<type> ret; \
1191        for (size_t n = 0; n < sizeof(a) / sizeof(type); n++) \
1192            ret[n] = -a[n]; \
1193        return ret; \
1194    } \
1195    \
1196    tprefix \
1197    inline type sqlength(tname<type> const &a) \
1198    { \
1199        type acc = 0; \
1200        for (size_t n = 0; n < sizeof(a) / sizeof(type); n++) \
1201            acc += a[n] * a[n]; \
1202        return acc; \
1203    } \
1204    \
1205    tprefix \
1206    inline double length(tname<type> const &a) \
1207    { \
1208        return sqrt((double)sqlength(a)); \
1209    } \
1210    \
1211    tprefix \
1212    inline tname<type> normalize(tname<type> const &a) \
1213    { \
1214        type norm = (type)length(a); \
1215        return norm ? a / norm : a * (type)0; \
1216    } \
1217    \
1218    tprefix \
1219    inline tname<type> abs(tname<type> const &a) \
1220    { \
1221        tname<type> ret; \
1222        for (size_t n = 0; n < sizeof(a) / sizeof(type); n++) \
1223            ret[n] = lol::abs(a[n]); \
1224        return ret; \
1225    }
1226
1227#define DECLARE_BINARY_NONVECTOR_COERCE_OPS(tname, tprefix, t1, t2, tf) \
1228    DECLARE_VECTOR_VECTOR_COERCE_OP(tname, -, tprefix, t1, t2, tf) \
1229    DECLARE_VECTOR_VECTOR_COERCE_OP(tname, +, tprefix, t1, t2, tf) \
1230    \
1231    DECLARE_VECTOR_SCALAR_COERCE_OP(tname, *, tprefix, t1, t2, tf) \
1232    DECLARE_VECTOR_SCALAR_COERCE_OP(tname, /, tprefix, t1, t2, tf) \
1233    DECLARE_SCALAR_VECTOR_COERCE_OP(tname, *, tprefix, t1, t2, tf) \
1234    \
1235    DECLARE_VECTOR_VECTOR_BOOLOP(tname, ==, ==, true, tprefix, t1, t2) \
1236    DECLARE_VECTOR_VECTOR_BOOLOP(tname, !=, ==, false, tprefix, t1, t2) \
1237    \
1238    tprefix \
1239    inline tf dot(tname<t1> const &a, tname<t2> const &b) \
1240    { \
1241        tf ret = 0; \
1242        for (size_t n = 0; n < sizeof(a) / sizeof(t1); n++) \
1243            ret += a[n] * b[n]; \
1244        return ret; \
1245    } \
1246    tprefix \
1247    inline tname<tf> lerp(tname<t1> const &a, tname<t2> const &b, tf x) \
1248    { \
1249        tname<tf> ret; \
1250        for (size_t n = 0; n < sizeof(a) / sizeof(t1); n++) \
1251            ret[n] = a[n] + (a[n] - b[n]) * x; \
1252        return ret; \
1253    }
1254
1255#define DECLARE_BINARY_VECTOR_COERCE_OPS(tname, tprefix, t1, t2, tf) \
1256    DECLARE_SCALAR_VECTOR_COERCE_OP(tname, /, tprefix, t1, t2, tf)
1257
1258#define DECLARE_VEC_3_COERCE_OPS(tname, tprefix, t1, t2, tf) \
1259    tprefix \
1260    inline tname<tf> cross(tname<t1> const &a, tname<t2> const &b) \
1261    { \
1262        return tname<tf>((tf)(a.y * b.z) - (tf)(a.z * b.y), \
1263                         (tf)(a.z * b.x) - (tf)(a.x * b.z), \
1264                         (tf)(a.x * b.y) - (tf)(a.y * b.x)); \
1265    }
1266
1267#define DECLARE_BINARY_NONVECTOR_OPS(tname, tprefix, type) \
1268    DECLARE_BINARY_NONVECTOR_COERCE_OPS(tname, tprefix, type, type, type) \
1269    \
1270    DECLARE_VECTOR_VECTOR_OP(tname, -, tprefix, type) \
1271    DECLARE_VECTOR_VECTOR_OP(tname, +, tprefix, type) \
1272    \
1273    DECLARE_VECTOR_SCALAR_OP(tname, *, tprefix, type) \
1274    DECLARE_VECTOR_SCALAR_OP(tname, /, tprefix, type)
1275
1276#define DECLARE_BINARY_VECTOR_OPS(tname, tprefix, type) \
1277    DECLARE_BINARY_VECTOR_COERCE_OPS(tname, tprefix, type, type, type) \
1278    \
1279    DECLARE_VECTOR_MINMAX_OP(tname, min, tprefix, type) \
1280    DECLARE_VECTOR_MINMAX_OP(tname, max, tprefix, type) \
1281    DECLARE_VECTOR_MINMAX_OP(tname, fmod, tprefix, type) \
1282    DECLARE_VECTOR_CLAMP_OP(tname, tprefix, type)
1283
1284#define DECLARE_VECTOR_COERCE_OPS(tname, tprefix, t1, t2, tf) \
1285    DECLARE_VECTOR_VECTOR_COERCE_OP(tname, *, tprefix, t1, t2, tf) \
1286    DECLARE_VECTOR_VECTOR_COERCE_OP(tname, /, tprefix, t1, t2, tf) \
1287    \
1288    DECLARE_VECTOR_VECTOR_BOOLOP(tname, <=, <=, true, tprefix, t1, t2) \
1289    DECLARE_VECTOR_VECTOR_BOOLOP(tname, >=, >=, true, tprefix, t1, t2) \
1290    DECLARE_VECTOR_VECTOR_BOOLOP(tname, <, <, true, tprefix, t1, t2) \
1291    DECLARE_VECTOR_VECTOR_BOOLOP(tname, >, >, true, tprefix, t1, t2)
1292
1293#define DECLARE_VECTOR_OPS(tname, tprefix, type) \
1294    DECLARE_VECTOR_COERCE_OPS(tname, static, type, type, type) \
1295    \
1296    DECLARE_VECTOR_VECTOR_OP(tname, *, tprefix, type) \
1297    DECLARE_VECTOR_VECTOR_OP(tname, /, tprefix, type)
1298
1299#define DECLARE_ALL_NONVECTOR_OPS(tname) \
1300    DECLARE_BINARY_NONVECTOR_OPS(tname, template<typename T> static, T) \
1301    DECLARE_UNARY_OPS(tname, template<typename T> static, T)
1302
1303#define DECLARE_ALL_VECTOR_OPS_INNER(tname, type) \
1304    DECLARE_BINARY_VECTOR_OPS(tname, static, type) \
1305    DECLARE_BINARY_NONVECTOR_OPS(tname, static, type) \
1306    DECLARE_UNARY_OPS(tname, static, type) \
1307    DECLARE_VECTOR_OPS(tname, static, type) \
1308
1309#define DECLARE_ALL_VECTOR_OPS(type) \
1310    namespace x##type \
1311    { \
1312        DECLARE_ALL_VECTOR_OPS_INNER(Vec2, type) \
1313        DECLARE_ALL_VECTOR_OPS_INNER(Vec3, type) \
1314        DECLARE_ALL_VECTOR_OPS_INNER(Vec4, type) \
1315        \
1316        DECLARE_VEC_3_COERCE_OPS(Vec3, static, type, type, type) \
1317    }
1318
1319#define DECLARE_VEC_ANY_COERCE_OPS(tname, tlow, thigh) \
1320    DECLARE_BINARY_NONVECTOR_COERCE_OPS(tname, static, tlow, thigh, thigh) \
1321    DECLARE_BINARY_NONVECTOR_COERCE_OPS(tname, static, thigh, tlow, thigh) \
1322    DECLARE_BINARY_VECTOR_COERCE_OPS(tname, static, tlow, thigh, thigh) \
1323    DECLARE_BINARY_VECTOR_COERCE_OPS(tname, static, thigh, tlow, thigh) \
1324    \
1325    DECLARE_VECTOR_COERCE_OPS(tname, static, tlow, thigh, thigh) \
1326    DECLARE_VECTOR_COERCE_OPS(tname, static, thigh, tlow, thigh)
1327
1328#define DECLARE_ALL_VECTOR_COERCE_OPS(tlow, thigh) \
1329    namespace x##tlow##thigh \
1330    { \
1331        DECLARE_VEC_ANY_COERCE_OPS(Vec2, tlow, thigh) \
1332        DECLARE_VEC_ANY_COERCE_OPS(Vec3, tlow, thigh) \
1333        DECLARE_VEC_ANY_COERCE_OPS(Vec4, tlow, thigh) \
1334    } \
1335    namespace y##tlow##thigh \
1336    { \
1337        DECLARE_VEC_3_COERCE_OPS(Vec3, static, tlow, thigh, thigh) \
1338        DECLARE_VEC_3_COERCE_OPS(Vec3, static, thigh, tlow, thigh) \
1339    } \
1340
1341DECLARE_ALL_NONVECTOR_OPS(Cmplx)
1342DECLARE_ALL_NONVECTOR_OPS(Quat)
1343
1344/* Disable warning about unary operator applied to unsigned type */
1345#if defined _MSC_VER
1346#   pragma warning(push)
1347#   pragma warning(disable: 4146)
1348#endif
1349
1350DECLARE_ALL_VECTOR_OPS(half)
1351DECLARE_ALL_VECTOR_OPS(float)
1352DECLARE_ALL_VECTOR_OPS(double)
1353DECLARE_ALL_VECTOR_OPS(ldouble)
1354DECLARE_ALL_VECTOR_OPS(int8_t)
1355DECLARE_ALL_VECTOR_OPS(uint8_t)
1356DECLARE_ALL_VECTOR_OPS(int16_t)
1357DECLARE_ALL_VECTOR_OPS(uint16_t)
1358DECLARE_ALL_VECTOR_OPS(int32_t)
1359DECLARE_ALL_VECTOR_OPS(uint32_t)
1360DECLARE_ALL_VECTOR_OPS(int64_t)
1361DECLARE_ALL_VECTOR_OPS(uint64_t)
1362
1363#if defined _MSC_VER
1364#   pragma warning(pop)
1365#endif
1366
1367/* Disable warnings in the >= > etc. operators about comparing signed and
1368 * unsigned. Ideally we would like to get these warnings only when the
1369 * inlined operators are actually used, but they seem to be triggered at
1370 * the code parsing step, so the compiler does not yet know whether they
1371 * will be used.
1372 * Also we do this for the whole block of declarations, because GCC prior
1373 * to 4.6.3 does not appear to support _Pragma() inside a macro. */
1374#if defined __GNUC__ && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6))
1375#   pragma GCC diagnostic push
1376#   pragma GCC diagnostic ignored "-Wsign-compare"
1377#elif defined _MSC_VER
1378#   pragma warning(push)
1379#   pragma warning(disable: 4018)
1380#endif
1381
1382/* Apply the same coercion rules as in the C++ standard. However, instead
1383 * of always promoting smaller types to int, we allow int8_t op int16_t to
1384 * return an int16_t. */
1385DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, uint8_t)
1386DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, int16_t)
1387DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, uint16_t)
1388DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, int32_t)
1389DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, uint32_t)
1390DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, int64_t)
1391DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, uint64_t)
1392DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, float)
1393DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, double)
1394DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, ldouble)
1395
1396DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, int16_t)
1397DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, uint16_t)
1398DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, int32_t)
1399DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, uint32_t)
1400DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, int64_t)
1401DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, uint64_t)
1402DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, float)
1403DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, double)
1404DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, ldouble)
1405
1406DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, uint16_t)
1407DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, int32_t)
1408DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, uint32_t)
1409DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, int64_t)
1410DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, uint64_t)
1411DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, float)
1412DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, double)
1413DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, ldouble)
1414
1415DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, int32_t)
1416DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, uint32_t)
1417DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, int64_t)
1418DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, uint64_t)
1419DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, float)
1420DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, double)
1421DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, ldouble)
1422
1423DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, uint32_t)
1424DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, int64_t)
1425DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, uint64_t)
1426DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, float)
1427DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, double)
1428DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, ldouble)
1429
1430DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, int64_t)
1431DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, uint64_t)
1432DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, float)
1433DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, double)
1434DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, ldouble)
1435
1436DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, uint64_t)
1437DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, float)
1438DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, double)
1439DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, ldouble)
1440
1441DECLARE_ALL_VECTOR_COERCE_OPS(uint64_t, float)
1442DECLARE_ALL_VECTOR_COERCE_OPS(uint64_t, double)
1443DECLARE_ALL_VECTOR_COERCE_OPS(uint64_t, ldouble)
1444
1445DECLARE_ALL_VECTOR_COERCE_OPS(float, double)
1446DECLARE_ALL_VECTOR_COERCE_OPS(float, ldouble)
1447
1448DECLARE_ALL_VECTOR_COERCE_OPS(double, ldouble)
1449
1450/* FIXME: vectors of "half" are deactivated for now, because they
1451 * induce extremely long compilation times (about 17 seconds per TU). */
1452
1453#if 0
1454/* All integer types are promoted to half; all floating point types
1455 * cause half to be promoted. */
1456DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, half)
1457DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, half)
1458DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, half)
1459DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, half)
1460DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, half)
1461DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, half)
1462DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, half)
1463DECLARE_ALL_VECTOR_COERCE_OPS(uint64_t, half)
1464
1465DECLARE_ALL_VECTOR_COERCE_OPS(half, float)
1466DECLARE_ALL_VECTOR_COERCE_OPS(half, double)
1467DECLARE_ALL_VECTOR_COERCE_OPS(half, ldouble)
1468#endif
1469
1470/* FIXME: vectors of "real" are deactivated for now, because we do
1471 * not implement all combinations of operators for these types yet. */
1472
1473#if 0
1474/* All types are promoted to real */
1475DECLARE_ALL_VECTOR_COERCE_OPS(int8_t, real)
1476DECLARE_ALL_VECTOR_COERCE_OPS(uint8_t, real)
1477DECLARE_ALL_VECTOR_COERCE_OPS(int16_t, real)
1478DECLARE_ALL_VECTOR_COERCE_OPS(uint16_t, real)
1479DECLARE_ALL_VECTOR_COERCE_OPS(int32_t, real)
1480DECLARE_ALL_VECTOR_COERCE_OPS(uint32_t, real)
1481DECLARE_ALL_VECTOR_COERCE_OPS(int64_t, real)
1482DECLARE_ALL_VECTOR_COERCE_OPS(uint64_t, real)
1483DECLARE_ALL_VECTOR_COERCE_OPS(half, real)
1484DECLARE_ALL_VECTOR_COERCE_OPS(float, real)
1485DECLARE_ALL_VECTOR_COERCE_OPS(double, real)
1486DECLARE_ALL_VECTOR_COERCE_OPS(ldouble, real)
1487#endif
1488
1489/* Activate all the namespaces that we created. Delaying this activation
1490 * reduces compilation times significantly. */
1491#define ACTIVATE_COERCE_NAMESPACES_INNER(tlow, thigh) \
1492    namespace x##tlow##thigh {} \
1493    namespace y##tlow##thigh {} \
1494    using namespace x##tlow##thigh; \
1495    using namespace y##tlow##thigh;
1496
1497#define ACTIVATE_COERCE_NAMESPACES(tlow) \
1498    namespace x##tlow {} \
1499    using namespace x##tlow; \
1500    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, int8_t) \
1501    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, uint8_t) \
1502    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, int16_t) \
1503    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, uint16_t) \
1504    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, int32_t) \
1505    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, uint32_t) \
1506    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, int64_t) \
1507    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, uint64_t) \
1508    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, half) \
1509    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, float) \
1510    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, double) \
1511    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, ldouble) \
1512    ACTIVATE_COERCE_NAMESPACES_INNER(tlow, real)
1513
1514ACTIVATE_COERCE_NAMESPACES(int8_t)
1515ACTIVATE_COERCE_NAMESPACES(uint8_t)
1516ACTIVATE_COERCE_NAMESPACES(int16_t)
1517ACTIVATE_COERCE_NAMESPACES(uint16_t)
1518ACTIVATE_COERCE_NAMESPACES(int32_t)
1519ACTIVATE_COERCE_NAMESPACES(uint32_t)
1520ACTIVATE_COERCE_NAMESPACES(int64_t)
1521ACTIVATE_COERCE_NAMESPACES(uint64_t)
1522ACTIVATE_COERCE_NAMESPACES(half)
1523ACTIVATE_COERCE_NAMESPACES(float)
1524ACTIVATE_COERCE_NAMESPACES(double)
1525ACTIVATE_COERCE_NAMESPACES(ldouble)
1526ACTIVATE_COERCE_NAMESPACES(real)
1527
1528#if defined __GNUC__ && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6))
1529#   pragma GCC diagnostic pop
1530#elif defined _MSC_VER
1531#   pragma warning(pop)
1532#endif
1533
1534#undef DECLARE_VECTOR_TYPEDEFS
1535#undef DECLARE_MEMBER_OPS
1536#undef DECLARE_VECTOR_VECTOR_OP
1537#undef DECLARE_VECTOR_VECTOR_BOOLOP
1538#undef DECLARE_VECTOR_SCALAR_OP
1539#undef DECLARE_BINARY_VECTOR_OPS
1540#undef DECLARE_BINARY_NONVECTOR_OPS
1541#undef DECLARE_UNARY_OPS
1542#undef DECLARE_ALL_NONVECTOR_OPS
1543#undef DECLARE_ALL_VECTOR_OPS_INNER
1544#undef DECLARE_ALL_VECTOR_OPS
1545
1546/*
1547 * Definition of additional functions requiring vector functions
1548 */
1549
1550template<typename T>
1551inline Quat<T> Quat<T>::operator *(Quat<T> const &val) const
1552{
1553    Quat<T> ret;
1554    Vec3<T> v1(x, y, z);
1555    Vec3<T> v2(val.x, val.y, val.z);
1556    Vec3<T> v3 = cross(v1, v2) + w * v2 + val.w * v1;
1557    return Quat<T>(w * val.w - dot(v1, v2), v3.x, v3.y, v3.z);
1558}
1559
1560/*
1561 * Magic vector swizzling (part 2/2)
1562 * Unfortunately these assignment operators cannot be used for now, because
1563 * we would also need to override the default copy assignment operator, and
1564 * in C++98 unions cannot contain such objects. This is why all the swizzling
1565 * magic objects are marked 'const' even those that could be lvalues.
1566 */
1567
1568template<typename T, int N>
1569inline Vec2<T> XVec2<T, N>::operator =(Vec2<T> const &that)
1570{
1571    for (int i = 0; i < 2; i++)
1572        *this[i] = that[i];
1573    return *this;
1574}
1575
1576template<typename T, int N>
1577inline Vec3<T> XVec3<T, N>::operator =(Vec3<T> const &that)
1578{
1579    for (int i = 0; i < 3; i++)
1580        *this[i] = that[i];
1581    return *this;
1582}
1583
1584template<typename T, int N>
1585inline Vec4<T> XVec4<T, N>::operator =(Vec4<T> const &that)
1586{
1587    for (int i = 0; i < 4; i++)
1588        *this[i] = that[i];
1589    return *this;
1590}
1591
1592/*
1593 * 2×2-element matrices
1594 */
1595
1596template <typename T> struct Mat2
1597{
1598    inline Mat2() {}
1599    inline Mat2(Vec2<T> V0, Vec2<T> V1)
1600      : v0(V0), v1(V1) {}
1601
1602    explicit inline Mat2(T val)
1603      : v0(val, (T)0),
1604        v1((T)0, val) {}
1605
1606    explicit inline Mat2(Mat4<T> const &mat)
1607      : v0(mat[0].xy),
1608        v1(mat[1].xy) {}
1609
1610    inline Vec2<T>& operator[](size_t n) { return (&v0)[n]; }
1611    inline Vec2<T> const& operator[](size_t n) const { return (&v0)[n]; }
1612
1613    /* Helpers for transformation matrices */
1614    static Mat2<T> rotate(T angle);
1615
1616    static inline Mat2<T> rotate(Mat2<T> mat, T angle)
1617    {
1618        return rotate(angle) * mat;
1619    }
1620
1621    void printf() const;
1622
1623    template<class U>
1624    friend std::ostream &operator<<(std::ostream &stream, Mat2<U> const &m);
1625
1626    inline Mat2<T> operator +(Mat2<T> const m) const
1627    {
1628        return Mat2<T>(v0 + m[0], v1 + m[1]);
1629    }
1630
1631    inline Mat2<T> operator +=(Mat2<T> const m)
1632    {
1633        return *this = *this + m;
1634    }
1635
1636    inline Mat2<T> operator -(Mat2<T> const m) const
1637    {
1638        return Mat2<T>(v0 - m[0], v1 - m[1]);
1639    }
1640
1641    inline Mat2<T> operator -=(Mat2<T> const m)
1642    {
1643        return *this = *this - m;
1644    }
1645
1646    inline Mat2<T> operator *(Mat2<T> const m) const
1647    {
1648        return Mat2<T>(*this * m[0], *this * m[1]);
1649    }
1650
1651    inline Mat2<T> operator *=(Mat2<T> const m)
1652    {
1653        return *this = *this * m;
1654    }
1655
1656    inline Vec2<T> operator *(Vec2<T> const m) const
1657    {
1658        Vec2<T> ret;
1659        for (int j = 0; j < 2; j++)
1660        {
1661            T tmp = 0;
1662            for (int k = 0; k < 2; k++)
1663                tmp += (*this)[k][j] * m[k];
1664            ret[j] = tmp;
1665        }
1666        return ret;
1667    }
1668
1669    Vec2<T> v0, v1;
1670};
1671
1672/*
1673 * 3×3-element matrices
1674 */
1675
1676template <typename T> struct Mat3
1677{
1678    inline Mat3() {}
1679    inline Mat3(Vec3<T> V0, Vec3<T> V1, Vec3<T> V2)
1680      : v0(V0), v1(V1), v2(V2) {}
1681
1682    explicit inline Mat3(T val)
1683      : v0(val, (T)0, (T)0),
1684        v1((T)0, val, (T)0),
1685        v2((T)0, (T)0, val) {}
1686
1687    explicit inline Mat3(Mat2<T> mat)
1688      : v0(mat[0], (T)0),
1689        v1(mat[1], (T)0),
1690        v2((T)0, (T)0, (T)0) {}
1691
1692    explicit inline Mat3(Mat2<T> mat, T val)
1693      : v0(Vec3<T>(mat[0], (T)0)),
1694        v1(Vec3<T>(mat[1], (T)0)),
1695        v2((T)0, (T)0, val) {}
1696
1697    explicit inline Mat3(Mat4<T> const &mat)
1698      : v0(mat[0].xyz),
1699        v1(mat[1].xyz),
1700        v2(mat[2].xyz) {}
1701
1702    explicit Mat3(Quat<T> const &q);
1703
1704    inline Vec3<T>& operator[](size_t n) { return (&v0)[n]; }
1705    inline Vec3<T> const& operator[](size_t n) const { return (&v0)[n]; }
1706
1707    /* Helpers for transformation matrices */
1708    static Mat3<T> scale(T x);
1709    static Mat3<T> scale(T x, T y, T z);
1710    static Mat3<T> scale(Vec3<T> v);
1711    static Mat3<T> rotate(T angle, T x, T y, T z);
1712    static Mat3<T> rotate(T angle, Vec3<T> v);
1713
1714    static Mat3<T> fromeuler_xyz(Vec3<T> const &v);
1715    static Mat3<T> fromeuler_xzy(Vec3<T> const &v);
1716    static Mat3<T> fromeuler_yxz(Vec3<T> const &v);
1717    static Mat3<T> fromeuler_yzx(Vec3<T> const &v);
1718    static Mat3<T> fromeuler_zxy(Vec3<T> const &v);
1719    static Mat3<T> fromeuler_zyx(Vec3<T> const &v);
1720    static Mat3<T> fromeuler_xyz(T phi, T theta, T psi);
1721    static Mat3<T> fromeuler_xzy(T phi, T theta, T psi);
1722    static Mat3<T> fromeuler_yxz(T phi, T theta, T psi);
1723    static Mat3<T> fromeuler_yzx(T phi, T theta, T psi);
1724    static Mat3<T> fromeuler_zxy(T phi, T theta, T psi);
1725    static Mat3<T> fromeuler_zyx(T phi, T theta, T psi);
1726
1727    static Mat3<T> fromeuler_xyx(Vec3<T> const &v);
1728    static Mat3<T> fromeuler_xzx(Vec3<T> const &v);
1729    static Mat3<T> fromeuler_yxy(Vec3<T> const &v);
1730    static Mat3<T> fromeuler_yzy(Vec3<T> const &v);
1731    static Mat3<T> fromeuler_zxz(Vec3<T> const &v);
1732    static Mat3<T> fromeuler_zyz(Vec3<T> const &v);
1733    static Mat3<T> fromeuler_xyx(T phi, T theta, T psi);
1734    static Mat3<T> fromeuler_xzx(T phi, T theta, T psi);
1735    static Mat3<T> fromeuler_yxy(T phi, T theta, T psi);
1736    static Mat3<T> fromeuler_yzy(T phi, T theta, T psi);
1737    static Mat3<T> fromeuler_zxz(T phi, T theta, T psi);
1738    static Mat3<T> fromeuler_zyz(T phi, T theta, T psi);
1739
1740    static inline Mat3<T> rotate(Mat3<T> mat, T angle, Vec3<T> v)
1741    {
1742        return rotate(angle, v) * mat;
1743    }
1744
1745    void printf() const;
1746
1747    template<class U>
1748    friend std::ostream &operator<<(std::ostream &stream, Mat3<U> const &m);
1749
1750    inline Mat3<T> operator +(Mat3<T> const m) const
1751    {
1752        return Mat3<T>(v0 + m[0], v1 + m[1], v2 + m[2]);
1753    }
1754
1755    inline Mat3<T> operator +=(Mat3<T> const m)
1756    {
1757        return *this = *this + m;
1758    }
1759
1760    inline Mat3<T> operator -(Mat3<T> const m) const
1761    {
1762        return Mat3<T>(v0 - m[0], v1 - m[1], v2 - m[2]);
1763    }
1764
1765    inline Mat3<T> operator -=(Mat3<T> const m)
1766    {
1767        return *this = *this - m;
1768    }
1769
1770    inline Mat3<T> operator *(Mat3<T> const m) const
1771    {
1772        return Mat3<T>(*this * m[0], *this * m[1], *this * m[2]);
1773    }
1774
1775    inline Mat3<T> operator *=(Mat3<T> const m)
1776    {
1777        return *this = *this * m;
1778    }
1779
1780    inline Vec3<T> operator *(Vec3<T> const m) const
1781    {
1782        Vec3<T> ret;
1783        for (int j = 0; j < 3; j++)
1784        {
1785            T tmp = 0;
1786            for (int k = 0; k < 3; k++)
1787                tmp += (*this)[k][j] * m[k];
1788            ret[j] = tmp;
1789        }
1790        return ret;
1791    }
1792
1793    Vec3<T> v0, v1, v2;
1794};
1795
1796/*
1797 * 4×4-element matrices
1798 */
1799
1800template <typename T> struct Mat4
1801{
1802    inline Mat4() {}
1803    inline Mat4(Vec4<T> V0, Vec4<T> V1, Vec4<T> V2, Vec4<T> V3)
1804      : v0(V0), v1(V1), v2(V2), v3(V3) {}
1805
1806    explicit inline Mat4(T val)
1807      : v0(val, (T)0, (T)0, (T)0),
1808        v1((T)0, val, (T)0, (T)0),
1809        v2((T)0, (T)0, val, (T)0),
1810        v3((T)0, (T)0, (T)0, val) {}
1811
1812    explicit inline Mat4(Mat2<T> mat)
1813      : v0(mat[0], (T)0, (T)0),
1814        v1(mat[1], (T)0, (T)0),
1815        v2((T)0, (T)0, (T)0, (T)0),
1816        v3((T)0, (T)0, (T)0, (T)0) {}
1817
1818    explicit inline Mat4(Mat2<T> mat, T val1, T val2)
1819      : v0(mat[0], (T)0, (T)0),
1820        v1(mat[1], (T)0, (T)0),
1821        v2((T)0, (T)0, val1, (T)0),
1822        v3((T)0, (T)0, (T)0, val2) {}
1823
1824    explicit inline Mat4(Mat3<T> mat)
1825      : v0(mat[0], (T)0),
1826        v1(mat[1], (T)0),
1827        v2(mat[2], (T)0),
1828        v3((T)0, (T)0, (T)0, (T)0) {}
1829
1830    explicit inline Mat4(Mat3<T> mat, T val)
1831      : v0(mat[0], (T)0),
1832        v1(mat[1], (T)0),
1833        v2(mat[2], (T)0),
1834        v3((T)0, (T)0, (T)0, val) {}
1835
1836    explicit Mat4(Quat<T> const &q);
1837
1838    inline Vec4<T>& operator[](size_t n) { return (&v0)[n]; }
1839    inline Vec4<T> const& operator[](size_t n) const { return (&v0)[n]; }
1840
1841    /* Helpers for transformation matrices */
1842    static Mat4<T> translate(T x, T y, T z);
1843    static Mat4<T> translate(Vec3<T> v);
1844
1845    static inline Mat4<T> scale(T x)
1846    {
1847        return Mat4<T>(Mat3<T>::scale(x), (T)1);
1848    }
1849
1850    static inline Mat4<T> scale(T x, T y, T z)
1851    {
1852        return Mat4<T>(Mat3<T>::scale(x, y, z), (T)1);
1853    }
1854
1855    static inline Mat4<T> scale(Vec3<T> v)
1856    {
1857        return Mat4<T>(Mat3<T>::scale(v), (T)1);
1858    }
1859
1860    static inline Mat4<T> translate(Mat4<T> const &mat, Vec3<T> v)
1861    {
1862        return translate(v) * mat;
1863    }
1864
1865    static inline Mat4<T> rotate(T angle, T x, T y, T z)
1866    {
1867        return Mat4<T>(Mat3<T>::rotate(angle, x, y, z), (T)1);
1868    }
1869
1870    static inline Mat4<T> rotate(T angle, Vec3<T> v)
1871    {
1872        return Mat4<T>(Mat3<T>::rotate(angle, v), (T)1);
1873    }
1874
1875    static inline Mat4<T> rotate(Mat4<T> &mat, T angle, Vec3<T> v)
1876    {
1877        return rotate(angle, v) * mat;
1878    }
1879
1880    static Mat4<T> fromeuler_xyz(Vec3<T> const &v);
1881    static Mat4<T> fromeuler_xzy(Vec3<T> const &v);
1882    static Mat4<T> fromeuler_yxz(Vec3<T> const &v);
1883    static Mat4<T> fromeuler_yzx(Vec3<T> const &v);
1884    static Mat4<T> fromeuler_zxy(Vec3<T> const &v);
1885    static Mat4<T> fromeuler_zyx(Vec3<T> const &v);
1886    static Mat4<T> fromeuler_xyz(T phi, T theta, T psi);
1887    static Mat4<T> fromeuler_xzy(T phi, T theta, T psi);
1888    static Mat4<T> fromeuler_yxz(T phi, T theta, T psi);
1889    static Mat4<T> fromeuler_yzx(T phi, T theta, T psi);
1890    static Mat4<T> fromeuler_zxy(T phi, T theta, T psi);
1891    static Mat4<T> fromeuler_zyx(T phi, T theta, T psi);
1892
1893    static Mat4<T> fromeuler_xyx(Vec3<T> const &v);
1894    static Mat4<T> fromeuler_xzx(Vec3<T> const &v);
1895    static Mat4<T> fromeuler_yxy(Vec3<T> const &v);
1896    static Mat4<T> fromeuler_yzy(Vec3<T> const &v);
1897    static Mat4<T> fromeuler_zxz(Vec3<T> const &v);
1898    static Mat4<T> fromeuler_zyz(Vec3<T> const &v);
1899    static Mat4<T> fromeuler_xyx(T phi, T theta, T psi);
1900    static Mat4<T> fromeuler_xzx(T phi, T theta, T psi);
1901    static Mat4<T> fromeuler_yxy(T phi, T theta, T psi);
1902    static Mat4<T> fromeuler_yzy(T phi, T theta, T psi);
1903    static Mat4<T> fromeuler_zxz(T phi, T theta, T psi);
1904    static Mat4<T> fromeuler_zyz(T phi, T theta, T psi);
1905
1906    /* Helpers for view matrices */
1907    static Mat4<T> lookat(Vec3<T> eye, Vec3<T> center, Vec3<T> up);
1908
1909    /* Helpers for projection matrices */
1910    static Mat4<T> ortho(T left, T right, T bottom, T top, T near, T far);
1911    static Mat4<T> ortho(T width, T height, T near, T far);
1912    static Mat4<T> frustum(T left, T right, T bottom, T top, T near, T far);
1913    static Mat4<T> perspective(T fov_y, T width, T height, T near, T far);
1914
1915    void printf() const;
1916
1917    template<class U>
1918    friend std::ostream &operator<<(std::ostream &stream, Mat4<U> const &m);
1919
1920    inline Mat4<T> operator +(Mat4<T> const &m) const
1921    {
1922        return Mat4<T>(v0 + m[0], v1 + m[1], v2 + m[2], v3 + m[3]);
1923    }
1924
1925    inline Mat4<T> operator +=(Mat4<T> const &m)
1926    {
1927        return *this = *this + m;
1928    }
1929
1930    inline Mat4<T> operator -(Mat4<T> const &m) const
1931    {
1932        return Mat4<T>(v0 - m[0], v1 - m[1], v2 - m[2], v3 - m[3]);
1933    }
1934
1935    inline Mat4<T> operator -=(Mat4<T> const &m)
1936    {
1937        return *this = *this - m;
1938    }
1939
1940    inline Mat4<T> operator *(Mat4<T> const &m) const
1941    {
1942        return Mat4<T>(*this * m[0], *this * m[1], *this * m[2], *this * m[3]);
1943    }
1944
1945    inline Mat4<T> operator *=(Mat4<T> const &m)
1946    {
1947        return *this = *this * m;
1948    }
1949
1950    inline Vec4<T> operator *(Vec4<T> const &m) const
1951    {
1952        Vec4<T> ret;
1953        for (int j = 0; j < 4; j++)
1954        {
1955            T tmp = 0;
1956            for (int k = 0; k < 4; k++)
1957                tmp += (*this)[k][j] * m[k];
1958            ret[j] = tmp;
1959        }
1960        return ret;
1961    }
1962
1963    Vec4<T> v0, v1, v2, v3;
1964};
1965
1966template<typename T> T determinant(Mat2<T> const &);
1967template<typename T> T determinant(Mat3<T> const &);
1968template<typename T> T determinant(Mat4<T> const &);
1969
1970template<typename T> Mat2<T> transpose(Mat2<T> const &);
1971template<typename T> Mat3<T> transpose(Mat3<T> const &);
1972template<typename T> Mat4<T> transpose(Mat4<T> const &);
1973
1974template<typename T> Mat2<T> inverse(Mat2<T> const &);
1975template<typename T> Mat3<T> inverse(Mat3<T> const &);
1976template<typename T> Mat4<T> inverse(Mat4<T> const &);
1977
1978/*
1979 * Arbitrarily-sized square matrices; for now this only supports
1980 * naive inversion and is used for the Remez inversion method.
1981 */
1982
1983template<int N, typename T> struct Mat
1984{
1985    inline Mat<N, T>() {}
1986
1987    Mat(T x)
1988    {
1989        for (int j = 0; j < N; j++)
1990            for (int i = 0; i < N; i++)
1991                if (i == j)
1992                    m[i][j] = x;
1993                else
1994                    m[i][j] = 0;
1995    }
1996
1997    /* Naive matrix inversion */
1998    Mat<N, T> inv() const
1999    {
2000        Mat a = *this, b((T)1);
2001
2002        /* Inversion method: iterate through all columns and make sure
2003         * all the terms are 1 on the diagonal and 0 everywhere else */
2004        for (int i = 0; i < N; i++)
2005        {
2006            /* If the expected coefficient is zero, add one of
2007             * the other lines. The first we meet will do. */
2008            if (!a.m[i][i])
2009            {
2010                for (int j = i + 1; j < N; j++)
2011                {
2012                    if (!a.m[i][j])
2013                        continue;
2014                    /* Add row j to row i */
2015                    for (int n = 0; n < N; n++)
2016                    {
2017                        a.m[n][i] += a.m[n][j];
2018                        b.m[n][i] += b.m[n][j];
2019                    }
2020                    break;
2021                }
2022            }
2023
2024            /* Now we know the diagonal term is non-zero. Get its inverse
2025             * and use that to nullify all other terms in the column */
2026            T x = (T)1 / a.m[i][i];
2027            for (int j = 0; j < N; j++)
2028            {
2029                if (j == i)
2030                    continue;
2031                T mul = x * a.m[i][j];
2032                for (int n = 0; n < N; n++)
2033                {
2034                    a.m[n][j] -= mul * a.m[n][i];
2035                    b.m[n][j] -= mul * b.m[n][i];
2036                }
2037            }
2038
2039            /* Finally, ensure the diagonal term is 1 */
2040            for (int n = 0; n < N; n++)
2041            {
2042                a.m[n][i] *= x;
2043                b.m[n][i] *= x;
2044            }
2045        }
2046
2047        return b;
2048    }
2049
2050    T m[N][N];
2051};
2052
2053} /* namespace lol */
2054
2055#endif // __LOL_MATH_VECTOR_H__
2056
Note: See TracBrowser for help on using the repository browser.