From 99b8cf1d8e75c788b88463d4c1df3f1415162324 Mon Sep 17 00:00:00 2001 From: Szymon Nowakowski Date: Thu, 24 Sep 2026 20:24:14 +0200 Subject: [PATCH] Increase cossin accuracy --- packages/vecmath/src/trig.zig | 54 +++++++++++++++++++++++------------ 1 file changed, 36 insertions(+), 18 deletions(-) diff --git a/packages/vecmath/src/trig.zig b/packages/vecmath/src/trig.zig index 66e816f..8c7212e 100644 --- a/packages/vecmath/src/trig.zig +++ b/packages/vecmath/src/trig.zig @@ -147,15 +147,24 @@ test sin_x8 { pub fn cossin(angle_turns: f32) vm.Complex { @setFloatMode(.optimized); + // zig fmt: off // 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 + const term_cos_0: f32 = 1.000000000e-0; // 1 + const term_cos_2: f32 = -1.233700550e-0; // -π² /8 + const term_cos_4: f32 = 2.536695079e-1; // π⁴ /384 + const term_cos_6: f32 = -2.086348076e-2; // -π⁶ /46080 + const term_cos_8: f32 = 9.192602748e-4; // π⁸ /10321920 + const term_cos_10: f32 = -2.520204237e-5; // -π¹⁰/3715891200 + const term_cos_12: f32 = 4.710874778e-7; // π¹²/1961990553600 // 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 term_sin_1: f32 = 1.570796326e-0; // π /2 + const term_sin_3: f32 = -6.459640975e-1; // -π³ /48 + const term_sin_5: f32 = 7.969262624e-2; // π⁵ /3840 + const term_sin_7: f32 = -4.681754135e-3; // -π⁷ /645120 + const term_sin_9: f32 = 1.604411847e-4; // π⁹ /185794560 + const term_sin_11: f32 = -3.598843235e-6; // -π¹¹/81749606400 + const term_sin_13: f32 = 5.692172921e-8; // π¹³/51011754393600 + // zig fmt: on const angle_01 = angle_turns - @floor(angle_turns); const angle_04 = 4.0 * angle_01; @@ -168,8 +177,8 @@ pub fn cossin(angle_turns: f32) vm.Complex { const x = angle_04 - @floor(angle_04); const x2 = x * x; - const c: u32 = @bitCast(((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0); - const s: u32 = @bitCast(((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x); + const c: u32 = @bitCast((((((term_cos_12 * x2 + term_cos_10) * x2 + term_cos_8) * x2 + term_cos_6) * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0); + const s: u32 = @bitCast(((((((term_sin_13 * x2 + term_sin_11) * x2 + term_sin_9) * x2 + term_sin_7) * x2 + term_sin_5) * x2 + term_sin_3) * x2 + term_sin_1) * x); const result_cos: f32 = @bitCast(((s & quadrant_odd) | (c & ~quadrant_odd)) ^ sign_mask_cos); const result_sin: f32 = @bitCast(((c & quadrant_odd) | (s & ~quadrant_odd)) ^ sign_mask_sin); @@ -190,15 +199,24 @@ test cossin { pub fn cossin_x8(angle_turns: vm.f32x8) vm.Complex_x8 { @setFloatMode(.optimized); + // zig fmt: off // 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 + const term_cos_0 = vm.ps( 1.000000000e-0); // 1 + const term_cos_2 = vm.ps(-1.233700550e-0); // -π² /8 + const term_cos_4 = vm.ps( 2.536695079e-1); // π⁴ /384 + const term_cos_6 = vm.ps(-2.086348076e-2); // -π⁶ /46080 + const term_cos_8 = vm.ps( 9.192602748e-4); // π⁸ /10321920 + const term_cos_10 = vm.ps(-2.520204237e-5); // -π¹⁰/3715891200 + const term_cos_12 = vm.ps( 4.710874778e-7); // π¹²/1961990553600 // 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 term_sin_1 = vm.ps( 1.570796326e-0); // π /2 + const term_sin_3 = vm.ps(-6.459640975e-1); // -π³ /48 + const term_sin_5 = vm.ps( 7.969262624e-2); // π⁵ /3840 + const term_sin_7 = vm.ps(-4.681754135e-3); // -π⁷ /645120 + const term_sin_9 = vm.ps( 1.604411847e-4); // π⁹ /185794560 + const term_sin_11 = vm.ps(-3.598843235e-6); // -π¹¹/81749606400 + const term_sin_13 = vm.ps( 5.692172921e-8); // π¹³/51011754393600 + // zig fmt: on const angle_01 = angle_turns - @floor(angle_turns); const angle_04 = vm.ps(4.0) * angle_01; @@ -211,8 +229,8 @@ pub fn cossin_x8(angle_turns: vm.f32x8) vm.Complex_x8 { const x = angle_04 - @floor(angle_04); const x2 = x * x; - const c: vm.u32x8 = @bitCast(((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0); - const s: vm.u32x8 = @bitCast(((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x); + const c: vm.u32x8 = @bitCast((((((term_cos_12 * x2 + term_cos_10) * x2 + term_cos_8) * x2 + term_cos_6) * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0); + const s: vm.u32x8 = @bitCast(((((((term_sin_13 * x2 + term_sin_11) * x2 + term_sin_9) * x2 + term_sin_7) * x2 + term_sin_5) * x2 + term_sin_3) * x2 + term_sin_1) * x); const result_cos: vm.f32x8 = @bitCast(@select(u32, quadrant_odd, s, c) ^ sign_mask_cos); const result_sin: vm.f32x8 = @bitCast(@select(u32, quadrant_odd, c, s) ^ sign_mask_sin);