jc/math/vec.jai
jesse 526e1c8e8e More Vector and Matrix math and their tests
added common.jai for common math procedures
Some common procedures for smaller fixed vector sizes were made more optimal. SIMD coming later
Added tests to the math/module.jai
2025-05-31 14:46:39 -07:00

645 lines
16 KiB
Text

/*
Vec is a generic set of N named values of type T (aka. a mathematical vector)
The values can be accessed via array index or their common component names:
- u, v, d
- x, y, z, w
- r, g, b, a
- x, y, width, height
- min_x, min_y, max_x, max_y
- c0, c1, c2, c3, ... cN
For most use cases, opt for the named variants of this type:
Vec2 : Vec(2, float)
Vec3 : Vec(3, float)
Vec4 : Vec(4, float)
Quat : Vec(4, float)
Sadly, Vec does *NOT* solve the problem of interfacing with
external vector types (including Jai's 'Math' module).
To fix this, create two helper macros:
to_jai :: (v: Vec4) -> jmath.Vector4 #expand {
return v.(jmath.Vector4,force);
}
to_jai :: (v: *Vec4) -> *jmath.Vector4 #expand {
return v.(*jmath.Vector4);
}
*/
Vec :: struct(N: int, T: Type)
#modify {
info := T.(*Type_Info);
return info.type == .INTEGER || info.type == .FLOAT, "Vec T must be a numeric type (int or float)";
} {
#assert (N > 0) "Vec N cannot be <= 0";
// Vecs are backed by an array internally. The #insert block below
// generates #place'd unions of named fields or cN fields when N is > 4.
#as components: [N]T;
#insert -> string {
b: basic.String_Builder;
basic.append(*b, "#place components;\n");
fields :: []string.[
string.[ "x", "r", "u", "min_x", "" ],
string.[ "y", "g", "v", "min_y", "" ],
string.[ "z", "b", "d", "max_x", "width" ],
string.[ "w", "a", "" , "max_y", "height" ],
];
for i: 0..N - 1 {
if i < fields.count {
basic.append(*b, "union {\n");
for field: fields[i] if field.count != 0 {
basic.print_to_builder(*b, "\t%: T = ---;\n", field);
}
basic.print_to_builder(*b, "\tc%: T = ---;\n", i);
basic.append(*b, "};\n");
}
else {
basic.print_to_builder(*b, "c%: T = ---;\n", i);
}
}
// compound accessors
basic.append(*b, "#place components;\n");
if N >= 4 {
basic.append(*b, "union {\n");
basic.append(*b, "\txy: Vec2 = ---;\n");
basic.append(*b, "\txyz: Vec3 = ---;\n");
basic.append(*b, "};");
}
else if N >= 3 {
basic.append(*b, "xy: T = ---;\n");
}
// Vec4 row accessors (Possibly becoming/used as SIMD lanes)
if N/4 > 1 {
basic.append(*b, "#place components;\n");
basic.print_to_builder(*b, "v4: [%]Vec4 = ---;\n", N/4);
}
// Matrix Accessors
if N == 4 { // Mat2
basic.append(*b, "#place components;\n");
for i: 0..3 {
basic.print_to_builder(*b, "_%0%: T = ---;\n", i/2, i%2);
}
}
if N == 16 { // Mat4
basic.append(*b, "#place components;\n");
for i: 0..15 {
basic.print_to_builder(*b, "_%0%: T = ---;\n", i/4, i%4);
}
}
return basic.builder_to_string(*b);
};
}
operator [] :: (v: Vec, $$idx: int) -> v.T #no_abc {
meta.check_bounds(idx, v.N);
return v.components[idx];
}
operator *[] :: (v: *Vec, $$idx: int) -> *v.T #no_abc {
meta.check_bounds(idx, v.N);
return *v.components[idx];
}
operator []= :: (v: *Vec, $$idx: int, value: v.T) #no_abc {
meta.check_bounds(idx, v.N);
v.components[idx] = value;
}
for_expansion :: (v: *Vec, body: Code, flags: For_Flags) #expand {
for i: 0..v.N - 1 {
`it_index := i;
#if flags & .POINTER {
`it := *v.components[i];
} else {
`it := v.components[i];
}
#insert,scope(body) body;
}
}
operator + :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x + r.x;
res.y = l.y + r.y;
#if l.N >= 3 then res.z = l.z + r.z;
#if l.N == 4 then res.w = l.w + r.w; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it + r[it_index];
}
return res;
}
operator - :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x - r.x;
res.y = l.y - r.y;
#if l.N >= 3 then res.z = l.z - r.z;
#if l.N == 4 then res.w = l.w - r.w; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it - r[it_index];
}
return res;
}
operator * :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x*r.x;
res.y = l.y*r.y;
#if l.N >= 3 then res.z = l.z*r.z;
#if l.N == 4 then res.w = l.w*r.w; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it * r[it_index];
}
return res;
}
operator / :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x/r.x;
res.y = l.y/r.y;
#if l.N >= 3 then res.z = l.z/r.z;
#if l.N == 4 then res.w = l.w/r.w; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it / r[it_index];
}
return res;
}
operator + :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc #symmetric
#modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x + r;
res.y = l.y + r;
#if l.N >= 3 then res.z = l.z + r;
#if l.N == 4 then res.w = l.w + r; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it + r;
}
return res;
}
operator - :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc
#modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x - r;
res.y = l.y - r;
#if l.N >= 3 then res.z = l.z - r;
#if l.N == 4 then res.w = l.w - r; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it - r;
}
return res;
}
operator - :: inline (l: $R, r: Vec) -> Vec(l.N, l.T) #no_abc
#modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l - r.x;
res.y = l - r.y;
#if r.N >= 3 then res.z = l - r.z;
#if r.N == 4 then res.w = l - r.w; // @todo(jesse): SIMD
}
else {
for l res[it_index] = r - it;
}
return res;
}
operator- :: inline(v: Vec) -> Vec(v.N, v.T) #no_abc {
res: Vec(v.N, v.T) = ---;
#if v.N <= 4 {
res.x = -v.x;
res.y = -v.y;
#if v.N >= 3 then res.z = -v.z;
#if v.N == 4 then res.w = -v.w; // @todo(jesse): SIMD
}
else {
for v res[it_index] = -it;
}
return res;
}
operator * :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc #symmetric
#modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x*r;
res.y = l.y*r;
#if l.N >= 3 then res.z = l.z*r;
#if l.N == 4 then res.w = l.w*r; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it*r;
}
return res;
}
operator / :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc
#modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = l.x/r;
res.y = l.y/r;
#if l.N >= 3 then res.z = l.z/r;
#if l.N == 4 then res.w = l.w/r; // @todo(jesse): SIMD
}
else {
for l res[it_index] = it/r;
}
return res;
}
operator == :: inline (l: Vec, r: Vec(l.N, l.T)) -> bool #no_abc {
#if l.N <= 4 {
res: bool = l.x == r.x && l.y == r.y;
#if l.N >= 3 then res &= l.z == r.z;
#if l.N == 4 then res &= l.w == r.w; // @todo(jesse): SIMD
return res;
}
else {
for l if it != r[it_index] return false;
return true;
}
}
min :: (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = ifx l.x < r.x then l.x else r.x;
res.y = ifx l.y < r.y then l.y else r.y;
#if l.N >= 3 then res.z = ifx l.z < r.z then l.z else r.z;
#if l.N == 4 then res.w = ifx l.w < r.w then l.w else r.w; // @todo(jesse): SIMD
}
else {
n := l.N - 1;
while n >= 0 {
if l[n] < r[n]
res[n] = l[n];
else
res[n] = r[n];
n -= 1;
}
}
return res;
}
max :: (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = ifx l.x > r.x then l.x else r.x;
res.y = ifx l.y > r.y then l.y else r.y;
#if l.N >= 3 then res.z = ifx l.z > r.z then l.z else r.z;
#if l.N == 4 then res.w = ifx l.w > r.w then l.w else r.w; // @todo(jesse): SIMD
}
else {
n := l.N - 1;
while n >= 0 {
if l[n] > r[n]
res[n] = l[n];
else
res[n] = r[n];
n -= 1;
}
}
return res;
}
min :: (l: Vec, r: l.T) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = ifx l.x > r then l.x else r;
res.y = ifx l.y > r then l.y else r;
#if l.N >= 3 then res.z = ifx l.z > r.z then l.z else r;
#if l.N == 4 then res.w = ifx l.w > r.w then l.w else r; // @todo(jesse): SIMD
}
else {
n := l.N - 1;
while n >= 0 {
if l[n] > r
res[n] = l[n];
else
res[n] = r;
n -= 1;
}
}
return res;
}
max :: (l: Vec, r: l.T) -> Vec(l.N, l.T) #no_abc {
res: Vec(l.N, l.T) = ---;
#if l.N <= 4 {
res.x = ifx l.x < r then l.x else r;
res.y = ifx l.y < r then l.y else r;
#if l.N >= 3 then res.z = ifx l.z < r then l.z else r;
#if l.N == 4 then res.w = ifx l.w < r then l.w else r; // @todo(jesse): SIMD
}
else {
n := l.N - 1;
while n >= 0 {
if l[n] < r
res[n] = l[n];
else
res[n] = r;
n -= 1;
}
}
return res;
}
// @todo(jesse): SIMD
ceil :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
r: Vec(l.N, l.T) = ---;
n := l.N - 1;
sign: float;
while n >= 0 {
int_part := l[n].(int);
add_part := ifx l[n] > int_part.(float) then 1.0 else 0.0;
r[n] = int_part + add_part;
n -= 1;
}
return r;
}
// @todo(jesse): SIMD
floor :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
r: Vec(l.N, l.T) = ---;
n := l.N - 1;
sign: float;
while n >= 0 {
int_part := l[n].(int);
sub_part := ifx l[n] < int_part.(float) then 1.0 else 0.0;
r[n] = int_part - sub_part;
n -= 1;
}
return r;
}
clamp :: (v: Vec, low: v.T, high: v.T) -> Vec(v.N, v.T) #no_abc {
r: Vec(v.N, v.T) = ---;
n := v.N - 1;
while n >= 0 {
if v[n] < low
r[n] = low;
else if v[n] > high
r[n] = high;
else
r[n] = v[n];
n -= 1;
}
return r;
}
dot :: (a: Vec, b: Vec(a.N, a.T)) -> a.T #no_abc {
sum: a.T;
n := a.N - 1;
while n >= 0 {
sum += a[n]*b[n];
n -= 1;
}
return sum;
}
length_squared :: (v: Vec) -> float #no_abc {
return inline dot(v, v);
}
length :: (v: Vec) -> float #no_abc {
return math.sqrt(inline dot(v, v));
}
abs :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
r: Vec(v.N, v.T) = ---;
n := v.N - 1;
while n >= 0 {
if v[n] < 0
r[n] = -v[n];
else
r[n] = v[n];
n -= 1;
}
return r;
}
// @todo(Jesse): SIMD
norm :: normalize;
normalize :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
inv_len := 1.0/length(v);
return inv_len*v;
}
lerp :: (a: Vec, b: Vec(a.N, a.T), t: float) -> Vec(a.N, a.T) #no_abc {
r: Vec(a.N, a.T) = ---;
n := a.N - 1;
while n >= 0 {
r[n] = a[n] + t*(b[n] - a[n]);
n -= 1;
}
return r;
}
// Note(Jesse): I don't think this is needed for bigger vectors
reflect :: (v: Vec3, p: Vec3) -> Vec3 #no_abc {
projection := p*dot(v, p)/length_squared(p);
return 2*projection - v;
}
reflect :: (v: Vec2, p: Vec2) -> Vec2 #no_abc {
projection := p*dot(v, p)/length_squared(p);
return 2*projection - v;
}
round :: (v: Vec($N, $T)) -> Vec(N, T) #no_abc
#modify { return meta.type_is_float(T), "Used non-float vector on round"; } {
r: Vec(N, T) = ---;
n := N - 1;
while n >= 0 {
if v[n] < 0
r[n] = (v[n] - 0.5).(int).(float);
else
r[n] = (v[n] + 0.5).(int).(float);
n -= 1;
}
return r;
}
// Concrete vector types (the usual cases)
Vec2 :: Vec(2, float);
Vec3 :: Vec(3, float);
Vec4 :: Vec(4, float);
Quat :: #type,distinct Vec4; // Note(Jesse): I Had to make this distinct, otherwise operators stomp on eachother
v2f :: (x: $T = 0, y: T = 0) -> Vec2
#modify { return meta.type_is_float(T), "use v2i for integer arguments"; }
#expand {
return .{ x = x, y = y };
}
v2i :: (x: $T = 0, y: T = 0) -> Vec(2, T)
#modify { return meta.type_is_integer(T), "use v2f for float arguments"; }
#expand {
return .{ x = x, y = y };
}
v3f :: (x: $T = 0, y: T = 0, z: T = 0) -> Vec3
#modify { return meta.type_is_float(T), "use v3i for integer arguments"; }
#expand {
return .{ x = x, y = y, z = z };
}
v3i :: (x: $T = 0, y: T = 0, z: T = 0) -> Vec(3, T)
#modify { return meta.type_is_integer(T), "use v3f for float arguments"; }
#expand {
return .{ x = x, y = y, z = z };
}
v4f :: (x: $T = 0, y: T = 0, z: T = 0, w: T = 0) -> Vec4
#modify { return meta.type_is_float(T), "use v4i for integer arguments"; }
#expand {
return .{ x = x, y = y, z = z, w = w };
}
v4f :: (v: Vec3, $$w: float) -> Vec4 #expand {
return .{ xyz=v, w=w };
}
v4i :: (x: $T = 0, y: T = 0, z: T = 0, w: T = 0) -> Vec(4, T)
#modify { return meta.type_is_integer(T), "use v4f for float arguments"; }
#expand {
return .{ x = x, y = y, z = z, w = w };
}
quat :: (x: float = 0, y: float = 0, z: float = 0, w: float = 1) -> Quat #expand {
return .{ x = x, y = y, z = z, w = w };
}
quat :: (xyz: Vec3, w: float) -> Quat #expand {
return .{xyz=xyz, w=w};
}
quat_identity :: Quat.{x=0, y=0, z=0, w=1};
operator* :: (a: Quat, b: Quat) -> Quat {
r: Quat = ---;
r.xyz = a.w*b.xyz + b.w*a.xyz + cross(a.xyz, b.xyz);
r.w = a.w*b.w - dot(a.xyz, b.xyz);
return r;
}
operator* :: (q: Quat, v: Vec3) -> Vec3 #symmetric {
r: Vec3 = 2*cross(q.xyz, v);
return v + q.w*r + cross(q.xyz, r);
}
operator* :: (q: Quat, r: $T) -> Quat #symmetric {
return .{x=q.x*r, y=q.y*r, z=q.z*r, w=q.w*r};
}
operator+ :: (a: Quat, b: Quat) -> Quat {
return .{x=a.x*b.x, y=a.y*b.y, z=a.z*b.z, w=a.w*b.w};
}
operator/ :: (l: Quat, r: $T) -> Quat {
return .{xyz=l.xyz/r, w=l.w/r};
}
operator== :: (l: Quat, r: Quat) -> bool {
return l.x == r.x && l.y == r.y && l.z == r.z && l.w == r.w;
}
conjugate :: (q: Quat) -> Quat #expand {
r: Quat = ---;
r.xyz = -q.xyz;
r.w = q.w;
return r;
}
// Note(Jesse): iff q is a unit vector
inverse_unit :: (q: Quat) -> Quat {
return inline conjugate(q);
}
inverse :: (q: Quat) -> Quat {
return conjugate(q)/length(q);
}
rotation_quat :: (theta: float, axis: Vec3) -> Quat {
q: Quat = ---;
a := normalize(axis);
q.xyz = sin(theta/2.0)*a;
q.w = cos(theta/2.0);
return q;
}
operator- :: (q: Quat) -> Quat {
r: Quat = ---;
r.x = -q.x;
r.y = -q.y;
r.z = -q.z;
r.w = -q.w;
return r;
}
length_squared :: (q: Quat) -> float {
return dot(q, q);
}
length :: (q: Quat) -> float {
return math.sqrt(dot(q, q));
}
dot :: (a: Quat, b: Quat) -> float {
return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
}
slerp :: (a: Quat, _b: Quat, t: float) -> Quat {
b := _b;
cos_a := dot(a, b);
if cos_a < 0.0 {
cos_a = -cos_a;
b = -b;
}
alpha := math.acos(cos_a);
sin_a := math.sin(alpha);
p1 := math.sin((1-t)*alpha)/sin_a;
p2 := math.sin(t*alpha)/sin_a;
return a*p1 + b*p2;
}
cross :: (a: Vec3, b: Vec3) -> Vec3 {
v: Vec3 = ---;
v.x = a.y*b.z - a.z*b.y;
v.y = a.z*b.x - a.x*b.z;
v.z = a.x*b.y - a.y*b.x;
return v;
}
#scope_file
meta :: #import "jc/meta";
math :: #import "Math"; // @future
basic :: #import "Basic"; // @future