diff --git a/load.c b/load.c
index b85a247c18f522..1c1fe1afa15109 100644
--- a/load.c
+++ b/load.c
@@ -447,6 +447,11 @@ get_loaded_features_index(vm_ns_t *vm_ns)
         VALUE previous_realpath_map = rb_hash_dup(realpath_map);
         rb_hash_clear(realpaths);
         rb_hash_clear(realpath_map);
+
+        /* We have to make a copy of features here because the StringValue call
+         * below could call a Ruby method, which could modify $LOADED_FEATURES
+         * and cause it to be corrupt. */
+        features = rb_ary_resurrect(features);
         for (i = 0; i < RARRAY_LEN(features); i++) {
             VALUE entry, as_str;
             as_str = entry = rb_ary_entry(features, i);
@@ -456,6 +461,10 @@ get_loaded_features_index(vm_ns_t *vm_ns)
                 rb_ary_store(features, i, as_str);
             features_index_add(vm_ns, as_str, INT2FIX(i));
         }
+        /* The user modified $LOADED_FEATURES, so we should restore the changes. */
+        if (!rb_ary_shared_with_p(features, CURRENT_NS_LOADED_FEATURES(vm_ns))) {
+            rb_ary_replace(CURRENT_NS_LOADED_FEATURES(vm_ns), features);
+        }
         reset_loaded_features_snapshot(vm_ns);
 
         features = CURRENT_NS_LOADED_FEATURES_SNAPSHOT(vm_ns);
diff --git a/missing/dtoa.c b/missing/dtoa.c
index 8859fcfa44f9fb..cbd6e6ebae23fb 100644
--- a/missing/dtoa.c
+++ b/missing/dtoa.c
@@ -210,6 +210,29 @@
 #include <locale.h>
 #endif
 
+#if defined(HAVE_STDCKDINT_H) || !defined(__has_include)
+#elif __has_include(<stdckdint.h>)
+# define HAVE_STDCKDINT_H 1
+#endif
+#ifdef HAVE_STDCKDINT_H
+# include <stdckdint.h>
+#endif
+
+#if !defined(ckd_add)
+static inline int /* bool */
+ckd_add(int *result, int x, int y)
+{
+    if (x < 0) {
+        if (y < INT_MIN - x) return 1;
+    }
+    else if (x > 0) {
+        if (y > INT_MAX - x) return 1;
+    }
+    *result = x + y;
+    return 0;
+}
+#endif
+
 #ifdef MALLOC
 extern void *MALLOC(size_t);
 #else
@@ -516,6 +539,7 @@ Balloc(int k)
 
     x = 1 << k;
     rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
+    if (!rv) return NULL;
     rv->k = k;
     rv->maxwds = x;
     rv->sign = rv->wds = 0;
@@ -572,6 +596,10 @@ multadd(Bigint *b, int m, int a)   /* multiply by m and add a */
     if (carry) {
         if (wds >= b->maxwds) {
             b1 = Balloc(b->k+1);
+            if (!b1) {
+                Bfree(b);
+                return NULL;
+            }
             Bcopy(b1, b);
             Bfree(b);
             b = b1;
@@ -593,10 +621,12 @@ s2b(const char *s, int nd0, int nd, ULong y9)
     for (k = 0, y = 1; x > y; y <<= 1, k++) ;
 #ifdef Pack_32
     b = Balloc(k);
+    if (!b) return NULL;
     b->x[0] = y9;
     b->wds = 1;
 #else
     b = Balloc(k+1);
+    if (!b) return NULL;
     b->x[0] = y9 & 0xffff;
     b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
 #endif
@@ -606,13 +636,16 @@ s2b(const char *s, int nd0, int nd, ULong y9)
         s += 9;
         do {
             b = multadd(b, 10, *s++ - '0');
+            if (!b) return NULL;
         } while (++i < nd0);
         s++;
     }
     else
         s += 10;
-    for (; i < nd; i++)
+    for (; i < nd; i++) {
         b = multadd(b, 10, *s++ - '0');
+        if (!b) return NULL;
+    }
     return b;
 }
 
@@ -694,11 +727,14 @@ i2b(int i)
     Bigint *b;
 
     b = Balloc(1);
+    if (!b) return NULL;
     b->x[0] = i;
     b->wds = 1;
     return b;
 }
 
