sleepy-day/Gears
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

remove gfm

Browse Source
This commit is contained in:
matthew 2025-07-28 08:23:22 +10:00
parent 516d6f9f3e
commit 0e0205ec6f
10 changed files with 8 additions and 4347 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
dub.json
View File

@ -9,8 +9,8 @@
"targetPath": "build",
"sourceFiles-linux": ["build/libvma.a", "build/libstb_image.a", "build/libm3d.a"],
"sourceFiles-windows": [],
"importPaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
"sourcePaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
"importPaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"sourcePaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"libs-linux": ["xcb", "X11", "X11-xcb", "vulkan", "stdc++"],
"libs-windows": [],
"preGenerateCommands-linux": ["./build.sh", "build/Packer"],
@ -22,8 +22,8 @@
"targetType": "executable",
"targetPath": "build",
"targetName": "Packer",
"importPaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
"sourcePaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
"importPaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"sourcePaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"sourceFiles-linux": ["build/libstb_image.a", "build/libm3d.a"],
"preGenerateCommands-linux": ["./build.sh"],
"postGenerateCommands-linux": [],
@ -35,8 +35,8 @@
"targetType": "executable",
"targetPath": "build",
"targetName": "Codegen",
"importPaths": ["src/codegen", "src/shared", "external/xxhash", "external/gfm", "external/inteli"],
"sourcePaths": ["src/codegen", "src/shared", "external/xxhash", "external/gfm", "external/inteli"],
"importPaths": ["src/codegen", "src/shared", "external/xxhash", "external/inteli"],
"sourcePaths": ["src/codegen", "src/shared", "external/xxhash", "external/inteli"],
"sourceFiles-linux": ["build/libstb_image.a"],
"preGenerateCommands-linux": ["./build.sh"],
"preGenerateCommands-windows": [],

659
external/gfm/math/box.d vendored
View File

@ -1,659 +0,0 @@
/**
This module implements a generic axis-aligned bounding box (AABB).
*/
module gfm.math.box;
import std.math,
std.traits,
std.conv,
std.string;
import gfm.math.vector,
gfm.math.funcs;
/// N-dimensional half-open interval [a, b[.
struct Box(T, int N)
{
static assert(N > 0);
public
{
alias bound_t = Vector!(T, N);
bound_t min; // not enforced, the box can have negative volume
bound_t max;
/// Construct a box which extends between 2 points.
/// Boundaries: min is inside the box, max is just outside.
@nogc this(bound_t min_, bound_t max_) pure nothrow
{
min = min_;
max = max_;
}
static if (N == 1)
{
@nogc this(T min_, T max_) pure nothrow
{
min.x = min_;
max.x = max_;
}
}
static if (N == 2)
{
@nogc this(T min_x, T min_y, T max_x, T max_y) pure nothrow
{
min = bound_t(min_x, min_y);
max = bound_t(max_x, max_y);
}
}
static if (N == 3)
{
@nogc this(T min_x, T min_y, T min_z, T max_x, T max_y, T max_z) pure nothrow
{
min = bound_t(min_x, min_y, min_z);
max = bound_t(max_x, max_y, max_z);
}
}
@property
{
/// Returns: Dimensions of the box.
@nogc bound_t size() pure const nothrow
{
return max - min;
}
/// Sets size of the box assuming min point is the pivot.
/// Returns: Dimensions of the box.
@nogc bound_t size(bound_t value) pure nothrow
{
max = min + value;
return value;
}
/// Returns: Center of the box.
@nogc bound_t center() pure const nothrow
{
return (min + max) / 2;
}
static if (N >= 1)
{
/// Returns: Width of the box, always applicable.
@nogc T width() pure const nothrow @property
{
return max.x - min.x;
}
/// Sets width of the box assuming min point is the pivot.
/// Returns: Width of the box, always applicable.
@nogc T width(T value) pure nothrow @property
{
max.x = min.x + value;
return value;
}
}
static if (N >= 2)
{
/// Returns: Height of the box, if applicable.
@nogc T height() pure const nothrow @property
{
return max.y - min.y;
}
/// Sets height of the box assuming min point is the pivot.
/// Returns: Height of the box, if applicable.
@nogc T height(T value) pure nothrow @property
{
max.y = min.y + value;
return value;
}
}
static if (N >= 3)
{
/// Returns: Depth of the box, if applicable.
@nogc T depth() pure const nothrow @property
{
return max.z - min.z;
}
/// Sets depth of the box assuming min point is the pivot.
/// Returns: Depth of the box, if applicable.
@nogc T depth(T value) pure nothrow @property
{
max.z = min.z + value;
return value;
}
}
/// Returns: Signed volume of the box.
@nogc T volume() pure const nothrow
{
T res = 1;
bound_t size = size();
for(int i = 0; i < N; ++i)
res *= size[i];
return res;
}
/// Returns: true if empty.
@nogc bool empty() pure const nothrow
{
bound_t size = size();
mixin(generateLoopCode!("if (min[@] == max[@]) return true;", N)());
return false;
}
}
/// Returns: true if it contains point.
@nogc bool contains(bound_t point) pure const nothrow
{
assert(isSorted());
for(int i = 0; i < N; ++i)
if ( !(point[i] >= min[i] && point[i] < max[i]) )
return false;
return true;
}
static if (N >= 2)
{
/// Returns: true if it contains point `x`, `y`.
@nogc bool contains(T x, T y) pure const nothrow
{
assert(isSorted());
if ( !(x >= min.x && x < max.x) )
return false;
if ( !(y >= min.y && y < max.y) )
return false;
return true;
}
}
static if (N >= 3)
{
/// Returns: true if it contains point `x`, `y`, `z`.
@nogc bool contains(T x, T y, T z) pure const nothrow
{
assert(isSorted());
if ( !(x >= min.x && x < max.x) )
return false;
if ( !(y >= min.y && y < max.y) )
return false;
if ( !(z >= min.z && z < max.z) )
return false;
return true;
}
}
/// Returns: true if it contains box other.
@nogc bool contains(Box other) pure const nothrow
{
assert(isSorted());
assert(other.isSorted());
mixin(generateLoopCode!("if ( (other.min[@] < min[@]) || (other.max[@] > max[@]) ) return false;", N)());
return true;
}
/// Euclidean squared distance from a point.
/// See_also: Numerical Recipes Third Edition (2007)
@nogc real squaredDistance(bound_t point) pure const nothrow
{
assert(isSorted());
real distanceSquared = 0;
for (int i = 0; i < N; ++i)
{
if (point[i] < min[i])
distanceSquared += (point[i] - min[i]) ^^ 2;
if (point[i] > max[i])
distanceSquared += (point[i] - max[i]) ^^ 2;
}
return distanceSquared;
}
/// Euclidean distance from a point.
/// See_also: squaredDistance.
@nogc real distance(bound_t point) pure const nothrow
{
return sqrt(squaredDistance(point));
}
/// Euclidean squared distance from another box.
/// See_also: Numerical Recipes Third Edition (2007)
@nogc real squaredDistance(Box o) pure const nothrow
{
assert(isSorted());
assert(o.isSorted());
real distanceSquared = 0;
for (int i = 0; i < N; ++i)
{
if (o.max[i] < min[i])
distanceSquared += (o.max[i] - min[i]) ^^ 2;
if (o.min[i] > max[i])
distanceSquared += (o.min[i] - max[i]) ^^ 2;
}
return distanceSquared;
}
/// Euclidean distance from another box.
/// See_also: squaredDistance.
@nogc real distance(Box o) pure const nothrow
{
return sqrt(squaredDistance(o));
}
/// Assumes sorted boxes.
/// This function deals with empty boxes correctly.
/// Returns: Intersection of two boxes.
@nogc Box intersection(Box o) pure const nothrow
{
assert(isSorted());
assert(o.isSorted());
// Return an empty box if one of the boxes is empty
if (empty())
return this;
if (o.empty())
return o;
Box result = void;
for (int i = 0; i < N; ++i)
{
T maxOfMins = (min.v[i] > o.min.v[i]) ? min.v[i] : o.min.v[i];
T minOfMaxs = (max.v[i] < o.max.v[i]) ? max.v[i] : o.max.v[i];
result.min.v[i] = maxOfMins;
result.max.v[i] = minOfMaxs >= maxOfMins ? minOfMaxs : maxOfMins;
}
return result;
}
/// Assumes sorted boxes.
/// This function deals with empty boxes correctly.
/// Returns: Intersection of two boxes.
@nogc bool intersects(Box other) pure const nothrow
{
Box inter = this.intersection(other);
return inter.isSorted() && !inter.empty();
}
/// Extends the area of this Box.
@nogc Box grow(bound_t space) pure const nothrow
{
Box res = this;
res.min -= space;
res.max += space;
return res;
}
/// Shrink the area of this Box. The box might became unsorted.
@nogc Box shrink(bound_t space) pure const nothrow
{
return grow(-space);
}
/// Extends the area of this Box.
@nogc Box grow(T space) pure const nothrow
{
return grow(bound_t(space));
}
/// Translate this Box.
@nogc Box translate(bound_t offset) pure const nothrow
{
return Box(min + offset, max + offset);
}
static if (N >= 2)
{
/// Translate this Box by `x`, `y`.
@nogc Box translate(T x, T y) pure const nothrow
{
Box res = this;
res.min.x += x;
res.min.y += y;
res.max.x += x;
res.max.y += y;
return res;
}
}
static if (N >= 3)
{
/// Translate this Box by `x`, `y`.
@nogc Box translate(T x, T y, T z) pure const nothrow
{
Box res = this;
res.min.x += x;
res.min.y += y;
res.min.z += z;
res.max.x += x;
res.max.y += y;
res.max.z += z;
return res;
}
}
/// Shrinks the area of this Box.
/// Returns: Shrinked box.
@nogc Box shrink(T space) pure const nothrow
{
return shrink(bound_t(space));
}
/// Expands the box to include point.
/// Returns: Expanded box.
@nogc Box expand(bound_t point) pure const nothrow
{
import vector = gfm.math.vector;
return Box(vector.minByElem(min, point), vector.maxByElem(max, point));
}
/// Expands the box to include another box.
/// This function deals with empty boxes correctly.
/// Returns: Expanded box.
@nogc Box expand(Box other) pure const nothrow
{
assert(isSorted());
assert(other.isSorted());
// handle empty boxes
if (empty())
return other;
if (other.empty())
return this;
Box result = void;
for (int i = 0; i < N; ++i)
{
T minOfMins = (min.v[i] < other.min.v[i]) ? min.v[i] : other.min.v[i];
T maxOfMaxs = (max.v[i] > other.max.v[i]) ? max.v[i] : other.max.v[i];
result.min.v[i] = minOfMins;
result.max.v[i] = maxOfMaxs;
}
return result;
}
/// Returns: true if each dimension of the box is >= 0.
@nogc bool isSorted() pure const nothrow
{
for(int i = 0; i < N; ++i)
{
if (min[i] > max[i])
return false;
}
return true;
}
/// Returns: Absolute value of the Box to ensure each dimension of the
/// box is >= 0.
@nogc Box abs() pure const nothrow
{
Box!(T, N) s = this;
for (int i = 0; i < N; ++i)
{
if (s.min.v[i] > s.max.v[i])
{
T tmp = s.min.v[i];
s.min.v[i] = s.max.v[i];
s.max.v[i] = tmp;
}
}
return s;
}
/// Assign with another box.
@nogc ref Box opAssign(U)(U x) nothrow if (isBox!U)
{
static if(is(U.element_t : T))
{
static if(U._size == _size)
{
min = x.min;
max = x.max;
}
else
{
static assert(false, "no conversion between boxes with different dimensions");
}
}
else
{
static assert(false, "no conversion from " ~ U.element_t.stringof ~ " to " ~ element_t.stringof);
}
return this;
}
/// Returns: true if comparing equal boxes.
@nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Box))
{
return (min == other.min) && (max == other.max);
}
/// Cast to other box types.
@nogc U opCast(U)() pure const nothrow if (isBox!U)
{
U b = void;
for(int i = 0; i < N; ++i)
{
b.min[i] = cast(U.element_t)(min[i]);
b.max[i] = cast(U.element_t)(max[i]);
}
return b; // return a box where each element has been casted
}
static if (N == 2)
{
/// Helper function to create rectangle with a given point, width and height.
static @nogc Box rectangle(T x, T y, T width, T height) pure nothrow
{
return Box(x, y, x + width, y + height);
}
}
}
private
{
enum _size = N;
alias T element_t;
}
}
/// Instanciate to use a 2D box.
template box2(T)
{
alias Box!(T, 2) box2;
}
/// Instanciate to use a 3D box.
template box3(T)
{
alias Box!(T, 3) box3;
}
alias box2!int box2i; /// 2D box with integer coordinates.
alias box3!int box3i; /// 3D box with integer coordinates.
alias box2!float box2f; /// 2D box with float coordinates.
alias box3!float box3f; /// 3D box with float coordinates.
alias box2!double box2d; /// 2D box with double coordinates.
alias box3!double box3d; /// 3D box with double coordinates.
/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
box2i rectangle(int x, int y, int width, int height) pure nothrow @nogc
{
return box2i(x, y, x + width, y + height);
}
/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
box2f rectanglef(float x, float y, float width, float height) pure nothrow @nogc
{
return box2f(x, y, x + width, y + height);
}
/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
box2d rectangled(double x, double y, double width, double height) pure nothrow @nogc
{
return box2d(x, y, x + width, y + height);
}
unittest
{
box2i a = box2i(1, 2, 3, 4);
assert(a.width == 2);
assert(a.height == 2);
assert(a.volume == 4);
box2i b = box2i(vec2i(1, 2), vec2i(3, 4));
assert(a == b);
box3i q = box3i(-3, -2, -1, 0, 1, 2);
q.bound_t s = q.bound_t(11, 17, 19);
q.bound_t q_min = q.min;
assert((q.size = s) == s);
assert(q.size == s);
assert(q.min == q_min);
assert(q.max == q.min + s);
assert(q.max - q.min == s);
assert((q.width = s.z) == s.z);
assert(q.width == s.z);
assert(q.min.x == q_min.x);
assert(q.max.x == q.min.x + s.z);
assert(q.max.x - q.min.x == s.z);
assert((q.height = s.y) == s.y);
assert(q.height == s.y);
assert(q.min.y == q_min.y);
assert(q.max.y == q.min.y + s.y);
assert(q.max.y - q.min.y == s.y);
assert((q.depth = s.x) == s.x);
assert(q.depth == s.x);
assert(q.min.z == q_min.z);
assert(q.max.z == q.min.z + s.x);
assert(q.max.z - q.min.z == s.x);
assert(q.size == s.zyx);
box3i n = box3i(2, 1, 0, -1, -2, -3);
assert(n.abs == box3i(-1, -2, -3, 2, 1, 0));
box2f bf = cast(box2f)b;
assert(bf == box2f(1.0f, 2.0f, 3.0f, 4.0f));
box3f qf = box3f(-0, 1f, 2.5f, 3.25f, 5.125f, 7.0625f);
qf.bound_t sf = qf.bound_t(-11.5f, -17.25f, -19.125f);
qf.bound_t qf_min = qf.min;
assert((qf.size = sf) == sf);
assert(qf.size == sf);
assert(qf.min == qf_min);
assert(qf.max == qf.min + sf);
assert(qf.max - qf.min == sf);
assert((qf.width = sf.z) == sf.z);
assert(qf.width == sf.z);
assert(qf.min.x == qf_min.x);
assert(qf.max.x == qf.min.x + sf.z);
assert(qf.max.x - qf.min.x == sf.z);
assert((qf.height = sf.y) == sf.y);
assert(qf.height == sf.y);
assert(qf.min.y == qf_min.y);
assert(qf.max.y == qf.min.y + sf.y);
assert(qf.max.y - qf.min.y == sf.y);
assert((qf.depth = sf.x) == sf.x);
assert(qf.depth == sf.x);
assert(qf.min.z == qf_min.z);
assert(qf.max.z == qf.min.z + sf.x);
assert(qf.max.z - qf.min.z == sf.x);
assert(qf.size == sf.zyx);
box2i c = box2i(0, 0, 1,1);
assert(c.translate(vec2i(3, 3)) == box2i(3, 3, 4, 4));
assert(c.translate(3, 3) == box2i(3, 3, 4, 4));
assert(c.contains(vec2i(0, 0)));
assert(c.contains(0, 0));
assert(!c.contains(vec2i(1, 1)));
assert(!c.contains(1, 1));
assert(b.contains(b));
box2i d = c.expand(vec2i(3, 3));
assert(d.contains(vec2i(2, 2)));
assert(d == d.expand(d));
assert(!box2i(0, 0, 4, 4).contains(box2i(2, 2, 6, 6)));
assert(box2f(0, 0, 0, 0).empty());
assert(!box2f(0, 2, 1, 1).empty());
assert(!box2f(0, 0, 1, 1).empty());
assert(box2i(260, 100, 360, 200).intersection(box2i(100, 100, 200, 200)).empty());
// union with empty box is identity
assert(a.expand(box2i(10, 4, 10, 6)) == a);
// intersection with empty box is empty
assert(a.intersection(box2i(10, 4, 10, 6)).empty);
assert(box2i.rectangle(1, 2, 3, 4) == box2i(1, 2, 4, 6));
assert(rectangle(1, 2, 3, 4) == box2i(1, 2, 4, 6));
assert(rectanglef(1, 2, 3, 4) == box2f(1, 2, 4, 6));
assert(rectangled(1, 2, 3, 4) == box2d(1, 2, 4, 6));
}
/// True if `T` is a kind of Box
enum isBox(T) = is(T : Box!U, U...);
unittest
{
static assert( isBox!box2f);
static assert( isBox!box3d);
static assert( isBox!(Box!(real, 2)));
static assert(!isBox!vec2f);
}
/// Get the numeric type used to measure a box's dimensions.
alias DimensionType(T : Box!U, U...) = U[0];
///
unittest
{
static assert(is(DimensionType!box2f == float));
static assert(is(DimensionType!box3d == double));
}
private
{
static string generateLoopCode(string formatString, int N)() pure nothrow
{
string result;
for (int i = 0; i < N; ++i)
{
string index = ctIntToString(i);
// replace all @ by indices
result ~= formatString.replace("@", index);
}
return result;
}
// Speed-up CTFE conversions
static string ctIntToString(int n) pure nothrow
{
static immutable string[16] table = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"];
if (n < 10)
return table[n];
else
return to!string(n);
}
}

