fully generic

Author: tgross35

src/math/generic/fma.rs

@@ -5,11 +5,15 @@ use super::super::{CastFrom, Float, Int, MinInt};
 
 const ZEROINFNAN: i32 = 0x7ff - 0x3ff - 52 - 1;
 
-type F = f64;
-
 /// Fused multiply-add.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fma(x: f64, y: f64, z: f64) -> f64 {
+pub fn fma<F>(x: F, y: F, z: F) -> F
+where
+    F: Float + Helper,
+    F: CastFrom<F::SignedInt>,
+    F: CastFrom<i8>,
+    F::Int: HInt,
+{
     // let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
     let one = IntTy::<F>::ONE;
     let zero = IntTy::<F>::ZERO;
@@ -32,8 +36,8 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
     }
 
     /* mul: r = x*y */
-    let zhi: u64;
-    let zlo: u64;
+    let zhi: F::Int;
+    let zlo: F::Int;
     let (mut rlo, mut rhi) = nx.m.widen_mul(ny.m).lo_hi();
 
     /* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
@@ -55,8 +59,9 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
             d -= sbits;
             if d == 0 {
             } else if d < sbits {
-                rlo =
-                    (rhi << (sbits - d)) | (rlo >> d) | IntTy::<F>::from((rlo << (sbits - d)) != 0);
+                rlo = (rhi << (sbits - d))
+                    | (rlo >> d)
+                    | IntTy::<F>::from((rlo << (sbits - d)) != zero);
                 rhi = rhi >> d;
             } else {
                 rlo = one;
@@ -69,7 +74,7 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
         if d == 0 {
             zlo = nz.m;
         } else if d < sbits {
-            zlo = (nz.m >> d) | IntTy::<F>::from((nz.m << (sbits - d)) != 0);
+            zlo = (nz.m >> d) | IntTy::<F>::from((nz.m << (sbits - d)) != zero);
         } else {
             zlo = one;
         }
@@ -88,25 +93,24 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
         let (res, borrow) = rlo.overflowing_sub(zlo);
         rlo = res;
         rhi = rhi.wrapping_sub(zhi.wrapping_add(IntTy::<F>::from(borrow)));
-        if (rhi >> (F::BITS - 1)) != 0 {
+        if (rhi >> (F::BITS - 1)) != zero {
             rlo = rlo.signed().wrapping_neg().unsigned();
-            rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::<F>::from(rlo != 0);
+            rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::<F>::from(rlo != zero);
             neg = !neg;
-            // sign = (sign == 0) as i32;
         }
-        nonzero = (rhi != 0) as i32;
+        nonzero = (rhi != zero) as i32;
     }
 
     /* set rhi to top 63bit of the result (last bit is sticky) */
     if nonzero != 0 {
         e += sbits;
         d = rhi.leading_zeros() as i32 - 1;
         /* note: d > 0 */
-        rhi = (rhi << d) | (rlo >> (sbits - d)) | IntTy::<F>::from((rlo << d) != 0);
-    } else if rlo != 0 {
+        rhi = (rhi << d) | (rlo >> (sbits - d)) | IntTy::<F>::from((rlo << d) != zero);
+    } else if rlo != zero {
         d = rlo.leading_zeros() as i32 - 1;
         if d < 0 {
-            rhi = (rlo >> 1) | (rlo & 1);
+            rhi = (rlo >> 1) | (rlo & one);
         } else {
             rhi = rlo << d;
         }
@@ -117,17 +121,17 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
     e -= d;
 
     /* convert to double */
-    let mut i: i64 = rhi.signed(); /* i is in [1<<62,(1<<63)-1] */
+    let mut i: F::SignedInt = rhi.signed(); /* i is in [1<<62,(1<<63)-1] */
     if neg {
         i = -i;
     }
 
-    let mut r: f64 = f64::cast_from_lossy(i); /* |r| is in [0x1p62,0x1p63] */
+    let mut r: F = F::cast_from_lossy(i); /* |r| is in [0x1p62,0x1p63] */
 
     if e < -(F::EXP_BIAS as i32 - 1) - (sbits - 2) {
         /* result is subnormal before rounding */
         if e == -(F::EXP_BIAS as i32 - 1) - (sbits - 1) {
-            let mut c: f64 = magic;
+            let mut c: F = magic;
             if neg {
                 c = -c;
             }
@@ -139,13 +143,14 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
             }
             /* one bit is lost when scaled, add another top bit to
              * only round once at conversion if it is inexact */
-            if (rhi << F::SIG_BITS) != 0 {
-                i = ((rhi >> 1) | (rhi & 1) | (1 << 62)).signed();
+            if (rhi << F::SIG_BITS) != zero {
+                let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << 62);
+                i = iu.signed();
                 if neg {
                     i = -i;
                 }
                 r = F::cast_from(i);
-                r = 2.0 * r - c; /* remove top bit */
+                r = F::cast_from(2i8) * r - c; /* remove top bit */
 
                 /* raise underflow portably, such that it
                  * cannot be optimized away */
@@ -154,11 +159,12 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
         } else {
             /* only round once when scaled */
             d = 10;
-            i = (((rhi >> d) | IntTy::<F>::from(rhi << (F::BITS as i32 - d) != 0)) << d).signed();
+            i = (((rhi >> d) | IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero)) << d)
+                .signed();
             if neg {
                 i = -i;
             }
-            r = f64::cast_from(i);
+            r = F::cast_from(i);
         }
     }
 
