vecmath file split mayhem

packages/vecmath/src/trig.zig

@@ -0,0 +1,157 @@
+const std = @import("std");
+const vm = @import("root");
+
+pub fn cos(angle_turns: f32) f32 {
+    return cossin(angle_turns).re;
+}
+
+test "cos" {
+    try std.testing.expectEqual(1, cos(-1));
+    try std.testing.expectEqual(0, cos(-0.75));
+    try std.testing.expectEqual(-1, cos(-0.5));
+    try std.testing.expectEqual(0, cos(-0.25));
+    try std.testing.expectEqual(1, cos(0));
+    try std.testing.expectEqual(0, cos(0.25));
+    try std.testing.expectEqual(-1, cos(0.5));
+    try std.testing.expectEqual(0, cos(0.75));
+}
+
+pub fn cos_x8(angle_turns: vm.f32x8) vm.f32x8 {
+    return cossin_x8(angle_turns).re;
+}
+
+test "cos_x8" {
+    try std.testing.expectEqual(
+        .{ 1, 0, -1, 0, 1, 0, -1, 0 },
+        cos_x8(.{ -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 }),
+    );
+}
+
+pub fn sin(angle_turns: f32) f32 {
+    return cossin(angle_turns).im;
+}
+
+test "sin" {
+    try std.testing.expectEqual(0, sin(-1));
+    try std.testing.expectEqual(1, sin(-0.75));
+    try std.testing.expectEqual(0, sin(-0.5));
+    try std.testing.expectEqual(-1, sin(-0.25));
+    try std.testing.expectEqual(0, sin(0));
+    try std.testing.expectEqual(1, sin(0.25));
+    try std.testing.expectEqual(0, sin(0.5));
+    try std.testing.expectEqual(-1, sin(0.75));
+}
+
+pub fn sin_x8(angle_turns: vm.f32x8) vm.f32x8 {
+    return cossin_x8(angle_turns).im;
+}
+
+test "sin_x8" {
+    try std.testing.expectEqual(
+        .{ 0, 1, 0, -1, 0, 1, 0, -1 },
+        sin_x8(.{ -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 }),
+    );
+}
+
+pub fn cossin(angle_turns: f32) vm.Complex {
+    @setFloatMode(.optimized);
+    // Taylor series expansion for f(x)=cos(xπ/2)
+    const term_cos_0: f32 = 1.0;
+    const term_cos_2: f32 = -1.23370055; // -π²/8
+    const term_cos_4: f32 = 0.253669508; // π⁴/384
+    const term_cos_6: f32 = -0.020863481; // -π⁶/46080
+    // Taylor series expansion for f(x)=sin(xπ/2)
+    const term_sin_1: f32 = 1.570796327; // π/2
+    const term_sin_3: f32 = -0.645964098; // -π³/48
+    const term_sin_5: f32 = 0.079692626; // π⁵/3840
+
+    const angle_01 = angle_turns - @floor(angle_turns);
+    const angle_04 = 4.0 * angle_01;
+
+    const quadrant: u32 = @intFromFloat(angle_04);
+
+    const x = angle_04 - @floor(angle_04);
+    const x2 = x * x;
+
+    const c = ((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0;
+    const s = ((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x;
+
+    return switch (quadrant) {
+        0 => .init(c, s),
+        1 => .init(-s, c),
+        2 => .init(-c, -s),
+        3 => .init(s, -c),
+        else => unreachable,
+    };
+}
+
+test "cossin" {
+    try std.testing.expectEqual(vm.Vector2.unit_x, cossin(-1));
+    try std.testing.expectEqual(vm.Vector2.unit_y, cossin(-0.75));
+    try std.testing.expectEqual(vm.Vector2.unit_nx, cossin(-0.5));
+    try std.testing.expectEqual(vm.Vector2.unit_ny, cossin(-0.25));
+    try std.testing.expectEqual(vm.Vector2.unit_x, cossin(0));
+    try std.testing.expectEqual(vm.Vector2.unit_y, cossin(0.25));
+    try std.testing.expectEqual(vm.Vector2.unit_nx, cossin(0.5));
+    try std.testing.expectEqual(vm.Vector2.unit_ny, cossin(0.75));
+}
+
+pub fn cossin_x8(angle_turns: vm.f32x8) vm.Complex_x8 {
+    @setFloatMode(.optimized);
+    // Taylor series expansion for f(x)=cos(xπ/2)
+    const term_cos_0 = vm.ps(1.0);
+    const term_cos_2 = vm.ps(-1.23370055); // -π²/8
+    const term_cos_4 = vm.ps(0.253669508); // π⁴/384
+    const term_cos_6 = vm.ps(-0.020863481); // -π⁶/46080
+    // Taylor series expansion for f(x)=sin(xπ/2)
+    const term_sin_1 = vm.ps(1.570796327); // π/2
+    const term_sin_3 = vm.ps(-0.645964098); // -π³/48
+    const term_sin_5 = vm.ps(0.079692626); // π⁵/3840
+
+    const angle_01 = angle_turns - @floor(angle_turns);
+    const angle_04 = vm.ps(4.0) * angle_01;
+
+    const quadrant: vm.u32x8 = @intFromFloat(angle_04);
+    const quadrant_odd = (quadrant & vm.epu32(1)) != vm.epu32(0);
+    const sign_mask_cos = ((quadrant + vm.epu32(1)) & vm.epu32(0b10)) << @splat(30);
+    const sign_mask_sin = (quadrant & vm.epu32(0b10)) << @splat(30);
+
+    const x = angle_04 - @floor(angle_04);
+    const x2 = x * x;
+
+    const c = ((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0;
+    const s = ((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x;
+
+    var result_cos = @select(f32, quadrant_odd, s, c);
+    var result_sin = @select(f32, quadrant_odd, c, s);
+
+    result_cos = @bitCast(@as(vm.u32x8, @bitCast(result_cos)) ^ sign_mask_cos);
+    result_sin = @bitCast(@as(vm.u32x8, @bitCast(result_sin)) ^ sign_mask_sin);
+
+    return .init(result_cos, result_sin);
+}
+
+test "cossin_x8" {
+    try std.testing.expectEqual(
+        vm.Vector2x8.initArrayOfVectors(.{
+            vm.Vector2.unit_x,
+            vm.Vector2.unit_y,
+            vm.Vector2.unit_nx,
+            vm.Vector2.unit_ny,
+            vm.Vector2.unit_x,
+            vm.Vector2.unit_y,
+            vm.Vector2.unit_nx,
+            vm.Vector2.unit_ny,
+        }),
+        cossin_x8(.{
+            -1,
+            -0.75,
+            -0.5,
+            -0.25,
+            0,
+            0.25,
+            0.5,
+            0.75,
+        }),
+    );
+}