486
external/gfm/math/funcs.d vendored
View File

@ -1,486 +0,0 @@
/**
Useful math functions and range-based statistic computations.
If you need real statistics, consider using the $(WEB github.com/dsimcha/dstats,Dstats) library.
*/
module gfm.math.funcs;
import std.math,
std.traits,
std.range,
std.math;
import inteli.xmmintrin;
import gfm.math.vector : Vector;
/// Convert from radians to degrees.
@nogc T degrees(T)(in T x) pure nothrow
if (isFloatingPoint!T || (is(T : Vector!(U, n), U, int n) && isFloatingPoint!U))
{
static if (is(T : Vector!(U, n), U, int n))
return x * U(180 / PI);
else
return x * T(180 / PI);
}
/// Convert from degrees to radians.
@nogc T radians(T)(in T x) pure nothrow
if (isFloatingPoint!T || (is(T : Vector!(U, n), U, int n) && isFloatingPoint!U))
{
static if (is(T : Vector!(U, n), U, int n))
return x * U(PI / 180);
else
return x * T(PI / 180);
}
/// Linear interpolation, akin to GLSL's mix.
@nogc S lerp(S, T)(S a, S b, T t) pure nothrow
if (is(typeof(t * b + (1 - t) * a) : S))
{
return t * b + (1 - t) * a;
}
/// Clamp x in [min, max], akin to GLSL's clamp.
@nogc T clamp(T)(T x, T min, T max) pure nothrow
{
if (x < min)
return min;
else if (x > max)
return max;
else
return x;
}
/// Integer truncation.
@nogc long ltrunc(real x) nothrow // may be pure but trunc isn't pure
{
return cast(long)(trunc(x));
}
/// Integer flooring.
@nogc long lfloor(real x) nothrow // may be pure but floor isn't pure
{
return cast(long)(floor(x));
}
/// Returns: Fractional part of x.
@nogc T fract(T)(real x) nothrow
{
return x - lfloor(x);
}
/// Safe asin: input clamped to [-1, 1]
@nogc T safeAsin(T)(T x) pure nothrow
{
return asin(clamp!T(x, -1, 1));
}
/// Safe acos: input clamped to [-1, 1]
@nogc T safeAcos(T)(T x) pure nothrow
{
return acos(clamp!T(x, -1, 1));
}
/// Same as GLSL step function.
/// 0.0 is returned if x < edge, and 1.0 is returned otherwise.
@nogc T step(T)(T edge, T x) pure nothrow
{
return (x < edge) ? 0 : 1;
}
/// Same as GLSL smoothstep function.
/// See: http://en.wikipedia.org/wiki/Smoothstep
@nogc T smoothStep(T)(T a, T b, T t) pure nothrow
{
if (t <= a)
return 0;
else if (t >= b)
return 1;
else
{
T x = (t - a) / (b - a);
return x * x * (3 - 2 * x);
}
}
/// Returns: true of i is a power of 2.
@nogc bool isPowerOf2(T)(T i) pure nothrow if (isIntegral!T)
{
assert(i >= 0);
return (i != 0) && ((i & (i - 1)) == 0);
}
/// Integer log2
@nogc int ilog2(T)(T i) nothrow if (isIntegral!T)
{
assert(i > 0);
assert(isPowerOf2(i));
import core.bitop : bsr;
return bsr(i);
}
/// Computes next power of 2.
@nogc int nextPowerOf2(int i) pure nothrow
{
int v = i - 1;
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
v++;
assert(isPowerOf2(v));
return v;
}
/// Computes next power of 2.
@nogc long nextPowerOf2(long i) pure nothrow
{
long v = i - 1;
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
v |= v >> 32;
v++;
assert(isPowerOf2(v));
return v;
}
/// Computes sin(x)/x accurately.
/// See_also: $(WEB www.plunk.org/~hatch/rightway.php)
@nogc T sinOverX(T)(T x) pure nothrow
{
if (1 + x * x == 1)
return 1;
else
return sin(x) / x;
}
/// Signed integer modulo a/b where the remainder is guaranteed to be in [0..b[,
/// even if a is negative. Only support positive dividers.
@nogc T moduloWrap(T)(T a, T b) pure nothrow if (isSigned!T)
{
assert(b > 0);
if (a >= 0)
a = a % b;
else
{
auto rem = a % b;
x = (rem == 0) ? 0 : (-rem + b);
}
assert(x >= 0 && x < b);
return x;
}
unittest
{
assert(nextPowerOf2(13) == 16);
}
/**
* Find the root of a linear polynomial a + b x = 0
* Returns: Number of roots.
*/
@nogc int solveLinear(T)(T a, T b, out T root) pure nothrow if (isFloatingPoint!T)
{
if (b == 0)
{
return 0;
}
else
{
root = -a / b;
return 1;
}
}
/**
* Finds the root roots of a quadratic polynomial a + b x + c x^2 = 0
* Params:
* a = Coefficient.
* b = Coefficient.
* c = Coefficient.
* outRoots = array of root results, should have room for at least 2 elements.
* Returns: Number of roots in outRoots.
*/
@nogc int solveQuadratic(T)(T a, T b, T c, T[] outRoots) pure nothrow if (isFloatingPoint!T)
{
assert(outRoots.length >= 2);
if (c == 0)
return solveLinear(a, b, outRoots[0]);
T delta = b * b - 4 * a * c;
if (delta < 0.0 )
return 0;
delta = sqrt(delta);
T oneOver2a = 0.5 / a;
outRoots[0] = oneOver2a * (-b - delta);
outRoots[1] = oneOver2a * (-b + delta);
return 2;
}
/**
* Finds the roots of a cubic polynomial a + b x + c x^2 + d x^3 = 0
* Params:
* a = Coefficient.
* b = Coefficient.
* c = Coefficient.
* d = Coefficient.
* outRoots = array of root results, should have room for at least 2 elements.
* Returns: Number of roots in outRoots.
* See_also: $(WEB www.codeguru.com/forum/archive/index.php/t-265551.html)
*/
@nogc int solveCubic(T)(T a, T b, T c, T d, T[] outRoots) pure nothrow if (isFloatingPoint!T)
{
assert(outRoots.length >= 3);
if (d == 0)
return solveQuadratic(a, b, c, outRoots);
// adjust coefficients
T a1 = c / d,
a2 = b / d,
a3 = a / d;
T Q = (a1 * a1 - 3 * a2) / 9,
R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) / 54;
T Qcubed = Q * Q * Q;
T d2 = Qcubed - R * R;
if (d2 >= 0)
{
// 3 real roots
if (Q < 0.0)
return 0;
T P = R / sqrt(Qcubed);
assert(-1 <= P && P <= 1);
T theta = acos(P);
T sqrtQ = sqrt(Q);
outRoots[0] = -2 * sqrtQ * cos(theta / 3) - a1 / 3;
outRoots[1] = -2 * sqrtQ * cos((theta + 2 * T(PI)) / 3) - a1 / 3;
outRoots[2] = -2 * sqrtQ * cos((theta + 4 * T(PI)) / 3) - a1 / 3;
return 3;
}
else
{
// 1 real root
T e = (sqrt(-d) + abs(R)) ^^ cast(T)(1.0 / 3.0);
if (R > 0)
e = -e;
outRoots[0] = e + Q / e - a1 / 3.0;
return 1;
}
}
/**
* Returns the roots of a quartic polynomial a + b x + c x^2 + d x^3 + e x^4 = 0
*
* Returns number of roots. roots slice should have room for up to 4 elements.
* Bugs: doesn't pass unit-test!
* See_also: $(WEB mathworld.wolfram.com/QuarticEquation.html)
*/
@nogc int solveQuartic(T)(T a, T b, T c, T d, T e, T[] roots) pure nothrow if (isFloatingPoint!T)
{
assert(roots.length >= 4);
if (e == 0)
return solveCubic(a, b, c, d, roots);
// Adjust coefficients
T a0 = a / e,
a1 = b / e,
a2 = c / e,
a3 = d / e;
// Find a root for the following cubic equation:
// y^3 - a2 y^2 + (a1 a3 - 4 a0) y + (4 a2 a0 - a1 ^2 - a3^2 a0) = 0
// aka Resolvent cubic
T b0 = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0;
T b1 = a1 * a3 - 4 * a0;
T b2 = -a2;
T[3] resolventCubicRoots;
int numRoots = solveCubic!T(b0, b1, b2, 1, resolventCubicRoots[]);
assert(numRoots == 3);
T y = resolventCubicRoots[0];
if (y < resolventCubicRoots[1]) y = resolventCubicRoots[1];
if (y < resolventCubicRoots[2]) y = resolventCubicRoots[2];
// Compute R, D & E
T R = 0.25f * a3 * a3 - a2 + y;
if (R < 0.0)
return 0;
R = sqrt(R);
T D = void,
E = void;
if (R == 0)
{
T d1 = 0.75f * a3 * a3 - 2 * a2;
T d2 = 2 * sqrt(y * y - 4 * a0);
D = sqrt(d1 + d2) * 0.5f;
E = sqrt(d1 - d2) * 0.5f;
}
else
{
T Rsquare = R * R;
T Rrec = 1 / R;
T d1 = 0.75f * a3 * a3 - Rsquare - 2 * a2;
T d2 = 0.25f * Rrec * (4 * a3 * a2 - 8 * a1 - a3 * a3 * a3);
D = sqrt(d1 + d2) * 0.5f;
E = sqrt(d1 - d2) * 0.5f;
}
// Compute the 4 roots
a3 *= -0.25f;
R *= 0.5f;
roots[0] = a3 + R + D;
roots[1] = a3 + R - D;
roots[2] = a3 - R + E;
roots[3] = a3 - R - E;
return 4;
}
unittest
{
bool arrayContainsRoot(double[] arr, double root)
{
foreach(e; arr)
if (abs(e - root) < 1e-7)
return true;
return false;
}
// test quadratic
{
double[3] roots;
int numRoots = solveCubic!double(-2, -3 / 2.0, 3 / 4.0, 1 / 4.0, roots[]);
assert(numRoots == 3);
assert(arrayContainsRoot(roots[], -4));
assert(arrayContainsRoot(roots[], -1));
assert(arrayContainsRoot(roots[], 2));
}
// test quartic
{
double[4] roots;
int numRoots = solveQuartic!double(0, -2, -1, 2, 1, roots[]);
assert(numRoots == 4);
assert(arrayContainsRoot(roots[], -2));
assert(arrayContainsRoot(roots[], -1));
assert(arrayContainsRoot(roots[], 0));
assert(arrayContainsRoot(roots[], 1));
}
}
/// Arithmetic mean.
double average(R)(R r) if (isInputRange!R)
{
if (r.empty)
return double.nan;
typeof(r.front()) sum = 0;
long count = 0;
foreach(e; r)
{
sum += e;
++count;
}
return sum / count;
}
/// Minimum of a range.
double minElement(R)(R r) if (isInputRange!R)
{
// do like Javascript for an empty range
if (r.empty)
return double.infinity;
return minmax!("<", R)(r);
}
/// Maximum of a range.
double maxElement(R)(R r) if (isInputRange!R)
{
// do like Javascript for an empty range
if (r.empty)
return -double.infinity;
return minmax!(">", R)(r);
}
/// Variance of a range.
double variance(R)(R r) if (isForwardRange!R)
{
if (r.empty)
return double.nan;
auto avg = average(r.save); // getting the average
typeof(avg) sum = 0;
long count = 0;
foreach(e; r)
{
sum += (e - avg) ^^ 2;
++count;
}
if (count <= 1)
return 0.0;
else
return (sum / (count - 1.0)); // using sample std deviation as estimator
}
/// Standard deviation of a range.
double standardDeviation(R)(R r) if (isForwardRange!R)
{
return sqrt(variance(r));
}
private
{
typeof(R.front()) minmax(string op, R)(R r) if (isInputRange!R)
{
assert(!r.empty);
auto best = r.front();
r.popFront();
foreach(e; r)
{
mixin("if (e " ~ op ~ " best) best = e;");
}
return best;
}
}
/// SSE approximation of reciprocal square root.
@nogc T inverseSqrt(T)(T x) pure nothrow if (isFloatingPoint!T)
{
static if (is(T == float))
{
__m128 V = _mm_set_ss(x);
V = _mm_rsqrt_ss(V);
return _mm_cvtss_f32(V);
}
else
{
return 1 / sqrt(x);
}
}
unittest
{
assert(abs( inverseSqrt!float(1) - 1) < 1e-3 );
assert(abs( inverseSqrt!double(1) - 1) < 1e-3 );
}

948
external/gfm/math/matrix.d vendored
View File