@@ -197,13 +203,13 @@ impl<F: Float> Norm<F> {
 }
 
 // Need to figure out how to do this better.
-trait RaiseUnderflow {
+pub trait Helper {
     fn raise_underflow(self) -> Self;
     fn raise_underflow2(self) -> Self;
     fn scalbn(self, n: i32) -> Self;
 }
 
-impl RaiseUnderflow for f64 {
+impl Helper for f64 {
     fn raise_underflow(self) -> Self {
         let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
         let fltmin: f32 = (x0_ffffff8p_63 * f32::MIN_POSITIVE as f64 * self) as f32;

src/math/support/float_traits.rs

@@ -23,7 +23,9 @@ pub trait Float:
     type Int: Int<OtherSign = Self::SignedInt, Unsigned = Self::Int>;
 
     /// A int of the same width as the float
-    type SignedInt: Int + MinInt<OtherSign = Self::Int, Unsigned = Self::Int>;
+    type SignedInt: Int
+        + MinInt<OtherSign = Self::Int, Unsigned = Self::Int>
+        + ops::Neg<Output = Self::SignedInt>;
 
     const ZERO: Self;
     const NEG_ZERO: Self;

src/math/support/int_traits.rs

@@ -52,6 +52,9 @@ pub trait Int:
     + ops::Sub<Output = Self>
     + ops::Mul<Output = Self>
     + ops::Div<Output = Self>
+    + ops::Shl<i32, Output = Self>
+    + ops::Shl<u32, Output = Self>
+    + ops::Shr<i32, Output = Self>
     + ops::Shr<u32, Output = Self>
     + ops::BitXor<Output = Self>
     + ops::BitAnd<Output = Self>
@@ -93,6 +96,7 @@ pub trait Int:
     fn wrapping_shr(self, other: u32) -> Self;
     fn rotate_left(self, other: u32) -> Self;
     fn overflowing_add(self, other: Self) -> (Self, bool);
+    fn overflowing_sub(self, other: Self) -> (Self, bool);
     fn leading_zeros(self) -> u32;
     fn ilog2(self) -> u32;
 }
@@ -151,6 +155,10 @@ macro_rules! int_impl_common {
             <Self>::overflowing_add(self, other)
         }
 
+        fn overflowing_sub(self, other: Self) -> (Self, bool) {
+            <Self>::overflowing_sub(self, other)
+        }
+
         fn leading_zeros(self) -> u32 {
             <Self>::leading_zeros(self)
         }
@@ -437,6 +445,7 @@ cast_into!(i64);
 cast_into!(u128);
 cast_into!(i128);
 
+cast_into_float!(i8);
 cast_into_float!(i16);
 cast_into_float!(i32);
 cast_into_float!(i64);