math: cbrt fix (#14395)

David 'Epper' Marshall 2022-05-14 14:06:38 -04:00 committed by Jef Roosens
parent 46ee0ff572
commit 153b095518
Signed by: Jef Roosens
GPG Key ID: B75D4F293C7052DB
2 changed files with 21 additions and 2 deletions

View File

@ -1,5 +1,17 @@
module math
// The vlang code is a modified version of the original C code from
// http://www.netlib.org/fdlibm/s_cbrt.c and came with this notice.
//
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
//
// Developed at SunSoft, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.
// ====================================================
// cbrt returns the cube root of a.
//
// special cases are:
@ -25,12 +37,12 @@ pub fn cbrt(a f64) f64 {
sign = true
}
// rough cbrt to 5 bits
mut t := f64_from_bits(f64_bits(x) / u64(3 + (u64(b1) << 32)))
mut t := f64_from_bits(f64_bits(x) / u64(3) + (u64(b1) << 32))
if x < smallest_normal {
// subnormal number
t = f64(u64(1) << 54) // set t= 2**54
t *= x
t = f64_from_bits(f64_bits(t) / u64(3 + (u64(b2) << 32)))
t = f64_from_bits(f64_bits(t) / u64(3) + (u64(b2) << 32))
}
// new cbrt to 23 bits
mut r := t * t / x

View File

@ -506,6 +506,13 @@ fn test_mod() {
assert (0.6447968302508578) == f
}
fn test_cbrt() {
cbrts := [2.0, 10, 56]
for idx, i in [8.0, 1000, 175_616] {
assert cbrt(i) == cbrts[idx]
}
}
fn test_exp() {
for i := 0; i < math.vf_.len; i++ {
f := exp(math.vf_[i])