@ -1,948 +0,0 @@
/// Custom sized 2-dimension Matrices
module gfm.math.matrix;
import std.math,
std.typetuple,
std.traits,
std.string,
std.typecons,
std.conv;
import gfm.math.vector,
gfm.math.shapes,
gfm.math.quaternion;
/// Generic non-resizeable matrix with R rows and C columns.
/// Intended for 3D use (size 3x3 and 4x4).
/// Important: <b>Matrices here are in row-major order whereas OpenGL is column-major.</b>
/// Params:
/// T = type of elements
/// R = number of rows
/// C = number of columns
struct Matrix(T, int R, int C)
{
public
{
static assert(R >= 1 && C >= 1);
alias Vector!(T, C) row_t;
alias Vector!(T, R) column_t;
enum bool isSquare = (R == C);
// fields definition
union
{
T[C*R] v; // all elements
row_t[R] rows; // all rows
T[C][R] c; // components
}
@nogc this(U...)(U values) pure nothrow
{
static if ((U.length == C*R) && allSatisfy!(isTAssignable, U))
{
// construct with components
foreach(int i, x; values)
v[i] = x;
}
else static if ((U.length == 1) && (isAssignable!(U[0])) && (!is(U[0] : Matrix)))
{
// construct with assignment
opAssign!(U[0])(values[0]);
}
else static assert(false, "cannot create a matrix from given arguments");
}
/// Construct a matrix from columns.
@nogc static Matrix fromColumns(column_t[] columns) pure nothrow
{
assert(columns.length == C);
Matrix res;
for (int i = 0; i < R; ++i)
for (int j = 0; j < C; ++j)
{
res.c[i][j] = columns[j][i];
}
return res;
}
/// Construct a matrix from rows.
@nogc static Matrix fromRows(row_t[] rows) pure nothrow
{
assert(rows.length == R);
Matrix res;
res.rows[] = rows[];
return res;
}
/// Construct matrix with a scalar.
@nogc this(U)(T x) pure nothrow
{
for (int i = 0; i < _N; ++i)
v[i] = x;
}
/// Assign with a scalar.
@nogc ref Matrix opAssign(U : T)(U x) pure nothrow
{
for (int i = 0; i < R * C; ++i)
v[i] = x;
return this;
}
/// Assign with a samey matrice.
@nogc ref Matrix opAssign(U : Matrix)(U x) pure nothrow
{
for (int i = 0; i < R * C; ++i)
v[i] = x.v[i];
return this;
}
/// Assign from other small matrices (same size, compatible type).
@nogc ref Matrix opAssign(U)(U x) pure nothrow
if (isMatrixInstantiation!U
&& is(U._T : _T)
&& (!is(U: Matrix))
&& (U._R == R) && (U._C == C))
{
for (int i = 0; i < R * C; ++i)
v[i] = x.v[i];
return this;
}
/// Assign with a static array of size R * C.
@nogc ref Matrix opAssign(U)(U x) pure nothrow
if ((isStaticArray!U)
&& is(typeof(x[0]) : T)
&& (U.length == R * C))
{
for (int i = 0; i < R * C; ++i)
v[i] = x[i];
return this;
}
/// Assign with a static array of shape (R, C).
@nogc ref Matrix opAssign(U)(U x) pure nothrow
if ((isStaticArray!U) && isStaticArray!(typeof(x[0]))
&& is(typeof(x[0][0]) : T)
&& (U.length == R)
&& (x[0].length == C))
{
foreach (i; 0..R)
rows[i] = x[i];
return this;
}
/// Assign with a dynamic array of size R * C.
@nogc ref Matrix opAssign(U)(U x) pure nothrow
if ((isDynamicArray!U)
&& is(typeof(x[0]) : T))
{
assert(x.length == R * C);
for (int i = 0; i < R * C; ++i)
v[i] = x[i];
return this;
}
/// Assign with a dynamic array of shape (R, C).
@nogc ref Matrix opAssign(U)(U x) pure nothrow
if ((isDynamicArray!U) && isDynamicArray!(typeof(x[0]))
&& is(typeof(x[0][0]) : T))
{
assert(x.length == R);
foreach (i; 0..R)
{
assert(x[i].length == C);
rows[i] = x[i];
}
return this;
}
/// Return a pointer to content.
@nogc inout(T)* ptr() pure inout nothrow @property
{
return v.ptr;
}
/// Returns a column as a vector
/// Returns: column j as a vector.
@nogc column_t column(int j) pure const nothrow
{
column_t res = void;
for (int i = 0; i < R; ++i)
res.v[i] = c[i][j];
return res;
}
/// Returns a row as a vector
/// Returns: row i as a vector.
@nogc row_t row(int i) pure const nothrow
{
return rows[i];
}
/// Covnerts to pretty string.
string toString() const nothrow
{
try
return format("%s", v);
catch (Exception e)
assert(false); // should not happen since format is right
}
/// Matrix * scalar multiplication.
@nogc Matrix opBinary(string op)(T factor) pure const nothrow if (op == "*")
{
Matrix result = void;
for (int i = 0; i < R; ++i)
{
for (int j = 0; j < C; ++j)
{
result.c[i][j] = c[i][j] * factor;
}
}
return result;
}
/// Matrix * vector multiplication.
@nogc column_t opBinary(string op)(row_t x) pure const nothrow if (op == "*")
{
column_t res = void;
for (int i = 0; i < R; ++i)
{
T sum = 0;
for (int j = 0; j < C; ++j)
{
sum += c[i][j] * x.v[j];
}
res.v[i] = sum;
}
return res;
}
/// Matrix * matrix multiplication.
@nogc auto opBinary(string op, U)(U x) pure const nothrow
if (isMatrixInstantiation!U && (U._R == C) && (op == "*"))
{
Matrix!(T, R, U._C) result = void;
for (int i = 0; i < R; ++i)
{
for (int j = 0; j < U._C; ++j)
{
T sum = 0;
for (int k = 0; k < C; ++k)
sum += c[i][k] * x.c[k][j];
result.c[i][j] = sum;
}
}
return result;
}
/// Matrix add and substraction.
@nogc Matrix opBinary(string op, U)(U other) pure const nothrow
if (is(U : Matrix) && (op == "+" || op == "-"))
{
Matrix result = void;
for (int i = 0; i < R; ++i)
{
for (int j = 0; j < C; ++j)
{
mixin("result.c[i][j] = c[i][j] " ~ op ~ " other.c[i][j];");
}
}
return result;
}
// matrix *= scalar
@nogc ref Matrix opOpAssign(string op, U : T)(U x) pure nothrow if (op == "*")
{
for (int i = 0; i < R * C; ++i)
v[i] *= x;
return this;
}
/// Assignment operator with another samey matrix.
@nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow
if (is(U : Matrix) && (op == "*" || op == "+" || op == "-"))
{
mixin("Matrix result = this " ~ op ~ " operand;");
return opAssign!Matrix(result);
}
/// Matrix += <something convertible to a Matrix>
/// Matrix -= <something convertible to a Matrix>
@nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow
if ((isConvertible!U) && (op == "*" || op == "+" || op == "-"))
{
Matrix conv = operand;
return opOpAssign!op(conv);
}
/// Cast to other matrix types.
/// If the size are different, the resulting matrix is truncated
/// and/or filled with identity coefficients.
@nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U)
{
U res = U.identity();
enum minR = R < U._R ? R : U._R;
enum minC = C < U._C ? C : U._C;
for (int i = 0; i < minR; ++i)
for (int j = 0; j < minC; ++j)
{
res.c[i][j] = cast(U._T)(c[i][j]);
}
return res;
}
@nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Matrix))
{
for (int i = 0; i < R * C; ++i)
if (v[i] != other.v[i])
return false;
return true;
}
@nogc bool opEquals(U)(U other) pure const nothrow
if ((isAssignable!U) && (!is(U: Matrix)))
{
Matrix conv = other;
return opEquals(conv);
}
// +matrix, -matrix, ~matrix, !matrix
@nogc Matrix opUnary(string op)() pure const nothrow if (op == "+" || op == "-" || op == "~" || op == "!")
{
Matrix res = void;
for (int i = 0; i < N; ++i)
mixin("res.v[i] = " ~ op ~ "v[i];");
return res;
}
/// Convert 3x3 rotation matrix to quaternion.
/// See_also: 3D Math Primer for Graphics and Game Development.
@nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U
&& is(U._T : _T)
&& (_R == 3) && (_C == 3))
{
T fourXSquaredMinus1 = c[0][0] - c[1][1] - c[2][2];
T fourYSquaredMinus1 = c[1][1] - c[0][0] - c[2][2];
T fourZSquaredMinus1 = c[2][2] - c[0][0] - c[1][1];
T fourWSquaredMinus1 = c[0][0] + c[1][1] + c[2][2];
int biggestIndex = 0;
T fourBiggestSquaredMinus1 = fourWSquaredMinus1;
if(fourXSquaredMinus1 > fourBiggestSquaredMinus1)
{
fourBiggestSquaredMinus1 = fourXSquaredMinus1;
biggestIndex = 1;
}
if(fourYSquaredMinus1 > fourBiggestSquaredMinus1)
{
fourBiggestSquaredMinus1 = fourYSquaredMinus1;
biggestIndex = 2;
}
if(fourZSquaredMinus1 > fourBiggestSquaredMinus1)
{
fourBiggestSquaredMinus1 = fourZSquaredMinus1;
biggestIndex = 3;
}
T biggestVal = sqrt(fourBiggestSquaredMinus1 + 1) / 2;
T mult = 1 / (biggestVal * 4);
U quat;
switch(biggestIndex)
{
case 1:
quat.w = (c[1][2] - c[2][1]) * mult;
quat.x = biggestVal;
quat.y = (c[0][1] + c[1][0]) * mult;
quat.z = (c[2][0] + c[0][2]) * mult;
break;
case 2:
quat.w = (c[2][0] - c[0][2]) * mult;
quat.x = (c[0][1] + c[1][0]) * mult;
quat.y = biggestVal;
quat.z = (c[1][2] + c[2][1]) * mult;
break;
case 3:
quat.w = (c[0][1] - c[1][0]) * mult;
quat.x = (c[2][0] + c[0][2]) * mult;
quat.y = (c[1][2] + c[2][1]) * mult;
quat.z = biggestVal;
break;
default: // biggestIndex == 0
quat.w = biggestVal;
quat.x = (c[1][2] - c[2][1]) * mult;
quat.y = (c[2][0] - c[0][2]) * mult;
quat.z = (c[0][1] - c[1][0]) * mult;
break;
}
return quat;
}
/// Converts a 4x4 rotation matrix to quaternion.
@nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U
&& is(U._T : _T)
&& (_R == 4) && (_C == 4))
{
auto m3 = cast(mat3!T)(this);
return cast(U)(m3);
}
static if (isSquare && isFloatingPoint!T && R == 1)
{
/// Returns an inverted copy of this matrix
/// Returns: inverse of matrix.
/// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
@nogc Matrix inverse() pure const nothrow
{
assert(c[0][0] != 0); // Programming error if matrix is not invertible.
return Matrix( 1 / c[0][0]);
}
}
static if (isSquare && isFloatingPoint!T && R == 2)
{
/// Returns an inverted copy of this matrix
/// Returns: inverse of matrix.
/// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
@nogc Matrix inverse() pure const nothrow
{
T det = (c[0][0] * c[1][1] - c[0][1] * c[1][0]);
assert(det != 0); // Programming error if matrix is not invertible.
T invDet = 1 / det;
return Matrix( c[1][1] * invDet, -c[0][1] * invDet,
-c[1][0] * invDet, c[0][0] * invDet);
}
}
static if (isSquare && isFloatingPoint!T && R == 3)
{
/// Returns an inverted copy of this matrix
/// Returns: inverse of matrix.
/// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
@nogc Matrix inverse() pure const nothrow
{
T det = c[0][0] * (c[1][1] * c[2][2] - c[2][1] * c[1][2])
- c[0][1] * (c[1][0] * c[2][2] - c[1][2] * c[2][0])
+ c[0][2] * (c[1][0] * c[2][1] - c[1][1] * c[2][0]);
assert(det != 0); // Programming error if matrix is not invertible.
T invDet = 1 / det;
Matrix res = void;
res.c[0][0] = (c[1][1] * c[2][2] - c[2][1] * c[1][2]) * invDet;
res.c[0][1] = -(c[0][1] * c[2][2] - c[0][2] * c[2][1]) * invDet;
res.c[0][2] = (c[0][1] * c[1][2] - c[0][2] * c[1][1]) * invDet;
res.c[1][0] = -(c[1][0] * c[2][2] - c[1][2] * c[2][0]) * invDet;
res.c[1][1] = (c[0][0] * c[2][2] - c[0][2] * c[2][0]) * invDet;
res.c[1][2] = -(c[0][0] * c[1][2] - c[1][0] * c[0][2]) * invDet;
res.c[2][0] = (c[1][0] * c[2][1] - c[2][0] * c[1][1]) * invDet;
res.c[2][1] = -(c[0][0] * c[2][1] - c[2][0] * c[0][1]) * invDet;
res.c[2][2] = (c[0][0] * c[1][1] - c[1][0] * c[0][1]) * invDet;
return res;
}
}
static if (isSquare && isFloatingPoint!T && R == 4)
{
/// Returns an inverted copy of this matrix
/// Returns: inverse of matrix.
/// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
@nogc Matrix inverse() pure const nothrow
{
T det2_01_01 = c[0][0] * c[1][1] - c[0][1] * c[1][0];
T det2_01_02 = c[0][0] * c[1][2] - c[0][2] * c[1][0];
T det2_01_03 = c[0][0] * c[1][3] - c[0][3] * c[1][0];
T det2_01_12 = c[0][1] * c[1][2] - c[0][2] * c[1][1];
T det2_01_13 = c[0][1] * c[1][3] - c[0][3] * c[1][1];
T det2_01_23 = c[0][2] * c[1][3] - c[0][3] * c[1][2];
T det3_201_012 = c[2][0] * det2_01_12 - c[2][1] * det2_01_02 + c[2][2] * det2_01_01;
T det3_201_013 = c[2][0] * det2_01_13 - c[2][1] * det2_01_03 + c[2][3] * det2_01_01;
T det3_201_023 = c[2][0] * det2_01_23 - c[2][2] * det2_01_03 + c[2][3] * det2_01_02;
T det3_201_123 = c[2][1] * det2_01_23 - c[2][2] * det2_01_13 + c[2][3] * det2_01_12;
T det = - det3_201_123 * c[3][0] + det3_201_023 * c[3][1] - det3_201_013 * c[3][2] + det3_201_012 * c[3][3];
assert(det != 0); // Programming error if matrix is not invertible.
T invDet = 1 / det;
T det2_03_01 = c[0][0] * c[3][1] - c[0][1] * c[3][0];
T det2_03_02 = c[0][0] * c[3][2] - c[0][2] * c[3][0];
T det2_03_03 = c[0][0] * c[3][3] - c[0][3] * c[3][0];
T det2_03_12 = c[0][1] * c[3][2] - c[0][2] * c[3][1];
T det2_03_13 = c[0][1] * c[3][3] - c[0][3] * c[3][1];
T det2_03_23 = c[0][2] * c[3][3] - c[0][3] * c[3][2];
T det2_13_01 = c[1][0] * c[3][1] - c[1][1] * c[3][0];
T det2_13_02 = c[1][0] * c[3][2] - c[1][2] * c[3][0];
T det2_13_03 = c[1][0] * c[3][3] - c[1][3] * c[3][0];
T det2_13_12 = c[1][1] * c[3][2] - c[1][2] * c[3][1];
T det2_13_13 = c[1][1] * c[3][3] - c[1][3] * c[3][1];
T det2_13_23 = c[1][2] * c[3][3] - c[1][3] * c[3][2];
T det3_203_012 = c[2][0] * det2_03_12 - c[2][1] * det2_03_02 + c[2][2] * det2_03_01;
T det3_203_013 = c[2][0] * det2_03_13 - c[2][1] * det2_03_03 + c[2][3] * det2_03_01;
T det3_203_023 = c[2][0] * det2_03_23 - c[2][2] * det2_03_03 + c[2][3] * det2_03_02;
T det3_203_123 = c[2][1] * det2_03_23 - c[2][2] * det2_03_13 + c[2][3] * det2_03_12;
T det3_213_012 = c[2][0] * det2_13_12 - c[2][1] * det2_13_02 + c[2][2] * det2_13_01;
T det3_213_013 = c[2][0] * det2_13_13 - c[2][1] * det2_13_03 + c[2][3] * det2_13_01;
T det3_213_023 = c[2][0] * det2_13_23 - c[2][2] * det2_13_03 + c[2][3] * det2_13_02;
T det3_213_123 = c[2][1] * det2_13_23 - c[2][2] * det2_13_13 + c[2][3] * det2_13_12;
T det3_301_012 = c[3][0] * det2_01_12 - c[3][1] * det2_01_02 + c[3][2] * det2_01_01;
T det3_301_013 = c[3][0] * det2_01_13 - c[3][1] * det2_01_03 + c[3][3] * det2_01_01;
T det3_301_023 = c[3][0] * det2_01_23 - c[3][2] * det2_01_03 + c[3][3] * det2_01_02;
T det3_301_123 = c[3][1] * det2_01_23 - c[3][2] * det2_01_13 + c[3][3] * det2_01_12;
Matrix res = void;
res.c[0][0] = - det3_213_123 * invDet;
res.c[1][0] = + det3_213_023 * invDet;
res.c[2][0] = - det3_213_013 * invDet;
res.c[3][0] = + det3_213_012 * invDet;
res.c[0][1] = + det3_203_123 * invDet;
res.c[1][1] = - det3_203_023 * invDet;
res.c[2][1] = + det3_203_013 * invDet;
res.c[3][1] = - det3_203_012 * invDet;
res.c[0][2] = + det3_301_123 * invDet;
res.c[1][2] = - det3_301_023 * invDet;
res.c[2][2] = + det3_301_013 * invDet;
res.c[3][2] = - det3_301_012 * invDet;
res.c[0][3] = - det3_201_123 * invDet;
res.c[1][3] = + det3_201_023 * invDet;
res.c[2][3] = - det3_201_013 * invDet;
res.c[3][3] = + det3_201_012 * invDet;
return res;
}
}
/// Returns a transposed copy of this matrix
/// Returns: transposed matrice.
@nogc Matrix!(T, C, R) transposed() pure const nothrow
{
Matrix!(T, C, R) res;
for (int i = 0; i < C; ++i)
for (int j = 0; j < R; ++j)
res.c[i][j] = c[j][i];
return res;
}
static if (isSquare && R > 1)
{
/// Makes a diagonal matrix from a vector.
@nogc static Matrix diag(Vector!(T, R) v) pure nothrow
{
Matrix res = void;
for (int i = 0; i < R; ++i)
for (int j = 0; j < C; ++j)
res.c[i][j] = (i == j) ? v.v[i] : 0;
return res;
}
/// In-place translate by (v, 1)
@nogc void translate(Vector!(T, R-1) v) pure nothrow
{
for (int i = 0; i < R; ++i)
{
T dot = 0;
for (int j = 0; j + 1 < C; ++j)
dot += v.v[j] * c[i][j];
c[i][C-1] += dot;
}
}
/// Make a translation matrix.
@nogc static Matrix translation(Vector!(T, R-1) v) pure nothrow
{
Matrix res = identity();
for (int i = 0; i + 1 < R; ++i)
res.c[i][C-1] += v.v[i];
return res;
}
/// In-place matrix scaling.
void scale(Vector!(T, R-1) v) pure nothrow
{
for (int i = 0; i < R; ++i)
for (int j = 0; j + 1 < C; ++j)
c[i][j] *= v.v[j];
}
/// Make a scaling matrix.
@nogc static Matrix scaling(Vector!(T, R-1) v) pure nothrow
{
Matrix res = identity();
for (int i = 0; i + 1 < R; ++i)
res.c[i][i] = v.v[i];
return res;
}
}
// rotations are implemented for 3x3 and 4x4 matrices.
static if (isSquare && (R == 3 || R == 4) && isFloatingPoint!T)
{
@nogc public static Matrix rotateAxis(int i, int j)(T angle) pure nothrow
{
Matrix res = identity();
const T cosa = cos(angle);
const T sina = sin(angle);
res.c[i][i] = cosa;
res.c[i][j] = -sina;
res.c[j][i] = sina;
res.c[j][j] = cosa;
return res;
}
/// Rotate along X axis
/// Returns: rotation matrix along axis X
alias rotateAxis!(1, 2) rotateX;
/// Rotate along Y axis
/// Returns: rotation matrix along axis Y
alias rotateAxis!(2, 0) rotateY;
/// Rotate along Z axis
/// Returns: rotation matrix along axis Z
alias rotateAxis!(0, 1) rotateZ;
/// Similar to the glRotate matrix, however the angle is expressed in radians
/// See_also: $(LINK http://www.cs.rutgers.edu/~decarlo/428/gl_man/rotate.html)
@nogc static Matrix rotation(T angle, vec3!T axis) pure nothrow
{
Matrix res = identity();
const T c = cos(angle);
const oneMinusC = 1 - c;
const T s = sin(angle);
axis = axis.normalized();
T x = axis.x,
y = axis.y,
z = axis.z;
T xy = x * y,
yz = y * z,
xz = x * z;
res.c[0][0] = x * x * oneMinusC + c;
res.c[0][1] = x * y * oneMinusC - z * s;
res.c[0][2] = x * z * oneMinusC + y * s;
res.c[1][0] = y * x * oneMinusC + z * s;
res.c[1][1] = y * y * oneMinusC + c;
res.c[1][2] = y * z * oneMinusC - x * s;
res.c[2][0] = z * x * oneMinusC - y * s;
res.c[2][1] = z * y * oneMinusC + x * s;
res.c[2][2] = z * z * oneMinusC + c;
return res;
}
}
// 4x4 specific transformations for 3D usage
static if (isSquare && R == 4 && isFloatingPoint!T)
{
/// Orthographic projection
/// Returns: orthographic projection.
@nogc static Matrix orthographic(T left, T right, T bottom, T top, T near, T far) pure nothrow
{
T dx = right - left,
dy = top - bottom,
dz = far - near;
T tx = -(right + left) / dx;
T ty = -(top + bottom) / dy;
T tz = -(far + near) / dz;
return Matrix(2 / dx, 0, 0, tx,
0, 2 / dy, 0, ty,
0, 0, -2 / dz, tz,
0, 0, 0, 1);
}
/// Perspective projection
/// Returns: perspective projection.
@nogc static Matrix perspective(T FOVInRadians, T aspect, T zNear, T zFar) pure nothrow
{
T f = 1 / tan(FOVInRadians / 2);
T d = 1 / (zNear - zFar);
return Matrix(f / aspect, 0, 0, 0,
0, f, 0, 0,
0, 0, (zFar + zNear) * d, 2 * d * zFar * zNear,
0, 0, -1, 0);
}
/// Look At projection
/// Returns: "lookAt" projection.
/// Thanks to vuaru for corrections.
@nogc static Matrix lookAt(vec3!T eye, vec3!T target, vec3!T up) pure nothrow
{
vec3!T Z = (eye - target).normalized();
vec3!T X = cross(-up, Z).normalized();
vec3!T Y = cross(Z, -X);
return Matrix(-X.x, -X.y, -X.z, dot(X, eye),
Y.x, Y.y, Y.z, -dot(Y, eye),
Z.x, Z.y, Z.z, -dot(Z, eye),
0, 0, 0, 1);
}
/// Extract frustum from a 4x4 matrice.
@nogc Frustum!T frustum() pure const nothrow
{
auto left = Plane!T(row(3) + row(0));
auto right = Plane!T(row(3) - row(0));
auto top = Plane!T(row(3) - row(1));
auto bottom = Plane!T(row(3) + row(1));
auto near = Plane!T(row(3) + row(2));
auto far = Plane!T(row(3) - row(2));
return Frustum!T(left, right, top, bottom, near, far);
}
}
}
package
{
alias T _T;
enum _R = R;
enum _C = C;
}
private
{
template isAssignable(T)
{
enum bool isAssignable = std.traits.isAssignable!(Matrix, T);
}
template isConvertible(T)
{
enum bool isConvertible = (!is(T : Matrix)) && isAssignable!T;
}
template isTAssignable(U)
{
enum bool isTAssignable = std.traits.isAssignable!(T, U);
}
template isRowConvertible(U)
{
enum bool isRowConvertible = is(U : row_t);
}
template isColumnConvertible(U)
{
enum bool isColumnConvertible = is(U : column_t);
}
}
public
{
/// Construct an identity matrix
/// Returns: an identity matrix.
/// Note: the identity matrix, while only meaningful for square matrices,
/// is also defined for non-square ones.
@nogc static Matrix identity() pure nothrow
{
Matrix res = void;
for (int i = 0; i < R; ++i)
for (int j = 0; j < C; ++j)
res.c[i][j] = (i == j) ? 1 : 0;
return res;
}
/// Construct an constant matrix
/// Returns: a constant matrice.
@nogc static Matrix constant(U)(U x) pure nothrow
{
Matrix res = void;
for (int i = 0; i < R * C; ++i)
res.v[i] = cast(T)x;
return res;
}
}
}
template isMatrixInstantiation(U)
{
private static void isMatrix(T, int R, int C)(Matrix!(T, R, C) x)
{
}
enum bool isMatrixInstantiation = is(typeof(isMatrix(U.init)));
}
// GLSL is a big inspiration here
// we defines types with more or less the same names
///
template mat2x2(T) { alias Matrix!(T, 2, 2) mat2x2; }
///
template mat3x3(T) { alias Matrix!(T, 3, 3) mat3x3; }
///
template mat4x4(T) { alias Matrix!(T, 4, 4) mat4x4; }
// WARNING: in GLSL, first number is _columns_, second is rows
// It is the opposite here: first number is rows, second is columns
// With this convention mat2x3 * mat3x4 -> mat2x4.
///
template mat2x3(T) { alias Matrix!(T, 2, 3) mat2x3; }
///
template mat2x4(T) { alias Matrix!(T, 2, 4) mat2x4; }
///
template mat3x2(T) { alias Matrix!(T, 3, 2) mat3x2; }
///
template mat3x4(T) { alias Matrix!(T, 3, 4) mat3x4; }
///
template mat4x2(T) { alias Matrix!(T, 4, 2) mat4x2; }
///
template mat4x3(T) { alias Matrix!(T, 4, 3) mat4x3; }
// shorter names for most common matrices
alias mat2x2 mat2;///
alias mat3x3 mat3;///
alias mat4x4 mat4;///
// Define a lot of type names
// Most useful are probably mat4f and mat4d
alias mat2!float mat2f;///
alias mat2!double mat2d;///
alias mat3!float mat3f;///
alias mat3!double mat3d;///
alias mat4!float mat4f;///
alias mat4!double mat4d;///
alias mat2x2!float mat2x2f;///
alias mat2x2!double mat2x2d;///
alias mat3x3!float mat3x3f;///
alias mat3x3!double mat3x3d;///
alias mat4x4!float mat4x4f;///
alias mat4x4!double mat4x4d;///
unittest
{
alias mat2i = mat2!int;
alias mat2x3f = mat2x3!float;
alias mat3x4f = mat3x4!float;
alias mat2x4f = mat2x4!float;
mat2i x = mat2i(0, 1,
2, 3);
assert(x.c[0][0] == 0 && x.c[0][1] == 1 && x.c[1][0] == 2 && x.c[1][1] == 3);
vec2i[2] cols = [vec2i(0, 2), vec2i(1, 3)];
mat2i y = mat2i.fromColumns(cols[]);
assert(y.c[0][0] == 0 && y.c[0][1] == 1 && y.c[1][0] == 2 && y.c[1][1] == 3);
y = mat2i.fromRows(cols[]);
assert(y.c[0][0] == 0 && y.c[1][0] == 1 && y.c[0][1] == 2 && y.c[1][1] == 3);
y = y.transposed();
assert(x == y);
x = [0, 1, 2, 3];
assert(x == y);
mat2i z = x * y;
assert(z == mat2i([2, 3, 6, 11]));
vec2i vz = z * vec2i(2, -1);
assert(vz == vec2i(1, 1));
mat2f a = z;
mat2d ad = a;
ad += a;
mat2f w = [4, 5, 6, 7];
z = cast(mat2i)w;
assert(w == z);
{
mat2x3f A;
mat3x4f B;
mat2x4f C = A * B;
}
assert(mat2i.diag(vec2i(1, 2)) == mat2i(1, 0,
0, 2));
// Construct with a single scalar
auto D = mat4f(1.0f);
assert(D.v[] == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ]);
{
double[4][3] starray = [
[ 0, 1, 2, 3],
[ 4, 5, 6, 7,],
[ 8, 9, 10, 11,],
];
// starray has the shape 3x4
assert(starray.length == 3);
assert(starray[0].length == 4);
auto m = mat3x4!double(starray);
assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]);
}
{
auto dynarray = [
[ 0, 1, 2, 3],
[ 4, 5, 6, 7,],
[ 8, 9, 10, 11,],
];
// dynarray has the shape 3x4
assert(dynarray.length == 3);
assert(dynarray[0].length == 4);
auto m = mat3x4!double(dynarray);
assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]);
}
}
// Issue #206 (matrix *= scalar) not yielding matrix * scalar but matrix * matrix(scalar)
unittest
{
mat4f mvp = mat4f.identity;
mvp *= 2;
assert(mvp == mat4f(2, 0, 0, 0,
0, 2, 0, 0,
0, 0, 2, 0,
0, 0, 0, 2));
mvp = mat4f.identity * 2;
assert(mvp == mat4f(2, 0, 0, 0,
0, 2, 0, 0,
0, 0, 2, 0,
0, 0, 0, 2));
mvp = mat4f(1) * mat4f(1);
assert(mvp == mat4f(4, 4, 4, 4,
4, 4, 4, 4,
4, 4, 4, 4,
4, 4, 4, 4));
mvp = mat4f(1);
mvp *= mat4f(1);
assert(mvp == mat4f(4, 4, 4, 4,
4, 4, 4, 4,
4, 4, 4, 4,
4, 4, 4, 4));
}

