diff --git a/lib/std/math/math.c3 b/lib/std/math/math.c3 index 37f72b72f..a7624cd17 100644 --- a/lib/std/math/math.c3 +++ b/lib/std/math/math.c3 @@ -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; -define Complex64 = Complex; +fault MatrixError +{ + MATRIX_INVERSE_DOESNT_EXIST, +} + +define Complexf = Complex; +define Complex = Complex; +define complexf_identity = complex::identity; +define complex_identity = complex::identity; + +define Quaternionf = Quaternion; +define Quaternion = Quaternion; +define quaternionf_identity = quaternion::identity; +define quaternion_identity = quaternion::identity; + +define Matrix2f = Matrix2x2; +define Matrix2 = Matrix2x2; +define Matrix3f = Matrix3x3; +define Matrix3 = Matrix3x3; +define Matrix4f = Matrix4x4; +define Matrix4 = Matrix4x4; +define matrix4_ortho = matrix::ortho; +define matrix4_perspective = matrix::perspective; +define matrix4f_ortho = matrix::ortho; +define matrix4f_perspective = matrix::perspective; /** * @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(); +} \ No newline at end of file diff --git a/lib/std/math/math.complex.c3 b/lib/std/math/math.complex.c3 index f43c4276b..ea2eda2ad 100644 --- a/lib/std/math/math.complex.c3 +++ b/lib/std/math/math.complex.c3 @@ -1,31 +1,27 @@ -module std::math::complex; +module std::math::complex; 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 }; -} \ No newline at end of file +} diff --git a/lib/std/math/math.matrix.c3 b/lib/std/math/math.matrix.c3 index eaaf9a5ee..8d9b3f0e3 100644 --- a/lib/std/math/math.matrix.c3 +++ b/lib/std/math/math.matrix.c3 @@ -1,9 +1,4 @@ -module std::math::matrix; - -fault MatrixError -{ - MATRIX_INVERSE_DOESNT_EXIST, -} +module std::math::matrix; 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 }; +} diff --git a/lib/std/math/math_quaternion.c3 b/lib/std/math/math_quaternion.c3 new file mode 100644 index 000000000..009110359 --- /dev/null +++ b/lib/std/math/math_quaternion.c3 @@ -0,0 +1,68 @@ +module std::math::quaternion; +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 }; +} + + diff --git a/lib/std/math/math_vector.c3 b/lib/std/math/math_vector.c3 new file mode 100644 index 000000000..e4f7cb500 --- /dev/null +++ b/lib/std/math/math_vector.c3 @@ -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; +} \ No newline at end of file diff --git a/src/compiler/enums.h b/src/compiler/enums.h index 08ea4e3d2..e9c00f6e7 100644 --- a/src/compiler/enums.h +++ b/src/compiler/enums.h @@ -100,6 +100,7 @@ typedef enum CAST_FPBOOL, CAST_INTBOOL, CAST_INTENUM, + CAST_NUMVEC, CAST_FPFP, CAST_FPSI, CAST_FPUI, diff --git a/src/compiler/expr.c b/src/compiler/expr.c index 692ade16f..abef2a05b 100644 --- a/src/compiler/expr.c +++ b/src/compiler/expr.c @@ -370,6 +370,7 @@ static inline bool expr_cast_is_constant_eval(Expr *expr, ConstantEvalKind eval_ case CAST_VECARR: case CAST_ARRVEC: case CAST_BOOLVECINT: + case CAST_NUMVEC: if (eval_kind != CONSTANT_EVAL_NO_SIDE_EFFECTS) return false; return exprid_is_constant_eval(expr->cast_expr.expr, eval_kind); case CAST_XIPTR: diff --git a/src/compiler/llvm_codegen_expr.c b/src/compiler/llvm_codegen_expr.c index 9a0e19339..5a55e1535 100644 --- a/src/compiler/llvm_codegen_expr.c +++ b/src/compiler/llvm_codegen_expr.c @@ -367,7 +367,6 @@ LLVMValueRef llvm_coerce_int_ptr(GenContext *c, LLVMValueRef value, LLVMTypeRef LLVMValueRef llvm_emit_coerce(GenContext *c, LLVMTypeRef coerced, BEValue *value, Type *original_type) { - assert(type_flatten_distinct(original_type) == type_flatten_distinct(value->type)); LLVMTypeRef llvm_source_type = llvm_get_type(c, value->type); // 1. If the types match then we're done, just load. @@ -1164,6 +1163,19 @@ void llvm_emit_array_to_vector_cast(GenContext *c, BEValue *value, Type *to_type llvm_value_set(value, vector, to_type); } +void llvm_emit_num_to_vec_cast(GenContext *c, BEValue *value, Type *to_type, Type *from_type) +{ + llvm_value_rvalue(c, value); + LLVMTypeRef type = llvm_get_type(c, to_type); + unsigned elements = LLVMGetVectorSize(type); + LLVMValueRef res = LLVMGetUndef(type); + for (unsigned i = 0; i < elements; i++) + { + res = LLVMBuildInsertElement(c->builder, res, value->value, llvm_const_int(c, type_usize, i), ""); + } + llvm_value_set(value, res, to_type); +} + void llvm_emit_bool_to_intvec_cast(GenContext *c, BEValue *value, Type *to_type, Type *from_type) { llvm_value_rvalue(c, value); @@ -1181,6 +1193,9 @@ void llvm_emit_cast(GenContext *c, CastKind cast_kind, Expr *expr, BEValue *valu switch (cast_kind) { + case CAST_NUMVEC: + llvm_emit_num_to_vec_cast(c, value, to_type, from_type); + return; case CAST_BOOLVECINT: llvm_emit_bool_to_intvec_cast(c, value, to_type, from_type); return; @@ -2394,7 +2409,7 @@ static void llvm_emit_unary_expr(GenContext *c, BEValue *value, Expr *expr) case UNARYOP_NEG: llvm_emit_expr(c, value, inner); llvm_value_rvalue(c, value); - if (type_is_float(type)) + if (type_is_floatlike(type)) { value->value = LLVMBuildFNeg(c->builder, value->value, "fneg"); return; @@ -2402,7 +2417,7 @@ static void llvm_emit_unary_expr(GenContext *c, BEValue *value, Expr *expr) assert(type->canonical != type_bool); llvm_emit_expr(c, value, expr->unary_expr.expr); llvm_value_rvalue(c, value); - if (active_target.feature.trap_on_wrap) + if (active_target.feature.trap_on_wrap && !type_flat_is_vector(value->type)) { LLVMValueRef zero = llvm_get_zero(c, expr->unary_expr.expr->type); LLVMTypeRef type_to_use = llvm_get_type(c, type->canonical); diff --git a/src/compiler/parse_expr.c b/src/compiler/parse_expr.c index 17514667c..b1af9d0c4 100644 --- a/src/compiler/parse_expr.c +++ b/src/compiler/parse_expr.c @@ -477,6 +477,10 @@ static Expr *parse_type_expr(ParseContext *c, Expr *left) assert(!left && "Unexpected left hand side"); Expr *expr = EXPR_NEW_TOKEN(EXPR_TYPEINFO); ASSIGN_TYPE_OR_RET(TypeInfo *type, parse_optional_type(c), poisoned_expr); + if (tok_is(c, TOKEN_LBRACE)) + { + return parse_type_compound_literal_expr_after_type(c, type); + } expr->span = type->span; expr->type_expr = type; if (tok_is(c, TOKEN_SCOPE)) diff --git a/src/compiler/parse_global.c b/src/compiler/parse_global.c index ee533626c..a62e5260d 100644 --- a/src/compiler/parse_global.c +++ b/src/compiler/parse_global.c @@ -835,6 +835,11 @@ Decl *parse_decl(ParseContext *c) Expr *parse_decl_or_expr(ParseContext *c, Decl **decl_ref) { + if (tok_is(c, TOKEN_VAR)) + { + ASSIGN_DECL_OR_RET(*decl_ref, parse_var_decl(c), poisoned_expr); + return NULL; + } TypeInfo *type_info; Expr *expr = parse_expr(c); if (expr->expr_kind != EXPR_TYPEINFO) return expr; diff --git a/src/compiler/sema_casts.c b/src/compiler/sema_casts.c index ee6ca92b2..055756157 100644 --- a/src/compiler/sema_casts.c +++ b/src/compiler/sema_casts.c @@ -301,6 +301,22 @@ static bool int_conversion(Expr *expr, CastKind kind, Type *canonical, Type *typ return true; } +static bool int_vector_conversion(Expr *expr, Type *canonical, Type *type) +{ + // Fold pointer casts if narrowing + Type *base = type_get_indexed_type(type); + cast(expr, base); + return insert_cast(expr, CAST_NUMVEC, type); +} + +static bool float_vector_conversion(Expr *expr, Type *canonical, Type *type) +{ + // Fold pointer casts if narrowing + Type *base = type_get_indexed_type(type); + cast(expr, base); + return insert_cast(expr, CAST_NUMVEC, type); +} + /** * Cast a signed or unsigned integer -> floating point @@ -1231,6 +1247,9 @@ RETRY: case TYPE_BOOL: if (is_explicit) goto CAST; goto REQUIRE_CAST; + case TYPE_VECTOR: + to = to->array.base->canonical; + goto RETRY; case TYPE_ENUM: if (from->type_kind == TYPE_ENUM) break; to = to->decl->enums.type_info->type->canonical; @@ -1313,6 +1332,9 @@ RETRY: case TYPE_BOOL: if (is_explicit) goto CAST; goto REQUIRE_CAST; + case TYPE_VECTOR: + to = to->array.base->canonical; + goto RETRY; case ALL_FLOATS: { if (is_explicit) goto CAST; @@ -1850,6 +1872,7 @@ static bool cast_inner(Expr *expr, Type *from_type, Type *to, Type *to_type) if (to == type_bool) return integer_to_bool(expr, to_type); if (to->type_kind == TYPE_POINTER) return int_to_pointer(expr, to_type); if (to->type_kind == TYPE_ENUM) return integer_to_enum(expr, to, to_type); + if (type_kind_is_any_vector(to->type_kind)) return int_vector_conversion(expr, to, to_type); break; case ALL_UNSIGNED_INTS: if (type_is_integer_unsigned(to)) return int_conversion(expr, CAST_UIUI, to, to_type); @@ -1858,11 +1881,13 @@ static bool cast_inner(Expr *expr, Type *from_type, Type *to, Type *to_type) if (to == type_bool) return integer_to_bool(expr, to_type); if (to->type_kind == TYPE_POINTER) return int_to_pointer(expr, to_type); if (to->type_kind == TYPE_ENUM) return integer_to_enum(expr, to, to_type); + if (type_kind_is_any_vector(to->type_kind)) return int_vector_conversion(expr, to, to_type); break; case ALL_FLOATS: if (type_is_integer(to)) return float_to_integer(expr, to, to_type); if (to == type_bool) return float_to_bool(expr, to_type); if (type_is_float(to)) return float_to_float(expr, to, to_type); + if (type_kind_is_any_vector(to->type_kind)) return float_vector_conversion(expr, to, to_type); break; case TYPE_TYPEID: case TYPE_POINTER: diff --git a/src/compiler/tilde_codegen_expr.c b/src/compiler/tilde_codegen_expr.c index 6880bba1b..802595dd8 100644 --- a/src/compiler/tilde_codegen_expr.c +++ b/src/compiler/tilde_codegen_expr.c @@ -1456,6 +1456,9 @@ void tilde_emit_cast(TildeContext *c, CastKind cast_kind, Expr *expr, TBEValue * switch (cast_kind) { + case CAST_NUMVEC: + TODO // llvm_emit_num_to_vec_cast(c, value, to_type, from_type); + return; case CAST_BOOLVECINT: TODO // llvm_emit_bool_to_intvec_cast(c, value, to_type, from_type); return; diff --git a/src/compiler/types.c b/src/compiler/types.c index 0018fbe26..55a3a565e 100644 --- a/src/compiler/types.c +++ b/src/compiler/types.c @@ -1904,9 +1904,11 @@ Type *type_find_max_type(Type *type, Type *other) case ALL_INTS: if (other->type_kind == TYPE_DISTINCT && type_underlying_is_numeric(other)) return other; if (other->type_kind == TYPE_ENUM) return type_find_max_type(type, other->decl->enums.type_info->type->canonical); + if (other->type_kind == TYPE_VECTOR) return other; return type_find_max_num_type(type, other); case ALL_FLOATS: if (other->type_kind == TYPE_DISTINCT && type_is_float(other->decl->distinct_decl.base_type)) return other; + if (other->type_kind == TYPE_VECTOR) return other; return type_find_max_num_type(type, other); case TYPE_POINTER: if (type->pointer->type_kind == TYPE_ARRAY) diff --git a/src/version.h b/src/version.h index 014bf095e..f743f7cb7 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define COMPILER_VERSION "0.4.15" \ No newline at end of file +#define COMPILER_VERSION "0.4.16" \ No newline at end of file