shishantbiswas/c3c
1
1
Fork 0
mirror of https://github.com/c3lang/c3c.git synced 2026-02-27 12:01:16 +00:00
Code Issues Packages Projects Releases Wiki Activity

Updated complex / matrix. Added quaternion math, vectors. Possible to add and mult scalar with vector. Fix where negating a float vector would be lowered incorrectly. Fix where $typeof(x) { ... } would not be valid compound literal. Fix where var would not be recognized as starting a declaration (e.g. in if (var x = ...)

Browse Source
This commit is contained in:
Christoffer Lerno
2023-01-20 21:27:36 +01:00
committed by Christoffer Lerno
parent e09628b664
commit 5151586450
14 changed files with 561 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files

83
lib/std/math/math.c3
View File

@@ -3,6 +3,8 @@
// a copy of which can be found in the LICENSE_STDLIB file.
module std::math;
import std::math::complex;
import std::math::matrix;
import std::math::quaternion;
// TODO Define these using quad precision.
const E = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466;
@@ -78,8 +80,31 @@ enum RoundingMode : int
TOWARD_NEG_INFINITY
}
define Complex32 = Complex<float>;
define Complex64 = Complex<double>;
fault MatrixError
{
MATRIX_INVERSE_DOESNT_EXIST,
}
define Complexf = Complex<float>;
define Complex = Complex<double>;
define complexf_identity = complex::identity<float>;
define complex_identity = complex::identity<double>;
define Quaternionf = Quaternion<float>;
define Quaternion = Quaternion<double>;
define quaternionf_identity = quaternion::identity<float>;
define quaternion_identity = quaternion::identity<double>;
define Matrix2f = Matrix2x2<float>;
define Matrix2 = Matrix2x2<double>;
define Matrix3f = Matrix3x3<float>;
define Matrix3 = Matrix3x3<double>;
define Matrix4f = Matrix4x4<float>;
define Matrix4 = Matrix4x4<double>;
define matrix4_ortho = matrix::ortho<double>;
define matrix4_perspective = matrix::perspective<double>;
define matrix4f_ortho = matrix::ortho<float>;
define matrix4f_perspective = matrix::perspective<float>;
/**
* @require types::is_numerical($typeof(x)) `The input must be a numerical value or numerical vector`
@@ -333,6 +358,19 @@ macro tanh(x) = (exp(2.0 * x) - 1.0) / (exp(2.0 * x) + 1.0);
**/
macro trunc(x) = $$trunc(x);
private macro lerp(x, y, amount) = x + (y - x) * amount;
private macro reflect(x, y)
{
var dot = x.dot(y);
return x - 2 * y * dot;
}
private macro normalize(x)
{
var len = x.length();
if (len == 0) return x;
return x * (1 / len);
}
macro float float.ceil(float x) = $$ceil(x);
macro float float.clamp(float x, float lower, float upper) = $$max(lower, $$min(x, upper));
macro float float.copysign(float mag, float sgn) = $$copysign(mag, sgn);
@@ -361,6 +399,20 @@ macro float[<*>] float[<*>].rint(float[<*>] x) = $$rint(x);
macro float[<*>] float[<*>].round(float[<*>] x) = $$round(x);
macro float[<*>] float[<*>].roundeven(float[<*>] x) = $$roundeven(x);
macro float[<*>] float[<*>].trunc(float[<*>] x) = $$trunc(x);
macro float float[<*>].dot(float[<*>] x, float[<*>] y) = (x * y).sum();
macro float float[<*>].length(float[<*>] x) = $$sqrt(x.dot(x));
macro float float[<*>].distance(float[<*>] x, float[<*>] y) = (x - y).length();
macro float[<*>] float[<*>].normalize(float[<*>] x) = normalize(x);
macro float[<*>] float[<*>].lerp(float[<*>] x, float[<*>] y, float amount) = lerp(x, y, amount);
macro float[<*>] float[<*>].reflect(float[<*>] x, float[<*>] y) = reflect(x, y);
macro bool float[<*>].equals(float[<*>] x, float[<*>] y) = equals_vec(x, y);
macro bool[<*>] float[<*>].comp_lt(float[<*>] x, float[<*>] y) = $$veccomplt(x, y);
macro bool[<*>] float[<*>].comp_le(float[<*>] x, float[<*>] y) = $$veccomple(x, y);
macro bool[<*>] float[<*>].comp_eq(float[<*>] x, float[<*>] y) = $$veccompeq(x, y);
macro bool[<*>] float[<*>].comp_gt(float[<*>] x, float[<*>] y) = $$veccompgt(x, y);
macro bool[<*>] float[<*>].comp_ge(float[<*>] x, float[<*>] y) = $$veccompge(x, y);
macro bool[<*>] float[<*>].comp_ne(float[<*>] x, float[<*>] y) = $$veccompne(x, y);
macro double double.ceil(double x) = $$ceil(x);
macro double double.clamp(double x, double lower, double upper) = $$max(lower, $$min(x, upper));
@@ -390,6 +442,20 @@ macro double[<*>] double[<*>].rint(double[<*>] x) = $$rint(x);
macro double[<*>] double[<*>].round(double[<*>] x) = $$round(x);
macro double[<*>] double[<*>].roundeven(double[<*>] x) = $$roundeven(x);
macro double[<*>] double[<*>].trunc(double[<*>] x) = $$trunc(x);
macro double double[<*>].dot(double[<*>] x, double[<*>] y) = (x * y).sum();
macro double double[<*>].length(double[<*>] x) = $$sqrt(x.dot(x));
macro double double[<*>].distance(double[<*>] x, double[<*>] y) = (x - y).length();
macro double[<*>] double[<*>].normalize(double[<*>] x) = normalize(x);
macro double[<*>] double[<*>].reflect(double[<*>] x, double[<*>] y) = reflect(x, y);
macro double[<*>] double[<*>].lerp(double[<*>] x, double[<*>] y, double amount) = lerp(x, y, amount);
macro bool double[<*>].equals(double[<*>] x, double[<*>] y) = equals_vec(x, y);
macro bool[<*>] double[<*>].comp_lt(double[<*>] x, double[<*>] y) = $$veccomplt(x, y);
macro bool[<*>] double[<*>].comp_le(double[<*>] x, double[<*>] y) = $$veccomple(x, y);
macro bool[<*>] double[<*>].comp_eq(double[<*>] x, double[<*>] y) = $$veccompeq(x, y);
macro bool[<*>] double[<*>].comp_gt(double[<*>] x, double[<*>] y) = $$veccompgt(x, y);
macro bool[<*>] double[<*>].comp_ge(double[<*>] x, double[<*>] y) = $$veccompge(x, y);
macro bool[<*>] double[<*>].comp_ne(double[<*>] x, double[<*>] y) = $$veccompne(x, y);
macro ichar ichar[<*>].sum(ichar[<*>] x) = $$reduce_add(x);
macro ichar ichar[<*>].product(ichar[<*>] x) = $$reduce_mul(x);
@@ -508,3 +574,16 @@ macro next_power_of_2(x)
while (y < x) y += y;
return y;
}
import std::io;
private macro equals_vec(v1, v2)
{
var $elements = v1.len;
var abs_diff = math::abs(v1 - v2);
io::printfn("diff %s", abs_diff);
var abs_v1 = math::abs(v1);
var abs_v2 = math::abs(v2);
io::printfn("abs %s", abs_v1);
$typeof(abs_v2) eps = 1;
return abs_diff.comp_le(FLOAT_EPSILON * math::max(abs_v1, abs_v2, eps)).and();
}