9
external/gfm/math/package.d vendored
View File

@ -1,9 +0,0 @@
module gfm.math;
public import gfm.math.funcs,
gfm.math.vector,
gfm.math.box,
gfm.math.matrix,
gfm.math.quaternion,
gfm.math.shapes,
gfm.math.simplerng;

339
external/gfm/math/quaternion.d vendored
View File

@ -1,339 +0,0 @@
///
module gfm.math.quaternion;
import std.math,
std.string;
import gfm.math.vector,
gfm.math.matrix,
funcs = gfm.math.funcs;
/// Quaternion implementation.
/// Holds a rotation + angle in a proper but wild space.
struct Quaternion(T)
{
public
{
union
{
Vector!(T, 4u) v;
struct
{
T x, y, z, w;
}
}
/// Construct a Quaternion from a value.
@nogc this(U)(U x) pure nothrow if (isAssignable!U)
{
opAssign!U(x);
}
/// Constructs a Quaternion from coordinates.
/// Warning: order of coordinates is different from storage.
@nogc this(T qw, T qx, T qy, T qz) pure nothrow
{
x = qx;
y = qy;
z = qz;
w = qw;
}
/// Constructs a Quaternion from axis + angle.
@nogc static Quaternion fromAxis(Vector!(T, 3) axis, T angle) pure nothrow
{
Quaternion q = void;
axis.normalize();
T cos_a = cos(angle / 2);
T sin_a = sin(angle / 2);
q.x = sin_a * axis.x;
q.y = sin_a * axis.y;
q.z = sin_a * axis.z;
q.w = cos_a;
return q; // should be normalized
}
/// Constructs a Quaternion from Euler angles.
/// All paramers given in radians.
/// Roll->X axis, Pitch->Y axis, Yaw->Z axis
/// See_also: $(LINK https://www.cs.princeton.edu/~gewang/projects/darth/stuff/quat_faq.html)
@nogc static Quaternion fromEulerAngles(T roll, T pitch, T yaw) pure nothrow
{
Quaternion q = void;
T sinPitch = sin(pitch / 2);
T cosPitch = cos(pitch / 2);
T sinYaw = sin(yaw / 2);
T cosYaw = cos(yaw / 2);
T sinRoll = sin(roll / 2);
T cosRoll = cos(roll / 2);
T cosPitchCosYaw = cosPitch * cosYaw;
T sinPitchSinYaw = sinPitch * sinYaw;
q.x = sinRoll * cosPitchCosYaw - cosRoll * sinPitchSinYaw;
q.y = cosRoll * sinPitch * cosYaw + sinRoll * cosPitch * sinYaw;
q.z = cosRoll * cosPitch * sinYaw - sinRoll * sinPitch * cosYaw;
q.w = cosRoll * cosPitchCosYaw + sinRoll * sinPitchSinYaw;
return q;
}
/// Converts a quaternion to Euler angles.
/// TODO: adds a EulerAngles type.
/// Returns: A vector which contains roll-pitch-yaw as x-y-z.
@nogc vec3!T toEulerAngles() pure const nothrow
{
mat3!T m = cast(mat3!T)(this);
T pitch, yaw, roll;
T s = sqrt(m.c[0][0] * m.c[0][0] + m.c[1][0] * m.c[1][0]);
if (s > T.epsilon)
{
pitch = - atan2(m.c[2][0], s);
yaw = atan2(m.c[1][0], m.c[0][0]);
roll = atan2(m.c[2][1], m.c[2][2]);
}
else
{
pitch = m.c[2][0] < 0.0f ? T(PI) /2 : -T(PI) / 2;
yaw = -atan2(m.c[0][1], m.c[1][1]);
roll = 0.0f;
}
return vec3!T(roll, pitch, yaw);
}
/// Assign from another Quaternion.
@nogc ref Quaternion opAssign(U)(U u) pure nothrow if (isQuaternionInstantiation!U && is(U._T : T))
{
v = u.v;
return this;
}
/// Assign from a vector of 4 elements.
@nogc ref Quaternion opAssign(U)(U u) pure nothrow if (is(U : Vector!(T, 4u)))
{
v = u;
return this;
}
/// Converts to a pretty string.
string toString() const nothrow
{
try
return format("%s", v);
catch (Exception e)
assert(false); // should not happen since format is right
}
/// Normalizes a quaternion.
@nogc void normalize() pure nothrow
{
v.normalize();
}
/// Returns: Normalized quaternion.
@nogc Quaternion normalized() pure const nothrow
{
Quaternion res = void;
res.v = v.normalized();
return res;
}
/// Inverses a quaternion in-place.
@nogc void inverse() pure nothrow
{
x = -x;
y = -y;
z = -z;
}
/// Returns: Inverse of quaternion.
@nogc Quaternion inversed() pure const nothrow
{
Quaternion res = void;
res.v = v;
res.inverse();
return res;
}
@nogc ref Quaternion opOpAssign(string op, U)(U q) pure nothrow
if (is(U : Quaternion) && (op == "*"))
{
T nx = w * q.x + x * q.w + y * q.z - z * q.y,
ny = w * q.y + y * q.w + z * q.x - x * q.z,
nz = w * q.z + z * q.w + x * q.y - y * q.x,
nw = w * q.w - x * q.x - y * q.y - z * q.z;
x = nx;
y = ny;
z = nz;
w = nw;
return this;
}
@nogc ref Quaternion opOpAssign(string op, U)(U operand) pure nothrow if (isConvertible!U)
{
Quaternion conv = operand;
return opOpAssign!op(conv);
}
@nogc Quaternion opBinary(string op, U)(U operand) pure const nothrow
if (is(U: Quaternion) || (isConvertible!U))
{
Quaternion temp = this;
return temp.opOpAssign!op(operand);
}
/// Compare two Quaternions.
bool opEquals(U)(U other) pure const if (is(U : Quaternion))
{
return v == other.v;
}
/// Compare Quaternion and other types.
bool opEquals(U)(U other) pure const nothrow if (isConvertible!U)
{
Quaternion conv = other;
return opEquals(conv);
}
/// Convert to a 3x3 rotation matrix.
/// TODO: check out why we can't do is(Unqual!U == mat3!T)
@nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U
&& is(U._T : _T)
&& (U._R == 3) && (U._C == 3))
{
// do not assume that quaternion is normalized
T norm = x*x + y*y + z*z + w*w;
T s = (norm > 0) ? 2 / norm : 0;
T xx = x * x * s, xy = x * y * s, xz = x * z * s, xw = x * w * s,
yy = y * y * s, yz = y * z * s, yw = y * w * s,
zz = z * z * s, zw = z * w * s;
return mat3!(U._T)
(
1 - (yy + zz) , (xy - zw) , (xz + yw) ,
(xy + zw) , 1 - (xx + zz) , (yz - xw) ,
(xz - yw) , (yz + xw) , 1 - (xx + yy)
);
}
/// Converts a to a 4x4 rotation matrix.
/// Bugs: check why we can't do is(Unqual!U == mat4!T)
@nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U
&& is(U._T : _T)
&& (U._R == 4) && (U._C == 4))
{
auto m3 = cast(mat3!(U._T))(this);
return cast(U)(m3);
}
/// Workaround Vector not being constructable through CTFE
@nogc static Quaternion identity() pure nothrow @property
{
Quaternion q;
q.x = q.y = q.z = 0;
q.w = 1;
return q;
}
}
private
{
alias T _T;
template isAssignable(T)
{
enum bool isAssignable =
is(typeof(
{
T x;
Quaternion v = x;
}()));
}
// define types that can be converted to Quaternion, but are not Quaternion
template isConvertible(T)
{
enum bool isConvertible = (!is(T : Quaternion)) && isAssignable!T;
}
}
}
template isQuaternionInstantiation(U)
{
private static void isQuaternion(T)(Quaternion!T x)
{
}
enum bool isQuaternionInstantiation = is(typeof(isQuaternion(U.init)));
}
alias Quaternion!float quatf;///
alias Quaternion!double quatd;///
/// Linear interpolation, for quaternions.
@nogc Quaternion!T lerp(T)(Quaternion!T a, Quaternion!T b, float t) pure nothrow
{
Quaternion!T res = void;
res.v = funcs.lerp(a.v, b.v, t);
return res;
}
/// Nlerp of quaternions
/// Returns: Nlerp of quaternions.
/// See_also: $(WEB keithmaggio.wordpress.com/2011/02/15/math-magician-lerp-slerp-and-nlerp/, Math Magician – Lerp, Slerp, and Nlerp)
@nogc Quaternion!T Nlerp(T)(Quaternion!T a, Quaternion!T b, float t) pure nothrow
{
assert(t >= 0 && t <= 1); // else probably doesn't make sense
Quaternion!T res = void;
res.v = funcs.lerp(a.v, b.v, t);
res.v.normalize();
return res;
}
/// Slerp of quaternions
/// Returns: Slerp of quaternions. Slerp is more expensive than Nlerp.
/// See_also: "Understanding Slerp, Then Not Using It"
@nogc Quaternion!T slerp(T)(Quaternion!T a, Quaternion!T b, T t) pure nothrow
{
assert(t >= 0 && t <= 1); // else probably doesn't make sense
Quaternion!T res = void;
T dotProduct = dot(a.v, b.v);
// spherical lerp always has 2 potential paths
// here we always take the shortest
if (dotProduct < 0)
{
b.v *= -1;
dotProduct = dot(a.v, b.v);
}
immutable T threshold = 10 * T.epsilon; // idMath uses 1e-6f for 32-bits float precision
if ((1 - dotProduct) > threshold) // if small difference, use lerp
return lerp(a, b, t);
T theta_0 = funcs.safeAcos(dotProduct); // angle between this and other
T theta = theta_0 * t; // angle between this and result
vec3!T v2 = dot(b.v, a.v * dotProduct);
v2.normalize();
res.v = dot(b.v, a.v * dotProduct);
res.v.normalize();
res.v = a.v * cos(theta) + res.v * sin(theta);
return res;
}
unittest
{
quatf a = quatf.fromAxis(vec3f(1, 0, 0), 1);
quatf b = quatf.fromAxis(vec3f(0, 1, 0), 0);
a = a * b;
quatf c = lerp(a, b, 0.5f);
quatf d = Nlerp(a, b, 0.1f);
quatf e = slerp(a, b, 0.0f);
quatd f = quatd(1.0, 4, 5.0, 6.0);
quatf g = quatf.fromEulerAngles(-0.1f, 1.2f, -0.3f);
vec3f ga = g.toEulerAngles();
}

