Merge pull request #20 from carlosga/19-bid128_ldexp

Closes #19
carlosga · Feb 4, 2024 · 81cdf2a · 81cdf2a
2 parents c3347b8 + 90b7ac9
commit 81cdf2a
Show file tree

Hide file tree

Showing 6 changed files with 627 additions and 12 deletions.
diff --git a/src/bid128_ldexp.rs b/src/bid128_ldexp.rs
@@ -0,0 +1,88 @@
+/* ----------------------------------------------------------------------------- */
+/* decimal128 type from Intel decimal math library port to Rust.                 */
+/* decmathlib-rs - Copyright (C) 2023-2024 Carlos Guzmán Álvarez                 */
+/* ----------------------------------------------------------------------------- */
+/* Intel® Decimal Floating-Point Math Library - Copyright (c) 2018, Intel Corp.  */
+/* ----------------------------------------------------------------------------- */
+
+#![allow(unused_assignments)]
+#![allow(non_camel_case_types)]
+#![allow(non_upper_case_globals)]
+#![allow(non_snake_case)]
+#![allow(dead_code)]
+#![allow(unused_mut)]
+
+use crate::bid_internal::{__add_128_128, __set_status_flags, bid_get_BID128, bid_get_BID128_very_fast, unpack_BID128_value};
+use crate::constants::{DECIMAL_MAX_EXPON_128, QUIET_MASK64, SNAN_MASK64};
+use crate::core::StatusFlags;
+use crate::d128::{_IDEC_flags, BID_SINT64, BID_UINT128, BID_UINT32, BID_UINT64};
+
+pub (crate) fn bid128_ldexp(x: &BID_UINT128, n: i32, rnd_mode: u32, pfpsf: &mut _IDEC_flags) -> BID_UINT128 {
+    let mut CX: BID_UINT128 = BID_UINT128::default();
+    let mut CX2: BID_UINT128 = BID_UINT128::default();
+    let mut CBID_X8: BID_UINT128 = BID_UINT128::default();
+    let mut res: BID_UINT128 = BID_UINT128::default();
+    let mut exp64: BID_SINT64;
+    let mut sign_x: BID_UINT64 = 0;
+    let mut exponent_x: i32 = 0;
+
+    // unpack arguments, check for NaN or Infinity
+    if unpack_BID128_value(&mut sign_x, &mut exponent_x, &mut CX, x) == 0 {
+        // x is Inf. or NaN or 0
+        if (x.w[1] & SNAN_MASK64) == SNAN_MASK64 { // y is sNaN
+            __set_status_flags(pfpsf, StatusFlags::BID_INVALID_EXCEPTION);
+        }
+        res.w[1] = CX.w[1] & QUIET_MASK64;
+        res.w[0] = CX.w[0];
+        if CX.w[1] == 0 {
+            exp64 = (exponent_x as BID_SINT64) + (n as BID_SINT64);
+            if exp64 < 0 {
+                exp64 = 0;
+            }
+            if exp64 > DECIMAL_MAX_EXPON_128 as i64 {
+                exp64 = DECIMAL_MAX_EXPON_128 as BID_SINT64;
+            };
+            exponent_x = exp64 as i32;
+            bid_get_BID128_very_fast(&mut res, sign_x, exponent_x, &CX);
+        }
+        return res;
+    }
+
+    exp64      = (exponent_x as BID_SINT64) + (n as BID_SINT64);
+    exponent_x = exp64 as i32;
+
+    if (exponent_x as BID_UINT32) <= DECIMAL_MAX_EXPON_128 as u32 {
+        bid_get_BID128_very_fast(&mut res, sign_x, exponent_x, &CX);
+        return res;
+    }
+    // check for overflow
+    if exp64 > DECIMAL_MAX_EXPON_128 as i64 {
+        if CX.w[1] < 0x314dc6448d93u64 {
+            // try to normalize coefficient
+            loop {
+                CBID_X8.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61);
+                CBID_X8.w[0] =  CX.w[0] << 3;
+                CX2.w[1]     = (CX.w[1] << 1) | (CX.w[0] >> 63);
+                CX2.w[0]     =  CX.w[0] << 1;
+                CX           = __add_128_128(&CX2, &CBID_X8);
+
+                exponent_x -= 1;
+                exp64      -= 1;
+
+                if ! (CX.w[1] < 0x314dc6448d93u64 && exp64 > DECIMAL_MAX_EXPON_128 as i64) {
+                    break;
+                }
+            }
+        }
+        if exp64 <= DECIMAL_MAX_EXPON_128 as i64 {
+            bid_get_BID128_very_fast(&mut res, sign_x, exponent_x, &CX);
+            return res;
+        } else {
+            exponent_x = 0x7fffffff; // overflow
+        }
+    }
+    // exponent < 0
+    // the BID pack routine will round the coefficient
+    res = bid_get_BID128(sign_x, exponent_x, &CX, rnd_mode, pfpsf);
+    return res;
+}
diff --git a/src/bid128_scalbn.rs b/src/bid128_scalbn.rs
@@ -12,13 +12,10 @@
 #![allow(overflowing_literals)]
 
 use crate::bid_internal::{__add_128_128, __set_status_flags, bid_get_BID128, bid_get_BID128_very_fast, unpack_BID128_value};