32
lib/std/math/math.complex.c3
View File

@@ -1,31 +1,27 @@
module std::math::complex<Type>;
module std::math::complex<Real>;
union Complex
{
struct
{
Type r, c;
Real r, c;
}
Type[<2>] v;
Real[<2>] v;
}
macro Complex Complex.add(Complex a, Complex b)
{
return Complex { .v = a.v + b.v };
}
macro Complex Complex.sub(Complex a, Complex b)
macro Complex identity() = { 1, 0 };
macro Complex Complex.add(Complex a, Complex b) = Complex { .v = a.v + b.v };
macro Complex Complex.add_each(Complex a, Real b) = Complex { .v = a.v + b };
macro Complex Complex.sub(Complex a, Complex b) = Complex { .v = a.v - b.v };
macro Complex Complex.sub_each(Complex a, Real b) = Complex { .v = a.v - b };
macro Complex Complex.scale(Complex a, Real s) = Complex { .v = a.v * s };
macro Complex Complex.mul(Complex a, Complex b)
{
return Complex { .v = a.v - b.v };
return { .r = a.r * b.r - a.c * b.c, .c = a.r * b.c + b.r * a.c };
}
macro Complex Complex.mult(Complex a, Complex b)
macro Complex Complex.div(Complex a, Complex b)
{
return Complex { .r = a.r * b.r - a.c * b.c, .c = a.r * b.c + b.r * a.c };
}
fn Complex Complex.div(Complex a, Complex b)
{
Type div = b.r * b.r + b.c * b.c;
Real div = b.r * b.r + b.c * b.c;
return Complex { .r = (a.r * b.r + a.c * b.c) / div, .c = (a.c * b.r - a.r * b.c) / div };
}
}

160
lib/std/math/math.matrix.c3
View File

