From f2ab89869d88a3350fa89571efc7ba969222656a Mon Sep 17 00:00:00 2001 From: Delyan Angelov Date: Fri, 21 Nov 2025 19:22:15 +0200 Subject: [PATCH] math: reduce the nested level of expressions, to reduce the chances of stack overflows on windows --- vlib/math/erf.v | 30 ++++++++++++++++-------------- vlib/math/gamma.v | 10 ++++++---- vlib/math/math.v | 6 ++++-- 3 files changed, 26 insertions(+), 20 deletions(-) diff --git a/vlib/math/erf.v b/vlib/math/erf.v index db08a084a..980c5c4cc 100644 --- a/vlib/math/erf.v +++ b/vlib/math/erf.v @@ -271,14 +271,15 @@ pub fn erf(a f64) f64 { mut r := 0.0 mut s := 0.0 if x < 1.0 / 0.35 { // |x| < 1 / 0.35 ~ 2.857143 - r = ra0 + s_ * (ra1 + s_ * (ra2 + s_ * (ra3 + s_ * (ra4 + s_ * (ra5 + s_ * (ra6 + - s_ * ra7)))))) - s = 1.0 + s_ * (sa1 + s_ * (sa2 + s_ * (sa3 + s_ * (sa4 + s_ * (sa5 + s_ * (sa6 + - s_ * (sa7 + s_ * sa8))))))) + tmp41 := s_ * (ra5 + s_ * (ra6 + s_ * ra7)) + r = ra0 + s_ * (ra1 + s_ * (ra2 + s_ * (ra3 + s_ * (ra4 + tmp41)))) + tmp42 := s_ * (sa5 + s_ * (sa6 + s_ * (sa7 + s_ * sa8))) + s = 1.0 + s_ * (sa1 + s_ * (sa2 + s_ * (sa3 + s_ * (sa4 + tmp42)))) } else { // |x| >= 1 / 0.35 ~ 2.857143 - r = rb0 + s_ * (rb1 + s_ * (rb2 + s_ * (rb3 + s_ * (rb4 + s_ * (rb5 + s_ * rb6))))) - s = 1.0 + s_ * (sb1 + s_ * (sb2 + s_ * (sb3 + s_ * (sb4 + s_ * (sb5 + s_ * (sb6 + - s_ * sb7)))))) + tmp31 := rb4 + s_ * (rb5 + s_ * rb6) + r = rb0 + s_ * (rb1 + s_ * (rb2 + s_ * (rb3 + s_ * tmp31))) + tmp32 := sb4 + s_ * (sb5 + s_ * (sb6 + s_ * sb7)) + s = 1.0 + s_ * (sb1 + s_ * (sb2 + s_ * (sb3 + s_ * tmp32))) } z := f64_from_bits(f64_bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x r_ := exp(-z * z - 0.5625) * exp((z - x) * (z + x) + r / s) @@ -346,17 +347,18 @@ pub fn erfc(a f64) f64 { mut r := 0.0 mut s := 0.0 if x < 1.0 / 0.35 { // |x| < 1 / 0.35 ~ 2.857143 - r = ra0 + s_ * (ra1 + s_ * (ra2 + s_ * (ra3 + s_ * (ra4 + s_ * (ra5 + s_ * (ra6 + - s_ * ra7)))))) - s = 1.0 + s_ * (sa1 + s_ * (sa2 + s_ * (sa3 + s_ * (sa4 + s_ * (sa5 + s_ * (sa6 + - s_ * (sa7 + s_ * sa8))))))) + tmp281 := ra4 + s_ * (ra5 + s_ * (ra6 + s_ * ra7)) + r = ra0 + s_ * (ra1 + s_ * (ra2 + s_ * (ra3 + s_ * tmp281))) + tmp282 := sa4 + s_ * (sa5 + s_ * (sa6 + s_ * (sa7 + s_ * sa8))) + s = 1.0 + s_ * (sa1 + s_ * (sa2 + s_ * (sa3 + s_ * tmp282))) } else { // |x| >= 1 / 0.35 ~ 2.857143 if sign && x > 6 { return 2.0 // x < -6 } - r = rb0 + s_ * (rb1 + s_ * (rb2 + s_ * (rb3 + s_ * (rb4 + s_ * (rb5 + s_ * rb6))))) - s = 1.0 + s_ * (sb1 + s_ * (sb2 + s_ * (sb3 + s_ * (sb4 + s_ * (sb5 + s_ * (sb6 + - s_ * sb7)))))) + tmp291 := rb3 + s_ * (rb4 + s_ * (rb5 + s_ * rb6)) + r = rb0 + s_ * (rb1 + s_ * (rb2 + s_ * tmp291)) + tmp292 := sb3 + s_ * (sb4 + s_ * (sb5 + s_ * (sb6 + s_ * sb7))) + s = 1.0 + s_ * (sb1 + s_ * (sb2 + s_ * tmp292)) } z := f64_from_bits(f64_bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x r_ := exp(-z * z - 0.5625) * exp((z - x) * (z + x) + r / s) diff --git a/vlib/math/gamma.v b/vlib/math/gamma.v index 7edec56d9..015c81a56 100644 --- a/vlib/math/gamma.v +++ b/vlib/math/gamma.v @@ -245,10 +245,12 @@ pub fn log_gamma_sign(a f64) (f64, int) { } else if x < 8 { // 2 <= x < 8 i := int(x) y := x - f64(i) - p := y * (lgamma_s[0] + y * (lgamma_s[1] + y * (lgamma_s[2] + y * (lgamma_s[3] + - y * (lgamma_s[4] + y * (lgamma_s[5] + y * lgamma_s[6])))))) - q := 1.0 + y * (lgamma_r[1] + y * (lgamma_r[2] + y * (lgamma_r[3] + y * (lgamma_r[4] + - y * (lgamma_r[5] + y * lgamma_r[6]))))) + tmp23456 := lgamma_s[2] + y * (lgamma_s[3] + y * (lgamma_s[4] + y * (lgamma_s[5] + + y * lgamma_s[6]))) + p := y * (lgamma_s[0] + y * (lgamma_s[1] + y * tmp23456)) + tmpr23456 := lgamma_r[2] + y * (lgamma_r[3] + y * (lgamma_r[4] + y * (lgamma_r[5] + + y * lgamma_r[6]))) + q := 1.0 + y * (lgamma_r[1] + y * tmpr23456) lgamma = 0.5 * y + p / q mut z := 1.0 // lgamma(1+s) = log(s) + lgamma(s) if i == 7 { diff --git a/vlib/math/math.v b/vlib/math/math.v index 23dd3391f..a639c63d1 100644 --- a/vlib/math/math.v +++ b/vlib/math/math.v @@ -13,7 +13,8 @@ pub fn aprox_sin(a f64) f64 { a5 := 2.08026600266304389e-2 a6 := -3.03996055049204407e-3 a7 := 1.38235642404333740e-4 - return a0 + a * (a1 + a * (a2 + a * (a3 + a * (a4 + a * (a5 + a * (a6 + a * a7)))))) + tmp := a4 + a * (a5 + a * (a6 + a * a7)) + return a0 + a * (a1 + a * (a2 + a * (a3 + a * tmp))) } // aprox_cos returns an approximation of cos(a) made using lolremez @@ -27,7 +28,8 @@ pub fn aprox_cos(a f64) f64 { a6 := -3.8510875386947414e-3 a7 := 4.7196604604366623e-4 a8 := -1.8776444013090451e-5 - return a0 + a * (a1 + a * (a2 + a * (a3 + a * (a4 + a * (a5 + a * (a6 + a * (a7 + a * a8))))))) + tmp := a4 + a * (a5 + a * (a6 + a * (a7 + a * a8))) + return a0 + a * (a1 + a * (a2 + a * (a3 + a * tmp))) } // copysign returns a value with the magnitude of x and the sign of y -- 2.39.5