-use crate::constants::{QUIET_MASK64, SNAN_MASK64};
+use crate::constants::{DECIMAL_MAX_EXPON_128, QUIET_MASK64, SNAN_MASK64};
 use crate::core::StatusFlags;
 use crate::d128::{_IDEC_flags, BID_SINT64, BID_UINT128, BID_UINT32, BID_UINT64};
 
-const DECIMAL_EXPONENT_BIAS_128: i64 = 6176;
-const MAX_DECIMAL_EXPONENT_128: i64 = 12287;
-
 pub (crate) fn bid128_scalbn(x: &BID_UINT128, n: i32, rnd_mode: u32, pfpsf: &mut _IDEC_flags) -> BID_UINT128 {
     let mut CX: BID_UINT128 = BID_UINT128::default();
     let mut CX2: BID_UINT128 = BID_UINT128::default();
@@ -41,8 +38,8 @@ pub (crate) fn bid128_scalbn(x: &BID_UINT128, n: i32, rnd_mode: u32, pfpsf: &mut
             if exp64 < 0 {
                 exp64 = 0;
             }
-            if exp64 > MAX_DECIMAL_EXPONENT_128 {
-                exp64 = MAX_DECIMAL_EXPONENT_128;
+            if exp64 > DECIMAL_MAX_EXPON_128 as i64 {
+                exp64 = DECIMAL_MAX_EXPON_128 as BID_SINT64;
             }
             exponent_x = exp64 as i32;
             bid_get_BID128_very_fast(&mut res, sign_x, exponent_x, &CX);
@@ -53,12 +50,12 @@ pub (crate) fn bid128_scalbn(x: &BID_UINT128, n: i32, rnd_mode: u32, pfpsf: &mut
     exp64      = (exponent_x as BID_SINT64) + (n as BID_SINT64);
     exponent_x = exp64 as i32;
 
-    if (exponent_x as BID_UINT32) <= (MAX_DECIMAL_EXPONENT_128 as BID_UINT32) {
+    if (exponent_x as BID_UINT32) <= (DECIMAL_MAX_EXPON_128 as BID_UINT32) {
       bid_get_BID128_very_fast(&mut res, sign_x, exponent_x, &CX);
       return res;
     }
     // check for overflow
-    if exp64 > MAX_DECIMAL_EXPONENT_128 {
+    if exp64 > DECIMAL_MAX_EXPON_128 as i64 {
       if CX.w[1] < 0x314dc6448d93u64 {
         // try to normalize coefficient
         loop {
@@ -71,12 +68,12 @@ pub (crate) fn bid128_scalbn(x: &BID_UINT128, n: i32, rnd_mode: u32, pfpsf: &mut
           exponent_x -= 1;
           exp64      -= 1;
 
-          if !(CX.w[1] < 0x314dc6448d93u64 && exp64 > MAX_DECIMAL_EXPONENT_128) {
+          if !(CX.w[1] < 0x314dc6448d93u64 && exp64 > DECIMAL_MAX_EXPON_128 as i64) {
             break;
           }
         }
       }
-      if exp64 <= MAX_DECIMAL_EXPONENT_128 {
+      if exp64 <= DECIMAL_MAX_EXPON_128 as i64 {
         bid_get_BID128_very_fast (&mut res, sign_x, exponent_x, &CX);
         return res;
       } else {

diff --git a/src/d128.rs b/src/d128.rs
@@ -15,6 +15,7 @@ use std::str::FromStr;
 use crate::bid128_add::{bid128_add, bid128_sub};
 use crate::bid128_compare::{bid128_quiet_equal, bid128_quiet_greater, bid128_quiet_greater_equal, bid128_quiet_less, bid128_quiet_less_equal, bid128_quiet_not_equal};
 use crate::bid128_div::bid128_div;
+use crate::bid128_ldexp::bid128_ldexp;
 use crate::bid128_mul::bid128_mul;
 
 use crate::bid128_noncomp::*;
@@ -104,7 +105,7 @@ pub struct BID_UINT128 {
 pub type d128 = BID_UINT128;
 
 #[macro_export]
-macro_rules! d128 {
+macro_rules! dec128 {
     ($t:tt) => {{
         $crate::d128::d128::from_str(stringify!($t)).expect("Invalid decimal number literal")
     }}
@@ -257,6 +258,11 @@ impl d128 {
         bid128_scalbln(self, n, rnd_mode.unwrap_or(DEFAULT_ROUNDING_MODE), pfpsf)
     }
 
+    /// multiply a 128-bit decimal floating-point value by an integral power of 2.
+    pub fn ldexp(&self, n: i32, rnd_mode: Option<u32>, pfpsf: &mut _IDEC_flags) -> Self {
+        bid128_ldexp(self, n, rnd_mode.unwrap_or(DEFAULT_ROUNDING_MODE), pfpsf)
+    }
+
     /// Convert 128-bit decimal floating-point value to 64-bit decimal floating-point format (binary encoding)
     pub fn to_decimal64(&self, rnd_mode: Option<u32>, status: &mut _IDEC_flags) -> decimal64 {
         bid128_to_bid64(self, rnd_mode.unwrap_or(DEFAULT_ROUNDING_MODE), status)

diff --git a/src/lib.rs b/src/lib.rs
@@ -29,4 +29,4 @@ mod bid_round;
 mod constants;
 mod convert;
 pub mod core;
-pub mod d128;mod bid128_scalbn;mod bid128_scalbln;
+pub mod d128;mod bid128_scalbn;mod bid128_scalbln;mod bid128_ldexp;