diff --git a/library/compiler-builtins/libm/src/math/cosf.rs b/library/compiler-builtins/libm/src/math/cosf.rs index ef46d4a15f2d..79df97e354a3 100644 --- a/library/compiler-builtins/libm/src/math/cosf.rs +++ b/library/compiler-builtins/libm/src/math/cosf.rs @@ -3,32 +3,36 @@ use super::{k_cosf, k_sinf, rem_pio2f}; use core::f64::consts::FRAC_PI_2; /* Small multiples of pi/2 rounded to double precision. */ -const C1_PIO2 : f64 = 1.*FRAC_PI_2; /* 0x3FF921FB, 0x54442D18 */ -const C2_PIO2 : f64 = 2.*FRAC_PI_2; /* 0x400921FB, 0x54442D18 */ -const C3_PIO2 : f64 = 3.*FRAC_PI_2; /* 0x4012D97C, 0x7F3321D2 */ -const C4_PIO2 : f64 = 4.*FRAC_PI_2; /* 0x401921FB, 0x54442D18 */ +const C1_PIO2: f64 = 1. * FRAC_PI_2; /* 0x3FF921FB, 0x54442D18 */ +const C2_PIO2: f64 = 2. * FRAC_PI_2; /* 0x400921FB, 0x54442D18 */ +const C3_PIO2: f64 = 3. * FRAC_PI_2; /* 0x4012D97C, 0x7F3321D2 */ +const C4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */ #[inline] pub fn cosf(x: f32) -> f32 { let x64 = x as f64; - + let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120 - + let mut ix = x.to_bits(); let sign = (ix >> 31) != 0; ix &= 0x7fffffff; - if ix <= 0x3f490fda { /* |x| ~<= pi/4 */ - if ix < 0x39800000 { /* |x| < 2**-12 */ + if ix <= 0x3f490fda { + /* |x| ~<= pi/4 */ + if ix < 0x39800000 { + /* |x| < 2**-12 */ /* raise inexact if x != 0 */ force_eval!(x + x1p120); return 1.; } return k_cosf(x64); } - if ix <= 0x407b53d1 { /* |x| ~<= 5*pi/4 */ - if ix > 0x4016cbe3 { /* |x| ~> 3*pi/4 */ - return -k_cosf(if sign { x64+C2_PIO2 } else { x64-C2_PIO2 }); + if ix <= 0x407b53d1 { + /* |x| ~<= 5*pi/4 */ + if ix > 0x4016cbe3 { + /* |x| ~> 3*pi/4 */ + return -k_cosf(if sign { x64 + C2_PIO2 } else { x64 - C2_PIO2 }); } else { if sign { return k_sinf(x64 + C1_PIO2); @@ -37,9 +41,11 @@ pub fn cosf(x: f32) -> f32 { } } } - if ix <= 0x40e231d5 { /* |x| ~<= 9*pi/4 */ - if ix > 0x40afeddf { /* |x| ~> 7*pi/4 */ - return k_cosf(if sign { x64+C4_PIO2 } else { x64-C4_PIO2 }); + if ix <= 0x40e231d5 { + /* |x| ~<= 9*pi/4 */ + if ix > 0x40afeddf { + /* |x| ~> 7*pi/4 */ + return k_cosf(if sign { x64 + C4_PIO2 } else { x64 - C4_PIO2 }); } else { if sign { return k_sinf(-x64 - C3_PIO2); @@ -51,15 +57,15 @@ pub fn cosf(x: f32) -> f32 { /* cos(Inf or NaN) is NaN */ if ix >= 0x7f800000 { - return x-x; + return x - x; } /* general argument reduction needed */ let (n, y) = rem_pio2f(x); - match n&3 { - 0 => { k_cosf( y) }, - 1 => { k_sinf(-y) }, - 2 => { -k_cosf( y) }, - _ => { k_sinf( y) }, + match n & 3 { + 0 => k_cosf(y), + 1 => k_sinf(-y), + 2 => -k_cosf(y), + _ => k_sinf(y), } } diff --git a/library/compiler-builtins/libm/src/math/k_cosf.rs b/library/compiler-builtins/libm/src/math/k_cosf.rs index 5d1ede7af922..83d13b2e998a 100644 --- a/library/compiler-builtins/libm/src/math/k_cosf.rs +++ b/library/compiler-builtins/libm/src/math/k_cosf.rs @@ -1,13 +1,13 @@ /* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */ -const C0 : f64 = -0.499999997251031003120; /* -0x1ffffffd0c5e81.0p-54 */ -const C1 : f64 = 0.0416666233237390631894; /* 0x155553e1053a42.0p-57 */ -const C2 : f64 = -0.00138867637746099294692; /* -0x16c087e80f1e27.0p-62 */ -const C3 : f64 = 0.0000243904487962774090654; /* 0x199342e0ee5069.0p-68 */ +const C0: f64 = -0.499999997251031003120; /* -0x1ffffffd0c5e81.0p-54 */ +const C1: f64 = 0.0416666233237390631894; /* 0x155553e1053a42.0p-57 */ +const C2: f64 = -0.00138867637746099294692; /* -0x16c087e80f1e27.0p-62 */ +const C3: f64 = 0.0000243904487962774090654; /* 0x199342e0ee5069.0p-68 */ #[inline] -pub(crate) fn k_cosf(x : f64) -> f32 { - let z = x*x; - let w = z*z; - let r = C2+z*C3; - (((1.0+z*C0) + w*C1) + (w*z)*r) as f32 +pub(crate) fn k_cosf(x: f64) -> f32 { + let z = x * x; + let w = z * z; + let r = C2 + z * C3; + (((1.0 + z * C0) + w * C1) + (w * z) * r) as f32 } diff --git a/library/compiler-builtins/libm/src/math/k_sinf.rs b/library/compiler-builtins/libm/src/math/k_sinf.rs index 5eedab3aa64a..bb2183afc7ca 100644 --- a/library/compiler-builtins/libm/src/math/k_sinf.rs +++ b/library/compiler-builtins/libm/src/math/k_sinf.rs @@ -1,14 +1,14 @@ /* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */ -const S1 : f64 = -0.166666666416265235595; /* -0x15555554cbac77.0p-55 */ -const S2 : f64 = 0.0083333293858894631756; /* 0x111110896efbb2.0p-59 */ -const S3 : f64 = -0.000198393348360966317347; /* -0x1a00f9e2cae774.0p-65 */ -const S4 : f64 = 0.0000027183114939898219064; /* 0x16cd878c3b46a7.0p-71 */ +const S1: f64 = -0.166666666416265235595; /* -0x15555554cbac77.0p-55 */ +const S2: f64 = 0.0083333293858894631756; /* 0x111110896efbb2.0p-59 */ +const S3: f64 = -0.000198393348360966317347; /* -0x1a00f9e2cae774.0p-65 */ +const S4: f64 = 0.0000027183114939898219064; /* 0x16cd878c3b46a7.0p-71 */ #[inline] -pub(crate) fn k_sinf(x : f64) -> f32 { - let z = x*x; - let w = z*z; - let r = S3 + z*S4; - let s = z*x; - ((x + s*(S1 + z*S2)) + s*w*r) as f32 +pub(crate) fn k_sinf(x: f64) -> f32 { + let z = x * x; + let w = z * z; + let r = S3 + z * S4; + let s = z * x; + ((x + s * (S1 + z * S2)) + s * w * r) as f32 } diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 00945b258c89..013e70870bbe 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -7,6 +7,7 @@ macro_rules! force_eval { } mod ceilf; +mod cosf; mod expf; mod fabs; mod fabsf; @@ -27,28 +28,22 @@ mod scalbn; mod scalbnf; mod sqrt; mod sqrtf; -mod cosf; mod trunc; mod truncf; pub use self::{ - ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmodf::fmodf, - hypot::hypot, hypotf::hypotf, log10::log10, log10f::log10f, log2::log2, log2f::log2f, - logf::logf, powf::powf, round::round, roundf::roundf, scalbn::scalbn, scalbnf::scalbnf, - sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf, cosf::cosf, + ceilf::ceilf, cosf::cosf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, + fmodf::fmodf, hypot::hypot, hypotf::hypotf, log10::log10, log10f::log10f, log2::log2, + log2f::log2f, logf::logf, powf::powf, round::round, roundf::roundf, scalbn::scalbn, + scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf, }; mod k_cosf; mod k_sinf; -mod rem_pio2f; mod rem_pio2_large; +mod rem_pio2f; -use self::{ - k_cosf::k_cosf, - k_sinf::k_sinf, - rem_pio2f::rem_pio2f, - rem_pio2_large::rem_pio2_large, -}; +use self::{k_cosf::k_cosf, k_sinf::k_sinf, rem_pio2_large::rem_pio2_large, rem_pio2f::rem_pio2f}; fn isnanf(x: f32) -> bool { x.to_bits() & 0x7fffffff > 0x7f800000 diff --git a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs index e1f5d1333905..52b47279cd91 100644 --- a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs +++ b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs @@ -1,8 +1,8 @@ -use super::scalbn; use super::floor; +use super::scalbn; // initial value for jk -const INIT_JK : [usize; 4] = [3,4,4,6]; +const INIT_JK: [usize; 4] = [3, 4, 4, 6]; // Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi // @@ -15,140 +15,99 @@ const INIT_JK : [usize; 4] = [3,4,4,6]; // NB: This table must have at least (e0-3)/24 + jk terms. // For quad precision (e0 <= 16360, jk = 6), this is 686. #[cfg(target_pointer_width = "32")] -const IPIO2 : [i32; 66] = [ - 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, - 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, - 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, - 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, - 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, - 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, - 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, - 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, - 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, - 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, - 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, +const IPIO2: [i32; 66] = [ + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, + 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, + 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, + 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, + 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292, + 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, + 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, + 0x73A8C9, 0x60E27B, 0xC08C6B, ]; #[cfg(target_pointer_width = "64")] -const IPIO2 : [i32; 690] = [ - 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, - 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, - 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, - 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, - 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, - 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, - 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, - 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, - 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, - 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, - 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, - 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, - 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, - 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35, - 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30, - 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, - 0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4, - 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770, - 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, - 0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19, - 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522, - 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, - 0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6, - 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E, - 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, - 0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3, - 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF, - 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, - 0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612, - 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929, - 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, - 0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B, - 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C, - 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, - 0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB, - 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC, - 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, - 0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F, - 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5, - 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, - 0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B, - 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA, - 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, - 0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3, - 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3, - 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, - 0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F, - 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61, - 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, - 0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51, - 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0, - 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, - 0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6, - 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC, - 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, - 0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328, - 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D, - 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, - 0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B, - 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4, - 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, - 0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F, - 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD, - 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, - 0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4, - 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761, - 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, - 0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30, - 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262, - 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, - 0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1, - 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C, - 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, - 0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08, - 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196, - 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, - 0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4, - 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC, - 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, - 0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0, - 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C, - 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, - 0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC, - 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22, - 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, - 0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7, - 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5, - 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, - 0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4, - 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF, - 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, - 0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2, - 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138, - 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, - 0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569, - 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34, - 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, - 0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D, - 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F, - 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, - 0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569, - 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B, - 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, - 0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41, - 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49, - 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, - 0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110, - 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8, - 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, - 0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A, - 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270, - 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, - 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616, - 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, +const IPIO2: [i32; 690] = [ + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, + 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, + 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, + 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, + 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292, + 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, + 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, + 0x73A8C9, 0x60E27B, 0xC08C6B, 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, + 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, 0xDE4F98, 0x327DBB, 0xC33D26, + 0xEF6B1E, 0x5EF89F, 0x3A1F35, 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30, + 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, 0x467D86, 0x2D71E3, 0x9AC69B, + 0x006233, 0x7CD2B4, 0x97A7B4, 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770, + 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, 0xCB2324, 0x778AD6, 0x23545A, + 0xB91F00, 0x1B0AF1, 0xDFCE19, 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522, + 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, 0xDE3B58, 0x929BDE, 0x2822D2, + 0xE88628, 0x4D58E2, 0x32CAC6, 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E, + 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, 0xD36710, 0xD8DDAA, 0x425FAE, + 0xCE616A, 0xA4280A, 0xB499D3, 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF, + 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, 0x36D9CA, 0xD2A828, 0x8D61C2, + 0x77C912, 0x142604, 0x9B4612, 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929, + 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, 0xC3E7B3, 0x28F8C7, 0x940593, + 0x3E71C1, 0xB3092E, 0xF3450B, 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C, + 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, 0x9794E8, 0x84E6E2, 0x973199, + 0x6BED88, 0x365F5F, 0x0EFDBB, 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC, + 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, 0x90AA47, 0x02E774, 0x24D6BD, + 0xA67DF7, 0x72486E, 0xEF169F, 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5, + 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, 0x10D86D, 0x324832, 0x754C5B, + 0xD4714E, 0x6E5445, 0xC1090B, 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA, + 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, 0x6AE290, 0x89D988, 0x50722C, + 0xBEA404, 0x940777, 0x7030F3, 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3, + 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, 0x3BDF08, 0x2B3715, 0xA0805C, + 0x93805A, 0x921110, 0xD8E80F, 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61, + 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, 0xAA140A, 0x2F2689, 0x768364, + 0x333B09, 0x1A940E, 0xAA3A51, 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0, + 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, 0x5BC3D8, 0xC492F5, 0x4BADC6, + 0xA5CA4E, 0xCD37A7, 0x36A9E6, 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC, + 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, 0x306529, 0xBF5657, 0x3AFF47, + 0xB9F96A, 0xF3BE75, 0xDF9328, 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D, + 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, 0xA8654F, 0xA5C1D2, 0x0F3F0B, + 0xCD785B, 0x76F923, 0x048B7B, 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4, + 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, 0xDA4886, 0xA05DF7, 0xF480C6, + 0x2FF0AC, 0x9AECDD, 0xBC5C3F, 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD, + 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, 0x2A1216, 0x2DB7DC, 0xFDE5FA, + 0xFEDB89, 0xFDBE89, 0x6C76E4, 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761, + 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, 0x48D784, 0x16DF30, 0x432DC7, + 0x356125, 0xCE70C9, 0xB8CB30, 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262, + 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, 0xC4F133, 0x5F6E13, 0xE4305D, + 0xA92E85, 0xC3B21D, 0x3632A1, 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C, + 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, 0xCBDA11, 0xD0BE7D, 0xC1DB9B, + 0xBD17AB, 0x81A2CA, 0x5C6A08, 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196, + 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, 0x4F6A68, 0xA82A4A, 0x5AC44F, + 0xBCF82D, 0x985AD7, 0x95C7F4, 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC, + 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, 0xD0C0B2, 0x485551, 0x0EFB1E, + 0xC37295, 0x3B06A3, 0x3540C0, 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C, + 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, 0x3C3ABA, 0x461846, 0x5F7555, + 0xF5BDD2, 0xC6926E, 0x5D2EAC, 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22, + 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, 0x745D7C, 0xB2AD6B, 0x9D6ECD, + 0x7B723E, 0x6A11C6, 0xA9CFF7, 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5, + 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, 0xBEFDFD, 0xEF4556, 0x367ED9, + 0x13D9EC, 0xB9BA8B, 0xFC97C4, 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF, + 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, 0x9C2A3E, 0xCC5F11, 0x4A0BFD, + 0xFBF4E1, 0x6D3B8E, 0x2C86E2, 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138, + 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, 0xCC2254, 0xDC552A, 0xD6C6C0, + 0x96190B, 0xB8701A, 0x649569, 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34, + 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, 0x9B5861, 0xBC57E1, 0xC68351, + 0x103ED8, 0x4871DD, 0xDD1C2D, 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F, + 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, 0x382682, 0x9BE7CA, 0xA40D51, + 0xB13399, 0x0ED7A9, 0x480569, 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B, + 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, 0x5FD45E, 0xA4677B, 0x7AACBA, + 0xA2F655, 0x23882B, 0x55BA41, 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49, + 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, 0xAE5ADB, 0x86C547, 0x624385, + 0x3B8621, 0x94792C, 0x876110, 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8, + 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, 0xB1933D, 0x0B7CBD, 0xDC51A4, + 0x63DD27, 0xDDE169, 0x19949A, 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270, + 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, 0x4D7E6F, 0x5119A5, 0xABF9B5, + 0xD6DF82, 0x61DD96, 0x023616, 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0, ]; -const PIO2 : [f64; 8] = [ +const PIO2: [f64; 8] = [ 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ @@ -251,92 +210,91 @@ const PIO2 : [f64; 8] = [ /// more accurately, = 0 mod 8 ). Thus the number of operations are /// independent of the exponent of the input. #[inline] -pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) -> i32 { - let x1p24 = f64::from_bits(0x4170000000000000); // 0x1p24 === 2 ^ 24 +pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32 { + let x1p24 = f64::from_bits(0x4170000000000000); // 0x1p24 === 2 ^ 24 let x1p_24 = f64::from_bits(0x3e70000000000000); // 0x1p_24 === 2 ^ (-24) - + #[cfg(target_pointer_width = "64")] assert!(e0 <= 16360); - + let nx = x.len(); - - let mut fw : f64; - let mut n : i32; - let mut ih : i32; - let mut z : f64; - let mut f : [f64;20] = [0.;20]; - let mut fq : [f64;20] = [0.;20]; - let mut q : [f64;20] = [0.;20]; - let mut iq : [i32;20] = [0;20]; - + + let mut fw: f64; + let mut n: i32; + let mut ih: i32; + let mut z: f64; + let mut f: [f64; 20] = [0.; 20]; + let mut fq: [f64; 20] = [0.; 20]; + let mut q: [f64; 20] = [0.; 20]; + let mut iq: [i32; 20] = [0; 20]; + /* initialize jk*/ let jk = INIT_JK[prec]; let jp = jk; /* determine jx,jv,q0, note that 3>q0 */ - let jx = nx-1; - let mut jv = (e0-3)/24; + let jx = nx - 1; + let mut jv = (e0 - 3) / 24; if jv < 0 { - jv=0; + jv = 0; } - let mut q0 = e0-24*(jv+1); + let mut q0 = e0 - 24 * (jv + 1); let jv = jv as usize; /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ - let mut j = (jv-jx) as i32; - let m = jx+jk; + let mut j = (jv - jx) as i32; + let m = jx + jk; for i in 0..=m { - f[i] = if j<0 { - 0. - } else { - IPIO2[j as usize] as f64 - }; + f[i] = if j < 0 { 0. } else { IPIO2[j as usize] as f64 }; j += 1 } - + /* compute q[0],q[1],...q[jk] */ for i in 0..=jk { fw = 0f64; for j in 0..=jx { - fw += x[j]*f[jx+i-j]; + fw += x[j] * f[jx + i - j]; } q[i] = fw; } - + let mut jz = jk; - + 'recompute: loop { /* distill q[] into iq[] reversingly */ let mut i = 0i32; z = q[jz]; for j in (1..=jz).rev() { - fw = (x1p_24*z) as i32 as f64; - iq[i as usize] = (z - x1p24*fw) as i32; - z = q[j-1]+fw; + fw = (x1p_24 * z) as i32 as f64; + iq[i as usize] = (z - x1p24 * fw) as i32; + z = q[j - 1] + fw; i += 1; } - + /* compute n */ - z = scalbn(z, q0); /* actual value of z */ - z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ - n = z as i32; + z = scalbn(z, q0); /* actual value of z */ + z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ + n = z as i32; z -= n as f64; ih = 0; - if q0 > 0 { /* need iq[jz-1] to determine n */ - i = iq[jz-1] >> (24-q0); + if q0 > 0 { + /* need iq[jz-1] to determine n */ + i = iq[jz - 1] >> (24 - q0); n += i; - iq[jz-1] -= i << (24-q0); - ih = iq[jz-1] >> (23-q0); + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); } else if q0 == 0 { - ih = iq[jz-1]>>23; + ih = iq[jz - 1] >> 23; } else if z >= 0.5 { ih = 2; } - - if ih > 0 { /* q > 0.5 */ + + if ih > 0 { + /* q > 0.5 */ n += 1; let mut carry = 0i32; - for i in 0..jz { /* compute 1-q */ + for i in 0..jz { + /* compute 1-q */ let j = iq[i]; if carry == 0 { if j != 0 { @@ -347,10 +305,15 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) iq[i] = 0xffffff - j; } } - if q0 > 0 { /* rare case: chance is 1 in 12 */ + if q0 > 0 { + /* rare case: chance is 1 in 12 */ match q0 { - 1 => { iq[jz-1] &= 0x7fffff; }, - 2 => { iq[jz-1] &= 0x3fffff; }, + 1 => { + iq[jz - 1] &= 0x7fffff; + } + 2 => { + iq[jz - 1] &= 0x3fffff; + } _ => {} } } @@ -361,24 +324,26 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) } } } - + /* check if recomputation is needed */ if z == 0. { let mut j = 0; - for i in (jk..=jz-1).rev() { + for i in (jk..=jz - 1).rev() { j |= iq[i]; } - if j == 0 { /* need recomputation */ + if j == 0 { + /* need recomputation */ let mut k = 1; - while iq[jk-k]==0 { - k += 1; /* k = no. of terms needed */ + while iq[jk - k] == 0 { + k += 1; /* k = no. of terms needed */ } - - for i in (jz+1)..=(jz+k) { /* add q[jz+1] to q[jz+k] */ - f[jx+i] = IPIO2[jv+i] as f64; + + for i in (jz + 1)..=(jz + k) { + /* add q[jz+1] to q[jz+k] */ + f[jx + i] = IPIO2[jv + i] as f64; fw = 0f64; for j in 0..=jx { - fw += x[j]*f[jx+i-j]; + fw += x[j] * f[jx + i - j]; } q[i] = fw; } @@ -389,7 +354,7 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) break; } - + /* chop off zero terms */ if z == 0. { jz -= 1; @@ -398,11 +363,12 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) jz -= 1; q0 -= 24; } - } else { /* break z into 24-bit if necessary */ + } else { + /* break z into 24-bit if necessary */ z = scalbn(z, -q0); if z >= x1p24 { - fw = (x1p_24*z) as i32 as f64; - iq[jz] = (z - x1p24*fw) as i32; + fw = (x1p_24 * z) as i32 as f64; + iq[jz] = (z - x1p24 * fw) as i32; jz += 1; q0 += 24; iq[jz] = fw as i32; @@ -410,25 +376,25 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) iq[jz] = z as i32; } } - + /* convert integer "bit" chunk to floating-point value */ fw = scalbn(1., q0); for i in (0..=jz).rev() { - q[i] = fw*(iq[i] as f64); + q[i] = fw * (iq[i] as f64); fw *= x1p_24; } - + /* compute PIo2[0,...,jp]*q[jz,...,0] */ for i in (0..=jz).rev() { fw = 0f64; let mut k = 0; - while (k <= jp) && (k <= jz-i) { - fw += PIO2[k]*q[i+k]; + while (k <= jp) && (k <= jz - i) { + fw += PIO2[k] * q[i + k]; k += 1; } - fq[jz-i] = fw; + fq[jz - i] = fw; } - + /* compress fq[] into y[] */ match prec { 0 => { @@ -437,7 +403,7 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) fw += fq[i]; } y[0] = if ih == 0 { fw } else { -fw }; - }, + } 1 | 2 => { fw = 0f64; for i in (0..=jz).rev() { @@ -446,38 +412,39 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) // TODO: drop excess precision here once double_t is used fw = fw as f64; y[0] = if ih == 0 { fw } else { -fw }; - fw = fq[0]-fw; + fw = fq[0] - fw; for i in 1..=jz { fw += fq[i]; } y[1] = if ih == 0 { fw } else { -fw }; - }, - 3 => { /* painful */ + } + 3 => { + /* painful */ for i in (1..=jz).rev() { - fw = fq[i-1]+fq[i]; - fq[i] += fq[i-1]-fw; - fq[i-1] = fw; + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; } for i in (2..=jz).rev() { - fw = fq[i-1]+fq[i]; - fq[i] += fq[i-1]-fw; - fq[i-1] = fw; + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; } fw = 0f64; for i in (2..=jz).rev() { fw += fq[i]; } - if ih==0 { - y[0] = fq[0]; - y[1] = fq[1]; - y[2] = fw; + if ih == 0 { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; } else { y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; } - }, - _ => { unreachable!() } + } + _ => unreachable!(), } n & 7 } diff --git a/library/compiler-builtins/libm/src/math/rem_pio2f.rs b/library/compiler-builtins/libm/src/math/rem_pio2f.rs index a908ccd9f751..73ec3775d3dc 100644 --- a/library/compiler-builtins/libm/src/math/rem_pio2f.rs +++ b/library/compiler-builtins/libm/src/math/rem_pio2f.rs @@ -2,40 +2,42 @@ use super::rem_pio2_large; use core::f64; -const TOINT : f64 = 1.5 / f64::EPSILON; +const TOINT: f64 = 1.5 / f64::EPSILON; /// 53 bits of 2/pi -const INV_PIO2 : f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */ +const INV_PIO2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */ /// first 25 bits of pi/2 -const PIO2_1 : f64 = 1.57079631090164184570e+00; /* 0x3FF921FB, 0x50000000 */ +const PIO2_1: f64 = 1.57079631090164184570e+00; /* 0x3FF921FB, 0x50000000 */ /// pi/2 - pio2_1 -const PIO2_1T : f64 = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ +const PIO2_1T: f64 = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ /// Return the remainder of x rem pi/2 in *y /// /// use double precision for everything except passing x /// use __rem_pio2_large() for large x #[inline] -pub(crate) fn rem_pio2f(x : f32) -> (i32, f64) { +pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) { let x64 = x as f64; - - let mut tx : [f64; 1] = [0.,]; - let mut ty : [f64; 1] = [0.,]; + + let mut tx: [f64; 1] = [0.]; + let mut ty: [f64; 1] = [0.]; let ix = x.to_bits() & 0x7fffffff; /* 25+53 bit pi is good enough for medium size */ - if ix < 0x4dc90fdb { /* |x| ~< 2^28*(pi/2), medium size */ + if ix < 0x4dc90fdb { + /* |x| ~< 2^28*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - let f_n = x64*INV_PIO2 + TOINT - TOINT; - return (f_n as i32, x64 - f_n*PIO2_1 - f_n*PIO2_1T); + let f_n = x64 * INV_PIO2 + TOINT - TOINT; + return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T); } - if ix>=0x7f800000 { /* x is inf or NaN */ - return (0, x64-x64); + if ix >= 0x7f800000 { + /* x is inf or NaN */ + return (0, x64 - x64); } /* scale x into [2^23, 2^24-1] */ let sign = (x.to_bits() >> 31) != 0; - let e0 = ((ix>>23) - (0x7f+23)) as i32; /* e0 = ilogb(|x|)-23, positive */ - tx[0] = f32::from_bits(ix - (e0<<23) as u32) as f64; + let e0 = ((ix >> 23) - (0x7f + 23)) as i32; /* e0 = ilogb(|x|)-23, positive */ + tx[0] = f32::from_bits(ix - (e0 << 23) as u32) as f64; let n = rem_pio2_large(&tx, &mut ty, e0, 0); if sign { return (-n, -ty[0]);