src/libcompiler_builtins/src/float/mul.rs

   1 use int::{CastInto, Int, WideInt};
   2 use float::Float;
   3
   4 fn mul<F: Float>(a: F, b: F) -> F
   5 where
   6     u32: CastInto<F::Int>,
   7     F::Int: CastInto<u32>,
   8     i32: CastInto<F::Int>,
   9     F::Int: CastInto<i32>,
  10     F::Int: WideInt,
  11 {
  12     let one = F::Int::ONE;
  13     let zero = F::Int::ZERO;
  14
  15     let bits = F::BITS;
  16     let significand_bits = F::SIGNIFICAND_BITS;
  17     let max_exponent = F::EXPONENT_MAX;
  18
  19     let exponent_bias = F::EXPONENT_BIAS;
  20
  21     let implicit_bit = F::IMPLICIT_BIT;
  22     let significand_mask = F::SIGNIFICAND_MASK;
  23     let sign_bit = F::SIGN_MASK as F::Int;
  24     let abs_mask = sign_bit - one;
  25     let exponent_mask = F::EXPONENT_MASK;
  26     let inf_rep = exponent_mask;
  27     let quiet_bit = implicit_bit >> 1;
  28     let qnan_rep = exponent_mask | quiet_bit;
  29     let exponent_bits = F::EXPONENT_BITS;
  30
  31     let a_rep = a.repr();
  32     let b_rep = b.repr();
  33
  34     let a_exponent = (a_rep >> significand_bits) & max_exponent.cast();
  35     let b_exponent = (b_rep >> significand_bits) & max_exponent.cast();
  36     let product_sign = (a_rep ^ b_rep) & sign_bit;
  37
  38     let mut a_significand = a_rep & significand_mask;
  39     let mut b_significand = b_rep & significand_mask;
  40     let mut scale = 0;
  41
  42     // Detect if a or b is zero, denormal, infinity, or NaN.
  43     if a_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
  44         || b_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
  45     {
  46         let a_abs = a_rep & abs_mask;
  47         let b_abs = b_rep & abs_mask;
  48
  49         // NaN + anything = qNaN
  50         if a_abs > inf_rep {
  51             return F::from_repr(a_rep | quiet_bit);
  52         }
  53         // anything + NaN = qNaN
  54         if b_abs > inf_rep {
  55             return F::from_repr(b_rep | quiet_bit);
  56         }
  57
  58         if a_abs == inf_rep {
  59             if b_abs != zero {
  60                 // infinity * non-zero = +/- infinity
  61                 return F::from_repr(a_abs | product_sign);
  62             } else {
  63                 // infinity * zero = NaN
  64                 return F::from_repr(qnan_rep);
  65             }
  66         }
  67
  68         if b_abs == inf_rep {
  69             if a_abs != zero {
  70                 // infinity * non-zero = +/- infinity
  71                 return F::from_repr(b_abs | product_sign);
  72             } else {
  73                 // infinity * zero = NaN
  74                 return F::from_repr(qnan_rep);
  75             }
  76         }
  77
  78         // zero * anything = +/- zero
  79         if a_abs == zero {
  80             return F::from_repr(product_sign);
  81         }
  82
  83         // anything * zero = +/- zero
  84         if b_abs == zero {
  85             return F::from_repr(product_sign);
  86         }
  87
  88         // one or both of a or b is denormal, the other (if applicable) is a
  89         // normal number.  Renormalize one or both of a and b, and set scale to
  90         // include the necessary exponent adjustment.
  91         if a_abs < implicit_bit {
  92             let (exponent, significand) = F::normalize(a_significand);
  93             scale += exponent;
  94             a_significand = significand;
  95         }
  96
  97         if b_abs < implicit_bit {
  98             let (exponent, significand) = F::normalize(b_significand);
  99             scale += exponent;
 100             b_significand = significand;
 101         }
 102     }
 103
 104     // Or in the implicit significand bit.  (If we fell through from the
 105     // denormal path it was already set by normalize( ), but setting it twice
 106     // won't hurt anything.)
 107     a_significand |= implicit_bit;
 108     b_significand |= implicit_bit;
 109
 110     // Get the significand of a*b.  Before multiplying the significands, shift
 111     // one of them left to left-align it in the field.  Thus, the product will
 112     // have (exponentBits + 2) integral digits, all but two of which must be
 113     // zero.  Normalizing this result is just a conditional left-shift by one
 114     // and bumping the exponent accordingly.
 115     let (mut product_high, mut product_low) =
 116         <F::Int as WideInt>::wide_mul(a_significand, b_significand << exponent_bits);
 117
 118     let a_exponent_i32: i32 = a_exponent.cast();
 119     let b_exponent_i32: i32 = b_exponent.cast();
 120     let mut product_exponent: i32 = a_exponent_i32
 121         .wrapping_add(b_exponent_i32)
 122         .wrapping_add(scale)
 123         .wrapping_sub(exponent_bias as i32);
 124
 125     // Normalize the significand, adjust exponent if needed.
 126     if (product_high & implicit_bit) != zero {
 127         product_exponent = product_exponent.wrapping_add(1);
 128     } else {
 129         <F::Int as WideInt>::wide_shift_left(&mut product_high, &mut product_low, 1);
 130     }
 131
 132     // If we have overflowed the type, return +/- infinity.
 133     if product_exponent >= max_exponent as i32 {
 134         return F::from_repr(inf_rep | product_sign);
 135     }
 136
 137     if product_exponent <= 0 {
 138         // Result is denormal before rounding
 139         //
 140         // If the result is so small that it just underflows to zero, return
 141         // a zero of the appropriate sign.  Mathematically there is no need to
 142         // handle this case separately, but we make it a special case to
 143         // simplify the shift logic.
 144         let shift = one.wrapping_sub(product_exponent.cast()).cast();
 145         if shift >= bits as i32 {
 146             return F::from_repr(product_sign);
 147         }
 148
 149         // Otherwise, shift the significand of the result so that the round
 150         // bit is the high bit of productLo.
 151         <F::Int as WideInt>::wide_shift_right_with_sticky(
 152             &mut product_high,
 153             &mut product_low,
 154             shift,
 155         )
 156     } else {
 157         // Result is normal before rounding; insert the exponent.
 158         product_high &= significand_mask;
 159         product_high |= product_exponent.cast() << significand_bits;
 160     }
 161
 162     // Insert the sign of the result:
 163     product_high |= product_sign;
 164
 165     // Final rounding.  The final result may overflow to infinity, or underflow
 166     // to zero, but those are the correct results in those cases.  We use the
 167     // default IEEE-754 round-to-nearest, ties-to-even rounding mode.
 168     if product_low > sign_bit {
 169         product_high += one;
 170     }
 171
 172     if product_low == sign_bit {
 173         product_high += product_high & one;
 174     }
 175
 176     return F::from_repr(product_high);
 177 }
 178
 179 intrinsics! {
 180     #[aapcs_on_arm]
 181     #[arm_aeabi_alias = __aeabi_fmul]
 182     pub extern "C" fn __mulsf3(a: f32, b: f32) -> f32 {
 183         mul(a, b)
 184     }
 185
 186     #[aapcs_on_arm]
 187     #[arm_aeabi_alias = __aeabi_dmul]
 188     pub extern "C" fn __muldf3(a: f64, b: f64) -> f64 {
 189         mul(a, b)
 190     }
 191
 192     #[cfg(target_arch = "arm")]
 193     pub extern "C" fn __mulsf3vfp(a: f32, b: f32) -> f32 {
 194         a * b
 195     }
 196
 197     #[cfg(target_arch = "arm")]
 198     pub extern "C" fn __muldf3vfp(a: f64, b: f64) -> f64 {
 199         a * b
 200     }
 201 }