+#define Bzero_p(b) (!(b)->x[0] && (b)->wds <= 1)
+
 static Bigint *
 mult(Bigint *a, Bigint *b)
 {
@@ -715,6 +751,14 @@ mult(Bigint *a, Bigint *b)
 #endif
 #endif
 
+    if (Bzero_p(a) || Bzero_p(b)) {
+        c = Balloc(0);
+        if (!c) return NULL;
+        c->wds = 1;
+        c->x[0] = 0;
+        return c;
+    }
+
     if (a->wds < b->wds) {
         c = a;
         a = b;
@@ -727,6 +771,7 @@ mult(Bigint *a, Bigint *b)
     if (wc > a->maxwds)
         k++;
     c = Balloc(k);
+    if (!c) return NULL;
     for (x = c->x, xa = x + wc; x < xa; x++)
         *x = 0;
     xa = a->x;
@@ -806,50 +851,47 @@ static Bigint *
 pow5mult(Bigint *b, int k)
 {
     Bigint *b1, *p5, *p51;
-    Bigint *p5tmp;
     int i;
     static const int p05[3] = { 5, 25, 125 };
 
-    if ((i = k & 3) != 0)
+    if ((i = k & 3) != 0) {
         b = multadd(b, p05[i-1], 0);
+        if (!b) return NULL;
+    }
+
+#define b_cache(var, addr, new_expr) \
+    if ((var = addr) != 0) {} else { \
+        Bigint *tmp = 0; \
+        ACQUIRE_DTOA_LOCK(1); \
+        if (!(var = addr) && (var = (new_expr)) != 0) { \
+            var->next = 0; \
+            tmp = ATOMIC_PTR_CAS(addr, NULL, var); \
+        } \
+        FREE_DTOA_LOCK(1); \
+        if (UNLIKELY(tmp)) { \
+            Bfree(var); \
+            var = tmp; \
+        } \
+        else if (!var) { \
+            Bfree(b); \
+            return NULL; \
+        } \
+    }
 
     if (!(k >>= 2))
         return b;
-    if (!(p5 = p5s)) {
-        /* first time */
-        ACQUIRE_DTOA_LOCK(1);
-        if (!(p5 = p5s)) {
-            p5 = i2b(625);
-            p5->next = 0;
-            p5tmp = ATOMIC_PTR_CAS(p5s, NULL, p5);
-            if (UNLIKELY(p5tmp)) {
-                Bfree(p5);
-                p5 = p5tmp;
-            }
-        }
-        FREE_DTOA_LOCK(1);
-    }
+    /* first time */
+    b_cache(p5, p5s, i2b(625));
     for (;;) {
         if (k & 1) {
             b1 = mult(b, p5);
             Bfree(b);
             b = b1;
+            if (!b) return NULL;
         }
         if (!(k >>= 1))
             break;
-        if (!(p51 = p5->next)) {
-            ACQUIRE_DTOA_LOCK(1);
-            if (!(p51 = p5->next)) {
-                p51 = mult(p5,p5);
-                p51->next = 0;
-                p5tmp = ATOMIC_PTR_CAS(p5->next, NULL, p51);
-                if (UNLIKELY(p5tmp)) {
-                    Bfree(p51);
-                    p51 = p5tmp;
-                }
-            }
-            FREE_DTOA_LOCK(1);
-        }
+        b_cache(p51, p5->next, mult(p5, p5));
         p5 = p51;
     }
     return b;
@@ -862,6 +904,8 @@ lshift(Bigint *b, int k)
     Bigint *b1;
     ULong *x, *x1, *xe, z;
 
+    if (!k || Bzero_p(b)) return b;
+
 #ifdef Pack_32
     n = k >> 5;
 #else
@@ -872,6 +916,10 @@ lshift(Bigint *b, int k)
     for (i = b->maxwds; n1 > i; i <<= 1)
         k1++;
     b1 = Balloc(k1);
+    if (!b1) {
+        Bfree(b);
+        return NULL;
+    }
     x1 = b1->x;
     for (i = 0; i < n; i++)
         *x1++ = 0;
@@ -957,6 +1005,7 @@ diff(Bigint *a, Bigint *b)
     i = cmp(a,b);
     if (!i) {
         c = Balloc(0);
+        if (!c) return NULL;
         c->wds = 1;
         c->x[0] = 0;
         return c;
@@ -970,6 +1019,7 @@ diff(Bigint *a, Bigint *b)
     else
         i = 0;
     c = Balloc(a->k);
+    if (!c) return NULL;
     c->sign = i;
     wa = a->wds;
     xa = a->x;
@@ -1155,6 +1205,7 @@ d2b(double d_, int *e, int *bits)
 #else
     b = Balloc(2);
 #endif
+    if (!b) return NULL;
     x = b->x;
 
     z = d0 & Frac_mask;
@@ -1905,12 +1956,16 @@ strtod(const char *s00, char **se)
     /* Put digits into bd: true value = bd * 10^e */
 
     bd0 = s2b(s0, nd0, nd, y);
+    if (!bd0) goto ret;
 
     for (;;) {
         bd = Balloc(bd0->k);
+        if (!bd) goto retfree;
         Bcopy(bd, bd0);
         bb = d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
+        if (!bb) goto retfree;
         bs = i2b(1);
+        if (!bs) goto retfree;
 
         if (e >= 0) {
             bb2 = bb5 = 0;
@@ -1967,19 +2022,30 @@ strtod(const char *s00, char **se)
         }
         if (bb5 > 0) {
             bs = pow5mult(bs, bb5);
+            if (!bs) goto retfree;
             bb1 = mult(bs, bb);
             Bfree(bb);
             bb = bb1;
+            if (!bb) goto retfree;
         }
-        if (bb2 > 0)
+        if (bb2 > 0) {
             bb = lshift(bb, bb2);
-        if (bd5 > 0)
+            if (!bb) goto retfree;
+        }
+        if (bd5 > 0) {
             bd = pow5mult(bd, bd5);
-        if (bd2 > 0)
+            if (!bd) goto retfree;
+        }
+        if (bd2 > 0) {
             bd = lshift(bd, bd2);
-        if (bs2 > 0)
+            if (!bd) goto retfree;
+        }
+        if (bs2 > 0) {
             bs = lshift(bs, bs2);
+            if (!bs) goto retfree;
+        }
         delta = diff(bb, bd);
+        if (!delta) goto retfree;
         dsign = delta->sign;
         delta->sign = 0;
         i = cmp(delta, bs);
@@ -2012,6 +2078,7 @@ strtod(const char *s00, char **se)
 #endif
                         {
                             delta = lshift(delta,Log2P);
+                            if (!delta) goto nomem;
                             if (cmp(delta, bs) <= 0)
                                 adj = -0.5;
                         }
@@ -2101,6 +2168,7 @@ strtod(const char *s00, char **se)
                 break;
             }
             delta = lshift(delta,Log2P);
+            if (!delta) goto retfree;
             if (cmp(delta, bs) > 0)
                 goto drop_down;
             break;
@@ -2503,6 +2571,7 @@ nrv_alloc(const char *s, char **rve, size_t n)
     char *rv, *t;
 
     t = rv = rv_alloc(n);
+    if (!rv) return NULL;
     while ((*t = *s++) != 0) t++;
     if (rve)
         *rve = t;
@@ -2675,6 +2744,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
 #endif
 
     b = d2b(dval(d), &be, &bbits);
+    if (!b) return NULL;
 #ifdef Sudden_Underflow
     i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 #else
@@ -2794,13 +2864,20 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         leftright = 0;
         /* no break */
       case 5:
-        i = ndigits + k + 1;
+        if (ckd_add(&i, ndigits, k + 1)) { /* k + 1 should be safe */
+            Bfree(b);
+            return NULL;
+        }
         ilim = i;
         ilim1 = i - 1;
         if (i <= 0)
             i = 1;
     }
     s = s0 = rv_alloc(i+1);
+    if (!s) {
+        Bfree(b);
+        return NULL;
+    }
 
 #ifdef Honor_FLT_ROUNDS
     if (mode > 1 && rounding != 1)
@@ -2981,6 +3058,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         b2 += i;
         s2 += i;
         mhi = i2b(1);
+        if (!mhi) goto nomem;
     }
     if (m2 > 0 && s2 > 0) {
         i = m2 < s2 ? m2 : s2;
@@ -2992,19 +3070,28 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         if (leftright) {
             if (m5 > 0) {
                 mhi = pow5mult(mhi, m5);
+                if (!mhi) goto nomem;
                 b1 = mult(mhi, b);
                 Bfree(b);
                 b = b1;
+                if (!b) goto nomem;
             }
-            if ((j = b5 - m5) != 0)
+            if ((j = b5 - m5) != 0) {
                 b = pow5mult(b, j);
+                if (!b) goto nomem;
+            }
         }
-        else
+        else {
             b = pow5mult(b, b5);
+            if (!b) goto nomem;
+        }
     }
     S = i2b(1);
-    if (s5 > 0)
+    if (!S) goto nomem;
+    if (s5 > 0) {
         S = pow5mult(S, s5);
+        if (!S) goto nomem;
+    }
 
     /* Check for special case that d is a normalized power of 2. */
 
@@ -3052,16 +3139,23 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         m2 += i;
         s2 += i;
     }
-    if (b2 > 0)
+    if (b2 > 0) {
         b = lshift(b, b2);
-    if (s2 > 0)
+        if (!b) goto nomem;
+    }
+    if (s2 > 0) {
         S = lshift(S, s2);
+        if (!S) goto nomem;
+    }
     if (k_check) {
         if (cmp(b,S) < 0) {
             k--;
             b = multadd(b, 10, 0);  /* we botched the k estimate */
-            if (leftright)
+            if (!b) goto nomem;
+            if (leftright) {
                 mhi = multadd(mhi, 10, 0);
+                if (!mhi) goto nomem;
+            }
             ilim = ilim1;
         }
     }
@@ -3078,8 +3172,10 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         goto ret;
     }
     if (leftright) {
-        if (m2 > 0)
+        if (m2 > 0) {
             mhi = lshift(mhi, m2);
+            if (!mhi) goto nomem;
+        }
 
         /* Compute mlo -- check for special case
          * that d is a normalized power of 2.
@@ -3088,8 +3184,10 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
         mlo = mhi;
         if (spec_case) {
             mhi = Balloc(mhi->k);
+            if (!mhi) goto nomem;
             Bcopy(mhi, mlo);
             mhi = lshift(mhi, Log2P);
+            if (!mhi) goto nomem;
         }
 
         for (i = 1;;i++) {
@@ -3099,6 +3197,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
              */
             j = cmp(b, mlo);
             delta = diff(S, mhi);
+            if (!delta) goto nomem;
             j1 = delta->sign ? 1 : cmp(b, delta);
             Bfree(delta);
 #ifndef ROUND_BIASED
@@ -3139,6 +3238,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
 #endif /*Honor_FLT_ROUNDS*/
                 if (j1 > 0) {
                     b = lshift(b, 1);
+                    if (!b) goto nomem;
                     j1 = cmp(b, S);
                     if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
                         goto round_9_up;
@@ -3167,11 +3267,16 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
             if (i == ilim)
                 break;
             b = multadd(b, 10, 0);
-            if (mlo == mhi)
+            if (!b) goto nomem;
+            if (mlo == mhi) {
                 mlo = mhi = multadd(mhi, 10, 0);
+                if (!mlo) goto nomem;
+            }
             else {
                 mlo = multadd(mlo, 10, 0);
+                if (!mlo) goto nomem;
                 mhi = multadd(mhi, 10, 0);
+                if (!mhi) goto nomem;
             }
         }
     }
