Skip to content
This repository was archived by the owner on Apr 28, 2025. It is now read-only.

Generic cbrt #516

Closed
wants to merge 1 commit into from
+261 −106
Closed

Generic cbrt #516

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion crates/libm-test/src/domain.rs
View file
Original file line number Diff line number Diff line change
@@ -171,7 +171,7 @@ impl<F: Float, I: Int> EitherPrim<Domain<F>, Domain<I>> {
Box::new((0..u8::MAX).map(|scale| {
let mut base = F::ZERO;
for _ in 0..scale {
base = base - F::ONE;
base -= F::ONE;
}
base
}))
7 changes: 7 additions & 0 deletions crates/libm-test/src/f8_impl.rs
View file
Original file line number Diff line number Diff line change
@@ -39,6 +39,13 @@ impl Float for f8 {
const NEG_PI: Self = Self::ZERO;
const FRAC_PI_2: Self = Self::ZERO;

/// `2^sig_bits`
const TWO_POW_SIG_BITS: Self =
Self(((Self::SIG_BITS + Self::EXP_BIAS) as Self::Int) << Self::SIG_BITS);
/// `2^-sig_bits`
const TWO_POW_NEG_SIG_BITS: Self =
Self(((-(Self::SIG_BITS as i32) + Self::EXP_BIAS as i32) as Self::Int) << Self::SIG_BITS);

const BITS: u32 = 8;
const SIG_BITS: u32 = 3;
const SIGN_MASK: Self::Int = 0b1_0000_000;
23 changes: 23 additions & 0 deletions etc/consts-cbrt.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@

using Printf
using Remez

function main()
run_one("f64", "hf64!", 53)
end

function run_one(name::String, hf::String, precision::Integer)
setprecision(precision)

println("Constants for ", name)

println("const ESCALE: [Self; 3] = [")
for n in 0:2
val = big(2) ^ (n / 3)
@printf " %s(\"%a\"),\n" hf val
end
print("];")

end

main()
4 changes: 4 additions & 0 deletions etc/update-api-list.py
View file
Original file line number Diff line number Diff line change
@@ -25,6 +25,8 @@

# These files do not trigger a retest.
IGNORED_SOURCES = ["src/libm_helper.rs"]
# Same as above, limited to specific functions
IGNORED_SOURCES_MAP = {"fma": ["src/math/cbrt.rs"]}

IndexTy: TypeAlias = dict[str, dict[str, Any]]
"""Type of the `index` item in rustdoc's JSON output"""
@@ -138,6 +140,8 @@ def _init_defs(self, index: IndexTy) -> None:

for src in IGNORED_SOURCES:
sources.discard(src)
for src in IGNORED_SOURCES_MAP.get(name, []):
sources.discard(src)

# Sort the set
self.defs = {k: sorted(v) for (k, v) in defs.items()}
316 changes: 211 additions & 105 deletions src/math/cbrt.rs
View file
Original file line number Diff line number Diff line change
@@ -5,208 +5,314 @@
*/

use super::Float;
use super::support::{FpResult, Round, cold_path};
use super::support::{CastFrom, FpResult, Int, MinInt, Round, cold_path};

/// Compute the cube root of the argument.
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn cbrt(x: f64) -> f64 {
cbrt_round(x, Round::Nearest).val
}

pub fn cbrt_round(x: f64, round: Round) -> FpResult<f64> {
const ESCALE: [f64; 3] = [
1.0,
hf64!("0x1.428a2f98d728bp+0"), /* 2^(1/3) */
hf64!("0x1.965fea53d6e3dp+0"), /* 2^(2/3) */
];

/* the polynomial c0+c1*x+c2*x^2+c3*x^3 approximates x^(1/3) on [1,2]
with maximal error < 9.2e-5 (attained at x=2) */
const C: [f64; 4] = [
hf64!("0x1.1b0babccfef9cp-1"),
hf64!("0x1.2c9a3e94d1da5p-1"),
hf64!("-0x1.4dc30b1a1ddbap-3"),
hf64!("0x1.7a8d3e4ec9b07p-6"),
];

let u0: f64 = hf64!("0x1.5555555555555p-2");
let u1: f64 = hf64!("0x1.c71c71c71c71cp-3");

let rsc = [1.0, -1.0, 0.5, -0.5, 0.25, -0.25];

let off = [hf64!("0x1p-53"), 0.0, 0.0, 0.0];

/* rm=0 for rounding to nearest, and other values for directed roundings */
let hx: u64 = x.to_bits();
let mut mant: u64 = hx & f64::SIG_MASK;
let sign: u64 = hx >> 63;

let mut e: u32 = (hx >> f64::SIG_BITS) as u32 & f64::EXP_SAT;

if ((e + 1) & f64::EXP_SAT) < 2 {
// /// Compute the cube root of the argument.
// #[cfg(f128_enabled)]
// #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
// pub fn cbrtf128(x: f128) -> f128 {
// cbrt_round(x, Round::Nearest).val
// }

/// Correctly rounded cube root.
///
/// Algorithm:
/// - Minimax initial approximation
/// - `F`-sized newton iteration
/// - `2xF`-sized newton iteration
pub fn cbrt_round<F: Float + CbrtHelper>(x: F, round: Round) -> FpResult<F>
where
F::Int: CastFrom<u64>,
F::Int: CastFrom<u32>,
F::Int: From<u8>,
{
let zero = F::Int::ZERO;
let one = F::Int::ONE;
let u0: F = F::U0;
let u1: F = F::U1;
let off = F::OFF;

let hx = x.to_bits();
let mut mant: F::Int = hx & F::SIG_MASK;
let sign: F::Int = hx & F::SIGN_MASK;
let neg = x.is_sign_negative();

let mut e: u32 = x.exp();

// Handle 0, infinity, NaN, and subnormals
if ((e + 1) & F::EXP_SAT) < 2 {
cold_path();

let ix: u64 = hx & !f64::SIGN_MASK;
let ix = hx & !F::SIGN_MASK;

/* 0, inf, nan: we return x + x instead of simply x,
to that for x a signaling NaN, it correctly triggers
the invalid exception. */
if e == f64::EXP_SAT || ix == 0 {
if e == F::EXP_SAT || ix == zero {
// 0, infinity, NaN; use x + x to trigger exceptions
return FpResult::ok(x + x);
}

let nz = ix.leading_zeros() - 11; /* subnormal */
// Normalize subnormals
let nz = ix.leading_zeros() - F::EXP_BITS;
mant <<= nz;
mant &= f64::SIG_MASK;
mant &= F::SIG_MASK;
e = e.wrapping_sub(nz - 1);
}

e = e.wrapping_add(3072);
let cvt1: u64 = mant | (0x3ffu64 << 52);
let mut cvt5: u64 = cvt1;
// Set the exponent to 0, z is now [1, 2)
let iz = mant | (F::Int::cast_from(F::EXP_BIAS) << F::SIG_BITS);

let et: u32 = e / 3;
let it: u32 = e % 3;

/* 2^(3k+it) <= x < 2^(3k+it+1), with 0 <= it <= 3 */
cvt5 += u64::from(it) << f64::SIG_BITS;
cvt5 |= sign << 63;
let zz: f64 = f64::from_bits(cvt5);
// 2^(3k+it) <= x < 2^(3k+it+1), with 0 <= it <= 3
// `zz` is `x` reduced to [1, 8)
let izz = (iz + (F::Int::cast_from(it) << F::SIG_BITS)) | sign;
let zz: F = F::from_bits(izz);

/* cbrt(x) = cbrt(zz)*2^(et-1365) where 1 <= zz < 8 */
let mut isc: u64 = ESCALE[it as usize].to_bits(); // todo: index
isc |= sign << 63;
let cvt2: u64 = isc;
let z: f64 = f64::from_bits(cvt1);
let isc = F::ESCALE[it as usize].to_bits() | sign;
let z: F = F::from_bits(iz);

/* cbrt(zz) = cbrt(z)*isc, where isc encodes 1, 2^(1/3) or 2^(2/3),
and 1 <= z < 2 */
let r: f64 = 1.0 / z;
let rr: f64 = r * rsc[((it as usize) << 1) | sign as usize];
let z2: f64 = z * z;
let c0: f64 = C[0] + z * C[1];
let c2: f64 = C[2] + z * C[3];
let mut y: f64 = c0 + z2 * c2;
let mut y2: f64 = y * y;
let r: F = F::ONE / z;
let rr: F = r * F::RSCALE[((it as usize) << 1) | neg as usize];
let z2: F = z * z;
let c0: F = F::C[0] + z * F::C[1];
let c2: F = F::C[2] + z * F::C[3];

/* y is an approximation of z^(1/3) */
let mut h: f64 = y2 * (y * r) - 1.0;
let mut y: F = c0 + z2 * c2;
let mut y2: F = y * y;

/* h determines the error between y and z^(1/3) */
y -= (h * y) * (u0 - u1 * h);
let mut h: F = y2 * (y * r) - F::ONE;

/* The correction y -= (h*y)*(u0 - u1*h) corresponds to a cubic variant
of Newton's method, with the function f(y) = 1-z/y^3. */
y *= f64::from_bits(cvt2);
y -= (h * y) * (u0 - u1 * h);

y *= F::from_bits(isc);

/* Now y is an approximation of zz^(1/3),
* and rr an approximation of 1/zz. We now perform another iteration of
* Newton-Raphson, this time with a linear approximation only. */
y2 = y * y;
let mut y2l: f64 = fmaf64(y, y, -y2);
let mut y2l: F = y.fma(y, -y2);

/* y2 + y2l = y^2 exactly */
let mut y3: f64 = y2 * y;
let mut y3l: f64 = fmaf64(y, y2, -y3) + y * y2l;
let mut y3: F = y2 * y;
let mut y3l: F = y.fma(y2, -y3) + y * y2l;

/* y3 + y3l approximates y^3 with about 106 bits of accuracy */
h = ((y3 - zz) + y3l) * rr;
let mut dy: f64 = h * (y * u0);
let mut dy: F = h * (y * u0);

/* the approximation of zz^(1/3) is y - dy */
let mut y1: f64 = y - dy;
let mut y1: F = y - dy;
dy = (y - y1) - dy;

/* the approximation of zz^(1/3) is now y1 + dy, where |dy| < 1/2 ulp(y)
* (for rounding to nearest) */
let mut ady: f64 = dy.abs();
let mut ady: F = dy.abs();

/* For directed roundings, ady0 is tiny when dy is tiny, or ady0 is near
* from ulp(1);
* for rounding to nearest, ady0 is tiny when dy is near from 1/2 ulp(1),
* or from 3/2 ulp(1). */
let mut ady0: f64 = (ady - off[round as usize]).abs();
let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();
let mut ady0: F = (ady - off[round as usize]).abs();
let mut ady1: F = (ady - (F::TWO_POW_NEG_SIG_BITS + off[round as usize])).abs();

if ady0 < hf64!("0x1p-75") || ady1 < hf64!("0x1p-75") {
let magic = F::from_parts(false, (-75 + F::EXP_BIAS as i32) as u32, zero);

if ady0 < magic || ady1 < magic {
cold_path();

y2 = y1 * y1;
y2l = fmaf64(y1, y1, -y2);
y2l = y1.fma(y1, -y2);
y3 = y2 * y1;
y3l = fmaf64(y1, y2, -y3) + y1 * y2l;
y3l = y1.fma(y2, -y3) + y1 * y2l;
h = ((y3 - zz) + y3l) * rr;
dy = h * (y1 * u0);
y = y1 - dy;
dy = (y1 - y) - dy;
y1 = y;
ady = dy.abs();
ady0 = (ady - off[round as usize]).abs();
ady1 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();
ady1 = (ady - (F::TWO_POW_NEG_SIG_BITS + off[round as usize])).abs();

if ady0 < hf64!("0x1p-98") || ady1 < hf64!("0x1p-98") {
let magic2 = F::from_parts(false, (-98 + F::EXP_BIAS as i32) as u32, zero);
if ady0 < magic2 || ady1 < magic2 {
cold_path();
let azz: f64 = zz.abs();
let azz: F = zz.abs();

// ~ 0x1.79d15d0e8d59b80000000000000ffc3dp+0
if azz == hf64!("0x1.9b78223aa307cp+1") {
y1 = hf64!("0x1.79d15d0e8d59cp+0").copysign(zz);
if azz == F::AZMAGIC1 {
y1 = F::AZMAGIC2.copysign(zz);
}

// ~ 0x1.de87aa837820e80000000000001c0f08p+0
if azz == hf64!("0x1.a202bfc89ddffp+2") {
y1 = hf64!("0x1.de87aa837820fp+0").copysign(zz);
if azz == F::AZMAGIC3 {
y1 = F::AZMAGIC4.copysign(zz);
}

if round != Round::Nearest {
let wlist = [
(hf64!("0x1.3a9ccd7f022dbp+0"), hf64!("0x1.1236160ba9b93p+0")), // ~ 0x1.1236160ba9b930000000000001e7e8fap+0
(hf64!("0x1.7845d2faac6fep+0"), hf64!("0x1.23115e657e49cp+0")), // ~ 0x1.23115e657e49c0000000000001d7a799p+0
(hf64!("0x1.d1ef81cbbbe71p+0"), hf64!("0x1.388fb44cdcf5ap+0")), // ~ 0x1.388fb44cdcf5a0000000000002202c55p+0
(hf64!("0x1.0a2014f62987cp+1"), hf64!("0x1.46bcbf47dc1e8p+0")), // ~ 0x1.46bcbf47dc1e8000000000000303aa2dp+0
(hf64!("0x1.fe18a044a5501p+1"), hf64!("0x1.95decfec9c904p+0")), // ~ 0x1.95decfec9c9040000000000000159e8ep+0
(hf64!("0x1.a6bb8c803147bp+2"), hf64!("0x1.e05335a6401dep+0")), // ~ 0x1.e05335a6401de00000000000027ca017p+0
(hf64!("0x1.ac8538a031cbdp+2"), hf64!("0x1.e281d87098de8p+0")), // ~ 0x1.e281d87098de80000000000000ee9314p+0
];

for (a, b) in wlist {
for (a, b) in F::WLIST {
if azz == a {
let tmp = if round as u64 + sign == 2 { hf64!("0x1p-52") } else { 0.0 };
let tmp = if F::Int::from(round as u8 + neg as u8) == F::Int::cast_from(2) {
F::TWO_POW_NEG_SIG_BITS
} else {
F::ZERO
};
y1 = (b + tmp).copysign(zz);
}
}
}
}
}

let mut cvt3: u64 = y1.to_bits();
cvt3 = cvt3.wrapping_add(((et.wrapping_sub(342).wrapping_sub(1023)) as u64) << 52);
let m0: u64 = cvt3 << 30;
let mut cvt3 = y1.to_bits();
cvt3 = cvt3.wrapping_add((F::Int::cast_from(et.wrapping_sub(342).wrapping_sub(1023))) << 52);
let m0 = cvt3 << 30;
let m1 = m0 >> 63;

if (m0 ^ m1) <= (1u64 << 30) {
if (m0 ^ m1) <= (one << 30) {
cold_path();

let mut cvt4: u64 = y1.to_bits();
cvt4 = (cvt4 + (164 << 15)) & 0xffffffffffff0000u64;
let mut cvt4 = y1.to_bits();
cvt4 = (cvt4 + (F::Int::cast_from(164) << 15)) & F::Int::cast_from(0xffffffffffff0000u64);

if ((f64::from_bits(cvt4) - y1) - dy).abs() < hf64!("0x1p-60") || (zz).abs() == 1.0 {
cvt3 = (cvt3 + (1u64 << 15)) & 0xffffffffffff0000u64;
let magic3 = F::from_parts(false, (-60 + F::EXP_BIAS as i32) as u32, zero);
if ((F::from_bits(cvt4) - y1) - dy).abs() < magic3 || (zz).abs() == F::ONE {
cvt3 = (cvt3 + (one << 15)) & F::Int::cast_from(0xffffffffffff0000u64);
}
}

FpResult::ok(f64::from_bits(cvt3))
FpResult::ok(F::from_bits(cvt3))
}

pub trait CbrtHelper: Float {
/// 2^(n / 3) for n = [0, 1, 2]
const ESCALE: [Self; 3];
/// The polynomial `c0+c1*x+c2*x^2+c3*x^3` approximates `x^(1/3)` on `[1,2]`
/// with maximal error < 9.2e-5 (attained at x=2)
const C: [Self; 4];
const U0: Self;
const U1: Self;
const RSCALE: [Self; 6];
const OFF: [Self; 4];
const WLIST: [(Self, Self); 7];
const AZMAGIC1: Self;
const AZMAGIC2: Self;
const AZMAGIC3: Self;
const AZMAGIC4: Self;
fn fma(self, y: Self, z: Self) -> Self;
}

fn fmaf64(x: f64, y: f64, z: f64) -> f64 {
#[cfg(intrinsics_enabled)]
{
return unsafe { core::intrinsics::fmaf64(x, y, z) };
impl CbrtHelper for f64 {
const ESCALE: [Self; 3] = [
1.0,
hf64!("0x1.428a2f98d728bp+0"), /* 2^(1/3) */
hf64!("0x1.965fea53d6e3dp+0"), /* 2^(2/3) */
];

const C: [Self; 4] = [
// hf64!("0x1.1b850259b99ddp-1"),
// hf64!("0x1.2b9762efeffecp-1"),
// hf64!("-0x1.4af8eb64ea1ecp-3"),
// hf64!("0x1.7590cccfad50bp-6"),
hf64!("0x1.1b0babccfef9cp-1"),
hf64!("0x1.2c9a3e94d1da5p-1"),
hf64!("-0x1.4dc30b1a1ddbap-3"),
hf64!("0x1.7a8d3e4ec9b07p-6"),
];

// 0.33333333...
const U0: Self = hf64!("0x1.5555555555555p-2");

// 0.22222222...
const U1: Self = hf64!("0x1.c71c71c71c71cp-3");

const RSCALE: [Self; 6] = [1.0, -1.0, 0.5, -0.5, 0.25, -0.25];

const OFF: [Self; 4] = [hf64!("0x1p-53"), 0.0, 0.0, 0.0];

const WLIST: [(Self, Self); 7] = [
(hf64!("0x1.3a9ccd7f022dbp+0"), hf64!("0x1.1236160ba9b93p+0")), // ~ 0x1.1236160ba9b930000000000001e7e8fap+0
(hf64!("0x1.7845d2faac6fep+0"), hf64!("0x1.23115e657e49cp+0")), // ~ 0x1.23115e657e49c0000000000001d7a799p+0
(hf64!("0x1.d1ef81cbbbe71p+0"), hf64!("0x1.388fb44cdcf5ap+0")), // ~ 0x1.388fb44cdcf5a0000000000002202c55p+0
(hf64!("0x1.0a2014f62987cp+1"), hf64!("0x1.46bcbf47dc1e8p+0")), // ~ 0x1.46bcbf47dc1e8000000000000303aa2dp+0
(hf64!("0x1.fe18a044a5501p+1"), hf64!("0x1.95decfec9c904p+0")), // ~ 0x1.95decfec9c9040000000000000159e8ep+0
(hf64!("0x1.a6bb8c803147bp+2"), hf64!("0x1.e05335a6401dep+0")), // ~ 0x1.e05335a6401de00000000000027ca017p+0
(hf64!("0x1.ac8538a031cbdp+2"), hf64!("0x1.e281d87098de8p+0")), // ~ 0x1.e281d87098de80000000000000ee9314p+0
];

const AZMAGIC1: Self = hf64!("0x1.9b78223aa307cp+1");
const AZMAGIC2: Self = hf64!("0x1.79d15d0e8d59cp+0");
const AZMAGIC3: Self = hf64!("0x1.a202bfc89ddffp+2");
const AZMAGIC4: Self = hf64!("0x1.de87aa837820fp+0");

fn fma(self, y: Self, z: Self) -> Self {
#[cfg(intrinsics_enabled)]
{
return unsafe { core::intrinsics::fmaf64(self, y, z) };
}

#[cfg(not(intrinsics_enabled))]
{
return super::fma(self, y, z);
}
}
}

#[cfg(not(intrinsics_enabled))]
{
return super::fma(x, y, z);
#[cfg(f128_enabled)]
impl CbrtHelper for f128 {
const ESCALE: [Self; 3] = [
1.0,
hf128!("0x1.428a2f98d728acf826cc8664b665p+0"), /* 2^(1/3) */
hf128!("0x1.965fea53d6e3c53be1ca3482bf3ap+0"), /* 2^(2/3) */
];

const C: [Self; 4] = [
hf128!("0x1.1b850223b8bf644fcef50feeced1p-1"),
hf128!("0x1.2b97635e9e17d5240965cb56dc73p-1"),
hf128!("-0x1.4af8ec964bbc3767a6cf8ac456cbp-3"),
hf128!("0x1.7590ceecbb8c4c40d8c5e8b64d6bp-6"),
];

const U0: Self = 0.3333333333333333333333333333333333333333;

const U1: Self = 0.2222222222222222222222222222222222222222;

const RSCALE: [Self; 6] = [1.0, -1.0, 0.5, -0.5, 0.25, -0.25];

const OFF: [Self; 4] = [hf128!("0x1p-53"), 0.0, 0.0, 0.0];

// Other rounding modes unsupported for f128
const WLIST: [(Self, Self); 7] =
[(0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0)];

const AZMAGIC1: Self = hf128!("0x1.9b78223aa307cp+1");
const AZMAGIC2: Self = hf128!("0x1.79d15d0e8d59cp+0");
const AZMAGIC3: Self = hf128!("0x1.a202bfc89ddffp+2");
const AZMAGIC4: Self = hf128!("0x1.de87aa837820fp+0");

fn fma(self, y: Self, z: Self) -> Self {
#[cfg(intrinsics_enabled)]
{
return unsafe { core::intrinsics::fmaf128(self, y, z) };
}

#[cfg(not(intrinsics_enabled))]
{
return super::fmaf128(self, y, z);
}
}
}

1 change: 1 addition & 0 deletions src/math/mod.rs
View file
Original file line number Diff line number Diff line change
@@ -390,6 +390,7 @@ cfg_if! {
mod truncf128;
// verify-sorted-end

// pub use self::cbrt::cbrtf128;
// verify-sorted-start
pub use self::ceilf128::ceilf128;
pub use self::copysignf128::copysignf128;
14 changes: 14 additions & 0 deletions src/math/support/float_traits.rs
View file
Original file line number Diff line number Diff line change
@@ -10,6 +10,7 @@ pub trait Float:
+ PartialEq
+ PartialOrd
+ ops::AddAssign
+ ops::SubAssign
+ ops::MulAssign
+ ops::Add<Output = Self>
+ ops::Sub<Output = Self>
@@ -43,6 +44,11 @@ pub trait Float:

const MIN_POSITIVE_NORMAL: Self;

/// `2^sig_bits`, e.g. `0x1p-52` for `f64`. Used for normalization.
const TWO_POW_SIG_BITS: Self;
/// `2^-sig_bits`, e.g. `0x1p-52` for `f64`. Used for normalization.
const TWO_POW_NEG_SIG_BITS: Self;

/// The bitwidth of the float type
const BITS: u32;

@@ -205,6 +211,14 @@ macro_rules! float_impl {
// Exponent is a 1 in the LSB
const MIN_POSITIVE_NORMAL: Self = $from_bits(1 << Self::SIG_BITS);

/// `2^sig_bits`
const TWO_POW_SIG_BITS: Self =
$from_bits(((Self::SIG_BITS + Self::EXP_BIAS) as Self::Int) << Self::SIG_BITS);
/// `2^-sig_bits`
const TWO_POW_NEG_SIG_BITS: Self = $from_bits(
((-(Self::SIG_BITS as i32) + Self::EXP_BIAS as i32) as Self::Int) << Self::SIG_BITS,
);

const PI: Self = core::$ty::consts::PI;
const NEG_PI: Self = -Self::PI;
const FRAC_PI_2: Self = core::$ty::consts::FRAC_PI_2;