diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3085994c4d..4abf2b69d4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,8 +12,10 @@ You may also find the [Upgrade Guide](https://rust-random.github.io/book/update.
 
 ### Fixes
 - Fix `WeightedIndex` panic when the sum of float weights is infinite; return `Error::Overflow` instead ([#1808])
+- Fix spurious `Error::NonFinite` from `Uniform::new_inclusive` on large finite float ranges such as `0.0..=f64::MAX` ([#1821])
 
 [#1808]: https://github.com/rust-random/rand/pull/1808
+[#1821]: https://github.com/rust-random/rand/pull/1821
 
 ## [0.10.2] — 2026-07-02
 
diff --git a/src/distr/uniform_float.rs b/src/distr/uniform_float.rs
index 3e91c0f3a0..8ac518f269 100644
--- a/src/distr/uniform_float.rs
+++ b/src/distr/uniform_float.rs
@@ -57,6 +57,9 @@ macro_rules! uniform_float_impl {
             /// Construct, reducing `scale` as required to ensure that rounding
             /// can never yield values greater than `high`.
             ///
+            /// Requirements: `low` and `high` must be finite. `scale` may be
+            /// infinite but must not be NaN.
+            ///
             /// Note: though it may be tempting to use a variant of this method
             /// to ensure that samples from `[low, high)` are always strictly
             /// less than `high`, this approach may be very slow where
@@ -66,11 +69,11 @@ macro_rules! uniform_float_impl {
                 let max_rand = <$ty>::splat(1.0 as $f_scalar - $f_scalar::EPSILON);
 
                 loop {
-                    let mask = (scale * max_rand + low).gt_mask(high);
-                    if !mask.any() {
+                    let mask = (scale * max_rand + low).le_mask(high);
+                    if mask.all() {
                         break;
                     }
-                    scale = scale.decrease_masked(mask);
+                    scale = scale.decrease_masked(!mask);
                 }
 
                 debug_assert!(<$ty>::splat(0.0).all_le(scale));
@@ -126,12 +129,14 @@ macro_rules! uniform_float_impl {
                     return Err(Error::EmptyRange);
                 }
 
-                let max_rand = <$ty>::splat(1.0 as $f_scalar - $f_scalar::EPSILON);
-                let scale = (high - low) / max_rand;
-                if !scale.all_finite() {
+                let range = high - low;
+                if !range.all_finite() {
                     return Err(Error::NonFinite);
                 }
 
+                let max_rand = <$ty>::splat(1.0 as $f_scalar - $f_scalar::EPSILON);
+                let scale = range / max_rand;
+
                 Ok(Self::new_bounded(low, high, scale))
             }
 
@@ -238,6 +243,8 @@ mod tests {
                     (-<$f_scalar>::from_bits(7), -0.0),
                     (0.1 * $f_scalar::MAX, $f_scalar::MAX),
                     (-$f_scalar::MAX * 0.2, $f_scalar::MAX * 0.7),
+                    (0.0, $f_scalar::MAX),
+                    (-$f_scalar::MAX, 0.0),
                 ];
                 for &(low_scalar, high_scalar) in v.iter() {
                     for lane in 0..<$ty>::LEN {
@@ -360,6 +367,10 @@ mod tests {
     #[test]
     fn test_float_overflow() {
         assert_eq!(Uniform::try_from(f64::MIN..f64::MAX), Err(Error::NonFinite));
+        assert_eq!(
+            Uniform::try_from(f64::MIN..=f64::MAX),
+            Err(Error::NonFinite)
+        );
     }
 
     #[test]
diff --git a/src/distr/utils.rs b/src/distr/utils.rs
index d3cfbf6a0f..4d0dca11ec 100644
--- a/src/distr/utils.rs
+++ b/src/distr/utils.rs
@@ -219,11 +219,11 @@ pub(crate) trait FloatSIMDUtils {
     fn all_finite(self) -> bool;
 
     type Mask;
-    fn gt_mask(self, other: Self) -> Self::Mask;
+    fn le_mask(self, other: Self) -> Self::Mask;
 
     // Decrease all lanes where the mask is `true` to the next lower value
     // representable by the floating-point type. At least one of the lanes
-    // must be set.
+    // must be set. Inputs may be non-finite but must not be NaN.
     fn decrease_masked(self, mask: Self::Mask) -> Self;
 
     // Convert from int value. Conversion is done while retaining the numerical
@@ -262,12 +262,12 @@ impl IntAsSIMD for u32 {}
 impl IntAsSIMD for u64 {}
 
 pub(crate) trait BoolAsSIMD: Sized {
-    fn any(self) -> bool;
+    fn all(self) -> bool;
 }
 
 impl BoolAsSIMD for bool {
     #[inline(always)]
-    fn any(self) -> bool {
+    fn all(self) -> bool {
         self
     }
 }
@@ -294,8 +294,8 @@ macro_rules! scalar_float_impl {
             }
 
             #[inline(always)]
-            fn gt_mask(self, other: Self) -> Self::Mask {
-                self > other
+            fn le_mask(self, other: Self) -> Self::Mask {
+                self <= other
             }
 
             #[inline(always)]
@@ -357,8 +357,8 @@ macro_rules! simd_impl {
             }
 
             #[inline(always)]
-            fn gt_mask(self, other: Self) -> Self::Mask {
-                self.simd_gt(other)
+            fn le_mask(self, other: Self) -> Self::Mask {
+                self.simd_le(other)
             }
 
             #[inline(always)]
@@ -401,3 +401,41 @@ macro_rules! simd_impl {
 simd_impl!(f32, u32);
 #[cfg(feature = "simd_support")]
 simd_impl!(f64, u64);
+
+#[cfg(test)]
+mod test {
+    use crate::distr::utils::FloatSIMDUtils;
+    #[cfg(feature = "simd_support")]
+    use std::simd::{Mask, Simd};
+
+    #[test]
+    fn decrease_masked() {
+        assert_eq!((-1.0 - f32::EPSILON).decrease_masked(true), -1.0);
+        assert_eq!((1.0 + f64::EPSILON).decrease_masked(true), 1.0);
+
+        #[cfg(feature = "simd_support")]
+        assert_eq!(
+            Simd::<f32, 4>::splat(1.0 + f32::EPSILON).decrease_masked(Mask::splat(true)),
+            Simd::splat(1.0)
+        );
+
+        #[cfg(feature = "simd_support")]
+        assert_eq!(
+            Simd::<f64, 2>::from_array([-1.0, -1.0 - f64::EPSILON])
+                .decrease_masked(Mask::from_array([false, true])),
+            Simd::splat(-1.0)
+        );
+    }
+
+    #[test]
+    fn decrease_masked_infinity() {
+        assert_eq!(f32::INFINITY.decrease_masked(true), f32::MAX);
+        assert_eq!((-f64::INFINITY).decrease_masked(true), -f64::MAX);
+
+        #[cfg(feature = "simd_support")]
+        assert_eq!(
+            Simd::<f32, 2>::splat(-f32::INFINITY).decrease_masked(Mask::splat(true)),
+            Simd::splat(-f32::MAX)
+        );
+    }
+}