diff --git a/src/runtime/Math.Mod b/src/runtime/Math.Mod
index b3ca4e6a..0d1d61f5 100644
--- a/src/runtime/Math.Mod
+++ b/src/runtime/Math.Mod
@@ -741,6 +741,32 @@ BEGIN t := ABS(x);
   RETURN arcsinh(x/sqrt(ONE - x * x))
 END arctanh;
 
+PROCEDURE fcmp* (x, y, epsilon: REAL): INTEGER;
+(*  fcmp: this procedure determines whether `x` and `y` are approximately equal
+    to a relative accuracy `epsilon`.
+    References:
+    The implementation is based on the GNU Scientific Library (GSL).
+    https://www.gnu.org/software/gsl/doc/html/math.html#approximate-comparison-of-floating-point-numbers
+*)
+VAR  max, exponent0, delta, difference: REAL;
+BEGIN
+  IF ABS(x) > ABS(y) THEN
+    max := x;
+  ELSE
+    max := y;
+  END;
+  exponent0 := exponent(max);
+  delta := 2.0*epsilon*power(2.0, exponent0);
+  difference := x - y;
+  IF difference > delta THEN
+    RETURN 1;
+  ELSIF difference < -delta THEN
+    RETURN -1;
+  ELSE
+    RETURN 0; (* approximately equal *)
+  END;
+END fcmp;
+
 PROCEDURE ToREAL(h: HUGEINT): REAL;
 BEGIN RETURN SYSTEM.VAL(REAL, h)
 END ToREAL;
diff --git a/src/runtime/MathL.Mod b/src/runtime/MathL.Mod
index bd17b490..b1448dc0 100644
--- a/src/runtime/MathL.Mod
+++ b/src/runtime/MathL.Mod
@@ -665,6 +665,32 @@ BEGIN t := ABS(x);
   RETURN arcsinh(x/sqrt(ONE-x*x))
 END arctanh;
 
+PROCEDURE fcmp* (x, y, epsilon: LONGREAL): INTEGER;
+(*  fcmp: this procedure determines whether `x` and `y` are approximately equal
+    to a relative accuracy `epsilon`.
+    References:
+    The implementation is based on the GNU Scientific Library (GSL).
+    https://www.gnu.org/software/gsl/doc/html/math.html#approximate-comparison-of-floating-point-numbers
+*)
+VAR  max, exponent0, delta, difference: LONGREAL;
+BEGIN
+  IF ABS(x) > ABS(y) THEN
+    max := x;
+  ELSE
+    max := y;
+  END;
+  exponent0 := exponent(max);
+  delta := 2.0D0*epsilon*power(2.0D0, exponent0);
+  difference := x - y;
+  IF difference > delta THEN
+    RETURN 1;
+  ELSIF difference < -delta THEN
+    RETURN -1;
+  ELSE
+    RETURN 0; (* approximately equal *)
+  END;
+END fcmp;
+
 PROCEDURE ToLONGREAL(h: HUGEINT): LONGREAL;
 BEGIN RETURN SYSTEM.VAL(LONGREAL, h)
 END ToLONGREAL;