fmod: Normalize subnormals to maximum precision

Discussed at [1], there was an off-by-one mistake when converting from
the loop routine to using `leading_zeros` for normalization.

Currently, using `EXP_BITS` has the effect that `ix` after the branch
has its MSB _one bit to the left_ of the implicit bit's position,
whereas a shift by `EXP_BITS + 1` ensures that the MSB is exactly at the
implicit bit's position, matching what is done for normals (where the
implicit bit is set to be explicit). This doesn't seem to have any
effect in our implementation since the failing test cases from [1]
appear to still have correct results.

Since the result of using `EXP_BITS + 1` is more consistent with what is
done for normals, apply this here.

[1]: #469 (comment)

Author: tgross35

src/math/generic/fmod.rs

@@ -26,7 +26,7 @@ pub fn fmod<F: Float>(x: F, y: F) -> F {
 
     /* normalize x and y */
     if ex == 0 {
-        let i = ix << F::EXP_BITS;
+        let i = ix << (F::EXP_BITS + 1);
         ex -= i.leading_zeros() as i32;
         ix <<= -ex + 1;
     } else {
@@ -35,7 +35,7 @@ pub fn fmod<F: Float>(x: F, y: F) -> F {
     }
 
     if ey == 0 {
-        let i = iy << F::EXP_BITS;
+        let i = iy << (F::EXP_BITS + 1);
         ey -= i.leading_zeros() as i32;
         iy <<= -ey + 1;
     } else {