From dde8131c1e8ecb9c55d1b5f6888ce70d55119cd3 Mon Sep 17 00:00:00 2001 From: C Jones Date: Sat, 14 Jul 2018 19:23:01 -0400 Subject: [PATCH 1/2] Implement sinh This also adds expo2 for the __expo2 function, and combine_words() to replace the INSERT_WORDS macro. Closes rust-lang/libm#35 --- library/compiler-builtins/libm/src/lib.rs | 2 - .../compiler-builtins/libm/src/math/expo2.rs | 13 +++++ .../compiler-builtins/libm/src/math/mod.rs | 9 ++++ .../compiler-builtins/libm/src/math/sinh.rs | 48 +++++++++++++++++++ .../libm/test-generator/src/main.rs | 2 +- 5 files changed, 71 insertions(+), 3 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/expo2.rs create mode 100644 library/compiler-builtins/libm/src/math/sinh.rs diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0b9efeeb338a..ee864a36cf2e 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -414,7 +414,6 @@ pub trait F64Ext: private::Sealed { fn ln_1p(self) -> Self; - #[cfg(todo)] fn sinh(self) -> Self; #[cfg(todo)] @@ -595,7 +594,6 @@ impl F64Ext for f64 { log1p(self) } - #[cfg(todo)] #[inline] fn sinh(self) -> Self { sinh(self) diff --git a/library/compiler-builtins/libm/src/math/expo2.rs b/library/compiler-builtins/libm/src/math/expo2.rs new file mode 100644 index 000000000000..8a121b0a287a --- /dev/null +++ b/library/compiler-builtins/libm/src/math/expo2.rs @@ -0,0 +1,13 @@ +use super::{combine_words, exp}; + +/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */ +pub(crate) fn expo2(x: f64) -> f64 { + /* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */ + const K: i32 = 2043; + let kln2 = f64::from_bits(0x40962066151add8b); + + /* note that k is odd and scale*scale overflows */ + let scale = combine_words(((0x3ff + K / 2) as u32) << 20, 0); + /* exp(x - k ln2) * 2**(k-1) */ + return exp(x - kln2) * scale * scale; +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 7b3d9abee697..19d8ad58d866 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -50,6 +50,7 @@ mod scalbn; mod scalbnf; mod sin; mod sinf; +mod sinh; mod sqrt; mod sqrtf; mod tanf; @@ -100,6 +101,7 @@ pub use self::scalbn::scalbn; pub use self::scalbnf::scalbnf; pub use self::sin::sin; pub use self::sinf::sinf; +pub use self::sinh::sinh; pub use self::sqrt::sqrt; pub use self::sqrtf::sqrtf; pub use self::tanf::tanf; @@ -107,6 +109,7 @@ pub use self::trunc::trunc; pub use self::truncf::truncf; // Private modules +mod expo2; mod k_cos; mod k_cosf; mod k_sin; @@ -117,6 +120,7 @@ mod rem_pio2_large; mod rem_pio2f; // Private re-imports +use self::expo2::expo2; use self::k_cos::k_cos; use self::k_cosf::k_cosf; use self::k_sin::k_sin; @@ -151,3 +155,8 @@ pub fn with_set_low_word(f: f64, lo: u32) -> f64 { tmp |= lo as u64; f64::from_bits(tmp) } + +#[inline] +fn combine_words(hi: u32, lo: u32) -> f64 { + f64::from_bits((hi as u64) << 32 | lo as u64) +} diff --git a/library/compiler-builtins/libm/src/math/sinh.rs b/library/compiler-builtins/libm/src/math/sinh.rs new file mode 100644 index 000000000000..0579871719a8 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/sinh.rs @@ -0,0 +1,48 @@ +use super::{expm1, expo2}; + +// sinh(x) = (exp(x) - 1/exp(x))/2 +// = (exp(x)-1 + (exp(x)-1)/exp(x))/2 +// = x + x^3/6 + o(x^5) +// +pub fn sinh(x: f64) -> f64 { + // union {double f; uint64_t i;} u = {.f = x}; + // uint32_t w; + // double t, h, absx; + + let mut uf: f64 = x; + let mut ui: u64 = f64::to_bits(uf); + let w: u32; + let t: f64; + let mut h: f64; + let absx: f64; + + h = 0.5; + if ui >> 63 != 0 { + h = -h; + } + /* |x| */ + ui &= !1 / 2; + uf = f64::from_bits(ui); + absx = uf; + w = (ui >> 32) as u32; + + /* |x| < log(DBL_MAX) */ + if w < 0x40862e42 { + t = expm1(absx); + if w < 0x3ff00000 { + if w < 0x3ff00000 - (26 << 20) { + /* note: inexact and underflow are raised by expm1 */ + /* note: this branch avoids spurious underflow */ + return x; + } + return h * (2.0 * t - t * t / (t + 1.0)); + } + /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */ + return h * (t + t / (t + 1.0)); + } + + /* |x| > log(DBL_MAX) or nan */ + /* note: the result is stored to handle overflow */ + t = 2.0 * h * expo2(absx); + return t; +} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index c8e35cea7b00..758db9f527fb 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -714,7 +714,7 @@ f64_f64! { log2, round, sin, - // sinh, + sinh, sqrt, // tan, // tanh, From 0be652de3a29a7dcb5817c17f746543616980c1f Mon Sep 17 00:00:00 2001 From: C Jones Date: Sat, 14 Jul 2018 16:08:41 -0400 Subject: [PATCH 2/2] Implement tan Also includes implementing the private k_tan function. Closes rust-lang/libm#36 --- library/compiler-builtins/libm/src/lib.rs | 2 - .../compiler-builtins/libm/src/math/k_tan.rs | 105 ++++++++++++++++++ .../compiler-builtins/libm/src/math/mod.rs | 4 + .../compiler-builtins/libm/src/math/tan.rs | 69 ++++++++++++ .../libm/test-generator/src/main.rs | 2 +- 5 files changed, 179 insertions(+), 3 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/k_tan.rs create mode 100644 library/compiler-builtins/libm/src/math/tan.rs diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0b9efeeb338a..e6f286a9edc6 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -390,7 +390,6 @@ pub trait F64Ext: private::Sealed { fn cos(self) -> Self; - #[cfg(todo)] fn tan(self) -> Self; #[cfg(todo)] @@ -556,7 +555,6 @@ impl F64Ext for f64 { cos(self) } - #[cfg(todo)] #[inline] fn tan(self) -> Self { tan(self) diff --git a/library/compiler-builtins/libm/src/math/k_tan.rs b/library/compiler-builtins/libm/src/math/k_tan.rs new file mode 100644 index 000000000000..b0dd317a2e72 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/k_tan.rs @@ -0,0 +1,105 @@ +// origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */ +// +// ==================================================== +// Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. +// +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 +// Input x is assumed to be bounded by ~pi/4 in magnitude. +// Input y is the tail of x. +// Input odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned. +// +// Algorithm +// 1. Since tan(-x) = -tan(x), we need only to consider positive x. +// 2. Callers must return tan(-0) = -0 without calling here since our +// odd polynomial is not evaluated in a way that preserves -0. +// Callers may do the optimization tan(x) ~ x for tiny x. +// 3. tan(x) is approximated by a odd polynomial of degree 27 on +// [0,0.67434] +// 3 27 +// tan(x) ~ x + T1*x + ... + T13*x +// where +// +// |tan(x) 2 4 26 | -59.2 +// |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 +// | x | +// +// Note: tan(x+y) = tan(x) + tan'(x)*y +// ~ tan(x) + (1+x*x)*y +// Therefore, for better accuracy in computing tan(x+y), let +// 3 2 2 2 2 +// r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) +// then +// 3 2 +// tan(x+y) = x + (T1*x + (x *(r+y)+y)) +// +// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then +// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) +// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) +static T: [f64; 13] = [ + 3.33333333333334091986e-01, /* 3FD55555, 55555563 */ + 1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */ + 5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */ + 2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */ + 8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */ + 3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */ + 1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */ + 5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */ + 2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */ + 7.81794442939557092300e-05, /* 3F147E88, A03792A6 */ + 7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */ + -1.85586374855275456654e-05, /* BEF375CB, DB605373 */ + 2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */ +]; +const PIO4: f64 = 7.85398163397448278999e-01; /* 3FE921FB, 54442D18 */ +const PIO4_LO: f64 = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */ + +pub(crate) fn k_tan(mut x: f64, mut y: f64, odd: i32) -> f64 { + let hx = (f64::to_bits(x) >> 32) as u32; + let big = (hx & 0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */ + if big { + let sign = hx >> 31; + if sign != 0 { + x = -x; + y = -y; + } + x = (PIO4 - x) + (PIO4_LO - y); + y = 0.0; + } + let z = x * x; + let w = z * z; + /* + * Break x^5*(T[1]+x^2*T[2]+...) into + * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + + * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) + */ + let r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11])))); + let v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12]))))); + let s = z * x; + let r = y + z * (s * (r + v) + y) + s * T[0]; + let w = x + r; + if big { + let sign = hx >> 31; + let s = 1.0 - 2.0 * odd as f64; + let v = s - 2.0 * (x + (r - w * w / (w + s))); + return if sign != 0 { -v } else { v }; + } + if odd == 0 { + return w; + } + /* -1.0/(x+r) has up to 2ulp error, so compute it accurately */ + let w0 = zero_low_word(w); + let v = r - (w0 - x); /* w0+v = r+x */ + let a = -1.0 / w; + let a0 = zero_low_word(a); + a0 + a * (1.0 + a0 * w0 + a0 * v) +} + +#[inline] +fn zero_low_word(x: f64) -> f64 { + f64::from_bits(f64::to_bits(x) & 0xFFFF_FFFF_0000_0000) +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 7b3d9abee697..a6ce52394cf4 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -52,6 +52,7 @@ mod sin; mod sinf; mod sqrt; mod sqrtf; +mod tan; mod tanf; mod trunc; mod truncf; @@ -102,6 +103,7 @@ pub use self::sin::sin; pub use self::sinf::sinf; pub use self::sqrt::sqrt; pub use self::sqrtf::sqrtf; +pub use self::tan::tan; pub use self::tanf::tanf; pub use self::trunc::trunc; pub use self::truncf::truncf; @@ -111,6 +113,7 @@ mod k_cos; mod k_cosf; mod k_sin; mod k_sinf; +mod k_tan; mod k_tanf; mod rem_pio2; mod rem_pio2_large; @@ -121,6 +124,7 @@ use self::k_cos::k_cos; use self::k_cosf::k_cosf; use self::k_sin::k_sin; use self::k_sinf::k_sinf; +use self::k_tan::k_tan; use self::k_tanf::k_tanf; use self::rem_pio2::rem_pio2; use self::rem_pio2_large::rem_pio2_large; diff --git a/library/compiler-builtins/libm/src/math/tan.rs b/library/compiler-builtins/libm/src/math/tan.rs new file mode 100644 index 000000000000..92bbb6221229 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/tan.rs @@ -0,0 +1,69 @@ +// origin: FreeBSD /usr/src/lib/msun/src/s_tan.c */ +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +use super::{k_tan, rem_pio2}; + +// tan(x) +// Return tangent function of x. +// +// kernel function: +// k_tan ... tangent function on [-pi/4,pi/4] +// rem_pio2 ... argument reduction routine +// +// Method. +// Let S,C and T denote the sin, cos and tan respectively on +// [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 +// in [-pi/4 , +pi/4], and let n = k mod 4. +// We have +// +// n sin(x) cos(x) tan(x) +// ---------------------------------------------------------- +// 0 S C T +// 1 C -S -1/T +// 2 -S -C T +// 3 -C S -1/T +// ---------------------------------------------------------- +// +// Special cases: +// Let trig be any of sin, cos, or tan. +// trig(+-INF) is NaN, with signals; +// trig(NaN) is that NaN; +// +// Accuracy: +// TRIG(x) returns trig(x) nearly rounded +pub fn tan(x: f64) -> f64 { + let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120 + + let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff; + /* |x| ~< pi/4 */ + if ix <= 0x3fe921fb { + if ix < 0x3e400000 { + /* |x| < 2**-27 */ + /* raise inexact if x!=0 and underflow if subnormal */ + force_eval!(if ix < 0x00100000 { + x / x1p120 as f64 + } else { + x + x1p120 as f64 + }); + return x; + } + return k_tan(x, 0.0, 0); + } + + /* tan(Inf or NaN) is NaN */ + if ix >= 0x7ff00000 { + return x - x; + } + + /* argument reduction */ + let (n, y0, y1) = rem_pio2(x); + k_tan(y0, y1, n & 1) +} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index c8e35cea7b00..5d47cb2a45b6 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -716,7 +716,7 @@ f64_f64! { sin, // sinh, sqrt, - // tan, + tan, // tanh, trunc, fabs,