@@ -3187,6 +3292,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
             if (i >= ilim)
                 break;
             b = multadd(b, 10, 0);
+            if (!b) goto nomem;
         }
 
     /* Round off last digit */
@@ -3198,6 +3304,7 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
     }
 #endif
     b = lshift(b, 1);
+    if (!b) goto nomem;
     j = cmp(b, S);
     if (j > 0 || (j == 0 && (dig & 1))) {
  roundoff:
@@ -3239,6 +3346,16 @@ dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
     if (rve)
         *rve = s;
     return s0;
+  nomem:
+    if (S) Bfree(S);
+    if (mhi) {
+        if (mlo && mlo != mhi)
+            Bfree(mlo);
+        Bfree(mhi);
+    }
+    if (b) Bfree(b);
+    FREE(s0);
+    return NULL;
 }
 
 /*-
@@ -3347,6 +3464,7 @@ hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, char **rv
 	 */
 	bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
 	s0 = rv_alloc(bufsize+1);
+        if (!s0) return NULL;
 
 	/* Round to the desired number of digits. */
 	if (SIGFIGS > ndigits && ndigits > 0) {
diff --git a/test/ruby/test_require.rb b/test/ruby/test_require.rb
index 13e707639123e9..3b6cc1178df221 100644
--- a/test/ruby/test_require.rb
+++ b/test/ruby/test_require.rb
@@ -840,6 +840,30 @@ def test_require_with_loaded_features_pop
         p :ok
       end;
     }
+
+    # [Bug #21567]
+    assert_separately(%w[-rtempfile], "#{<<~"begin;"}\n#{<<~"end;"}")
+    begin;
+      class MyString
+        def initialize(path)
+          @path = path
+        end
+
+        def to_str
+          $LOADED_FEATURES.clear
+          @path
+        end
+
+        def to_path = @path
+      end
+
+      def create_ruby_file = Tempfile.create(["test", ".rb"]).path
+
+      require MyString.new(create_ruby_file)
+      $LOADED_FEATURES.unshift(create_ruby_file)
+      $LOADED_FEATURES << MyString.new(create_ruby_file)
+      require create_ruby_file
+    end;
   end
 
   def test_loading_fifo_threading_raise