diff --git a/benchmarks/src/random_integer_generation.cpp b/benchmarks/src/random_integer_generation.cpp
index 9c1615a5bc1..7a65f2b1f2e 100644
--- a/benchmarks/src/random_integer_generation.cpp
+++ b/benchmarks/src/random_integer_generation.cpp
@@ -31,6 +31,23 @@ void BM_lcg(benchmark::State& state) {
 }
 BENCHMARK(BM_lcg);
 
+/// Test discard()
+
+template <class Engine>
+void BM_discard(benchmark::State& state) {
+    Engine gen;
+    const auto n = static_cast<unsigned long long>(state.range(0));
+    for (auto _ : state) {
+        gen.discard(n);
+        benchmark::DoNotOptimize(gen());
+    }
+}
+BENCHMARK(BM_discard<std::mt19937>)->Range(0, 1 << 18);
+BENCHMARK(BM_discard<std::mt19937_64>)->Range(0, 1 << 18);
+BENCHMARK(BM_discard<std::minstd_rand>)->Range(0, 1 << 18);
+
+/// Support machinery for testing _Rng_from_urng and _Rng_from_urng_v2
+
 std::uint32_t GetMax() {
     std::mt19937 gen;
     std::uniform_int_distribution<std::uint32_t> dist(10'000'000, 20'000'000);
diff --git a/stl/inc/random b/stl/inc/random
index fa401b630cd..bc50e0c9f0e 100644
--- a/stl/inc/random
+++ b/stl/inc/random
@@ -617,12 +617,38 @@ public:
 
     void discard(unsigned long long _Nskip) noexcept /* strengthened */ {
         // discard _Nskip elements
-        auto _Temp = _Prev;
-        for (; 0 < _Nskip; --_Nskip) {
-            _Temp = _Next_linear_congruential_value<_Uint, _Ax, _Cx, _Mx>(_Temp);
-        }
+        if constexpr (_Cx == 0 && _Mx == 2147483647) {
+            // for minstd_rand and minstd_rand0 we can improve performance by
+            // performing fast exponentiation and avoiding constant divisions
+            auto _Temp = static_cast<unsigned long long>(_Prev);
+            auto _Mul  = static_cast<unsigned long long>(_Ax);
+
+            for (;;) {
+                if (_Nskip & 1) {
+                    _Temp = _Temp * _Mul;
+                    _Temp = (_Temp >> 31) + (_Temp & _Mx);
+                    _Temp = _Temp < _Mx ? _Temp : _Temp - _Mx;
+                }
 
-        _Prev = _Temp;
+                if (_Nskip >>= 1) {
+                    _Mul = _Mul * _Mul;
+                    _Mul = (_Mul >> 31) + (_Mul & _Mx);
+                    _Mul = _Mul < _Mx ? _Mul : _Mul - _Mx;
+                } else {
+                    break;
+                }
+            }
+
+            _Prev = static_cast<_Uint>(_Temp);
+        } else {
+            auto _Temp = _Prev;
+
+            for (; 0 < _Nskip; --_Nskip) {
+                _Temp = _Next_linear_congruential_value<_Uint, _Ax, _Cx, _Mx>(_Temp);
+            }
+
+            _Prev = _Temp;
+        }
     }
 
     _NODISCARD friend bool operator==(