645
external/gfm/math/shapes.d vendored
View File

@ -1,645 +0,0 @@
/**
This module implements some abstract geometric shapes:
$(UL
$(LI Line segments.)
$(LI Triangle.)
$(LI Circles/spheres.)
$(LI Rays)
$(LI Planes)
$(LI Frustum)
)
*/
module gfm.math.shapes;
import std.math,
std.traits;
import gfm.math.vector,
gfm.math.box;
/// A Segment is 2 points.
/// When considered like a vector, it represents the arrow from a to b.
struct Segment(T, int N)
{
public
{
alias Vector!(T, N) point_t;
point_t a, b;
static if (N == 3 && isFloatingPoint!T)
{
/// Segment vs plane intersection.
/// See_also: "Geometry for Computer Graphics: Formulae, Examples
/// and Proofs", Vince (2005), p. 62
/// Returns: $(D true) if the segment crosses the plane exactly
/// once, or $(D false) if the segment is parallel to or does
/// not reach the plane
/// Params:
/// plane = The plane to intersect.
/// intersection = Point of intersection.
/// progress = Set to the point's progress between endpoints
/// $(D a) at 0.0 and $(D b) at 1.0, or T.infinity
/// if parallel to the plane
@nogc bool intersect(Plane!T plane, out point_t intersection, out T progress) pure const nothrow
{
// direction vector from a to b
point_t dir = b-a;
// dot product will be 0 if angle to plane is 0
T dp = dot(plane.n, dir);
if (abs(dp) < T.epsilon)
{
progress = T.infinity;
return false; // parallel to plane
}
// relative distance along segment where intersection happens
progress = -(dot(plane.n, a) - plane.d) / dp;
intersection = progress*dir + a;
return progress >= 0 && progress <= 1;
}
}
}
}
alias Segment!(float, 2) seg2f; /// 2D float segment.
alias Segment!(float, 3) seg3f; /// 3D float segment.
alias Segment!(double, 2) seg2d; /// 2D double segment.
alias Segment!(double, 3) seg3d; /// 3D double segment.
alias Segment!(int, 2) seg2i; /// 2D integer segment.
alias Segment!(int, 3) seg3i; /// 3D integer segment.
/// A Triangle is 3 points.
struct Triangle(T, int N)
{
public
{
alias Vector!(T, N) point_t;
point_t a, b, c;
static if (N == 2)
{
/// Returns: Area of a 2D triangle.
@nogc T area() pure const nothrow
{
return abs(signedArea());
}
/// Returns: Signed area of a 2D triangle.
@nogc T signedArea() pure const nothrow
{
return ((b.x * a.y - a.x * b.y)
+ (c.x * b.y - b.x * c.y)
+ (a.x * c.y - c.x * a.y)) / 2;
}
}
static if (N == 3 && isFloatingPoint!T)
{
/// Returns: Triangle normal.
@nogc Vector!(T, 3) computeNormal() pure const nothrow
{
return cross(b - a, c - a).normalized();
}
}
}
}
alias Triangle!(float, 2) triangle2f; /// 2D float triangle.
alias Triangle!(float, 3) triangle3f; /// 3D float triangle.
alias Triangle!(double, 2) triangle2d; /// 2D double triangle.
alias Triangle!(double, 3) triangle3d; /// 3D double triangle.
/// A Sphere is a point + a radius.
struct Sphere(T, int N)
{
public nothrow
{
alias Vector!(T, N) point_t;
point_t center;
T radius;
/// Creates a sphere from a point and a radius.
@nogc this(in point_t center_, T radius_) pure nothrow
{
center = center_;
radius = radius_;
}
/// Sphere contains point test.
/// Returns: true if the point is inside the sphere.
@nogc bool contains(in Sphere s) pure const nothrow
{
if (s.radius > radius)
return false;
T innerRadius = radius - s.radius;
return squaredDistanceTo(s.center) < innerRadius * innerRadius;
}
/// Sphere vs point Euclidean distance squared.
@nogc T squaredDistanceTo(point_t p) pure const nothrow
{
return center.squaredDistanceTo(p);
}
/// Sphere vs sphere intersection.
/// Returns: true if the spheres intersect.
@nogc bool intersects(Sphere s) pure const nothrow
{
T outerRadius = radius + s.radius;
return squaredDistanceTo(s.center) < outerRadius * outerRadius;
}
static if (isFloatingPoint!T)
{
/// Sphere vs point Euclidean distance.
@nogc T distanceTo(point_t p) pure const nothrow
{
return center.distanceTo(p);
}
static if(N == 2)
{
/// Returns: Circle area.
@nogc T area() pure const nothrow
{
return T(PI) * (radius * radius);
}
}
}
}
}
alias Sphere!(float, 2) sphere2f; /// 2D float sphere (ie. a circle).
alias Sphere!(float, 3) sphere3f; /// 3D float sphere.
alias Sphere!(double, 2) sphere2d; /// 2D double sphere (ie. a circle).
alias Sphere!(double, 3) sphere3d; /// 3D double sphere (ie. a circle).
/// A Ray ir a point + a direction.
struct Ray(T, int N)
{
nothrow:
public
{
alias Vector!(T, N) point_t;
point_t orig;
point_t dir;
/// Returns: A point further along the ray direction.
@nogc point_t progress(T t) pure const nothrow
{
return orig + dir * t;
}
static if (N == 3 && isFloatingPoint!T)
{
/// Ray vs triangle intersection.
/// See_also: "Fast, Minimum Storage Ray/Triangle intersection", Mommer & Trumbore (1997)
/// Returns: Barycentric coordinates, the intersection point is at $(D (1 - u - v) * A + u * B + v * C).
@nogc bool intersect(Triangle!(T, 3) triangle, out T t, out T u, out T v) pure const nothrow
{
point_t edge1 = triangle.b - triangle.a;
point_t edge2 = triangle.c - triangle.a;
point_t pvec = cross(dir, edge2);
T det = dot(edge1, pvec);
if (abs(det) < T.epsilon)
return false; // no intersection
T invDet = 1 / det;
// calculate distance from triangle.a to ray origin
point_t tvec = orig - triangle.a;
// calculate U parameter and test bounds
u = dot(tvec, pvec) * invDet;
if (u < 0 || u > 1)
return false;
// prepare to test V parameter
point_t qvec = cross(tvec, edge1);
// calculate V parameter and test bounds
v = dot(dir, qvec) * invDet;
if (v < 0.0 || u + v > 1.0)
return false;
// calculate t, ray intersects triangle
t = dot(edge2, qvec) * invDet;
return true;
}
/// Ray vs plane intersection.
/// See_also: "Geometry for Computer Graphics: Formulae, Examples
/// and Proofs", Vince (2005), p. 62
/// Returns: $(D true) if the ray crosses the plane exactly once,
/// or $(D false) if the ray is parallel to or points away from
/// the plane
/// Params:
/// plane = Plane to intersect.
/// intersection = Point of intersection.
/// distance = set to the point's distance along the ray, or
/// T.infinity if parallel to the plane
@nogc bool intersect(Plane!T plane, out point_t intersection, out T distance) pure const nothrow
{
// dot product will be 0 if angle to plane is 0
T dp = dot(plane.n, dir);
if (abs(dp) < T.epsilon)
{
distance = T.infinity;
return false; // parallel to plane
}
// distance along ray where intersection happens
distance = -(dot(plane.n, orig) - plane.d) / dp;
intersection = distance*dir + orig;
return distance >= 0;
}
}
}
}
alias Ray!(float, 2) ray2f; /// 2D float ray.
alias Ray!(float, 3) ray3f; /// 3D float ray.
alias Ray!(double, 2) ray2d; /// 2D double ray.
alias Ray!(double, 3) ray3d; /// 3D double ray.
/// 3D plane.
/// See_also: Flipcode article by Nate Miller $(WEB www.flipcode.com/archives/Plane_Class.shtml).
struct Plane(T) if (isFloatingPoint!T)
{
public
{
vec3!T n; /// Normal (always stored normalized).
T d;
/// Create from four coordinates.
@nogc this(vec4!T abcd) pure nothrow
{
n = vec3!T(abcd.x, abcd.y, abcd.z).normalized();
d = abcd.w;
}
/// Create from a point and a normal.
@nogc this(vec3!T origin, vec3!T normal) pure nothrow
{
n = normal.normalized();
d = -dot(origin, n);
}
/// Create from 3 non-aligned points.
@nogc this(vec3!T A, vec3!T B, vec3!T C) pure nothrow
{
this(C, cross(B - A, C - A));
}
/// Assign a plane with another plane.
@nogc ref Plane opAssign(Plane other) pure nothrow
{
n = other.n;
d = other.d;
return this;
}
/// Returns: signed distance between a point and the plane.
@nogc T signedDistanceTo(vec3!T point) pure const nothrow
{
return dot(n, point) + d;
}
/// Returns: absolute distance between a point and the plane.
@nogc T distanceTo(vec3!T point) pure const nothrow
{
return abs(signedDistanceTo(point));
}
/// Returns: true if the point is in front of the plane.
@nogc bool isFront(vec3!T point) pure const nothrow
{
return signedDistanceTo(point) >= 0;
}
/// Returns: true if the point is in the back of the plane.
@nogc bool isBack(vec3!T point) pure const nothrow
{
return signedDistanceTo(point) < 0;
}
/// Returns: true if the point is on the plane, with a given epsilon.
@nogc bool isOn(vec3!T point, T epsilon) pure const nothrow
{
T sd = signedDistanceTo(point);
return (-epsilon < sd) && (sd < epsilon);
}
}
}
alias Plane!float planef; /// 3D float plane.
alias Plane!double planed; /// 3D double plane.
unittest
{
auto p = planed(vec4d(1.0, 2.0, 3.0, 4.0));
auto p2 = planed(vec3d(1.0, 0.0, 0.0),
vec3d(0.0, 1.0, 0.0),
vec3d(0.0, 0.0, 1.0));
assert(p2.isOn(vec3d(1.0, 0.0, 0.0), 1e-7));
assert(p2.isFront(vec3d(1.0, 1.0, 1.0)));
assert(p2.isBack(vec3d(0.0, 0.0, 0.0)));
}
/// Plane intersection tests
@nogc pure nothrow unittest
{
void testR(planed p, ray3d r, bool shouldIntersect, double expectedDistance, vec3d expectedPoint = vec3d.init) pure nothrow @nogc
{
vec3d point;
double distance;
assert(r.intersect(p, point, distance) == shouldIntersect);
assert(approxEqual(distance, expectedDistance));
if (shouldIntersect)
assert(approxEqual(point.v[], expectedPoint.v[]));
}
// ray facing plane
testR(planed(vec4d(1.0, 0.0, 0.0, 1.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(-1.0, 0.0, 0.0)),
true, 1.0, vec3d(1.0, 3.0, 4.0));
testR(planed(vec4d(1.0, 1.0, 1.0, -sqrt(3.0))), ray3d(vec3d(1.0, 1.0, 1.0), vec3d(-1.0, -1.0, -1.0).normalized()),
true, 2.0*sqrt(3.0), vec3d(-1.0, -1.0, -1.0));
// ray facing away
testR(planed(vec4d(1.0, 0.0, 0.0, 1.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(1.0, 0.0, 0.0)),
false, -1.0);
testR(planed(vec4d(1.0, 0.0, 0.0, 5.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(-1.0, 0.0, 0.0)),
false, -3.0);
// ray parallel
testR(planed(vec4d(0.0, 0.0, 1.0, 1.0)), ray3d(vec3d(1.0, 2.0, 3.0), vec3d(1.0, 0.0, 0.0)),
false, double.infinity);
void testS(planed p, seg3d s, bool shouldIntersect, double expectedProgress, vec3d expectedPoint = vec3d.init) pure nothrow @nogc
{
vec3d point;
double progress;
assert(s.intersect(p, point, progress) == shouldIntersect);
assert(approxEqual(progress, expectedProgress));
if (shouldIntersect)
assert(approxEqual(point.v[], expectedPoint.v[]));
}
// segment crossing plane
testS(planed(vec4d(1.0, 0.0, 0.0, 2.0)), seg3d(vec3d(1.0, 2.0, 3.0), vec3d(3.0, 4.0, 5.0)),
true, 0.5, vec3d(2.0, 3.0, 4.0));
// segment too short
testS(planed(vec4d(1.0, 0.0, 0.0, 0.0)), seg3d(vec3d(1.0, 2.0, 3.0), vec3d(3.0, 4.0, 5.0)),
false, -0.5);
}
/// 3D frustum.
/// See_also: Flipcode article by Dion Picco $(WEB www.flipcode.com/archives/Frustum_Culling.shtml).
/// Bugs: verify proper signedness of half-spaces
struct Frustum(T) if (isFloatingPoint!T)
{
public
{
enum int LEFT = 0,
RIGHT = 1,
TOP = 2,
BOTTOM = 3,
NEAR = 4,
FAR = 5;
Plane!T[6] planes;
/// Create a frustum from 6 planes.
@nogc this(Plane!T left, Plane!T right, Plane!T top, Plane!T bottom, Plane!T near, Plane!T far) pure nothrow
{
planes[LEFT] = left;
planes[RIGHT] = right;
planes[TOP] = top;
planes[BOTTOM] = bottom;
planes[NEAR] = near;
planes[FAR] = far;
}
enum : int
{
OUTSIDE, /// object is outside the frustum
INTERSECT, /// object intersects with the frustum
INSIDE /// object is inside the frustum
}
/// Point vs frustum intersection.
@nogc bool contains(vec3!T point) pure const nothrow
{
for(int i = 0; i < 6; ++i)
{
T distance = planes[i].signedDistanceTo(point);
if(distance < 0)
return false;
}
return true;
}
/// Sphere vs frustum intersection.
/// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE.
@nogc int contains(Sphere!(T, 3) sphere) pure const nothrow
{
// calculate our distances to each of the planes
for(int i = 0; i < 6; ++i)
{
// find the distance to this plane
T distance = planes[i].signedDistanceTo(sphere.center);
if(distance < -sphere.radius)
return OUTSIDE;
else if (distance < sphere.radius)
return INTERSECT;
}
// otherwise we are fully in view
return INSIDE;
}
/// AABB vs frustum intersection.
/// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE.
@nogc int contains(box3!T box) pure const nothrow
{
vec3!T[8] corners;
int totalIn = 0;
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k)
{
auto x = i == 0 ? box.min.x : box.max.x;
auto y = j == 0 ? box.min.y : box.max.y;
auto z = k == 0 ? box.min.z : box.max.z;
corners[i*4 + j*2 + k] = vec3!T(x, y, z);
}
// test all 8 corners against the 6 sides
// if all points are behind 1 specific plane, we are out
// if we are in with all points, then we are fully in
for(int p = 0; p < 6; ++p)
{
int inCount = 8;
int ptIn = 1;
for(int i = 0; i < 8; ++i)
{
// test this point against the planes
if (planes[p].isBack(corners[i]))
{
ptIn = 0;
--inCount;
}
}
// were all the points outside of plane p?
if (inCount == 0)
return OUTSIDE;
// check if they were all on the right side of the plane
totalIn += ptIn;
}
// so if totalIn is 6, then all are inside the view
if(totalIn == 6)
return INSIDE;
// we must be partly in then otherwise
return INTERSECT;
}
}
}
unittest
{
seg2f se;
triangle3f tr;
Frustum!double frust;
planed pl;
}
/// True if `T` is a kind of Segment
enum isSegment(T) = is(T : Segment!U, U...);
/// True if `T` is a kind of Triangle
enum isTriangle(T) = is(T : Triangle!U, U...);
/// True if `T` is a kind of Sphere
enum isSphere(T) = is(T : Sphere!U, U...);
/// True if `T` is a kind of Ray
enum isRay(T) = is(T : Ray!U, U...);
/// True if `T` is a kind of Plane
enum isPlane(T) = is(T : Plane!U, U);
/// True if `T` is a kind of Frustum
enum isFrustum(T) = is(T : Frustum!U, U);
/// True if `T` is a kind of 2 dimensional Segment
enum isSegment2D(T) = is(T : Segment!(U, 2), U);
/// True if `T` is a kind of 2 dimensional Triangle
enum isTriangle2D(T) = is(T : Triangle!(U, 2), U);
/// True if `T` is a kind of 2 dimensional Sphere
enum isSphere2D(T) = is(T : Sphere!(U, 2), U);
/// True if `T` is a kind of 2 dimensional Ray
enum isRay2D(T) = is(T : Ray!(U, 2), U);
/// True if `T` is a kind of 3 dimensional Segment
enum isSegment3D(T) = is(T : Segment!(U, 3), U);
/// True if `T` is a kind of 3 dimensional Triangle
enum isTriangle3D(T) = is(T : Triangle!(U, 3), U);
/// True if `T` is a kind of 3 dimensional Sphere
enum isSphere3D(T) = is(T : Sphere!(U, 3), U);
/// True if `T` is a kind of 3 dimensional Ray
enum isRay3D(T) = is(T : Ray!(U, 3), U);
unittest
{
enum ShapeType
{
segment,
triangle,
sphere,
ray,
plane,
frustum
}
void test(T, ShapeType type, int dim)()
{
with (ShapeType)
{
static assert(isSegment!T == (type == segment ));
static assert(isTriangle!T == (type == triangle));
static assert(isSphere!T == (type == sphere ));
static assert(isRay!T == (type == ray ));
static assert(isPlane!T == (type == plane ));
static assert(isFrustum!T == (type == frustum ));
static assert(isSegment2D!T == (type == segment && dim == 2));
static assert(isTriangle2D!T == (type == triangle && dim == 2));
static assert(isSphere2D!T == (type == sphere && dim == 2));
static assert(isRay2D!T == (type == ray && dim == 2));
static assert(isSegment3D!T == (type == segment && dim == 3));
static assert(isTriangle3D!T == (type == triangle && dim == 3));
static assert(isSphere3D!T == (type == sphere && dim == 3));
static assert(isRay3D!T == (type == ray && dim == 3));
}
}
with (ShapeType)
{
// test case type #dimensions
test!(seg2f , segment , 2);
test!(seg3d , segment , 3);
test!(triangle2d , triangle, 2);
test!(triangle3f , triangle, 3);
test!(sphere2d , sphere , 2);
test!(Sphere!(uint, 3), sphere , 3);
test!(ray2f , ray , 2);
test!(Ray!(real, 3) , ray , 3);
test!(planed , plane , 0); // ignore dimension (always 3D)
test!(Plane!float , plane , 0);
test!(Frustum!double , frustum , 0);
}
}
/// Get the numeric type used to measure a shape's dimensions.
alias DimensionType(T : Segment!U, U...) = U[0];
/// ditto
alias DimensionType(T : Triangle!U, U...) = U[0];
/// ditto
alias DimensionType(T : Sphere!U, U...) = U[0];
/// ditto
alias DimensionType(T : Ray!U, U...) = U[0];
/// ditto
alias DimensionType(T : Plane!U, U) = U;
/// ditto
alias DimensionType(T : Frustum!U, U) = U;
///
unittest
{
static assert(is(DimensionType!seg2i == int));
static assert(is(DimensionType!triangle3d == double));
static assert(is(DimensionType!sphere2d == double));
static assert(is(DimensionType!ray3f == float));
static assert(is(DimensionType!planed == double));
static assert(is(DimensionType!(Frustum!real) == real));
}

485
external/gfm/math/simplerng.d vendored
View File

@ -1,485 +0,0 @@
/**
Translation of SimpleRNG.
Removed the builtin RNG to use std.random, but kept the distribution functions.
John D. Cook confirmed this code as public domain.
Authors: John D. Cook.
See_also: $(WEB www.johndcook.com/cpp_random_number_generation.html)
*/
module gfm.math.simplerng;
public import std.random;
import std.math;
/// Returns: Normal (Gaussian) random sample.
/// See_also: Box-Muller algorithm.
double randNormal(RNG)(ref RNG rng, double mean = 0.0, double standardDeviation = 1.0)
{
assert(standardDeviation > 0);
double u1;
do
{
u1 = uniform01(rng);
} while (u1 == 0); // u1 must not be zero
double u2 = uniform01(rng);
double r = sqrt(-2.0 * log(u1));
double theta = 2.0 * double(PI) * u2;
return mean + standardDeviation * r * sin(theta);
}
/// Returns: Exponential random sample with specified mean.
double randExponential(RNG)(ref RNG rng, double mean = 1.0)
{
assert(mean > 0);
return -mean*log(uniform01(rng));
}
/// Returns: Gamma random sample.
/// See_also: "A Simple Method for Generating Gamma Variables"
/// by George Marsaglia and Wai Wan Tsang. ACM Transactions on Mathematical Software
/// Vol 26, No 3, September 2000, pages 363-372.
double randGamma(RNG)(ref RNG rng, double shape, double scale)
{
double d, c, x, xsquared, v, u;
if (shape >= 1.0)
{
d = shape - 1.0/3.0;
c = 1.0/sqrt(9.0*d);
for (;;)
{
do
{
x = randNormal(rng);
v = 1.0 + c*x;
}
while (v <= 0.0);
v = v*v*v;
u = uniform01(rng);
xsquared = x*x;
if (u < 1.0 -.0331*xsquared*xsquared || log(u) < 0.5*xsquared + d*(1.0 - v + log(v)))
return scale*d*v;
}
}
else
{
assert(shape > 0);
double g = randGamma(rng, shape+1.0, 1.0);
double w = uniform01(rng);
return scale*g*pow(w, 1.0/shape);
}
}
/// Returns: Chi-square sample.
double randChiSquare(RNG)(ref RNG rng, double degreesOfFreedom)
{
// A chi squared distribution with n degrees of freedom
// is a gamma distribution with shape n/2 and scale 2.
return randGamma(rng, 0.5 * degreesOfFreedom, 2.0);
}
/// Returns: Inverse-gamma sample.
double randInverseGamma(RNG)(ref RNG rng, double shape, double scale)
{
// If X is gamma(shape, scale) then
// 1/Y is inverse gamma(shape, 1/scale)
return 1.0 / randGamma(rng, shape, 1.0 / scale);
}
/// Returns: Weibull sample.
double randWeibull(RNG)(ref RNG rng, double shape, double scale)
{
assert(shape > 0 && scale > 0);
return scale * pow(-log(uniform01(rng)), 1.0 / shape);
}
/// Returns: Cauchy sample.
double randCauchy(RNG)(ref RNG rng, double median, double scale)
{
assert(scale > 0);
double p = uniform01(rng);
// Apply inverse of the Cauchy distribution function to a uniform
return median + scale*tan(double(PI)*(p - 0.5));
}
/// Returns: Student-t sample.
/// See_also: Seminumerical Algorithms by Knuth.
double randStudentT(RNG)(ref RNG rng, double degreesOfFreedom)
{
assert(degreesOfFreedom > 0);
double y1 = getNormal(rng);
double y2 = getChiSquare(rng, degreesOfFreedom);
return y1 / sqrt(y2 / degreesOfFreedom);
}
/// Returns: Laplace distribution random sample (also known as the double exponential distribution).
double randLaplace(RNG)(ref RNG rng, double mean, double scale)
{
double u = uniform01(rng);
return (u < 0.5) ? (mean + scale*log(2.0*u))
: (mean - scale*log(2*(1-u)));
}
/// Returns: Log-normal sample.
double randLogNormal(RNG)(ref RNG rng, double mu, double sigma)
{
return exp(getNormal(rng, mu, sigma));
}
/// Returns: Beta sample.
double randBeta(RNG)(ref RNG rng, double a, double b)
{
assert(a > 0 && b > 0);
// There are more efficient methods for generating beta samples.
// However such methods are a little more efficient and much more complicated.
// For an explanation of why the following method works, see
// http://www.johndcook.com/distribution_chart.html#gamma_beta
double u = getGamma(rng, a, 1.0);
double v = getGamma(rng, b, 1.0);
return u / (u + v);
}
/// Returns: Poisson sample.
int randPoisson(RNG)(ref RNG rng, double lambda)
{
return (lambda < 30.0) ? poissonSmall(rng, lambda) : poissonLarge(rng, lambda);
}
private
{
int poissonSmall(RNG)(ref RNG rng, double lambda)
{
// Algorithm due to Donald Knuth, 1969.
double p = 1.0, L = exp(-lambda);
int k = 0;
do
{
k++;
p *= uniform01(rng);
}
while (p > L);
return k - 1;
}
int poissonLarge(RNG)(ref RNG rng, double lambda)
{
// "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson
// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
// The article is on pages 29-35. The algorithm given here is on page 32.
double c = 0.767 - 3.36/lambda;
double beta = double(PI)/sqrt(3.0*lambda);
double alpha = beta*lambda;
double k = log(c) - lambda - log(beta);
for(;;)
{
double u = uniform01(rng);
double x = (alpha - log((1.0 - u)/u))/beta;
int n = cast(int)(floor(x + 0.5));
if (n < 0)
continue;
double v = uniform01(rng);
double y = alpha - beta*x;
double temp = 1.0 + exp(y);
double lhs = y + log(v/(temp*temp));
double rhs = k + n*log(lambda) - logFactorial(n);
if (lhs <= rhs)
return n;
}
}
double logFactorial(int n) nothrow
{
assert(n >= 0);
if (n > 254)
{
double x = n + 1;
return (x - 0.5)*log(x) - x + 0.5*log(2*double(PI)) + 1.0/(12.0*x);
}
else
{
return LOG_FACTORIAL[n];
}
}
}
private static immutable double[255] LOG_FACTORIAL =
[
0.000000000000000,
0.000000000000000,
0.693147180559945,
1.791759469228055,
3.178053830347946,
4.787491742782046,
6.579251212010101,
8.525161361065415,
10.604602902745251,
12.801827480081469,
15.104412573075516,
17.502307845873887,
19.987214495661885,
22.552163853123421,
25.191221182738683,
27.899271383840894,
30.671860106080675,
33.505073450136891,
36.395445208033053,
39.339884187199495,
42.335616460753485,
45.380138898476908,
48.471181351835227,
51.606675567764377,
54.784729398112319,
58.003605222980518,
61.261701761002001,
64.557538627006323,
67.889743137181526,
71.257038967168000,
74.658236348830158,
78.092223553315307,
81.557959456115029,
85.054467017581516,
88.580827542197682,
92.136175603687079,
95.719694542143202,
99.330612454787428,
102.968198614513810,
106.631760260643450,
110.320639714757390,
114.034211781461690,
117.771881399745060,
121.533081515438640,
125.317271149356880,
129.123933639127240,
132.952575035616290,
136.802722637326350,
140.673923648234250,
144.565743946344900,
148.477766951773020,
152.409592584497350,
156.360836303078800,
160.331128216630930,
164.320112263195170,
168.327445448427650,
172.352797139162820,
176.395848406997370,
180.456291417543780,
184.533828861449510,
188.628173423671600,
192.739047287844900,
196.866181672889980,
201.009316399281570,
205.168199482641200,
209.342586752536820,
213.532241494563270,
217.736934113954250,
221.956441819130360,
226.190548323727570,
230.439043565776930,
234.701723442818260,
238.978389561834350,
243.268849002982730,
247.572914096186910,
251.890402209723190,
256.221135550009480,
260.564940971863220,
264.921649798552780,
269.291097651019810,
273.673124285693690,
278.067573440366120,
282.474292687630400,
286.893133295426990,
291.323950094270290,
295.766601350760600,
300.220948647014100,
304.686856765668720,
309.164193580146900,
313.652829949878990,
318.152639620209300,
322.663499126726210,
327.185287703775200,
331.717887196928470,
336.261181979198450,
340.815058870798960,
345.379407062266860,
349.954118040770250,
354.539085519440790,
359.134205369575340,
363.739375555563470,
368.354496072404690,
372.979468885689020,
377.614197873918670,
382.258588773060010,
386.912549123217560,
391.575988217329610,
396.248817051791490,
400.930948278915760,
405.622296161144900,
410.322776526937280,
415.032306728249580,
419.750805599544780,
424.478193418257090,
429.214391866651570,
433.959323995014870,
438.712914186121170,
443.475088120918940,
448.245772745384610,
453.024896238496130,
457.812387981278110,
462.608178526874890,
467.412199571608080,
472.224383926980520,
477.044665492585580,
481.872979229887900,
486.709261136839360,
491.553448223298010,
496.405478487217580,
501.265290891579240,
506.132825342034830,
511.008022665236070,
515.890824587822520,
520.781173716044240,
525.679013515995050,
530.584288294433580,
535.496943180169520,
540.416924105997740,
545.344177791154950,
550.278651724285620,
555.220294146894960,
560.169054037273100,
565.124881094874350,
570.087725725134190,
575.057539024710200,
580.034272767130800,
585.017879388839220,
590.008311975617860,
595.005524249382010,
600.009470555327430,
605.020105849423770,
610.037385686238740,
615.061266207084940,
620.091704128477430,
625.128656730891070,
630.172081847810200,
635.221937855059760,
640.278183660408100,
645.340778693435030,
650.409682895655240,
655.484856710889060,
660.566261075873510,
665.653857411105950,
670.747607611912710,
675.847474039736880,
680.953419513637530,
686.065407301994010,
691.183401114410800,
696.307365093814040,
701.437263808737160,
706.573062245787470,
711.714725802289990,
716.862220279103440,
722.015511873601330,
727.174567172815840,
732.339353146739310,
737.509837141777440,
742.685986874351220,
747.867770424643370,
753.055156230484160,
758.248113081374300,
763.446610112640200,
768.650616799717000,
773.860102952558460,
779.075038710167410,
784.295394535245690,
789.521141208958970,
794.752249825813460,
799.988691788643450,
805.230438803703120,
810.477462875863580,
815.729736303910160,
820.987231675937890,
826.249921864842800,
831.517780023906310,
836.790779582469900,
842.068894241700490,
847.352097970438420,
852.640365001133090,
857.933669825857460,
863.231987192405430,
868.535292100464630,
873.843559797865740,
879.156765776907600,
884.474885770751830,
889.797895749890240,
895.125771918679900,
900.458490711945270,
905.796028791646340,
911.138363043611210,
916.485470574328820,
921.837328707804890,
927.193914982476710,
932.555207148186240,
937.921183163208070,
943.291821191335660,
948.667099599019820,
954.046996952560450,
959.431492015349480,
964.820563745165940,
970.214191291518320,
975.612353993036210,
981.015031374908400,
986.422203146368590,
991.833849198223450,
997.249949600427840,
1002.670484599700300,
1008.095434617181700,
1013.524780246136200,
1018.958502249690200,
1024.396581558613400,
1029.838999269135500,
1035.285736640801600,
1040.736775094367400,
1046.192096209724900,
1051.651681723869200,
1057.115513528895000,
1062.583573670030100,
1068.055844343701400,
1073.532307895632800,
1079.012946818975000,
1084.497743752465600,
1089.986681478622400,
1095.479742921962700,
1100.976911147256000,
1106.478169357800900,
1111.983500893733000,
1117.492889230361000,
1123.006317976526100,
1128.523770872990800,
1134.045231790853000,
1139.570684729984800,
1145.100113817496100,
1150.633503306223700,
1156.170837573242400,
];
unittest
{
Xorshift32 rng;
rng.seed(unpredictableSeed());
double x = randNormal!Xorshift32(rng, 0.0, 1.0);
x = randExponential!Xorshift32(rng);
x = randGamma!Xorshift32(rng, 1.2, 1.0);
x = randGamma!Xorshift32(rng, 0.8, 2.0);
x = randChiSquare!Xorshift32(rng, 2.0);
x = randInverseGamma!Xorshift32(rng, 1.1, 0.7);
x = randWeibull!Xorshift32(rng, 3.0, 0.7);
x = randCauchy!Xorshift32(rng, 5.0, 1.4);
}

770
external/gfm/math/vector.d vendored
View File

@ -1,770 +0,0 @@
/// N-dimension vector mathematical object
module gfm.math.vector;
import std.traits,
std.math,
std.conv,
std.array,
std.string;
import gfm.math.funcs;
/**
* Generic 1D small vector.
* Params:
* N = number of elements
* T = type of elements
*/
struct Vector(T, int N)
{
nothrow:
public
{
static assert(N >= 1);
// fields definition
union
{
T[N] v;
struct
{
static if (N >= 1)
{
T x;
alias x r;
}
static if (N >= 2)
{
T y;
alias y g;
}
static if (N >= 3)
{
T z;
alias z b;
}
static if (N >= 4)
{
T w;
alias w a;
}
}
}
/// Construct a Vector with a `T[]` or the values as arguments
@nogc this(Args...)(Args args) pure nothrow
{
static if (args.length == 1)
{
// Construct a Vector from a single value.
opAssign!(Args[0])(args[0]);
}
else
{
// validate the total argument count across scalars and vectors
template argCount(T...) {
static if(T.length == 0)
enum argCount = 0; // done recursing
else static if(isVector!(T[0]))
enum argCount = T[0]._N + argCount!(T[1..$]);
else
enum argCount = 1 + argCount!(T[1..$]);
}
static assert(argCount!Args <= N, "Too many arguments in vector constructor");
int index = 0;
foreach(arg; args)
{
static if (isAssignable!(T, typeof(arg)))
{
v[index] = arg;
index++; // has to be on its own line (DMD 2.068)
}
else static if (isVector!(typeof(arg)) && isAssignable!(T, arg._T))
{
mixin(generateLoopCode!("v[index + @] = arg[@];", arg._N)());
index += arg._N;
}
else
static assert(false, "Unrecognized argument in Vector constructor");
}
assert(index == N, "Bad arguments in Vector constructor");
}
}
size_t toHash() const nothrow @safe
{
size_t hash = 0;
foreach (elem; v) {
hash = elem.hashOf(hash);
}
return hash;
}
/// Assign a Vector from a compatible type.
@nogc ref Vector opAssign(U)(U x) pure nothrow if (isAssignable!(T, U))
{
mixin(generateLoopCode!("v[@] = x;", N)()); // copy to each component
return this;
}
/// Assign a Vector with a static array type.
@nogc ref Vector opAssign(U)(U arr) pure nothrow if ((isStaticArray!(U) && isAssignable!(T, typeof(arr[0])) && (arr.length == N)))
{
mixin(generateLoopCode!("v[@] = arr[@];", N)());
return this;
}
/// Assign with a dynamic array.
/// Size is checked in debug-mode.
@nogc ref Vector opAssign(U)(U arr) pure nothrow if (isDynamicArray!(U) && isAssignable!(T, typeof(arr[0])))
{
assert(arr.length == N);
mixin(generateLoopCode!("v[@] = arr[@];", N)());
return this;
}
/// Assign from a samey Vector.
@nogc ref Vector opAssign(U)(U u) pure nothrow if (is(U : Vector))
{
v[] = u.v[];
return this;
}
/// Assign from other vectors types (same size, compatible type).
@nogc ref Vector opAssign(U)(U x) pure nothrow if (isVector!U
&& isAssignable!(T, U._T)
&& (!is(U: Vector))
&& (U._N == _N))
{
mixin(generateLoopCode!("v[@] = x.v[@];", N)());
return this;
}
/// Returns: a pointer to content.
@nogc inout(T)* ptr() pure inout nothrow @property
{
return v.ptr;
}
/// Converts to a pretty string.
string toString() const nothrow
{
try
return format("%s", v);
catch (Exception e)
assert(false); // should not happen since format is right
}
@nogc bool opEquals(U)(U other) pure const nothrow
if (is(U : Vector))
{
for (int i = 0; i < N; ++i)
{
if (v[i] != other.v[i])
{
return false;
}
}
return true;
}
@nogc bool opEquals(U)(U other) pure const nothrow
if (isConvertible!U)
{
Vector conv = other;
return opEquals(conv);
}
@nogc Vector opUnary(string op)() pure const nothrow
if (op == "+" || op == "-" || op == "~" || op == "!")
{
Vector res = void;
mixin(generateLoopCode!("res.v[@] = " ~ op ~ " v[@];", N)());
return res;
}
@nogc ref Vector opOpAssign(string op, U)(U operand) pure nothrow
if (is(U : Vector))
{
mixin(generateLoopCode!("v[@] " ~ op ~ "= operand.v[@];", N)());
return this;
}
@nogc ref Vector opOpAssign(string op, U)(U operand) pure nothrow if (isConvertible!U)
{
Vector conv = operand;
return opOpAssign!op(conv);
}
@nogc Vector opBinary(string op, U)(U operand) pure const nothrow
if (is(U: Vector) || (isConvertible!U))
{
Vector result = void;
static if (is(U: T))
mixin(generateLoopCode!("result.v[@] = cast(T)(v[@] " ~ op ~ " operand);", N)());
else
{
Vector other = operand;
mixin(generateLoopCode!("result.v[@] = cast(T)(v[@] " ~ op ~ " other.v[@]);", N)());
}
return result;
}
@nogc Vector opBinaryRight(string op, U)(U operand) pure const nothrow if (isConvertible!U)
{
Vector result = void;
static if (is(U: T))
mixin(generateLoopCode!("result.v[@] = cast(T)(operand " ~ op ~ " v[@]);", N)());
else
{
Vector other = operand;
mixin(generateLoopCode!("result.v[@] = cast(T)(other.v[@] " ~ op ~ " v[@]);", N)());
}
return result;
}
@nogc ref T opIndex(size_t i) pure nothrow
{
return v[i];
}
@nogc ref const(T) opIndex(size_t i) pure const nothrow
{
return v[i];
}
@nogc T opIndexAssign(U : T)(U x, size_t i) pure nothrow
{
return v[i] = x;
}
/// Implements swizzling.
///
/// Example:
/// ---
/// vec4i vi = [4, 1, 83, 10];
/// assert(vi.zxxyw == [83, 4, 4, 1, 10]);
/// ---
@nogc @property auto opDispatch(string op, U = void)() pure const nothrow if (isValidSwizzle!(op))
{
alias Vector!(T, op.length) returnType;
returnType res = void;
enum indexTuple = swizzleTuple!op;
foreach(i, index; indexTuple)
res.v[i] = v[index];
return res;
}
/// Support swizzling assignment like in shader languages.
///
/// Example:
/// ---
/// vec3f v = [0, 1, 2];
/// v.yz = v.zx;
/// assert(v == [0, 2, 0]);
/// ---
@nogc @property void opDispatch(string op, U)(U x) pure
if ((op.length >= 2)
&& (isValidSwizzleUnique!op) // v.xyy will be rejected
&& is(typeof(Vector!(T, op.length)(x)))) // can be converted to a small vector of the right size
{
Vector!(T, op.length) conv = x;
enum indexTuple = swizzleTuple!op;
foreach(i, index; indexTuple)
v[index] = conv[i];
}
/// Casting to small vectors of the same size.
/// Example:
/// ---
/// vec4f vf;
/// vec4d vd = cast!(vec4d)vf;
/// ---
@nogc U opCast(U)() pure const nothrow if (isVector!U && (U._N == _N))
{
U res = void;
mixin(generateLoopCode!("res.v[@] = cast(U._T)v[@];", N)());
return res;
}
/// Implement slices operator overloading.
/// Allows to go back to slice world.
/// Returns: length.
@nogc int opDollar() pure const nothrow
{
return N;
}
/// Slice containing vector values
/// Returns: a slice which covers the whole Vector.
@nogc T[] opSlice() pure nothrow
{
return v[];
}
/// vec[a..b]
@nogc T[] opSlice(int a, int b) pure nothrow
{
return v[a..b];
}
/// Squared Euclidean length of the Vector
/// Returns: squared length.
@nogc T squaredMagnitude() pure const nothrow
{
T sumSquares = 0;
mixin(generateLoopCode!("sumSquares += v[@] * v[@];", N)());
return sumSquares;
}
/// Squared Euclidean distance between this vector and another one
/// Returns: squared Euclidean distance.
@nogc T squaredDistanceTo(Vector v) pure const nothrow
{
return (v - this).squaredMagnitude();
}
static if (isFloatingPoint!T)
{
/// Euclidean length of the vector
/// Returns: Euclidean length
@nogc T magnitude() pure const nothrow
{
return sqrt(squaredMagnitude());
}
/// Inverse Euclidean length of the vector
/// Returns: Inverse of Euclidean length.
@nogc T inverseMagnitude() pure const nothrow
{
return 1 / sqrt(squaredMagnitude());
}
alias fastInverseLength = fastInverseMagnitude;
/// Faster but less accurate inverse of Euclidean length.
/// Returns: Inverse of Euclidean length.
@nogc T fastInverseMagnitude() pure const nothrow
{
return inverseSqrt(squaredMagnitude());
}
/// Euclidean distance between this vector and another one
/// Returns: Euclidean distance between this and other.
@nogc T distanceTo(Vector other) pure const nothrow
{
return (other - this).magnitude();
}
/// In-place normalization.
@nogc void normalize() pure nothrow
{
auto invMag = inverseMagnitude();
mixin(generateLoopCode!("v[@] *= invMag;", N)());
}
/// Returns a normalized copy of this Vector
/// Returns: Normalized vector.
@nogc Vector normalized() pure const nothrow
{
Vector res = this;
res.normalize();
return res;
}
/// Faster but less accurate in-place normalization.
@nogc void fastNormalize() pure nothrow
{
auto invLength = fastInverseMagnitude();
mixin(generateLoopCode!("v[@] *= invLength;", N)());
}
/// Faster but less accurate vector normalization.
/// Returns: Normalized vector.
@nogc Vector fastNormalized() pure const nothrow
{
Vector res = this;
res.fastNormalize();
return res;
}
static if (N == 3)
{
/// Gets an orthogonal vector from a 3-dimensional vector.
/// Doesn’t normalize the output.
/// Authors: Sam Hocevar
/// See_also: Source at $(WEB lolengine.net/blog/2013/09/21/picking-orthogonal-vector-combing-coconuts).
@nogc Vector getOrthogonalVector() pure const nothrow
{
return abs(x) > abs(z) ? Vector(-y, x, 0.0) : Vector(0.0, -z, y);
}
}
}
}
private
{
enum _N = N;
alias T _T;
// define types that can be converted to this, but are not the same type
template isConvertible(T)
{
enum bool isConvertible = (!is(T : Vector))
&& is(typeof(
{
T x;
Vector v = x;
}()));
}
// define types that can't be converted to this
template isForeign(T)
{
enum bool isForeign = (!isConvertible!T) && (!is(T: Vector));
}
template isValidSwizzle(string op, int lastSwizzleClass = -1)
{
static if (op.length == 0)
enum bool isValidSwizzle = true;
else
{
enum len = op.length;
enum int swizzleClass = swizzleClassify!(op[0]);
enum bool swizzleClassValid = (lastSwizzleClass == -1 || (swizzleClass == lastSwizzleClass));
enum bool isValidSwizzle = (swizzleIndex!(op[0]) != -1)
&& swizzleClassValid
&& isValidSwizzle!(op[1..len], swizzleClass);
}
}
template searchElement(char c, string s)
{
static if (s.length == 0)
{
enum bool result = false;
}
else
{
enum string tail = s[1..s.length];
enum bool result = (s[0] == c) || searchElement!(c, tail).result;
}
}
template hasNoDuplicates(string s)
{
static if (s.length == 1)
{
enum bool result = true;
}
else
{
enum tail = s[1..s.length];
enum bool result = !(searchElement!(s[0], tail).result) && hasNoDuplicates!(tail).result;
}
}
// true if the swizzle has at the maximum one time each letter
template isValidSwizzleUnique(string op)
{
static if (isValidSwizzle!op)
enum isValidSwizzleUnique = hasNoDuplicates!op.result;
else
enum bool isValidSwizzleUnique = false;
}
template swizzleIndex(char c)
{
static if((c == 'x' || c == 'r') && N >= 1)
enum swizzleIndex = 0;
else static if((c == 'y' || c == 'g') && N >= 2)
enum swizzleIndex = 1;
else static if((c == 'z' || c == 'b') && N >= 3)
enum swizzleIndex = 2;
else static if ((c == 'w' || c == 'a') && N >= 4)
enum swizzleIndex = 3;
else
enum swizzleIndex = -1;
}
template swizzleClassify(char c)
{
static if(c == 'x' || c == 'y' || c == 'z' || c == 'w')
enum swizzleClassify = 0;
else static if(c == 'r' || c == 'g' || c == 'b' || c == 'a')
enum swizzleClassify = 1;
else
enum swizzleClassify = -1;
}
template swizzleTuple(string op)
{
enum opLength = op.length;
static if (op.length == 0)
enum swizzleTuple = [];
else
enum swizzleTuple = [ swizzleIndex!(op[0]) ] ~ swizzleTuple!(op[1..op.length]);
}
}
}
/// True if `T` is some kind of `Vector`
enum isVector(T) = is(T : Vector!U, U...);
///
unittest
{
static assert(isVector!vec2f);
static assert(isVector!vec3d);
static assert(isVector!(vec4!real));
static assert(!isVector!float);
}
/// Get the numeric type used to measure a vectors's coordinates.
alias DimensionType(T : Vector!U, U...) = U[0];
///
unittest
{
static assert(is(DimensionType!vec2f == float));
static assert(is(DimensionType!vec3d == double));
}
///
template vec2(T) { alias Vector!(T, 2) vec2; }
///
template vec3(T) { alias Vector!(T, 3) vec3; }
///
template vec4(T) { alias Vector!(T, 4) vec4; }
alias vec2!int vec2i; ///
alias vec2!float vec2f; ///
alias vec2!double vec2d; ///
alias vec3!int vec3i; ///
alias vec3!float vec3f; ///
alias vec3!double vec3d; ///
alias vec4!int vec4i; ///
alias vec4!float vec4f; ///
alias vec4!double vec4d; ///
private
{
static string generateLoopCode(string formatString, int N)() pure nothrow
{
string result;
for (int i = 0; i < N; ++i)
{
string index = ctIntToString(i);
// replace all @ by indices
result ~= formatString.replace("@", index);
}
return result;
}
// Speed-up CTFE conversions
static string ctIntToString(int n) pure nothrow
{
static immutable string[16] table = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"];
if (n < 10)
return table[n];
else
return to!string(n);
}
}
/// Element-wise minimum.
@nogc Vector!(T, N) minByElem(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
{
import std.algorithm: min;
Vector!(T, N) res = void;
mixin(generateLoopCode!("res.v[@] = min(a.v[@], b.v[@]);", N)());
return res;
}
/// Element-wise maximum.
@nogc Vector!(T, N) maxByElem(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
{
import std.algorithm: max;
Vector!(T, N) res = void;
mixin(generateLoopCode!("res.v[@] = max(a.v[@], b.v[@]);", N)());
return res;
}
/// Element-wise absolute value.
@nogc Vector!(T, N) absByElem(T, int N)(const Vector!(T, N) a) pure nothrow
{
Vector!(T, N) res = void;
mixin(generateLoopCode!("res.v[@] = abs(a.v[@]);", N)());
return res;
}
/// Dot product of two vectors
/// Returns: Dot product.
@nogc T dot(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
{
T sum = 0;
mixin(generateLoopCode!("sum += a.v[@] * b.v[@];", N)());
return sum;
}
/// Cross product of two 3D vectors
/// Returns: 3D cross product.
/// Thanks to vuaru for corrections.
@nogc Vector!(T, 3) cross(T)(const Vector!(T, 3) a, const Vector!(T, 3) b) pure nothrow
{
return Vector!(T, 3)(a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x);
}
/// 3D reflect, like the GLSL function.
/// Returns: a reflected by normal b.
@nogc Vector!(T, N) reflect(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
{
return a - (2 * dot(b, a)) * b;
}
///
@nogc unittest
{
// reflect a 2D vector across the x axis (the normal points along the y axis)
assert(vec2f(1,1).reflect(vec2f(0,1)) == vec2f(1,-1));
assert(vec2f(1,1).reflect(vec2f(0,-1)) == vec2f(1,-1));
// note that the normal must be, well, normalized:
assert(vec2f(1,1).reflect(vec2f(0,20)) != vec2f(1,-1));
// think of this like a ball hitting a flat floor at an angle.
// the x and y components remain unchanged, and the z inverts
assert(vec3f(2,3,-0.5).reflect(vec3f(0,0,1)) == vec3f(2,3,0.5));
}
/// Angle between two vectors
/// Returns: angle between vectors.
/// See_also: "The Right Way to Calculate Stuff" at $(WEB www.plunk.org/~hatch/rightway.php)
@nogc T angleBetween(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
{
auto aN = a.normalized();
auto bN = b.normalized();
auto dp = dot(aN, bN);
if (dp < 0)
return T(PI) - 2 * asin((-bN-aN).magnitude / 2);
else
return 2 * asin((bN-aN).magnitude / 2);
}
static assert(vec2f.sizeof == 8);
static assert(vec3d.sizeof == 24);
static assert(vec4i.sizeof == 16);
unittest
{
static assert(vec2i.isValidSwizzle!"xyx");
static assert(!vec2i.isValidSwizzle!"xyz");
static assert(vec4i.isValidSwizzle!"brra");
static assert(!vec4i.isValidSwizzle!"rgyz");
static assert(vec2i.isValidSwizzleUnique!"xy");
static assert(vec2i.isValidSwizzleUnique!"yx");
static assert(!vec2i.isValidSwizzleUnique!"xx");
alias vec2l = vec2!long;
alias vec3ui = vec3!uint;
alias vec4ub = vec4!ubyte;
assert(vec2l(0, 1) == vec2i(0, 1));
int[2] arr = [0, 1];
int[] arr2 = new int[2];
arr2[] = arr[];
vec2i a = vec2i([0, 1]);
vec2i a2 = vec2i(0, 1);
immutable vec2i b = vec2i(0);
assert(b[0] == 0 && b[1] == 0);
vec2i c = arr;
vec2l d = arr2;
assert(a == a2);
assert(a == c);
assert(vec2l(a) == vec2l(a));
assert(vec2l(a) == d);
int[vec2i] hashMap;
hashMap[a] = (c - a).squaredMagnitude;
assert(hashMap[a] == (c - a).squaredMagnitude);
vec4i x = [4, 5, 6, 7];
assert(x == x);
--x[0];
assert(x[0] == 3);
++x[0];
assert(x[0] == 4);
x[1] &= 1;
x[2] = 77 + x[2];
x[3] += 3;
assert(x == [4, 1, 83, 10]);
assert(x.xxywz == [4, 4, 1, 10, 83]);
assert(x.xxxxxxx == [4, 4, 4, 4, 4, 4, 4]);
assert(x.abgr == [10, 83, 1, 4]);
assert(a != b);
x = vec4i(x.xyz, 166);
assert(x == [4, 1, 83, 166]);
vec2l e = a;
vec2l f = a + b;
assert(f == vec2l(a));
vec3ui g = vec3i(78,9,4);
g ^= vec3i(78,9,4);
assert(g == vec3ui(0));
//g[0..2] = 1u;
//assert(g == [2, 1, 0]);
assert(vec2i(4, 5) + 1 == vec2i(5,6));
assert(vec2i(4, 5) - 1 == vec2i(3,4));
assert(1 + vec2i(4, 5) == vec2i(5,6));
assert(vec3f(1,1,1) * 0 == 0);
assert(1.0 * vec3d(4,5,6) == vec3f(4,5.0f,6.0));
auto dx = vec2i(1,2);
auto dy = vec2i(4,5);
auto dp = dot(dx, dy);
assert(dp == 14 );
vec3i h = cast(vec3i)(vec3d(0.5, 1.1, -2.2));
assert(h == [0, 1, -2]);
assert(h[] == [0, 1, -2]);
assert(h[1..3] == [1, -2]);
assert(h.zyx == [-2, 1, 0]);
h.yx = vec2i(5, 2); // swizzle assignment
assert(h.xy == [2, 5]);
assert(-h[1] == -5);
assert(++h[0] == 3);
//assert(h == [-2, 1, 0]);
assert(!__traits(compiles, h.xx = h.yy));
vec4ub j;
assert(lerp(vec2f(-10, -1), vec2f(10, 1), 0.5) == vec2f(0, 0));
// larger vectors
alias Vector!(float, 5) vec5f;
vec5f l = vec5f(1, 2.0f, 3.0, 4u, 5.0L);
l = vec5f(l.xyz, vec2i(1, 2));
// the ctor should not compile if given too many arguments
static assert(!is(typeof(vec2f(1, 2, 3))));
static assert(!is(typeof(vec2f(vec2f(1, 2), 3))));
static assert( is(typeof(vec3f(vec2f(1, 2), 3))));
static assert( is(typeof(vec3f(1, 2, 3))));
assert(absByElem(vec3i(-1, 0, 2)) == vec3i(1, 0, 2));
}

2
src/shared/math.d
View File

@ -542,6 +542,8 @@ struct Matrix(T, int D)
auto l = &this;
auto r = &x;
// TODO: Test assembly on windows, may need to push registers onto the stack, decrement by 192 and move registers 1 by 1
asm pure @trusted @nogc nothrow
{
mov R8, l;