comments

Author: tgross35

src/math/generic/scalbn.rs

@@ -1,5 +1,9 @@
+use super::super::support::Hexf;
 use super::super::{CastFrom, CastInto, Float, IntTy, MinInt};
 
+extern crate std;
+use std::dbg;
+
 /// Scale the exponent.
 ///
 /// From N3220:
@@ -40,76 +44,80 @@ where
     // 2 ^ sig_total_bits, representation of what can be accounted for with subnormals
     let f_exp_subnorm = F::from_parts(false, sig_total_bits + F::EXP_BIAS, zero);
 
+    dbg!((Hexf(f_exp_max), Hexf(f_exp_min), Hexf(f_exp_subnorm)));
+
+    dbg!(Hexf(x), n);
+
+    // The goal is to multiply `x` by a scale factor that applies `n`. However, there are cases
+    // where `2^n` is not representable by `F` but the result should be, e.g. `x = 2^Emin` with
+    // `n = -EMin + 2`. To get around this, reduce the magnitude of the final scale operation by
+    // prescaling by the max/min power representable by `F`.
+
     if n > exp_max {
+        // Worse case positive `n`: `x`  is the minimum subnormal value, the result is `F::MAX`.
+        // This can be reached by three scaling multiplications (two here and one final).
+        debug_assert!(-exp_min + F::SIG_BITS as i32 + exp_max <= exp_max * 3);
+
         x *= f_exp_max;
         n -= exp_max;
         if n > exp_max {
             x *= f_exp_max;
             n -= exp_max;
-
-            if F::BITS < 32 && n > exp_max {
-                x *= f_exp_max;
-                n -= exp_max;
-
-                if n > exp_max {
-                    x *= f_exp_max;
-                    n -= exp_max;
-                    if n > exp_max {
-                        x *= f_exp_max;
-                        n -= exp_max;
-                        if n > exp_max {
-                            x *= f_exp_max;
-                            n -= exp_max;
-                            if n > exp_max {
-                                n = exp_max;
-                            }
-                        }
-                    }
-                }
-            } else if n > exp_max {
+            if n > exp_max {
                 n = exp_max;
             }
         }
     } else if n < exp_min {
         let mul = f_exp_min * f_exp_subnorm;
         let add = (exp_max - 1) - sig_total_bits as i32;
+        dbg!(Hexf(mul), add);
 
         x *= mul;
         n += add;
+        dbg!(Hexf(x), n);
         if n < exp_min {
             x *= mul;
             n += add;
+            dbg!(Hexf(x), n);
             if F::BITS < 32 {
                 if n < exp_min {
                     x *= mul;
                     n += add;
 
+                    dbg!(Hexf(x), n);
                     if n < exp_min {
                         x *= mul;
                         n += add;
 
+                        dbg!(Hexf(x), n);
                         if n < exp_min {
                             x *= mul;
                             n += add;
 
+                            dbg!(Hexf(x), n);
                             if n < exp_min {
                                 x *= mul;
                                 n += add;
 
+                                dbg!(Hexf(x), n);
                                 if n < exp_min {
                                     x *= mul;
                                     n += add;
 
+                                    dbg!(Hexf(x), n);
                                     if n < exp_min {
                                         x *= mul;
                                         n += add;
 
+                                        dbg!(Hexf(x), n);
                                         if n < exp_min {
                                             x *= mul;
                                             n += add;
 
+                                            dbg!(Hexf(x), n);
                                             if n < exp_min {
                                                 n = exp_min;
+                                                dbg!(Hexf(x), n);
                                             }
                                         }
                                     }
@@ -120,11 +128,19 @@ where
                 }
             } else if n < exp_min {
                 n = exp_min;
+                dbg!(Hexf(x), n);
             }
         }
     }
 
-    x * F::from_parts(false, (F::EXP_BIAS as i32 + n) as u32, zero)
+    dbg!(Hexf(x), n);
+    let scale = F::from_parts(false, (F::EXP_BIAS as i32 + n) as u32, zero);
+    let ret = x * scale;
+    dbg!(Hexf(scale), Hexf(ret));
+    ret
+
+    // let ret = dbg!(x) * dbg!(F::from_parts(false, (F::EXP_BIAS as i32 + n) as u32, zero));
+    // dbg!(ret)
 }
 
 #[cfg(test)]
@@ -191,4 +207,28 @@ mod tests {
     fn spec_test_f128() {
         spec_test::<f128>();
     }
+
+    // #[test]
+    // fn foobar32() {
+    //     let x = hf32!("0x1.fffffep+127");
+    //     let n = -2147483639;
+    //     scalbn(x, n);
+    //     std::eprintln!();
+    //     let x = hf32!("0x1.fffffep-126");
+    //     let n = 2147483639;
+    //     scalbn(x, n);
+    //     panic!();
+    // }
+
+    // #[test]
+    // fn foobar16() {
+    //     let x = hf16!("0x1.ffp+15");
+    //     let n = -2147483639;
+    //     scalbn(x, n);
+    //     std::eprintln!();
+    //     let x = hf16!("0x1.ffp-15");
+    //     let n = 2147483639;
+    //     scalbn(x, n);
+    //     panic!();
+    // }
 }