@@ -1,9 +1,4 @@
module std::math::matrix;
fault MatrixError
{
MATRIX_INVERSE_DOESNT_EXIST,
}
module std::math::matrix<Real>;
struct Matrix2x2
{
@@ -11,10 +6,10 @@ struct Matrix2x2
{
struct
{
float m00, m01;
float m10, m11;
Real m00, m01;
Real m10, m11;
}
float[4] m;
Real[4] m;
}
}
@@ -24,11 +19,11 @@ struct Matrix3x3
{
struct
{
float m00; float m01; float m02;
float m10; float m11; float m12;
float m20; float m21; float m22;
Real m00, m01, m02;
Real m10, m11, m12;
Real m20, m21, m22;
}
float[9] m;
Real[9] m;
}
}
@@ -38,35 +33,35 @@ struct Matrix4x4
{
struct
{
float m00, m01, m02, m03;
float m10, m11, m12, m13;
float m20, m21, m22, m23;
float m30, m31, m32, m33;
Real m00, m01, m02, m03;
Real m10, m11, m12, m13;
Real m20, m21, m22, m23;
Real m30, m31, m32, m33;
}
float[16] m;
Real[16] m;
}
}
fn float[<2>] Matrix2x2.apply(Matrix2x2* mat, float[<2>] vec)
fn Real[<2>] Matrix2x2.apply(Matrix2x2* mat, Real[<2>] vec)
{
return float[<2>] {
return Real[<2>] {
mat.m00 * vec[0] + mat.m01 * vec[1],
mat.m10 * vec[0] + mat.m11 * vec[1],
};
}
fn float[<3>] Matrix3x3.apply(Matrix3x3* mat, float[<3>] vec)
fn Real[<3>] Matrix3x3.apply(Matrix3x3* mat, Real[<3>] vec)
{
return float[<3>] {
return Real[<3>] {
mat.m00 * vec[0] + mat.m01 * vec[1] + mat.m02 * vec[2],
mat.m10 * vec[0] + mat.m11 * vec[1] + mat.m12 * vec[2],
mat.m20 * vec[0] + mat.m21 * vec[1] + mat.m22 * vec[2],
};
}
fn float[<4>] Matrix4x4.apply(Matrix4x4* mat, float[<4>] vec)
fn Real[<4>] Matrix4x4.apply(Matrix4x4* mat, Real[<4>] vec)
{
return float[<4>] {
return Real[<4>] {
mat.m00 * vec[0] + mat.m01 * vec[1] + mat.m02 * vec[2] + mat.m03 * vec[3],
mat.m10 * vec[0] + mat.m11 * vec[1] + mat.m12 * vec[2] + mat.m13 * vec[3],
mat.m20 * vec[0] + mat.m21 * vec[1] + mat.m22 * vec[2] + mat.m23 * vec[3],
@@ -125,33 +120,18 @@ fn Matrix4x4 Matrix4x4.mul(Matrix4x4* a, Matrix4x4 b)
};
}
fn Matrix2x2 Matrix2x2.component_mul(Matrix2x2* mat, Real s) = matrix_component_mul(mat, s);
fn Matrix3x3 Matrix3x3.component_mul(Matrix3x3* mat, Real s) = matrix_component_mul(mat, s);
fn Matrix4x4 Matrix4x4.component_mul(Matrix4x4* mat, Real s) = matrix_component_mul(mat, s);
fn Matrix2x2 Matrix2x2.component_mul(Matrix2x2* mat, float s)
{
return Matrix2x2 {
mat.m00 * s, mat.m01 * s,
mat.m10 * s, mat.m11 * s,
};
}
fn Matrix2x2 Matrix2x2.add(Matrix2x2* mat, Matrix2x2 mat2) = matrix_add(mat, mat2);
fn Matrix3x3 Matrix3x3.add(Matrix3x3* mat, Matrix3x3 mat2) = matrix_add(mat, mat2);
fn Matrix4x4 Matrix4x4.add(Matrix4x4* mat, Matrix4x4 mat2) = matrix_add(mat, mat2);
fn Matrix3x3 Matrix3x3.component_mul(Matrix3x3* mat, float s)
{
return Matrix3x3 {
mat.m00 * s, mat.m01 * s, mat.m02 * s,
mat.m10 * s, mat.m11 * s, mat.m12 * s,
mat.m20 * s, mat.m21 * s, mat.m22 * s,
};
}
fn Matrix2x2 Matrix2x2.sub(Matrix2x2* mat, Matrix2x2 mat2) = matrix_sub(mat, mat2);
fn Matrix3x3 Matrix3x3.sub(Matrix3x3* mat, Matrix3x3 mat2) = matrix_sub(mat, mat2);
fn Matrix4x4 Matrix4x4.sub(Matrix4x4* mat, Matrix4x4 mat2) = matrix_sub(mat, mat2);
fn Matrix4x4 Matrix4x4.component_mul(Matrix4x4* mat, float s)
{
return Matrix4x4 {
mat.m00 * s, mat.m01 * s, mat.m02 * s, mat.m03 * s,
mat.m10 * s, mat.m11 * s, mat.m12 * s, mat.m13 * s,
mat.m20 * s, mat.m21 * s, mat.m22 * s, mat.m23 * s,
mat.m30 * s, mat.m31 * s, mat.m32 * s, mat.m33 * s,
};
}
fn Matrix2x2 Matrix2x2.transpose(Matrix2x2* mat)
@@ -182,12 +162,12 @@ fn Matrix4x4 Matrix4x4.transpose(Matrix4x4* mat)
}
fn float Matrix2x2.determinant(Matrix2x2* mat)
fn Real Matrix2x2.determinant(Matrix2x2* mat)
{
return mat.m00 * mat.m11 - mat.m01 * mat.m10;
}
fn float Matrix3x3.determinant(Matrix3x3* mat)
fn Real Matrix3x3.determinant(Matrix3x3* mat)
{
return
mat.m00 * (mat.m11 * mat.m22 - mat.m21 * mat.m12) -
@@ -195,7 +175,7 @@ fn float Matrix3x3.determinant(Matrix3x3* mat)
mat.m02 * (mat.m10 * mat.m21 - mat.m20 * mat.m11);
}
fn float Matrix4x4.determinant(Matrix4x4* mat)
fn Real Matrix4x4.determinant(Matrix4x4* mat)
{
return
mat.m00 * (mat.m11 * (mat.m22 * mat.m33 - mat.m32 * mat.m23) -
@@ -295,7 +275,7 @@ fn Matrix4x4 Matrix4x4.adjoint(Matrix4x4* mat)
fn Matrix2x2! Matrix2x2.inverse(Matrix2x2* m)
{
float det = m.determinant();
Real det = m.determinant();
if (det == 0) return MatrixError.MATRIX_INVERSE_DOESNT_EXIST!;
Matrix2x2 adj = m.adjoint();
return adj.component_mul(1 / det).transpose();
@@ -303,7 +283,7 @@ fn Matrix2x2! Matrix2x2.inverse(Matrix2x2* m)
fn Matrix3x3! Matrix3x3.inverse(Matrix3x3* m)
{
float det = m.determinant();
Real det = m.determinant();
if (det == 0) return MatrixError.MATRIX_INVERSE_DOESNT_EXIST!;
Matrix3x3 adj = m.adjoint();
return adj.component_mul(1 / det).transpose();
@@ -311,14 +291,14 @@ fn Matrix3x3! Matrix3x3.inverse(Matrix3x3* m)
fn Matrix4x4! Matrix4x4.inverse(Matrix4x4* m)
{
float det = m.determinant();
Real det = m.determinant();
if (det == 0) return MatrixError.MATRIX_INVERSE_DOESNT_EXIST!;
Matrix4x4 adj = m.adjoint();
return adj.component_mul(1 / det).transpose();
}
fn Matrix3x3 Matrix3x3.translate(Matrix3x3* m, float[<2>] v)
fn Matrix3x3 Matrix3x3.translate(Matrix3x3* m, Real[<2>] v)
{
return m.mul(Matrix3x3 {
1, 0, v[0],
@@ -327,7 +307,7 @@ fn Matrix3x3 Matrix3x3.translate(Matrix3x3* m, float[<2>] v)
});
}
fn Matrix4x4 Matrix4x4.translate(Matrix4x4* m, float[<3>] v)
fn Matrix4x4 Matrix4x4.translate(Matrix4x4* m, Real[<3>] v)
{
return m.mul(Matrix4x4 {
1, 0, 0, v[0],
@@ -338,7 +318,7 @@ fn Matrix4x4 Matrix4x4.translate(Matrix4x4* m, float[<3>] v)
}
// r in radians
fn Matrix3x3 Matrix3x3.rotate(Matrix3x3* m, float r)
fn Matrix3x3 Matrix3x3.rotate(Matrix3x3* m, Real r)
{
return m.mul(Matrix3x3 {
math::cos(r), -math::sin(r), 0,
@@ -348,7 +328,7 @@ fn Matrix3x3 Matrix3x3.rotate(Matrix3x3* m, float r)
}
// r in radians
fn Matrix4x4 Matrix4x4.rotate_z(Matrix4x4* m, float r)
fn Matrix4x4 Matrix4x4.rotate_z(Matrix4x4* m, Real r)
{
return m.mul(Matrix4x4 {
math::cos(r), -math::sin(r), 0, 0,
@@ -359,7 +339,7 @@ fn Matrix4x4 Matrix4x4.rotate_z(Matrix4x4* m, float r)
}
// r in radians
fn Matrix4x4 Matrix4x4.rotate_y(Matrix4x4* m, float r)
fn Matrix4x4 Matrix4x4.rotate_y(Matrix4x4* m, Real r)
{
return m.mul(Matrix4x4 {
math::cos(r), 0, -math::sin(r), 0,
@@ -370,7 +350,7 @@ fn Matrix4x4 Matrix4x4.rotate_y(Matrix4x4* m, float r)
}
// r in radians
fn Matrix4x4 Matrix4x4.rotate_x(Matrix4x4* m, float r)
fn Matrix4x4 Matrix4x4.rotate_x(Matrix4x4* m, Real r)
{
return m.mul(Matrix4x4 {
1, 0, 0, 0,
@@ -381,7 +361,7 @@ fn Matrix4x4 Matrix4x4.rotate_x(Matrix4x4* m, float r)
}
fn Matrix3x3 Matrix3x3.scale(Matrix3x3* m, float[<2>] v)
fn Matrix3x3 Matrix3x3.scale(Matrix3x3* m, Real[<2>] v)
{
return m.mul(Matrix3x3 {
v[0], 0, 0,
@@ -390,7 +370,11 @@ fn Matrix3x3 Matrix3x3.scale(Matrix3x3* m, float[<2>] v)
});
}
fn Matrix4x4 Matrix4x4.scale(Matrix4x4* m, float[<3>] v)
fn Real Matrix2x2.trace(Matrix2x2* m) = m.m00 + m.m11;
fn Real Matrix3x3.trace(Matrix3x3* m) = m.m00 + m.m11 + m.m22;
fn Real Matrix4x4.trace(Matrix4x4* m) = m.m00 + m.m11 + m.m22 + m.m33;
fn Matrix4x4 Matrix4x4.scale(Matrix4x4* m, Real[<3>] v)
{
return m.mul(Matrix4x4 {
v[0], 0, 0, 0,
@@ -400,32 +384,48 @@ fn Matrix4x4 Matrix4x4.scale(Matrix4x4* m, float[<3>] v)
});
}
fn Matrix4x4 ortho(float left, float right, float top, float bottom, float near, float far)
fn Matrix4x4 ortho(Real left, Real right, Real top, Real bottom, Real near, Real far)
{
float width = right - left;
float height = top - bottom;
float depth = far - near;
Real width = right - left;
Real height = top - bottom;
Real depth = far - near;
return Matrix4x4 {
2 / width, 0, 0, 0,
0, 2 / height, 0, 0,
0, 0, -2 / depth, 0,
-(right + left) / width, -(top + bottom) / height, -(far + near) / depth, 1,
2 / width, 0, 0, 0,
0, 2 / height, 0, 0,
0, 0, -2 / depth, 0,
-(right + left) / width, -(top + bottom) / height, -(far + near) / depth, 1
};
}
// fov in radians
fn Matrix4x4 perspective(float fov, float aspect_ratio, float near, float far)
fn Matrix4x4 perspective(Real fov, Real aspect_ratio, Real near, Real far)
{
float top = ((float)math::sin(fov / 2) / (float)math::cos(fov / 2)) * near;
float right = top * aspect_ratio;
float depth = far - near;
Real top = (math::sin(fov / 2) / math::cos(fov / 2)) * near;
Real right = top * aspect_ratio;
Real depth = far - near;
return Matrix4x4 {
1 / right, 0, 0, 0,
0, 1 / top, 0, 0,
0, 0, -2 / depth, 0,
0, 0, - (far + near) / depth, 1,
1 / right, 0, 0, 0,
0, 1 / top, 0, 0,
0, 0, -2 / depth, 0,
0, 0, -(far + near) / depth, 1,
};
}
private macro matrix_component_mul(mat, val)
{
var $Type = Real[<$typeof(mat.m).len>];
return $typeof(*mat) { .m = val * ($Type)mat.m };
}
private macro matrix_add(mat, mat2)
{
var $Type = Real[<$typeof(mat.m).len>];
return $typeof(*mat) { .m = ($Type)mat.m + ($Type)mat2.m };
}
private macro matrix_sub(mat, mat2)
{
var $Type = Real[<$typeof(mat.m).len>];
return $typeof(*mat) { .m = ($Type)mat.m - ($Type)mat2.m };
}

68
lib/std/math/math_quaternion.c3 Normal file
View File

@@ -0,0 +1,68 @@
module std::math::quaternion<Real>;
import std::math::vector;
union Quaternion
{
struct
{
Real i, j, k, l;
}
Real[<4>] v;
}
macro Quaternion identity() = { 0, 0, 0, 1 };
macro Quaternion Quaternion.add(Quaternion a, Quaternion b) = Quaternion { .v = a.v + b.v };
macro Quaternion Quaternion.add_each(Quaternion a, Real b) = Quaternion { .v = a.v + b };
macro Quaternion Quaternion.sub(Quaternion a, Quaternion b) = Quaternion { .v = a.v - b.v };
macro Quaternion Quaternion.sub_each(Quaternion a, Real b) = Quaternion { .v = a.v - b };
macro Quaternion Quaternion.scale(Quaternion a, Real s) = Quaternion { .v = a.v * s };
macro Quaternion Quaternion.normalize(Quaternion q) = { .v = q.v.normalize() };
macro Real Quaternion.length(Quaternion q) = q.v.length();
macro Quaternion Quaternion.lerp(Quaternion q1, Quaternion q2, Real amount) = { .v = q1.v.lerp(q2.v, amount) };
fn Quaternion Quaternion.nlerp(Quaternion q1, Quaternion q2, Real amount) = { .v = q1.v.lerp(q2.v, amount).normalize() };
fn Quaternion Quaternion.invert(Quaternion q)
{
Real length_sq = q.v.dot(q.v);
if (length_sq <= 0) return q;
Real inv_length = 1 / length_sq;
return { q.v[0] * -inv_length, q.v[1] * -inv_length, q.v[2] * -inv_length, q.v[3] * inv_length };
}
fn Quaternion Quaternion.slerp(Quaternion q1, Quaternion q2, Real amount)
{
Quaternion result = {};
Real[<4>] q2v = q2.v;
Real cos_half_theta = q1.v.dot(q2v);
if (cos_half_theta < 0)
{
q2v = -q2v;
cos_half_theta = -cos_half_theta;
}
if (cos_half_theta >= 1) return q1;
Real[<4>] q1v = q1.v;
if (cos_half_theta > 0.95f) return { .v = q1v.lerp(q2v, amount) };
Real half_theta = math::cos(cos_half_theta);
Real sin_half_theta = math::sqrt(1 - cos_half_theta * cos_half_theta);
if (math::abs(sin_half_theta) < 0.001f)
{
return { .v = (q1v + q2v) * 0.5f };
}
Real ratio_a = math::sin((1 - amount) * half_theta) / sin_half_theta;
Real ratio_b = math::sin(amount * half_theta) / sin_half_theta;
return { .v = q1v * ratio_a + q2v * ratio_b };
}
fn Quaternion Quaternion.mul(Quaternion a, Quaternion b)
{
return { a.i * b.l + a.l * b.i + a.j * b.k - a.k * b.j,
a.j * b.l + a.l * b.j + a.k * b.i - a.i * b.k,
a.k * b.l + a.l * b.k + a.i * b.j - a.j * b.i,
a.l * b.l - a.i * b.i - a.j * a.j - a.k * a.k };
}

258
lib/std/math/math_vector.c3 Normal file
View File

@@ -0,0 +1,258 @@
module std::math::vector;
import std::math;
define Vec2f = float[<2>];
define Vec3f = float[<3>];
define Vec4f = float[<4>];
define Vec2 = double[<2>];
define Vec3 = double[<3>];
define Vec4 = double[<4>];
macro Vec2f.length_sq(Vec2f v) = v.dot(v);
macro Vec3f.length_sq(Vec3f v) = v.dot(v);
macro Vec4f.length_sq(Vec4f v) = v.dot(v);
macro Vec2.length_sq(Vec2 v) = v.dot(v);
macro Vec3.length_sq(Vec3 v) = v.dot(v);
macro Vec4.length_sq(Vec4 v) = v.dot(v);
macro Vec2f.distance_sq(Vec2f v1, Vec2f v2) = (v1 - v2).length_sq();
macro Vec3f.distance_sq(Vec3f v1, Vec3f v2) = (v1 - v2).length_sq();
macro Vec4f.distance_sq(Vec4f v1, Vec4f v2) = (v1 - v2).length_sq();
macro Vec2.distance_sq(Vec2 v1, Vec2 v2) = (v1 - v2).length_sq();
macro Vec3.distance_sq(Vec3 v1, Vec3 v2) = (v1 - v2).length_sq();
macro Vec4.distance_sq(Vec4 v1, Vec4 v2) = (v1 - v2).length_sq();
macro Vec2f.transform(Vec2f v, Matrix4f mat) = transform2(v, mat);
macro Vec2f.rotate(Vec2f v, float angle) = rotate(v, angle);
macro Vec2f.angle(Vec2f v1, Vec2f v2) = math::atan2(v2[1], v2[0]) - math::atan2(v1[1], v2[0]);
macro Vec2.transform(Vec2 v, Matrix4 mat) = transform2(v, mat);
macro Vec2.rotate(Vec2 v, double angle) = rotate(v, angle);
macro Vec2.angle(Vec2 v1, Vec2 v2) = math::atan2(v2[1], v2[0]) - math::atan2(v1[1], v2[0]);
macro Vec2f.clamp_mag(Vec2f v, float min, float max) = clamp_magnitude(v, min, max);
macro Vec3f.clamp_mag(Vec3f v, float min, float max) = clamp_magnitude(v, min, max);
macro Vec4f.clamp_mag(Vec4f v, float min, float max) = clamp_magnitude(v, min, max);
macro Vec2.clamp_mag(Vec2 v, double min, double max) = clamp_magnitude(v, min, max);
macro Vec3.clamp_mag(Vec3 v, double min, double max) = clamp_magnitude(v, min, max);
macro Vec4.clamp_mag(Vec4 v, double min, double max) = clamp_magnitude(v, min, max);
fn Vec2f Vec2f.towards(Vec2f v, Vec2f target, float max_distance) = towards(v, target, max_distance);
fn Vec3f Vec3f.towards(Vec3f v, Vec3f target, float max_distance) = towards(v, target, max_distance);
fn Vec4f Vec4f.towards(Vec4f v, Vec4f target, float max_distance) = towards(v, target, max_distance);
fn Vec2 Vec2.towards(Vec2 v, Vec2 target, double max_distance) = towards(v, target, max_distance);
fn Vec3 Vec3.towards(Vec3 v, Vec3 target, double max_distance) = towards(v, target, max_distance);
fn Vec4 Vec4.towards(Vec4 v, Vec4 target, double max_distance) = towards(v, target, max_distance);
fn Vec3f Vec3f.cross(Vec3f v1, Vec3f v2) = cross3(v1, v2);
fn Vec3 Vec3.cross(Vec3 v1, Vec3 v2) = cross3(v1, v2);
fn Vec3f Vec3f.perpendicular(Vec3f v) = perpendicular3(v);
fn Vec3 Vec3.perpendicular(Vec3 v) = perpendicular3(v);
fn Vec3f Vec3f.barycenter(Vec3f p, Vec3f a, Vec3f b, Vec3f c) = barycenter3(p, a, b, c);
fn Vec3 Vec3.barycenter(Vec3 p, Vec3 a, Vec3 b, Vec3 c) = barycenter3(p, a, b, c);
fn Vec3f Vec3f.transform(Vec3f v, Matrix4f mat) = transform3(v, mat);
fn Vec3 Vec3.transform(Vec3 v, Matrix4 mat) = transform3(v, mat);
fn float Vec3f.angle(Vec3f v1, Vec3f v2) = angle3(v1, v2);
fn double Vec3.angle(Vec3 v1, Vec3 v2) = angle3(v1, v2);
fn Vec3f Vec3f.refract(Vec3f v, Vec3f n, float r) = refract3(v, n, r);
fn Vec3 Vec3.refract(Vec3 v, Vec3 n, double r) = refract3(v, n, r);
fn void ortho_normalize(Vec3f* v1, Vec3f* v2) = ortho_normalize3(v1, v2);
fn void ortho_normalized(Vec3* v1, Vec3* v2) = ortho_normalize3(v1, v2);
fn Matrix4f matrix4f_look_at(Vec3f eye, Vec3f target, Vec3f up) = matrix_look_at(Matrix4f, eye, target, up);
fn Matrix4 matrix4_look_at(Vec3 eye, Vec3 target, Vec3 up) = matrix_look_at(Matrix4, eye, target, up);
fn Vec3f Vec3f.rotate_quat(Vec3f v, Quaternionf q) = rotate_by_quat3(v, q);
fn Vec3 Vec3.rotate_quat(Vec3 v, Quaternion q) = rotate_by_quat3(v, q);
fn Vec3f Vec3f.rotate_axis(Vec3f v, Vec3f axis, float angle) = rotate_axis_angle(v, axis, angle);
fn Vec3 Vec3.rotate_axis(Vec3 v, Vec3 axis, double angle) = rotate_axis_angle(v, axis, angle);
fn Vec3f Vec3f.unproject(Vec3f v, Matrix4f projection, Matrix4f view) = unproject3(v, projection, view);
fn Vec3 Vec3.unproject(Vec3 v, Matrix4 projection, Matrix4 view) = unproject3(v, projection, view);
private macro towards(v, target, max_distance)
{
var delta = target - v;
var square = delta.length_sq();
if (square == 0 || (max_distance >= 0 && (square <= max_distance * max_distance))) return target;
var dist = math::sqrt(square);
return v + delta * max_distance / dist;
}
private macro clamp_magnitude(v, min, max)
{
var length = v.dot(v);
if (length > 0)
{
length = math::sqrt(length);
if (length < min) return v * (min / length);
if (length > max) return v * (max / length);
}
return v;
}
private macro rotate(v, angle)
{
var c = math::cos(angle);
var s = math::sin(angle);
return $typeof(v) { v[0] * c - v[1] * s, v[0] * s + v[1] * c };
}
private macro perpendicular3(v)
{
var min = math::abs(v[0]);
$typeof(v) cardinal_axis = { 1, 0, 0 };
if (var vy = math::abs(v[1]), vy < min)
{
min = vy;
cardinal_axis = { 0, 1, 0 };
}
if (var vz = math::abs(v[2]), vz < min)
{
cardinal_axis = { 0, 0, 1 };
}
return cross3(v, cardinal_axis);
}
private macro cross3(v1, v2)
{
var a = $$shufflevector(v1, v1, { 1, 2, 0 }) * $$shufflevector(v2, v2, { 2, 0, 1 });
var b = $$shufflevector(v1, v1, { 2, 0, 1 }) * $$shufflevector(v2, v2, { 1, 2, 0 });
return a - b;
}
private macro transform2(v, mat)
{
return $typeof(v) { mat.m00 * v[0] + mat.m10 * v[1] + mat.30,
mat.m01 * v[0] + mar.m11 * v[1] + mat.31 };
}
private macro transform3(v, mat)
{
return $typeof(v) {
mat.m00 * v[0] + mat.m10 * v[1] + mat.m20 * v[2] + mat.m30,
mat.m01 * v[0] + mat.m11 * v[1] + mat.m21 * v[2] + mat.m31,
mat.m02 * v[0] + mat.m12 * v[1] + mat.m22 * v[2] + mat.m32
};
}
private macro angle3(v1, v2)
{
var len = v1.cross(v2).length();
var dot = v1.dot(v2);
return math::atan2(len, dot);
}
private macro void ortho_normalize3(v1, v2)
{
var v1n = *v1 = v1.normalize();
var vn1 = v1n.cross(*v2).normalize();
*v2 = v1n.cross(vn1);
}
private macro rotate_by_quat3(v, q)
{
return $typeof(v) {
v[0] * (q.i * q.i + q.l * q.l - q.j * q.j - q.k * q.k)
+ v[1] * (2 * q.i * q.j - 2 * q.l * q.k)
+ v[2] * (2 * q.i * q.k - 2 * q.l * q.j),
v[0] * (2 * q.l * q.k + 2 * q.i * q.j)
+ v[1] * (q.l * q.l - q.i * q.i + q.j * q.j - q.k * q.k)
+ v[2] * (-2 * q.l * q.i + 2 * q.j * q.k),
v[0] * (-2 * q.l * q.j + 2 * q.i * q.k)
+ v[1] * (2 * q.l * q.i + 2 * q.j * q.k)
+ v[2] * (q.l * q.l - q.i * q.i - q.j * q.j + q.k * q.k)
};
}
private macro rotate_axis_angle(v, axis, angle)
{
axis = axis.normalize();
angle /= 2;
var w = axis * math::sin(angle);
var wv = w.cross(v);
var wwv = w.cross(wv);
wv *= math::cos(angle) * 2;
wwv *= 2;
return v + wv + wwv;
}
private macro matrix_look_at($Type, eye, target, up)
{
var vz = (eye - target).normalize();
var vx = up.cross(vz).normalize();
var vy = vz.cross(vx);
return $Type {
vx[0], vy[0], vz[0], 0,
vx[1], vy[1], vz[1], 0,
vx[2], vy[2], vz[2], 0,
- vx.dot(eye), - vy.dot(eye), - vz.dot(eye), 1
};
}
private macro unproject3(v, m1, m2)
{
return v;
/*
var view_proj = m1.mul(m2);
var invert = view_proj.invert();
// Create quaternion from source point
$if ($typeof(v[0]).typeid == float.typeid):
Quaternionf quat = { v.x, v.y, v.z, 1 };
$else:
Quaternion quat = { v.x, v.y, v.z, 1 };
$endif;
// Multiply quat point by unproject matrix
var qtransformed = quat.transform(invert);
// Normalized world points in vectors
return {
qtransformed.i / qtransformed.l,
qtransformed.j / qtransformed.l,
qtransformed.k / qtransformed.l
};*/
}
private macro barycenter3(p, a, b, c)
{
var v0 = b - a;
var v1 = c - a;
var v2 = p - a;
var d00 = v0.dot(v0);
var d01 = v0.dot(v1);
var d11 = v1.dot(v1);
var d20 = v2.dot(v0);
var d21 = v2.dot(v1);
var denom = d00 * d11 - d01 * d01;
var y = (d11 * d20 - d01 * d21) / denom;
var z = (d00 * d21 - d01 * d20) / denom;
return $typeof(p) { 1 - y - z, y, z };
}
private macro refract3(v, n, r)
{
var dot = v.dot(n);
var d = 1 - r * r * (1 - dot * dot);
return d < 0 ? v : r * v - (r * dot + math::sqrt(d)) * n;
}