renati/castle
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: renati:af5f035830394d083fab9c6da4445e048bbc7808
Branches Tags
renati:main
...
compare: renati:99b8cf1d8e75c788b88463d4c1df3f1415162324
Branches Tags
renati:main

1 Commits
af5f035830 ... 99b8cf1d8e

Author SHA1 Message Date
Szymon Nowakowski
99b8cf1d8e Increase cossin accuracy 2026-09-24 20:24:14 +02:00
1 changed files with 36 additions and 18 deletions
Expand all files Collapse all files

54
packages/vecmath/src/trig.zig
View File

@@ -147,15 +147,24 @@ test sin_x8 {
pub fn cossin(angle_turns: f32) vm.Complex { pub fn cossin(angle_turns: f32) vm.Complex {
@setFloatMode(.optimized); @setFloatMode(.optimized);
// zig fmt: off
// Taylor series expansion for f(x)=cos(xπ/2) // Taylor series expansion for f(x)=cos(xπ/2)
const term_cos_0: f32 = 1.0; const term_cos_0: f32 = 1.000000000e-0; // 1
const term_cos_2: f32 = -1.23370055; // -π²/8 const term_cos_2: f32 = -1.233700550e-0; // -π² /8
const term_cos_4: f32 = 0.253669508; // π⁴/384 const term_cos_4: f32 = 2.536695079e-1; // π⁴ /384
const term_cos_6: f32 = -0.020863481; // -π⁶/46080 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) // Taylor series expansion for f(x)=sin(xπ/2)
const term_sin_1: f32 = 1.570796327; // π/2 const term_sin_1: f32 = 1.570796326e-0; // π /2
const term_sin_3: f32 = -0.645964098; // -π³/48 const term_sin_3: f32 = -6.459640975e-1; // -π³ /48
const term_sin_5: f32 = 0.079692626; // π⁵/3840 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_01 = angle_turns - @floor(angle_turns);
const angle_04 = 4.0 * angle_01; 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 x = angle_04 - @floor(angle_04);
const x2 = x * x; const x2 = x * x;
const c: u32 = @bitCast(((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0); 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_5 * x2 + term_sin_3) * x2 + term_sin_1) * x); 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_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); 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 { pub fn cossin_x8(angle_turns: vm.f32x8) vm.Complex_x8 {
@setFloatMode(.optimized); @setFloatMode(.optimized);
// zig fmt: off
// Taylor series expansion for f(x)=cos(xπ/2) // Taylor series expansion for f(x)=cos(xπ/2)
const term_cos_0 = vm.ps(1.0); const term_cos_0 = vm.ps( 1.000000000e-0); // 1
const term_cos_2 = vm.ps(-1.23370055); // -π²/8 const term_cos_2 = vm.ps(-1.233700550e-0); // -π² /8
const term_cos_4 = vm.ps(0.253669508); // π⁴/384 const term_cos_4 = vm.ps( 2.536695079e-1); // π⁴ /384
const term_cos_6 = vm.ps(-0.020863481); // -π⁶/46080 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) // Taylor series expansion for f(x)=sin(xπ/2)
const term_sin_1 = vm.ps(1.570796327); // π/2 const term_sin_1 = vm.ps( 1.570796326e-0); // π /2
const term_sin_3 = vm.ps(-0.645964098); // -π³/48 const term_sin_3 = vm.ps(-6.459640975e-1); // -π³ /48
const term_sin_5 = vm.ps(0.079692626); // π⁵/3840 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_01 = angle_turns - @floor(angle_turns);
const angle_04 = vm.ps(4.0) * angle_01; 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 x = angle_04 - @floor(angle_04);
const x2 = x * x; 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 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_5 * x2 + term_sin_3) * x2 + term_sin_1) * x); 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_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); const result_sin: vm.f32x8 = @bitCast(@select(u32, quadrant_odd, c, s) ^ sign_mask_sin);