master
   1/**
   2 * This file has no copyright assigned and is placed in the Public Domain.
   3 * This file is part of the mingw-w64 runtime package.
   4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
   5 */
   6
   7/**
   8 * Normal users should include `windowsnumerics.h` instead of this header.
   9 * However, the cppwinrt headers set `_WINDOWS_NUMERICS_NAMESPACE_`,
  10 * `_WINDOWS_NUMERICS_BEGIN_NAMESPACE_` and `_WINDOWS_NUMERICS_END_NAMESPACE_`
  11 * to custom values and include `windowsnumerics.impl.h`. Therefore this shall
  12 * be considered a public header, and these macros are public API.
  13*/
  14
  15
  16#ifdef min
  17#  pragma push_macro("min")
  18#  undef min
  19#  define _WINDOWS_NUMERICS_IMPL_PUSHED_MIN_
  20#endif
  21
  22#ifdef max
  23#  pragma push_macro("max")
  24#  undef max
  25#  define _WINDOWS_NUMERICS_IMPL_PUSHED_MAX_
  26#endif
  27
  28#include <algorithm>
  29#include <cmath>
  30
  31#include "directxmath.h"
  32
  33
  34// === Internal macros ===
  35#define _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(_ty1, _op, _ty2) \
  36  inline _ty1 &operator _op ## =(_ty1 &val1, _ty2 val2) { \
  37    val1 = operator _op (val1, val2); \
  38    return val1; \
  39  }
  40
  41
  42// === Internal functions ===
  43_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
  44  namespace _impl {
  45
  46#if 0 && defined(__cpp_lib_clamp)
  47    using std::clamp;
  48#else
  49    constexpr const float &clamp(const float &val, const float &min, const float &max) {
  50      return val < min ? min : (val > max ? max : val);
  51    }
  52#endif
  53
  54#if 0 && defined(__cpp_lib_interpolate)
  55    using std::lerp;
  56#else
  57    constexpr float lerp(float val1, float val2, float amount) {
  58      // Don't do (val2 - val1) * amount + val1 as it has worse precision.
  59      return val2 * amount + val1 * (1.0f - amount);
  60    }
  61#endif
  62
  63  }
  64} _WINDOWS_NUMERICS_END_NAMESPACE_
  65
  66
  67// === Forward decls ===
  68_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
  69
  70  struct float2;
  71  struct float3;
  72  struct float4;
  73  struct float3x2;
  74  struct float4x4;
  75  struct plane;
  76  struct quaternion;
  77
  78} _WINDOWS_NUMERICS_END_NAMESPACE_
  79
  80
  81// === float2: Struct and function defs ===
  82_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
  83
  84  struct float2 {
  85    float2() = default;
  86    constexpr float2(float x, float y)
  87      : x(x), y(y)
  88    {}
  89    constexpr explicit float2(float val)
  90      : x(val), y(val)
  91    {}
  92
  93    static constexpr float2 zero() {
  94      return float2(0.0f);
  95    }
  96    static constexpr float2 one() {
  97      return float2(1.0f);
  98    }
  99    static constexpr float2 unit_x() {
 100      return { 1.0f, 0.0f };
 101    }
 102    static constexpr float2 unit_y() {
 103      return { 0.0f, 1.0f };
 104    }
 105
 106    float x;
 107    float y;
 108  };
 109
 110  // Forward decl functions
 111  inline float length(const float2 &val);
 112  inline float length_squared(const float2 &val);
 113  inline float distance(const float2 &val1, const float2 &val2);
 114  inline float distance_squared(const float2 &val1, const float2 &val2);
 115  inline float dot(const float2 &val1, const float2 &val2);
 116  inline float2 normalize(const float2 &val);
 117  inline float2 reflect(const float2 &vec, const float2 &norm);
 118  inline float2 min(const float2 &val1, const float2 &val2);
 119  inline float2 max(const float2 &val1, const float2 &val2);
 120  inline float2 clamp(const float2 &val, const float2 &min, const float2 &max);
 121  inline float2 lerp(const float2 &val1, const float2 &val2, float amount);
 122  inline float2 transform(const float2 &pos, const float3x2 &mat);
 123  inline float2 transform(const float2 &pos, const float4x4 &mat);
 124  inline float2 transform_normal(const float2 &norm, const float3x2 &mat);
 125  inline float2 transform_normal(const float2 &norm, const float4x4 &mat);
 126  inline float2 transform(const float2 &val, const quaternion &rot);
 127
 128  // Define operators
 129#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 130  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 131    return { val1.x _op val2.x, val1.y _op val2.y }; \
 132  }
 133  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float2, +)
 134  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float2, -)
 135  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float2, *)
 136  inline float2 operator*(const float2 &val1, float val2) {
 137    return { val1.x * val2, val1.y * val2 };
 138  }
 139  inline float2 operator*(float val1, const float2 &val2) {
 140    return operator*(val2, val1);
 141  }
 142  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float2, /)
 143  inline float2 operator/(const float2 &val1, float val2) {
 144    return operator*(val1, 1.0f / val2);
 145  }
 146  inline float2 operator-(const float2 &val) {
 147    return { -val.x, -val.y };
 148  }
 149#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 150  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, +, const float2 &)
 151  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, -, const float2 &)
 152  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, *, const float2 &)
 153  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, *, float)
 154  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, /, const float2 &)
 155  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float2, /, float)
 156  inline bool operator==(const float2 &val1, const float2 &val2) {
 157    return val1.x == val2.x && val1.y == val2.y;
 158  }
 159  inline bool operator!=(const float2 &val1, const float2 &val2) {
 160    return !operator==(val1, val2);
 161  }
 162
 163} _WINDOWS_NUMERICS_END_NAMESPACE_
 164
 165
 166// === float3: Struct and function defs ===
 167_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 168
 169  struct float3 {
 170    float3() = default;
 171    constexpr float3(float x, float y, float z)
 172      : x(x), y(y), z(z)
 173    {}
 174    constexpr float3(float2 val, float z)
 175      : x(val.x), y(val.y), z(z)
 176    {}
 177    constexpr explicit float3(float val)
 178      : x(val), y(val), z(val)
 179    {}
 180
 181    static constexpr float3 zero() {
 182      return float3(0.0);
 183    }
 184    static constexpr float3 one() {
 185      return float3(1.0);
 186    }
 187    static constexpr float3 unit_x() {
 188      return { 1.0f, 0.0f, 0.0f };
 189    }
 190    static constexpr float3 unit_y() {
 191      return { 0.0f, 1.0f, 0.0f };
 192    }
 193    static constexpr float3 unit_z() {
 194      return { 0.0f, 0.0f, 1.0f };
 195    }
 196
 197    float x;
 198    float y;
 199    float z;
 200  };
 201
 202  // Forward decl functions
 203  inline float length(const float3 &val);
 204  inline float length_squared(const float3 &val);
 205  inline float distance(const float3 &val1, const float3 &val2);
 206  inline float distance_squared(const float3 &val1, const float3 &val2);
 207  inline float dot(const float3 &val1, const float3 &val2);
 208  inline float3 normalize(const float3 &val);
 209  inline float3 cross(const float3 &val1, const float3 &val2);
 210  inline float3 reflect(const float3 &vec, const float3 &norm);
 211  inline float3 min(const float3 &val1, const float3 &val2);
 212  inline float3 max(const float3 &val1, const float3 &val2);
 213  inline float3 clamp(const float3 &val, const float3 &min, const float3 &max);
 214  inline float3 lerp(const float3 &val1, const float3 &val2, float amount);
 215  inline float3 transform(const float3 &pos, const float4x4 &mat);
 216  inline float3 transform_normal(const float3 &norm, const float4x4 &mat);
 217  inline float3 transform(const float3 &val, const quaternion &rot);
 218
 219  // Define operators
 220#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 221  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 222    return { val1.x _op val2.x, val1.y _op val2.y, val1.z _op val2.z }; \
 223  }
 224  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3, +)
 225  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3, -)
 226  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3, *)
 227  inline float3 operator*(const float3 &val1, float val2) {
 228    return { val1.x * val2, val1.y * val2, val1.z * val2 };
 229  }
 230  inline float3 operator*(float val1, const float3 &val2) {
 231    return operator*(val2, val1);
 232  }
 233  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3, /)
 234  inline float3 operator/(const float3 &val1, float val2) {
 235    return operator*(val1, 1.0f / val2);
 236  }
 237  inline float3 operator-(const float3 &val) {
 238    return { -val.x, -val.y, -val.z };
 239  }
 240#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 241  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, +, const float3 &)
 242  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, -, const float3 &)
 243  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, *, const float3 &)
 244  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, *, float)
 245  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, /, const float3 &)
 246  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3, /, float)
 247  inline bool operator==(const float3 &val1, const float3 &val2) {
 248    return val1.x == val2.x && val1.y == val2.y && val1.z == val2.z;
 249  }
 250  inline bool operator!=(const float3 &val1, const float3 &val2) {
 251    return !operator==(val1, val2);
 252  }
 253
 254} _WINDOWS_NUMERICS_END_NAMESPACE_
 255
 256
 257// === float4: Struct and function defs ===
 258_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 259
 260  struct float4 {
 261    float4() = default;
 262    constexpr float4(float x, float y, float z, float w)
 263      : x(x), y(y), z(z), w(w)
 264    {}
 265    constexpr float4(float2 val, float z, float w)
 266      : x(val.x), y(val.y), z(z), w(w)
 267    {}
 268    constexpr float4(float3 val, float w)
 269      : x(val.x), y(val.y), z(val.z), w(w)
 270    {}
 271    constexpr explicit float4(float val)
 272      : x(val), y(val), z(val), w(val)
 273    {}
 274
 275    static constexpr float4 zero() {
 276      return float4(0.0);
 277    }
 278    static constexpr float4 one() {
 279      return float4(1.0);
 280    }
 281    static constexpr float4 unit_x() {
 282      return { 1.0f, 0.0f, 0.0f, 0.0f };
 283    }
 284    static constexpr float4 unit_y() {
 285      return { 0.0f, 1.0f, 0.0f, 0.0f };
 286    }
 287    static constexpr float4 unit_z() {
 288      return { 0.0f, 0.0f, 1.0f, 0.0f };
 289    }
 290    static constexpr float4 unit_w() {
 291      return { 0.0f, 0.0f, 0.0f, 1.0f };
 292    }
 293
 294    float x;
 295    float y;
 296    float z;
 297    float w;
 298  };
 299
 300  // Forward decl functions
 301  inline float length(const float4 &val);
 302  inline float length_squared(const float4 &val);
 303  inline float distance(const float4 &val1, const float4 &val2);
 304  inline float distance_squared(const float4 &val1, const float4 &val2);
 305  inline float dot(const float4 &val1, const float4 &val2);
 306  inline float4 normalize(const float4 &val);
 307  inline float4 min(const float4 &val1, const float4 &val2);
 308  inline float4 max(const float4 &val1, const float4 &val2);
 309  inline float4 clamp(const float4 &val, const float4 &min, const float4 &max);
 310  inline float4 lerp(const float4 &val1, const float4 &val2, float amount);
 311  inline float4 transform(const float4 &pos, const float4x4 &mat);
 312  inline float4 transform4(const float3 &pos, const float4x4 &mat);
 313  inline float4 transform4(const float2 &pos, const float4x4 &mat);
 314  inline float4 transform(const float4 &val, const quaternion &rot);
 315  inline float4 transform4(const float3 &val, const quaternion &rot);
 316  inline float4 transform4(const float2 &val, const quaternion &rot);
 317
 318  // Define operators
 319#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 320  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 321    return { val1.x _op val2.x, val1.y _op val2.y, val1.z _op val2.z, val1.w _op val2.w }; \
 322  }
 323  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4, +)
 324  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4, -)
 325  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4, *)
 326  inline float4 operator*(const float4 &val1, float val2) {
 327    return { val1.x * val2, val1.y * val2, val1.z * val2, val1.w * val2 };
 328  }
 329  inline float4 operator*(float val1, const float4 &val2) {
 330    return operator*(val2, val1);
 331  }
 332  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4, /)
 333  inline float4 operator/(const float4 &val1, float val2) {
 334    return operator*(val1, 1.0f / val2);
 335  }
 336  inline float4 operator-(const float4 &val) {
 337    return { -val.x, -val.y, -val.z, -val.w };
 338  }
 339#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 340  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, +, const float4 &)
 341  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, -, const float4 &)
 342  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, *, const float4 &)
 343  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, *, float)
 344  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, /, const float4 &)
 345  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4, /, float)
 346  inline bool operator==(const float4 &val1, const float4 &val2) {
 347    return val1.x == val2.x && val1.y == val2.y && val1.z == val2.z && val2.w == val2.w;
 348  }
 349  inline bool operator!=(const float4 &val1, const float4 &val2) {
 350    return !operator==(val1, val2);
 351  }
 352
 353} _WINDOWS_NUMERICS_END_NAMESPACE_
 354
 355
 356// === float3x2: Struct and function defs ===
 357_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 358
 359  struct float3x2 {
 360    float3x2() = default;
 361    constexpr float3x2(
 362      float m11, float m12,
 363      float m21, float m22,
 364      float m31, float m32
 365    )
 366      : m11(m11), m12(m12)
 367      , m21(m21), m22(m22)
 368      , m31(m31), m32(m32)
 369    {}
 370
 371    static constexpr float3x2 identity() {
 372      return {
 373        1.0f, 0.0f,
 374        0.0f, 1.0f,
 375        0.0f, 0.0f
 376      };
 377    }
 378
 379    float m11; float m12;
 380    float m21; float m22;
 381    float m31; float m32;
 382  };
 383
 384  // Forward decl functions
 385  inline float3x2 make_float3x2_translation(const float2 &pos);
 386  inline float3x2 make_float3x2_translation(float xpos, float ypos);
 387  inline float3x2 make_float3x2_scale(float xscale, float yscale);
 388  inline float3x2 make_float3x2_scale(float xscale, float yscale, const float2 &center);
 389  inline float3x2 make_float3x2_scale(const float2 &xyscale);
 390  inline float3x2 make_float3x2_scale(const float2 &xyscale, const float2 &center);
 391  inline float3x2 make_float3x2_scale(float scale);
 392  inline float3x2 make_float3x2_scale(float scale, const float2 &center);
 393  inline float3x2 make_float3x2_skew(float xrad, float yrad);
 394  inline float3x2 make_float3x2_skew(float xrad, float yrad, const float2 &center);
 395  inline float3x2 make_float3x2_rotation(float rad);
 396  inline float3x2 make_float3x2_rotation(float rad, const float2 &center);
 397  inline bool is_identity(const float3x2 &val);
 398  inline float determinant(const float3x2 &val);
 399  inline float2 translation(const float3x2 &val);
 400  inline bool invert(const float3x2 &val, float3x2 *out);
 401  inline float3x2 lerp(const float3x2 &mat1, const float3x2 &mat2, float amount);
 402
 403  // Define operators
 404#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 405  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 406    return { \
 407      val1.m11 _op val2.m11, val1.m12 _op val2.m12, \
 408      val1.m21 _op val2.m21, val1.m22 _op val2.m22, \
 409      val1.m31 _op val2.m31, val1.m32 _op val2.m32, \
 410    }; \
 411  }
 412  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3x2, +)
 413  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float3x2, -)
 414  inline float3x2 operator*(const float3x2 &val1, const float3x2 &val2) {
 415    // 2D transformation matrix has an implied 3rd column with (0, 0, 1)
 416    const float3 v1r1(val1.m11, val1.m12, 0.0f);
 417    const float3 v1r2(val1.m21, val1.m22, 0.0f);
 418    const float3 v1r3(val1.m31, val1.m32, 1.0f);
 419    const float3 v2c1(val2.m11, val2.m21, val2.m31);
 420    const float3 v2c2(val2.m12, val2.m22, val2.m32);
 421    // const float3 v2c3(0.0f, 0.0f, 1.0f);
 422    return {
 423      dot(v1r1, v2c1), dot(v1r1, v2c2),
 424      dot(v1r2, v2c1), dot(v1r2, v2c2),
 425      dot(v1r3, v2c1), dot(v1r3, v2c2)
 426    };
 427  }
 428  inline float3x2 operator*(const float3x2 &val1, float val2) {
 429    return {
 430      val1.m11 * val2, val1.m12 * val2,
 431      val1.m21 * val2, val1.m22 * val2,
 432      val1.m31 * val2, val1.m32 * val2
 433    };
 434  }
 435  inline float3x2 operator-(const float3x2 &val) {
 436    return {
 437      -val.m11, -val.m12,
 438      -val.m21, -val.m22,
 439      -val.m31, -val.m32
 440    };
 441  }
 442#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 443  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3x2, +, const float3x2 &)
 444  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3x2, -, const float3x2 &)
 445  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3x2, *, const float3x2 &)
 446  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float3x2, *, float)
 447  inline bool operator==(const float3x2 &val1, const float3x2 &val2) {
 448    return
 449      val1.m11 == val2.m11 && val1.m12 == val2.m12 &&
 450      val1.m21 == val2.m21 && val1.m22 == val2.m22 &&
 451      val1.m31 == val2.m31 && val1.m32 == val2.m32;
 452  }
 453  inline bool operator!=(const float3x2 &val1, const float3x2 &val2) {
 454    return !operator==(val1, val2);
 455  }
 456
 457} _WINDOWS_NUMERICS_END_NAMESPACE_
 458
 459
 460// === float4x4: Struct and function defs ===
 461_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 462
 463  struct float4x4 {
 464    float4x4() = default;
 465    constexpr float4x4(
 466      float m11, float m12, float m13, float m14,
 467      float m21, float m22, float m23, float m24,
 468      float m31, float m32, float m33, float m34,
 469      float m41, float m42, float m43, float m44
 470    )
 471      : m11(m11), m12(m12), m13(m13), m14(m14)
 472      , m21(m21), m22(m22), m23(m23), m24(m24)
 473      , m31(m31), m32(m32), m33(m33), m34(m34)
 474      , m41(m41), m42(m42), m43(m43), m44(m44)
 475    {}
 476
 477    static constexpr float4x4 identity() {
 478      return {
 479        1.0f, 0.0f, 0.0f, 0.0f,
 480        0.0f, 1.0f, 0.0f, 0.0f,
 481        0.0f, 0.0f, 1.0f, 0.0f,
 482        0.0f, 0.0f, 0.0f, 1.0f
 483      };
 484    }
 485
 486    float m11; float m12; float m13; float m14;
 487    float m21; float m22; float m23; float m24;
 488    float m31; float m32; float m33; float m34;
 489    float m41; float m42; float m43; float m44;
 490  };
 491
 492  // Forward decl functions
 493  inline float4x4 make_float4x4_billboard(const float3 &objpos, const float3 &camerapos, const float3 &cameraup, const float3 &camerafwd);
 494  inline float4x4 make_float4x4_constrained_billboard(const float3 &objpos, const float3 &camerapos, const float3 &rotateaxis, const float3 &camerafwd, const float3 &objfwd);
 495  inline float4x4 make_float4x4_translation(const float3 &pos);
 496  inline float4x4 make_float4x4_translation(float xpos, float ypos, float zpos);
 497  inline float4x4 make_float4x4_scale(float xscale, float yscale, float zscale);
 498  inline float4x4 make_float4x4_scale(float xscale, float yscale, float zscale, const float3 &center);
 499  inline float4x4 make_float4x4_scale(const float3 &xyzscale);
 500  inline float4x4 make_float4x4_scale(const float3 &xyzscale, const float3 &center);
 501  inline float4x4 make_float4x4_scale(float scale);
 502  inline float4x4 make_float4x4_scale(float scale, const float3 &center);
 503  inline float4x4 make_float4x4_rotation_x(float rad);
 504  inline float4x4 make_float4x4_rotation_x(float rad, const float3 &center);
 505  inline float4x4 make_float4x4_rotation_y(float rad);
 506  inline float4x4 make_float4x4_rotation_y(float rad, const float3 &center);
 507  inline float4x4 make_float4x4_rotation_z(float rad);
 508  inline float4x4 make_float4x4_rotation_z(float rad, const float3 &center);
 509  inline float4x4 make_float4x4_from_axis_angle(const float3 &axis, float angle);
 510  inline float4x4 make_float4x4_perspective_field_of_view(float fov, float aspect, float nearplane, float farplane);
 511  inline float4x4 make_float4x4_perspective(float w, float h, float nearplane, float farplane);
 512  inline float4x4 make_float4x4_perspective_off_center(float left, float right, float bottom, float top, float nearplane, float farplane);
 513  inline float4x4 make_float4x4_orthographic(float w, float h, float znearplane, float zfarplane);
 514  inline float4x4 make_float4x4_orthographic_off_center(float left, float right, float bottom, float top, float znearplane, float zfarplane);
 515  inline float4x4 make_float4x4_look_at(const float3 &camerapos, const float3 &target, const float3 &cameraup);
 516  inline float4x4 make_float4x4_world(const float3 &pos, const float3 &fwd, const float3 &up);
 517  inline float4x4 make_float4x4_from_quaternion(const quaternion &quat);
 518  inline float4x4 make_float4x4_from_yaw_pitch_roll(float yaw, float pitch, float roll);
 519  inline float4x4 make_float4x4_shadow(const float3 &lightdir, const plane &plane);
 520  inline float4x4 make_float4x4_reflection(const plane &plane);
 521  inline bool is_identity(const float4x4 &val);
 522  inline float determinant(const float4x4 &val);
 523  inline float3 translation(const float4x4 &val);
 524  inline bool invert(const float4x4 &mat, float4x4 *out);
 525  inline bool decompose(const float4x4 &mat, float3 *out_scale, quaternion *out_rot, float3 *out_translate);
 526  inline float4x4 transform(const float4x4 &val, const quaternion &rot);
 527  inline float4x4 transpose(const float4x4 &val);
 528  inline float4x4 lerp(const float4x4 &val1, const float4x4 &val2, float amount);
 529
 530  // Define operators
 531#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 532  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 533    return { \
 534      val1.m11 _op val2.m11, val1.m12 _op val2.m12, val1.m13 _op val2.m13, val1.m14 _op val2.m14, \
 535      val1.m21 _op val2.m21, val1.m22 _op val2.m22, val1.m23 _op val2.m23, val1.m24 _op val2.m24, \
 536      val1.m31 _op val2.m31, val1.m32 _op val2.m32, val1.m33 _op val2.m33, val1.m34 _op val2.m34, \
 537      val1.m41 _op val2.m41, val1.m42 _op val2.m42, val1.m43 _op val2.m43, val1.m44 _op val2.m44, \
 538    }; \
 539  }
 540  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4x4, +)
 541  _WINDOWS_NUMERICS_IMPL_BINARY_OP(float4x4, -)
 542  inline float4x4 operator*(const float4x4 &val1, const float4x4 &val2) {
 543    const float4 v1r1(val1.m11, val1.m12, val1.m13, val1.m14);
 544    const float4 v1r2(val1.m21, val1.m22, val1.m23, val1.m24);
 545    const float4 v1r3(val1.m31, val1.m32, val1.m33, val1.m34);
 546    const float4 v1r4(val1.m41, val1.m42, val1.m43, val1.m44);
 547    const float4 v2c1(val2.m11, val2.m21, val2.m31, val2.m41);
 548    const float4 v2c2(val2.m12, val2.m22, val2.m32, val2.m42);
 549    const float4 v2c3(val2.m13, val2.m23, val2.m33, val2.m43);
 550    const float4 v2c4(val2.m14, val2.m24, val2.m34, val2.m44);
 551    return {
 552      dot(v1r1, v2c1), dot(v1r1, v2c2), dot(v1r1, v2c3), dot(v1r1, v2c4),
 553      dot(v1r2, v2c1), dot(v1r2, v2c2), dot(v1r2, v2c3), dot(v1r2, v2c4),
 554      dot(v1r3, v2c1), dot(v1r3, v2c2), dot(v1r3, v2c3), dot(v1r3, v2c4),
 555      dot(v1r4, v2c1), dot(v1r4, v2c2), dot(v1r4, v2c3), dot(v1r4, v2c4)
 556    };
 557  }
 558  inline float4x4 operator*(const float4x4 &val1, float val2) {
 559    return {
 560      val1.m11 * val2, val1.m12 * val2, val1.m13 * val2, val1.m14 * val2,
 561      val1.m21 * val2, val1.m22 * val2, val1.m23 * val2, val1.m24 * val2,
 562      val1.m31 * val2, val1.m32 * val2, val1.m33 * val2, val1.m34 * val2,
 563      val1.m41 * val2, val1.m42 * val2, val1.m43 * val2, val1.m44 * val2
 564    };
 565  }
 566  inline float4x4 operator-(const float4x4 &val) {
 567    return {
 568      -val.m11, -val.m12, -val.m13, -val.m14,
 569      -val.m21, -val.m22, -val.m23, -val.m24,
 570      -val.m31, -val.m32, -val.m33, -val.m34,
 571      -val.m41, -val.m42, -val.m43, -val.m44
 572    };
 573  }
 574#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 575  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4x4, +, const float4x4 &)
 576  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4x4, -, const float4x4 &)
 577  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4x4, *, const float4x4 &)
 578  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(float4x4, *, float)
 579  inline bool operator==(const float4x4 &val1, const float4x4 &val2) {
 580    return
 581      val1.m11 == val2.m11 && val1.m12 == val2.m12 && val1.m13 == val2.m13 && val1.m14 == val2.m14 &&
 582      val1.m21 == val2.m21 && val1.m22 == val2.m22 && val1.m23 == val2.m23 && val1.m24 == val2.m24 &&
 583      val1.m31 == val2.m31 && val1.m32 == val2.m32 && val1.m33 == val2.m33 && val1.m34 == val2.m34 &&
 584      val1.m41 == val2.m41 && val1.m42 == val2.m42 && val1.m43 == val2.m43 && val1.m44 == val2.m44;
 585  }
 586  inline bool operator!=(const float4x4 &val1, const float4x4 &val2) {
 587    return !operator==(val1, val2);
 588  }
 589
 590} _WINDOWS_NUMERICS_END_NAMESPACE_
 591
 592
 593// === plane: Struct and function defs ===
 594_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 595
 596  struct plane {
 597    plane() = default;
 598    constexpr plane(float x, float y, float z, float d)
 599      : normal(float3(x, y, z)), d(d)
 600    {}
 601    constexpr plane(float3 normal, float d)
 602      : normal(normal), d(d)
 603    {}
 604    constexpr explicit plane(float4 val)
 605      : normal(float3(val.x, val.y, val.z)), d(val.w)
 606    {}
 607
 608    float3 normal;
 609    float d;
 610  };
 611
 612  // Forward decl functions
 613  inline plane make_plane_from_vertices(const float3 &pt1, const float3 &pt2, const float3 &pt3);
 614  inline plane normalize(const plane &val);
 615  inline plane transform(const plane &plane, const float4x4 &mat);
 616  inline plane transform(const plane &plane, const quaternion &rot);
 617  inline float dot(const plane &plane, const float4 &val);
 618  inline float dot_coordinate(const plane &plane, const float3 &val);
 619  inline float dot_normal(const plane &plane, const float3 &val);
 620
 621  // Define operators
 622  inline bool operator==(const plane &val1, const plane &val2) {
 623    return val1.normal == val2.normal && val1.d == val2.d;
 624  }
 625  inline bool operator!=(const plane &val1, const plane &val2) {
 626    return !operator==(val1, val2);
 627  }
 628
 629} _WINDOWS_NUMERICS_END_NAMESPACE_
 630
 631
 632// === quaternion: Struct and function defs ===
 633_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 634
 635  struct quaternion {
 636    quaternion() = default;
 637    constexpr quaternion(float x, float y, float z, float w)
 638      : x(x), y(y), z(z), w(w)
 639    {}
 640    constexpr quaternion(float3 vecPart, float scalarPart)
 641      : x(vecPart.x), y(vecPart.y), z(vecPart.z), w(scalarPart)
 642    {}
 643
 644    static constexpr quaternion identity() {
 645      return { 0.0f, 0.0f, 0.0f, 1.0f };
 646    }
 647
 648    float x;
 649    float y;
 650    float z;
 651    float w;
 652  };
 653
 654  // Forward decl functions
 655  inline quaternion make_quaternion_from_axis_angle(const float3 &axis, float angle);
 656  inline quaternion make_quaternion_from_yaw_pitch_roll(float yaw, float pitch, float roll);
 657  inline quaternion make_quaternion_from_rotation_matrix(const float4x4 &mat);
 658  inline bool is_identity(const quaternion &val);
 659  inline float length(const quaternion &val);
 660  inline float length_squared(const quaternion &val);
 661  inline float dot(const quaternion &val1, const quaternion &val2);
 662  inline quaternion normalize(const quaternion &val);
 663  inline quaternion conjugate(const quaternion &val);
 664  inline quaternion inverse(const quaternion &val);
 665  inline quaternion slerp(const quaternion &val1, const quaternion &val2, float amount);
 666  inline quaternion lerp(const quaternion &val1, const quaternion &val2, float amount);
 667  inline quaternion concatenate(const quaternion &val1, const quaternion &val2);
 668
 669  // Define operators
 670#define _WINDOWS_NUMERICS_IMPL_BINARY_OP(_ty, _op) \
 671  inline _ty operator _op(const _ty &val1, const _ty &val2) { \
 672    return { val1.x _op val2.x, val1.y _op val2.y, val1.z _op val2.z, val1.w _op val2.w }; \
 673  }
 674  _WINDOWS_NUMERICS_IMPL_BINARY_OP(quaternion, +)
 675  _WINDOWS_NUMERICS_IMPL_BINARY_OP(quaternion, -)
 676  inline quaternion operator*(const quaternion &val1, const quaternion &val2) {
 677    return {
 678      val1.w * val2.x + val1.x * val2.w + val1.y * val2.z - val1.z * val2.y,
 679      val1.w * val2.y - val1.x * val2.z + val1.y * val2.w + val1.z * val2.x,
 680      val1.w * val2.z + val1.x * val2.y - val1.y * val2.x + val1.z * val2.w,
 681      val1.w * val2.w - val1.x * val2.x - val1.y * val2.y - val1.z * val2.z
 682    }; }
 683  inline quaternion operator*(const quaternion &val1, float val2) {
 684    return { val1.x * val2, val1.y * val2, val1.z * val2, val1.w * val2 };
 685  }
 686  inline quaternion operator/(const quaternion &val1, const quaternion &val2) {
 687    return operator*(val1, inverse(val2));
 688  }
 689  inline quaternion operator-(const quaternion &val) {
 690    return { -val.x, -val.y, -val.z, -val.w };
 691  }
 692#undef _WINDOWS_NUMERICS_IMPL_BINARY_OP
 693  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(quaternion, +, const quaternion &)
 694  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(quaternion, -, const quaternion &)
 695  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(quaternion, *, const quaternion &)
 696  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(quaternion, *, float)
 697  _WINDOWS_NUMERICS_IMPL_ASSIGN_OP(quaternion, /, const quaternion &)
 698  inline bool operator==(const quaternion &val1, const quaternion &val2) {
 699    return val1.x == val2.x && val1.y == val2.y && val1.z == val2.z && val2.w == val2.w;
 700  }
 701  inline bool operator!=(const quaternion &val1, const quaternion &val2) {
 702    return !operator==(val1, val2);
 703  }
 704
 705} _WINDOWS_NUMERICS_END_NAMESPACE_
 706
 707
 708// === Function definitions ===
 709
 710// Define float2 functions
 711_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 712
 713  inline float length(const float2 &val) {
 714    return ::std::sqrt(length_squared(val));
 715  }
 716  inline float length_squared(const float2 &val) {
 717    return val.x * val.x + val.y * val.y;
 718  }
 719  inline float distance(const float2 &val1, const float2 &val2) {
 720    return length(val2 - val1);
 721  }
 722  inline float distance_squared(const float2 &val1, const float2 &val2) {
 723    return length_squared(val2 - val1);
 724  }
 725  inline float dot(const float2 &val1, const float2 &val2) {
 726    return val1.x * val2.x + val1.y * val2.y;
 727  }
 728  inline float2 normalize(const float2 &val) {
 729    return val / length(val);
 730  }
 731  inline float2 reflect(const float2 &vec, const float2 &norm) {
 732    // norm is assumed to be normalized.
 733    return vec - 2.0f * dot(vec, norm) * norm;
 734  }
 735  inline float2 min(const float2 &val1, const float2 &val2) {
 736    return { ::std::min(val1.x, val2.x), ::std::min(val1.y, val2.y) };
 737  }
 738  inline float2 max(const float2 &val1, const float2 &val2) {
 739    return { ::std::max(val1.x, val2.x), ::std::max(val1.y, val2.y) };
 740  }
 741  inline float2 clamp(const float2 &val, const float2 &min, const float2 &max) {
 742    return { _impl::clamp(val.x, min.x, max.x), _impl::clamp(val.y, min.y, max.y) };
 743  }
 744  inline float2 lerp(const float2 &val1, const float2 &val2, float amount) {
 745    return { _impl::lerp(val1.x, val2.x, amount), _impl::lerp(val1.y, val2.y, amount) };
 746  }
 747  // TODO: impl
 748  inline float2 transform(const float2 &pos, const float3x2 &mat);
 749  inline float2 transform(const float2 &pos, const float4x4 &mat);
 750  inline float2 transform_normal(const float2 &norm, const float3x2 &mat);
 751  inline float2 transform_normal(const float2 &norm, const float4x4 &mat);
 752  inline float2 transform(const float2 &val, const quaternion &rot) {
 753    // See comments in the float3 transform function.
 754    quaternion result = rot * quaternion(float3(val.x, val.y, 0.0f), 0.0f) * inverse(rot);
 755    return { result.x, result.y };
 756  }
 757
 758} _WINDOWS_NUMERICS_END_NAMESPACE_
 759
 760
 761// Define float3 functions
 762_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 763
 764  inline float length(const float3 &val) {
 765    return ::std::sqrt(length_squared(val));
 766  }
 767  inline float length_squared(const float3 &val) {
 768    return val.x * val.x + val.y * val.y + val.z * val.z;
 769  }
 770  inline float distance(const float3 &val1, const float3 &val2) {
 771    return length(val2 - val1);
 772  }
 773  inline float distance_squared(const float3 &val1, const float3 &val2) {
 774    return length_squared(val2 - val1);
 775  }
 776  inline float dot(const float3 &val1, const float3 &val2) {
 777    return val1.x * val2.x + val1.y * val2.y + val1.z * val2.z;
 778  }
 779  inline float3 normalize(const float3 &val) {
 780    return val / length(val);
 781  }
 782  inline float3 cross(const float3 &val1, const float3 &val2) {
 783    return {
 784      val1.y * val2.z - val2.y * val1.z,
 785      val1.z * val2.x - val2.z * val1.x,
 786      val1.x * val2.y - val2.x * val1.y
 787    };
 788  }
 789  inline float3 reflect(const float3 &vec, const float3 &norm) {
 790    // norm is assumed to be normalized.
 791    return vec - 2.0f * dot(vec, norm) * norm;
 792  }
 793  inline float3 min(const float3 &val1, const float3 &val2) {
 794    return { ::std::min(val1.x, val2.x), ::std::min(val1.y, val2.y), ::std::min(val1.z, val2.z) };
 795  }
 796  inline float3 max(const float3 &val1, const float3 &val2) {
 797    return { ::std::max(val1.x, val2.x), ::std::max(val1.y, val2.y), ::std::max(val1.z, val2.z) };
 798  }
 799  inline float3 clamp(const float3 &val, const float3 &min, const float3 &max) {
 800    return { _impl::clamp(val.x, min.x, max.x), _impl::clamp(val.y, min.y, max.y), _impl::clamp(val.z, min.z, max.z) };
 801  }
 802  inline float3 lerp(const float3 &val1, const float3 &val2, float amount) {
 803    return { _impl::lerp(val1.x, val2.x, amount), _impl::lerp(val1.y, val2.y, amount), _impl::lerp(val1.z, val2.z, amount) };
 804  }
 805  // TODO: impl
 806  inline float3 transform(const float3 &pos, const float4x4 &mat);
 807  inline float3 transform_normal(const float3 &norm, const float4x4 &mat);
 808  inline float3 transform(const float3 &val, const quaternion &rot) {
 809    // https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Using_quaternions_as_rotations
 810    // If assuming rot is a unit quaternion, this could use
 811    // conjugate() instead of inverse() too.
 812    // This can be expressed as a matrix operation too, with
 813    // https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix
 814    // (see make_float4x4_from_quaternion).
 815    quaternion result = rot * quaternion(val, 0.0f) * inverse(rot);
 816    return { result.x, result.y, result.z };
 817  }
 818
 819} _WINDOWS_NUMERICS_END_NAMESPACE_
 820
 821
 822// Define float4 functions
 823_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 824
 825  inline float length(const float4 &val) {
 826    return ::std::sqrt(length_squared(val));
 827  }
 828  inline float length_squared(const float4 &val) {
 829    return val.x * val.x + val.y * val.y + val.z * val.z + val.w * val.w;
 830  }
 831  inline float distance(const float4 &val1, const float4 &val2) {
 832    return length(val2 - val1);
 833  }
 834  inline float distance_squared(const float4 &val1, const float4 &val2) {
 835    return length_squared(val2 - val1);
 836  }
 837  inline float dot(const float4 &val1, const float4 &val2) {
 838    return val1.x * val2.x + val1.y * val2.y + val1.z * val2.z + val1.w * val2.w;
 839  }
 840  inline float4 normalize(const float4 &val) {
 841    return val / length(val);
 842  }
 843  inline float4 min(const float4 &val1, const float4 &val2) {
 844    return {
 845      ::std::min(val1.x, val2.x),
 846      ::std::min(val1.y, val2.y),
 847      ::std::min(val1.z, val2.z),
 848      ::std::min(val1.w, val2.w)
 849    };
 850  }
 851  inline float4 max(const float4 &val1, const float4 &val2) {
 852    return {
 853      ::std::max(val1.x, val2.x),
 854      ::std::max(val1.y, val2.y),
 855      ::std::max(val1.z, val2.z),
 856      ::std::max(val1.w, val2.w)
 857    };
 858  }
 859  inline float4 clamp(const float4 &val, const float4 &min, const float4 &max) {
 860    return {
 861      _impl::clamp(val.x, min.x, max.x),
 862      _impl::clamp(val.y, min.y, max.y),
 863      _impl::clamp(val.z, min.z, max.z),
 864      _impl::clamp(val.w, min.w, max.w)
 865    };
 866  }
 867  inline float4 lerp(const float4 &val1, const float4 &val2, float amount) {
 868    return {
 869      _impl::lerp(val1.x, val2.x, amount),
 870      _impl::lerp(val1.y, val2.y, amount),
 871      _impl::lerp(val1.z, val2.z, amount),
 872      _impl::lerp(val1.w, val2.w, amount)
 873    };
 874  }
 875  // TODO: impl
 876  inline float4 transform(const float4 &pos, const float4x4 &mat);
 877  inline float4 transform4(const float3 &pos, const float4x4 &mat);
 878  inline float4 transform4(const float2 &pos, const float4x4 &mat);
 879  inline float4 transform(const float4 &val, const quaternion &rot) {
 880    // See comments in the float3 transform function.
 881    quaternion result = rot * quaternion(float3(val.x, val.y, val.z), 0.0f) * inverse(rot);
 882    return { result.x, result.y, result.z, val.w };
 883  }
 884  inline float4 transform4(const float3 &val, const quaternion &rot) {
 885    quaternion result = rot * quaternion(val, 0.0f) * inverse(rot);
 886    return { result.x, result.y, result.z, 1.0f };
 887  }
 888  inline float4 transform4(const float2 &val, const quaternion &rot) {
 889    quaternion result = rot * quaternion(float3(val.x, val.y, 0.0f), 0.0f) * inverse(rot);
 890    return { result.x, result.y, result.z, 1.0f };
 891  }
 892
 893} _WINDOWS_NUMERICS_END_NAMESPACE_
 894
 895
 896// Define float3x2 functions
 897_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 898
 899  // TODO: impl
 900  inline float3x2 make_float3x2_translation(const float2 &pos);
 901  inline float3x2 make_float3x2_translation(float xpos, float ypos);
 902  inline float3x2 make_float3x2_scale(float xscale, float yscale);
 903  inline float3x2 make_float3x2_scale(float xscale, float yscale, const float2 &center);
 904  inline float3x2 make_float3x2_scale(const float2 &xyscale);
 905  inline float3x2 make_float3x2_scale(const float2 &xyscale, const float2 &center);
 906  inline float3x2 make_float3x2_scale(float scale);
 907  inline float3x2 make_float3x2_scale(float scale, const float2 &center);
 908  inline float3x2 make_float3x2_skew(float xrad, float yrad);
 909  inline float3x2 make_float3x2_skew(float xrad, float yrad, const float2 &center);
 910  inline float3x2 make_float3x2_rotation(float rad);
 911  inline float3x2 make_float3x2_rotation(float rad, const float2 &center);
 912  inline bool is_identity(const float3x2 &val) {
 913    return val == float3x2::identity();
 914  }
 915  inline float determinant(const float3x2 &val) {
 916    // 2D transformation matrix has an implied 3rd column with (0, 0, 1)
 917    // det = m11 * m22 * m33 + m12 * m23 * m31 + m13 * m21 * m32
 918    //     - m31 * m22 * m13 - m21 * m12 * m33 - m11 * m32 * m23;
 919    return val.m11 * val.m22 - val.m21 * val.m12;
 920  }
 921  inline float2 translation(const float3x2 &val) {
 922    return { val.m31, val.m32 };
 923  }
 924  inline bool invert(const float3x2 &val, float3x2 *out);
 925  inline float3x2 lerp(const float3x2 &mat1, const float3x2 &mat2, float amount);
 926
 927} _WINDOWS_NUMERICS_END_NAMESPACE_
 928
 929
 930// Define float4x4 functions
 931_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
 932
 933  // TODO: impl
 934  inline float4x4 make_float4x4_billboard(const float3 &objpos, const float3 &camerapos, const float3 &cameraup, const float3 &camerafwd);
 935  inline float4x4 make_float4x4_constrained_billboard(const float3 &objpos, const float3 &camerapos, const float3 &rotateaxis, const float3 &camerafwd, const float3 &objfwd);
 936  inline float4x4 make_float4x4_translation(const float3 &pos) {
 937    return {
 938      1.0f,  0.0f,  0.0f,  0.0f,
 939      0.0f,  1.0f,  0.0f,  0.0f,
 940      0.0f,  0.0f,  1.0f,  0.0f,
 941      pos.x, pos.y, pos.z, 1.0f
 942    };
 943  }
 944  inline float4x4 make_float4x4_translation(float xpos, float ypos, float zpos);
 945  inline float4x4 make_float4x4_scale(float xscale, float yscale, float zscale);
 946  inline float4x4 make_float4x4_scale(float xscale, float yscale, float zscale, const float3 &center);
 947  inline float4x4 make_float4x4_scale(const float3 &xyzscale) {
 948    return {
 949      xyzscale.x, 0.0f,       0.0f,       0.0f,
 950      0.0f,       xyzscale.y, 0.0f,       0.0f,
 951      0.0f,       0.0f,       xyzscale.z, 0.0f,
 952      0.0f,       0.0f,       0.0f,       1.0f
 953    };
 954  }
 955  inline float4x4 make_float4x4_scale(const float3 &xyzscale, const float3 &center);
 956  inline float4x4 make_float4x4_scale(float scale);
 957  inline float4x4 make_float4x4_scale(float scale, const float3 &center);
 958  inline float4x4 make_float4x4_rotation_x(float rad);
 959  inline float4x4 make_float4x4_rotation_x(float rad, const float3 &center);
 960  inline float4x4 make_float4x4_rotation_y(float rad);
 961  inline float4x4 make_float4x4_rotation_y(float rad, const float3 &center);
 962  inline float4x4 make_float4x4_rotation_z(float rad);
 963  inline float4x4 make_float4x4_rotation_z(float rad, const float3 &center);
 964  inline float4x4 make_float4x4_from_axis_angle(const float3 &axis, float angle);
 965  inline float4x4 make_float4x4_perspective_field_of_view(float fov, float aspect, float nearplane, float farplane);
 966  inline float4x4 make_float4x4_perspective(float w, float h, float nearplane, float farplane);
 967  inline float4x4 make_float4x4_perspective_off_center(float left, float right, float bottom, float top, float nearplane, float farplane);
 968  inline float4x4 make_float4x4_orthographic(float w, float h, float znearplane, float zfarplane);
 969  inline float4x4 make_float4x4_orthographic_off_center(float left, float right, float bottom, float top, float znearplane, float zfarplane);
 970  inline float4x4 make_float4x4_look_at(const float3 &camerapos, const float3 &target, const float3 &cameraup);
 971  inline float4x4 make_float4x4_world(const float3 &pos, const float3 &fwd, const float3 &up);
 972  inline float4x4 make_float4x4_from_quaternion(const quaternion &quat) {
 973    // https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
 974    float xx = quat.x * quat.x;
 975    float yy = quat.y * quat.y;
 976    float zz = quat.z * quat.z;
 977    float xy = quat.x * quat.y;
 978    float xz = quat.x * quat.z;
 979    float xw = quat.x * quat.w;
 980    float yz = quat.y * quat.z;
 981    float yw = quat.y * quat.w;
 982    float zw = quat.z * quat.w;
 983    return {
 984      1.0f - 2.0f*yy - 2.0f*zz, 2.0f*xy + 2.0f*zw,        2.0f*xz - 2.0f*yw,        0.0f,
 985      2.0f*xy - 2.0f*zw,        1.0f - 2.0f*xx - 2.0f*zz, 2.0f*yz + 2.0f*xw,        0.0f,
 986      2.0f*xz + 2.0f*yw,        2.0f*yz - 2.0f*xw,        1.0f - 2.0f*xx - 2.0f*yy, 0.0f,
 987      0.0f,                     0.0f,                     0.0f,                     1.0f
 988    };
 989  }
 990  inline float4x4 make_float4x4_from_yaw_pitch_roll(float yaw, float pitch, float roll);
 991  inline float4x4 make_float4x4_shadow(const float3 &lightdir, const plane &plane);
 992  inline float4x4 make_float4x4_reflection(const plane &plane);
 993  inline bool is_identity(const float4x4 &val) {
 994    return val == float4x4::identity();
 995  }
 996  inline float determinant(const float4x4 &val) {
 997    const float det_33_44 = (val.m33 * val.m44 - val.m34 * val.m43);
 998    const float det_32_44 = (val.m32 * val.m44 - val.m34 * val.m42);
 999    const float det_32_43 = (val.m32 * val.m43 - val.m33 * val.m42);
1000    const float det_31_44 = (val.m31 * val.m44 - val.m34 * val.m41);
1001    const float det_31_43 = (val.m31 * val.m43 - val.m33 * val.m41);
1002    const float det_31_42 = (val.m31 * val.m42 - val.m32 * val.m41);
1003    return val.m11 * (val.m22 * det_33_44 - val.m23 * det_32_44 + val.m24 * det_32_43)
1004      - val.m12 * (val.m21 * det_33_44 - val.m23 * det_31_44 + val.m24 * det_31_43)
1005      + val.m13 * (val.m21 * det_32_44 - val.m22 * det_31_44 + val.m24 * det_31_42)
1006      - val.m14 * (val.m21 * det_32_43 - val.m22 * det_31_43 + val.m23 * det_31_42);
1007  }
1008  inline float3 translation(const float4x4 &val) {
1009    return { val.m41, val.m42, val.m43 };
1010  }
1011  inline bool invert(const float4x4 &mat, float4x4 *out);
1012  inline bool decompose(const float4x4 &mat, float3 *out_scale, quaternion *out_rot, float3 *out_translate);
1013  inline float4x4 transform(const float4x4 &val, const quaternion &rot) {
1014    return val * make_float4x4_from_quaternion(rot);
1015  }
1016  inline float4x4 transpose(const float4x4 &val) {
1017    return  {
1018      val.m11, val.m21, val.m31, val.m41,
1019      val.m12, val.m22, val.m32, val.m42,
1020      val.m13, val.m23, val.m33, val.m43,
1021      val.m14, val.m24, val.m34, val.m44
1022    };
1023  }
1024  inline float4x4 lerp(const float4x4 &val1, const float4x4 &val2, float amount);
1025
1026} _WINDOWS_NUMERICS_END_NAMESPACE_
1027
1028
1029// Define plane functions
1030_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
1031
1032  // TODO: impl
1033  inline plane make_plane_from_vertices(const float3 &pt1, const float3 &pt2, const float3 &pt3);
1034  inline plane normalize(const plane &val) {
1035    const float invlen = 1.0f / length(val.normal);
1036    return { val.normal * invlen, val.d * invlen };
1037  }
1038  inline plane transform(const plane &plane, const float4x4 &mat);
1039  inline plane transform(const plane &plane, const quaternion &rot) {
1040    quaternion result = rot * quaternion(plane.normal, 0.0f) * inverse(rot);
1041    return { float3(result.x, result.y, result.z), plane.d };
1042  }
1043  inline float dot(const plane &plane, const float4 &val);
1044  inline float dot_coordinate(const plane &plane, const float3 &val);
1045  inline float dot_normal(const plane &plane, const float3 &val);
1046
1047} _WINDOWS_NUMERICS_END_NAMESPACE_
1048
1049
1050// Define quaternion functions
1051_WINDOWS_NUMERICS_BEGIN_NAMESPACE_ {
1052
1053  inline quaternion make_quaternion_from_axis_angle(const float3 &axis, float angle) {
1054    return quaternion(axis * ::std::sin(angle * 0.5f), ::std::cos(angle * 0.5f));
1055  }
1056  inline quaternion make_quaternion_from_yaw_pitch_roll(float yaw, float pitch, float roll) {
1057    quaternion yq = make_quaternion_from_axis_angle(float3(0.0f, 1.0f, 0.0f), yaw);
1058    quaternion pq = make_quaternion_from_axis_angle(float3(1.0f, 0.0f, 0.0f), pitch);
1059    quaternion rq = make_quaternion_from_axis_angle(float3(0.0f, 0.0f, 1.0f), roll);
1060    return concatenate(concatenate(rq, pq), yq);
1061  }
1062  inline quaternion make_quaternion_from_rotation_matrix(const float4x4 &mat) {
1063    // https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
1064    float t = mat.m11 + mat.m22 + mat.m33;
1065    if (t > 0) {
1066      float r = ::std::sqrt(1.0f + t);
1067      float s = 1.0f / (2.0f * r);
1068      return { (mat.m23 - mat.m32) * s, (mat.m31 - mat.m13) * s,
1069               (mat.m12 - mat.m21) * s, r * 0.5f };
1070    } else if (mat.m11 >= mat.m22 && mat.m11 >= mat.m33) {
1071      float r = ::std::sqrt(1.0f + mat.m11 - mat.m22 - mat.m33);
1072      float s = 1.0f / (2.0f * r);
1073      return { r * 0.5f, (mat.m21 + mat.m12) * s,
1074               (mat.m13 + mat.m31) * s, (mat.m23 - mat.m32) * s };
1075    } else if (mat.m22 >= mat.m11 && mat.m22 >= mat.m33) {
1076      float r = ::std::sqrt(1.0f + mat.m22 - mat.m11 - mat.m33);
1077      float s = 1.0f / (2.0f * r);
1078      return { (mat.m21 + mat.m12) * s, r * 0.5f,
1079               (mat.m32 + mat.m23) * s, (mat.m31 - mat.m13) * s };
1080    } else {
1081      float r = ::std::sqrt(1.0f + mat.m33 - mat.m11 - mat.m22);
1082      float s = 1.0f / (2.0f * r);
1083      return { (mat.m13 + mat.m31) * s, (mat.m32 + mat.m23) * s,
1084               r * 0.5f, (mat.m12 - mat.m21) * s };
1085    }
1086  }
1087  inline bool is_identity(const quaternion &val) {
1088    return val == quaternion::identity();
1089  }
1090  inline float length(const quaternion &val) {
1091    return ::std::sqrt(length_squared(val));
1092  }
1093  inline float length_squared(const quaternion &val) {
1094    return dot(val, val);
1095  }
1096  inline float dot(const quaternion &val1, const quaternion &val2) {
1097    return val1.x * val2.x + val1.y * val2.y + val1.z * val2.z + val1.w * val2.w;
1098  }
1099  inline quaternion normalize(const quaternion &val) {
1100    return operator*(val, 1.0f / length(val));
1101  }
1102  inline quaternion conjugate(const quaternion &val) {
1103    return { -val.x, -val.y, -val.z, val.w};
1104  }
1105  inline quaternion inverse(const quaternion &val) {
1106    return operator*(conjugate(val), 1.0f / length_squared(val));
1107  }
1108  inline quaternion slerp(const quaternion &val1, const quaternion &val2, float amount) {
1109    // https://en.wikipedia.org/wiki/Slerp#Geometric_Slerp
1110    float cosangle = dot(val1, val2);
1111    quaternion end = val2;
1112    if (cosangle < 0.0f) {
1113      end = -val2;
1114      cosangle = -cosangle;
1115    }
1116    float fact1, fact2;
1117    const float epsilon = 1.0e-6f;
1118    if (cosangle > 1.0f - epsilon) {
1119      // Very small rotation range, or non-normalized input quaternions.
1120      fact1 = (1.0f - amount);
1121      fact2 = amount;
1122    } else {
1123      float angle = ::std::acos(cosangle);
1124      float sinangle = ::std::sin(angle);
1125      fact1 = ::std::sin((1.0f - amount) * angle) / sinangle;
1126      fact2 = ::std::sin(amount * angle) / sinangle;
1127    }
1128    return val1 * fact1 + end * fact2;
1129  }
1130  inline quaternion lerp(const quaternion &val1, const quaternion &val2, float amount) {
1131    quaternion end = val2;
1132    if (dot(val1, val2) < 0.0f)
1133      end = -val2;
1134    return normalize(quaternion(
1135      _impl::lerp(val1.x, end.x, amount), _impl::lerp(val1.y, end.y, amount),
1136      _impl::lerp(val1.z, end.z, amount), _impl::lerp(val1.w, end.w, amount)
1137    ));
1138  }
1139  inline quaternion concatenate(const quaternion &val1, const quaternion &val2) {
1140    return val2 * val1;
1141  }
1142
1143} _WINDOWS_NUMERICS_END_NAMESPACE_
1144
1145
1146/**
1147 * FIXME: Implement interop functions with DirectXMath.
1148 * This is where we are supposed to define the functions to convert between
1149 * Windows::Foundation::Numerics types and XMVECTOR / XMMATRIX. But our
1150 * directxmath.h does not contain the definitions for these types...
1151 */
1152#if 0
1153// === DirectXMath Interop ===
1154namespace DirectX {
1155
1156  // TODO: impl
1157  XMVECTOR XMLoadFloat2(const _WINDOWS_NUMERICS_NAMESPACE_::float2 *src);
1158  XMVECTOR XMLoadFloat3(const _WINDOWS_NUMERICS_NAMESPACE_::float3 *src);
1159  XMVECTOR XMLoadFloat4(const _WINDOWS_NUMERICS_NAMESPACE_::float4 *src);
1160  XMMATRIX XMLoadFloat3x2(const _WINDOWS_NUMERICS_NAMESPACE_::float3x2 *src);
1161  XMMATRIX XMLoadFloat4x4(const _WINDOWS_NUMERICS_NAMESPACE_::float4x4 *src);
1162  XMVECTOR XMLoadPlane(const _WINDOWS_NUMERICS_NAMESPACE_::plane *src);
1163  XMVECTOR XMLoadQuaternion(const _WINDOWS_NUMERICS_NAMESPACE_::quaternion *src);
1164  void XMStoreFloat2(_WINDOWS_NUMERICS_NAMESPACE_::float2 *out, FXMVECTOR in);
1165  void XMStoreFloat3(_WINDOWS_NUMERICS_NAMESPACE_::float3 *out, FXMVECTOR in);
1166  void XMStoreFloat4(_WINDOWS_NUMERICS_NAMESPACE_::float4 *out, FXMVECTOR in);
1167  void XMStoreFloat3x2(_WINDOWS_NUMERICS_NAMESPACE_::float3x2 *out, FXMMATRIX in);
1168  void XMStoreFloat4x4(_WINDOWS_NUMERICS_NAMESPACE_::float4x4 *out, FXMMATRIX in);
1169  void XMStorePlane(_WINDOWS_NUMERICS_NAMESPACE_::plane *out, FXMVECTOR in);
1170  void XMStoreQuaternion(_WINDOWS_NUMERICS_NAMESPACE_::quaternion *out, FXMVECTOR in);
1171
1172} /* namespace DirectX */
1173#endif
1174
1175
1176#undef _WINDOWS_NUMERICS_IMPL_ASSIGN_OP
1177
1178#ifdef _WINDOWS_NUMERICS_IMPL_PUSHED_MIN_
1179#  undef _WINDOWS_NUMERICS_IMPL_PUSHED_MIN_
1180#  pragma pop_macro("min")
1181#endif
1182
1183#ifdef _WINDOWS_NUMERICS_IMPL_PUSHED_MAX_
1184#  undef _WINDOWS_NUMERICS_IMPL_PUSHED_MAX_
1185#  pragma pop_macro("max")
1186#endif