diff --git a/src/ansi-c/c_typecheck_expr.cpp b/src/ansi-c/c_typecheck_expr.cpp
index 35bc45cb257..578fbd38fda 100644
--- a/src/ansi-c/c_typecheck_expr.cpp
+++ b/src/ansi-c/c_typecheck_expr.cpp
@@ -3237,6 +3237,23 @@ exprt c_typecheck_baset::do_special_functions(
 
     return std::move(infl_expr);
   }
+  else if(identifier == CPROVER_PREFIX "round_to_integralf" ||
+          identifier == CPROVER_PREFIX "round_to_integrald" ||
+          identifier == CPROVER_PREFIX "round_to_integralld" )
+  {
+    if(expr.arguments().size() != 2)
+    {
+      error().source_location = f_op.source_location();
+      error() << identifier << " expects two arguments" << eom;
+      throw 0;
+    }
+
+    auto round_to_integral_expr =
+      floatbv_round_to_integral_exprt{expr.arguments()[0], expr.arguments()[1]};
+    round_to_integral_expr.add_source_location() = source_location;
+
+    return std::move(round_to_integral_expr);
+  }
   else if(
     identifier == CPROVER_PREFIX "abs" || identifier == CPROVER_PREFIX "labs" ||
     identifier == CPROVER_PREFIX "llabs" ||
diff --git a/src/ansi-c/cprover_builtin_headers.h b/src/ansi-c/cprover_builtin_headers.h
index f0ae5e74004..998e81f0445 100644
--- a/src/ansi-c/cprover_builtin_headers.h
+++ b/src/ansi-c/cprover_builtin_headers.h
@@ -93,6 +93,9 @@ int __CPROVER_islessgreaterf(float f, float g);
 int __CPROVER_islessgreaterd(double f, double g);
 int __CPROVER_isunorderedf(float f, float g);
 int __CPROVER_isunorderedd(double f, double g);
+float __CPROVER_round_to_integralf(float, int);
+double __CPROVER_round_to_integrald(double, int);
+long double __CPROVER_round_to_integralld(long double, int);
 
 // absolute value
 int __CPROVER_abs(int x);
diff --git a/src/ansi-c/library/math.c b/src/ansi-c/library/math.c
index 6d4fc796b24..225fee085d8 100644
--- a/src/ansi-c/library/math.c
+++ b/src/ansi-c/library/math.c
@@ -1244,126 +1244,6 @@ float fdimf(float f, float g) { return ((f > g) ? f - g : +0.0f); }
 
 long double fdiml(long double f, long double g) { return ((f > g) ? f - g : +0.0); }
 
-
-
-/* FUNCTION: __sort_of_CPROVER_round_to_integral */
-// TODO : Should be a real __CPROVER function to convert to SMT-LIB
-
-#ifndef __CPROVER_MATH_H_INCLUDED
-#include <math.h>
-#define __CPROVER_MATH_H_INCLUDED
-#endif
-
-#ifndef __CPROVER_FENV_H_INCLUDED
-#include <fenv.h>
-#define __CPROVER_FENV_H_INCLUDED
-#endif
-
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d)
-{
-  double magicConst = 0x1.0p+52;
-  double return_value;
-  int saved_rounding_mode = fegetround();
-  fesetround(rounding_mode);
-  
-  if (fabs(d) >= magicConst || d == 0.0)
-  {
-    return_value = d;
-  }
-  else
-  {
-    if (!signbit(d)) {
-      double tmp = d + magicConst;
-      return_value = tmp - magicConst;
-    } else {
-      double tmp = d - magicConst;
-      return_value = tmp + magicConst;    
-    }
-  }
-
-  fesetround(saved_rounding_mode);
-  return return_value;
-}
-
-/* FUNCTION: __sort_of_CPROVER_round_to_integralf */
-// TODO : Should be a real __CPROVER function to convert to SMT-LIB
-
-#ifndef __CPROVER_MATH_H_INCLUDED
-#include <math.h>
-#define __CPROVER_MATH_H_INCLUDED
-#endif
-
-#ifndef __CPROVER_FENV_H_INCLUDED
-#include <fenv.h>
-#define __CPROVER_FENV_H_INCLUDED
-#endif
-
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d)
-{
-  float magicConst = 0x1.0p+23f;  // 23 is significant
-  float return_value;
-  int saved_rounding_mode = fegetround();
-  fesetround(rounding_mode);
-  
-  if (fabsf(d) >= magicConst || d == 0.0)
-  {
-    return_value = d;
-  }
-  else
-  {
-    if (!signbit(d)) {
-      float tmp = d + magicConst;
-      return_value = tmp - magicConst;    
-    } else {
-      float tmp = d - magicConst;
-      return_value = tmp + magicConst;    
-    }
-  }
-
-  fesetround(saved_rounding_mode);
-  return return_value;
-}
-
-
-/* FUNCTION: __sort_of_CPROVER_round_to_integrall */
-// TODO : Should be a real __CPROVER function to convert to SMT-LIB
-
-#ifndef __CPROVER_MATH_H_INCLUDED
-#include <math.h>
-#define __CPROVER_MATH_H_INCLUDED
-#endif
-
-#ifndef __CPROVER_FENV_H_INCLUDED
-#include <fenv.h>
-#define __CPROVER_FENV_H_INCLUDED
-#endif
-
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d)
-{
-  long double magicConst = 0x1.0p+64;
-  long double return_value;
-  int saved_rounding_mode = fegetround();
-  fesetround(rounding_mode);
-  
-  if (fabsl(d) >= magicConst || d == 0.0)
-  {
-    return_value = d;
-  }
-  else
-  {
-    if (!signbit(d)) {
-      long double tmp = d + magicConst;
-      return_value = tmp - magicConst;    
-    } else {
-      long double tmp = d - magicConst;
-      return_value = tmp + magicConst;    
-    }
-  }
-
-  fesetround(saved_rounding_mode);
-  return return_value;
-}
-
 /* ISO 9899:2011
  *
  * The ceil functions compute the smallest integer value not less than
@@ -1382,11 +1262,9 @@ long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double ceil(double x)
 {
-  return __sort_of_CPROVER_round_to_integral(FE_UPWARD, x);
+  return __CPROVER_round_to_integral(x, FE_UPWARD);
 }
 
 /* FUNCTION: ceilf */
@@ -1401,11 +1279,9 @@ double ceil(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float ceilf(float x)
 {
-  return __sort_of_CPROVER_round_to_integralf(FE_UPWARD, x);
+  return __CPROVER_round_to_integral(x, FE_UPWARD);
 }
 
 
@@ -1421,11 +1297,9 @@ float ceilf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double ceill(long double x)
 {
-  return __sort_of_CPROVER_round_to_integrall(FE_UPWARD, x);
+  return __CPROVER_round_to_integral(x, FE_UPWARD);
 }
 
 
@@ -1446,11 +1320,9 @@ long double ceill(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double floor(double x)
 {
-  return __sort_of_CPROVER_round_to_integral(FE_DOWNWARD, x);
+  return __CPROVER_round_to_integral(x, FE_DOWNWARD);
 }
 
 /* FUNCTION: floorf */
@@ -1465,11 +1337,9 @@ double floor(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float floorf(float x)
 {
-  return __sort_of_CPROVER_round_to_integralf(FE_DOWNWARD, x);
+  return __CPROVER_round_to_integralf(x, FE_DOWNWARD);
 }
 
 
@@ -1485,11 +1355,9 @@ float floorf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double floorl(long double x)
 {
-  return __sort_of_CPROVER_round_to_integrall(FE_DOWNWARD, x);
+  return __CPROVER_round_to_integral(x, FE_DOWNWARD);
 }
 
 
@@ -1511,11 +1379,9 @@ long double floorl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double trunc(double x)
 {
-  return __sort_of_CPROVER_round_to_integral(FE_TOWARDZERO, x);
+  return __CPROVER_round_to_integral(x, FE_TOWARDZERO);
 }
 
 /* FUNCTION: truncf */
@@ -1530,11 +1396,9 @@ double trunc(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float truncf(float x)
 {
-  return __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, x);
+  return __CPROVER_round_to_integral(x, FE_TOWARDZERO);
 }
 
 
@@ -1550,11 +1414,9 @@ float truncf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double truncl(long double x)
 {
-  return __sort_of_CPROVER_round_to_integrall(FE_TOWARDZERO, x);
+  return __CPROVER_round_to_integral(x, FE_TOWARDZERO);
 }
 
 
@@ -1576,12 +1438,10 @@ long double truncl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double round(double x)
 {
   // Tempting but RNE not RNA
-  // return __sort_of_CPROVER_round_to_integral(FE_TONEAREST, x);
+  // return __CPROVER_round_to_integral(x, FE_TONEAREST);
 
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
@@ -1596,8 +1456,8 @@ double round(double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  return __sort_of_CPROVER_round_to_integral(FE_TOWARDZERO, xp);
+
+  return __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
 }
 
 /* FUNCTION: roundf */
@@ -1612,12 +1472,10 @@ double round(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float roundf(float x)
 {
   // Tempting but RNE not RNA
-  // return __sort_of_CPROVER_round_to_integralf(FE_TONEAREST, x);
+  // return __CPROVER_round_to_integralf(x, FE_TONEAREST);
 
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
@@ -1632,8 +1490,8 @@ float roundf(float x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  return __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, xp);
+
+  return __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
 }
 
 
@@ -1649,12 +1507,10 @@ float roundf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double roundl(long double x)
 {
   // Tempting but RNE not RNA
-  // return __sort_of_CPROVER_round_to_integrall(FE_TONEAREST, x);
+  // return __CPROVER_round_to_integrall(x, FE_TONEAREST);
 
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
@@ -1669,8 +1525,8 @@ long double roundl(long double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  return __sort_of_CPROVER_round_to_integrall(FE_TOWARDZERO, xp);
+
+  return __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
 }
 
 
@@ -1694,11 +1550,9 @@ long double roundl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double nearbyint(double x)
 {
-  return __sort_of_CPROVER_round_to_integral(fegetround(), x);
+  return __CPROVER_round_to_integral(x, fegetround());
 }
 
 /* FUNCTION: nearbyintf */
@@ -1713,11 +1567,9 @@ double nearbyint(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float nearbyintf(float x)
 {
-  return __sort_of_CPROVER_round_to_integralf(fegetround(), x);
+  return __CPROVER_round_to_integral(x, fegetround());
 }
 
 
@@ -1733,11 +1585,9 @@ float nearbyintf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double nearbyintl(long double x)
 {
-  return __sort_of_CPROVER_round_to_integrall(fegetround(), x);
+  return __CPROVER_round_to_integrall(x, fegetround());
 }
 
 
@@ -1761,11 +1611,9 @@ long double nearbyintl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double rint(double x)
 {
-  return __sort_of_CPROVER_round_to_integral(fegetround(), x);
+  return __CPROVER_round_to_integral(x, fegetround());
 }
 
 /* FUNCTION: rintf */
@@ -1780,11 +1628,9 @@ double rint(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float rintf(float x)
 {
-  return __sort_of_CPROVER_round_to_integralf(fegetround(), x);
+  return __CPROVER_round_to_integral(x, fegetround());
 }
 
 /* FUNCTION: rintl */
@@ -1799,11 +1645,9 @@ float rintf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double rintl(long double x)
 {
-  return __sort_of_CPROVER_round_to_integrall(fegetround(), x);
+  return __CPROVER_round_to_integral(x, fegetround());
 }
 
 
@@ -1829,13 +1673,9 @@ long double rintl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 long int lrint(double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  double rti = __sort_of_CPROVER_round_to_integral(fegetround(), x);
+  double rti = __CPROVER_round_to_integral(x, fegetround());
   return (long int)rti;
 }
 
@@ -1851,13 +1691,9 @@ long int lrint(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 long int lrintf(float x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  float rti = __sort_of_CPROVER_round_to_integralf(fegetround(), x);
+  float rti = __CPROVER_round_to_integral(x, fegetround());
   return (long int)rti;
 }
 
@@ -1874,13 +1710,9 @@ long int lrintf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long int lrintl(long double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  long double rti = __sort_of_CPROVER_round_to_integrall(fegetround(), x);
+  long double rti = __CPROVER_round_to_integral(x, fegetround());
   return (long int)rti;
 }
 
@@ -1897,13 +1729,9 @@ long int lrintl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 long long int llrint(double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  double rti = __sort_of_CPROVER_round_to_integral(fegetround(), x);
+  double rti = __CPROVER_round_to_integral(x, fegetround());
   return (long long int)rti;
 }
 
@@ -1919,13 +1747,9 @@ long long int llrint(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 long long int llrintf(float x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  float rti = __sort_of_CPROVER_round_to_integralf(fegetround(), x);
+  float rti = __CPROVER_round_to_integralf(x, fegetround());
   return (long long int)rti;
 }
 
@@ -1942,13 +1766,9 @@ long long int llrintf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long long int llrintl(long double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT
-  long double rti = __sort_of_CPROVER_round_to_integrall(fegetround(), x);
+  long double rti = __CPROVER_round_to_integrall(x, fegetround());
   return (long long int)rti;
 }
 
@@ -1974,12 +1794,9 @@ long long int llrintl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 long int lround(double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT, plus should use RNA
+  // TODO: should use RNA
 
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
@@ -1994,8 +1811,8 @@ long int lround(double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  double rti = __sort_of_CPROVER_round_to_integral(FE_TOWARDZERO, xp);
+
+  double rti = __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
   return (long int)rti;
 }
 
@@ -2011,12 +1828,9 @@ long int lround(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 long int lroundf(float x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT, plus should use RNA
+  // should use RNA
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
 
@@ -2030,8 +1844,8 @@ long int lroundf(float x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  float rti = __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, xp);
+
+  float rti = __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
   return (long int)rti;
 }
 
@@ -2048,15 +1862,12 @@ long int lroundf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long int lroundl(long double x)
 {
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
 
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT, plus should use RNA
+  // TODO: should use RNA
   long double xp;
   if (x > 0.0) {
     xp = x + 0.5;
@@ -2067,8 +1878,8 @@ long int lroundl(long double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  long double rti = __sort_of_CPROVER_round_to_integrall(FE_TOWARDZERO, xp);
+
+  long double rti = __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
   return (long int)rti;
 }
 
@@ -2085,8 +1896,6 @@ long int lroundl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 long long int llround(double x)
 {
   // TODO : should be an all-in-one __CPROVER function to allow
@@ -2104,8 +1913,8 @@ long long int llround(double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  double rti = __sort_of_CPROVER_round_to_integral(FE_TOWARDZERO, xp);
+
+  double rti = __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
   return (long long int)rti;
 }
 
@@ -2121,12 +1930,9 @@ long long int llround(double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 long long int llroundf(float x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT, plus should use RNA
+  // TODO: should use RNA
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
 
@@ -2140,8 +1946,8 @@ long long int llroundf(float x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  float rti = __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, xp);
+
+  float rti = __round_to_integral(xp, FE_TOWARDZERO);
   return (long long int)rti;
 }
 
@@ -2158,12 +1964,9 @@ long long int llroundf(float x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long long int llroundl(long double x)
 {
-  // TODO : should be an all-in-one __CPROVER function to allow
-  // conversion to SMT, plus should use RNA
+  // TODO: should use RNA
   int saved_rounding_mode = fegetround();
   fesetround(FE_TOWARDZERO);
 
@@ -2177,8 +1980,8 @@ long long int llroundl(long double x)
   }
 
   fesetround(saved_rounding_mode);
-  
-  long double rti = __sort_of_CPROVER_round_to_integrall(FE_TOWARDZERO, xp);
+
+  long double rti = __CPROVER_round_to_integral(xp, FE_TOWARDZERO);
   return (long long int)rti;
 }
 
@@ -2203,11 +2006,9 @@ long long int llroundl(long double x)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-
 double modf(double x, double *iptr)
 {
-  *iptr = __sort_of_CPROVER_round_to_integral(FE_TOWARDZERO, x);
+  *iptr = __round_to_integral(x, FE_TOWARDZERO);
   return (x - *iptr);
 }
 
@@ -2223,11 +2024,9 @@ double modf(double x, double *iptr)
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
-  float modff(float x, float *iptr)
+float modff(float x, float *iptr)
 {
-  *iptr = __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, x);
+  *iptr = __CPROVER_round_to_integral(x, FE_TOWARDZERO);
   return (x - *iptr);
 }
 
@@ -2244,11 +2043,9 @@ float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
 #define __CPROVER_FENV_H_INCLUDED
 #endif
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
-  long double modfl(long double x, long double *iptr)
+long double modfl(long double x, long double *iptr)
 {
-  *iptr = __sort_of_CPROVER_round_to_integralf(FE_TOWARDZERO, x);
+  *iptr = __CPROVER_round_to_integral(x, FE_TOWARDZERO);
   return (x - *iptr);
 }
 
@@ -2256,8 +2053,7 @@ long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double
 
 /* FUNCTION: __sort_of_CPROVER_remainder */
 // TODO : Should be a real __CPROVER function to convert to SMT-LIB
-double __sort_of_CPROVER_round_to_integral (int rounding_mode, double d);
-  
+
 double __sort_of_CPROVER_remainder (int rounding_mode, double x, double y)
 {
   if (x == 0.0 || __CPROVER_isinfd(y))
@@ -2265,7 +2061,7 @@ double __sort_of_CPROVER_remainder (int rounding_mode, double x, double y)
 
   // Extended precision helps... a bit...
   long double div = x/y;
-  long double n = __sort_of_CPROVER_round_to_integral(rounding_mode,div);
+  long double n = __CPROVER_round_to_integral(rounding_mode, div);
   long double res = (-y * n) + x;   // TODO : FMA would be an improvement
   return res;
 }
@@ -2273,8 +2069,6 @@ double __sort_of_CPROVER_remainder (int rounding_mode, double x, double y)
 /* FUNCTION: __sort_of_CPROVER_remainderf */
 // TODO : Should be a real __CPROVER function to convert to SMT-LIB
 
-float __sort_of_CPROVER_round_to_integralf (int rounding_mode, float d);
-
 float __sort_of_CPROVER_remainderf (int rounding_mode, float x, float y)
 {
   if (x == 0.0f || __CPROVER_isinff(y))
@@ -2282,7 +2076,7 @@ float __sort_of_CPROVER_remainderf (int rounding_mode, float x, float y)
 
   // Extended precision helps... a bit...
   long double div = x/y;
-  long double n = __sort_of_CPROVER_round_to_integral(rounding_mode,div);
+  long double n = __CPROVER_round_to_integral(rounding_mode, div);
   long double res = (-y * n) + x;   // TODO : FMA would be an improvement
   return res;
 }
@@ -2290,8 +2084,6 @@ float __sort_of_CPROVER_remainderf (int rounding_mode, float x, float y)
 /* FUNCTION: __sort_of_CPROVER_remainderl */
 // TODO : Should be a real __CPROVER function to convert to SMT-LIB
 
-long double __sort_of_CPROVER_round_to_integrall (int rounding_mode, long double d);
-
 long double __sort_of_CPROVER_remainderl (int rounding_mode, long double x, long double y)
 {
   if (x == 0.0 || __CPROVER_isinfld(y))
@@ -2299,7 +2091,7 @@ long double __sort_of_CPROVER_remainderl (int rounding_mode, long double x, long
 
   // Extended precision helps... a bit...
   long double div = x/y;
-  long double n = __sort_of_CPROVER_round_to_integral(rounding_mode,div);
+  long double n = __CPROVER_round_to_integral(rounding_mode, div);
   long double res = (-y * n) + x;   // TODO : FMA would be an improvement
   return res;
 }
diff --git a/src/solvers/flattening/boolbv.h b/src/solvers/flattening/boolbv.h
index 78987b734ec..62d2703ee10 100644
--- a/src/solvers/flattening/boolbv.h
+++ b/src/solvers/flattening/boolbv.h
@@ -33,6 +33,7 @@ class byte_update_exprt;
 class concatenation_exprt;
 class extractbit_exprt;
 class extractbits_exprt;
+class floatbv_round_to_integral_exprt;
 class floatbv_typecast_exprt;
 class ieee_float_op_exprt;
 class overflow_result_exprt;
@@ -176,6 +177,8 @@ class boolbvt:public arrayst
   virtual bvt convert_floatbv_op(const ieee_float_op_exprt &);
   virtual bvt convert_floatbv_mod_rem(const binary_exprt &);
   virtual bvt convert_floatbv_typecast(const floatbv_typecast_exprt &expr);
+  virtual bvt
+  convert_floatbv_round_to_integral(const floatbv_round_to_integral_exprt &);
   virtual bvt convert_member(const member_exprt &expr);
   virtual bvt convert_with(const with_exprt &expr);
   virtual bvt convert_update(const update_exprt &);
diff --git a/src/solvers/flattening/boolbv_floatbv_op.cpp b/src/solvers/flattening/boolbv_floatbv_op.cpp
index 87dd3a8cd31..eb6d7eec40d 100644
--- a/src/solvers/flattening/boolbv_floatbv_op.cpp
+++ b/src/solvers/flattening/boolbv_floatbv_op.cpp
@@ -80,6 +80,24 @@ bvt boolbvt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr)
     return conversion_failed(expr);
 }
 
+bvt boolbvt::convert_floatbv_round_to_integral(
+  const floatbv_round_to_integral_exprt &expr)
+{
+  if(expr.op().type().id() == ID_floatbv)
+  {
+    float_utilst float_utils(prop);
+
+    float_utils.set_rounding_mode(convert_bv(expr.rounding_mode()));
+    float_utils.spec = ieee_float_spect{to_floatbv_type(expr.op().type())};
+
+    auto op_bv = convert_bv(expr.op());
+
+    return float_utils.round_to_integral(op_bv);
+  }
+  else
+    return conversion_failed(expr);
+}
+
 bvt boolbvt::convert_floatbv_op(const ieee_float_op_exprt &expr)
 {
   const exprt &lhs = expr.lhs();
diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp
index 13706f06bd5..4826cc27b88 100644
--- a/src/solvers/floatbv/float_utils.cpp
+++ b/src/solvers/floatbv/float_utils.cpp
@@ -136,6 +136,13 @@ bvt float_utilst::to_integer(
   return result;
 }
 
+bvt float_utilst::round_to_integral(const bvt &src)
+{
+  PRECONDITION(src.size() == spec.width());
+
+  return bv_utils.build_constant(0, spec.width());
+}
+
 bvt float_utilst::build_constant(const ieee_floatt &src)
 {
   unbiased_floatt result;
diff --git a/src/solvers/floatbv/float_utils.h b/src/solvers/floatbv/float_utils.h
index 989c8d53377..3be526ce318 100644
--- a/src/solvers/floatbv/float_utils.h
+++ b/src/solvers/floatbv/float_utils.h
@@ -133,6 +133,7 @@ class float_utilst
   bvt to_signed_integer(const bvt &src, std::size_t int_width);
   bvt to_unsigned_integer(const bvt &src, std::size_t int_width);
   bvt conversion(const bvt &src, const ieee_float_spect &dest_spec);
+  bvt round_to_integral(const bvt &);
 
   // relations
   enum class relt { LT, LE, EQ, GT, GE };
diff --git a/src/solvers/smt2/smt2_conv.cpp b/src/solvers/smt2/smt2_conv.cpp
index 3662b8f5e1e..f65331c9bd0 100644
--- a/src/solvers/smt2/smt2_conv.cpp
+++ b/src/solvers/smt2/smt2_conv.cpp
@@ -1199,6 +1199,10 @@ void smt2_convt::convert_expr(const exprt &expr)
   {
     convert_floatbv_typecast(to_floatbv_typecast_expr(expr));
   }
+  else if(expr.id() == ID_floatbv_round_to_integral)
+  {
+    convert_floatbv_round_to_integral(to_floatbv_round_to_integral_expr(expr));
+  }
   else if(expr.id()==ID_struct)
   {
     convert_struct(to_struct_expr(expr));
@@ -3291,6 +3295,23 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr)
   }
 }
 
+void smt2_convt::convert_floatbv_round_to_integral(
+  const floatbv_round_to_integral_exprt &expr)
+{
+  PRECONDITION(expr.type().id() == ID_floatbv);
+
+  if(use_FPA_theory)
+  {
+    out << "(fp.roundToIntegral ";
+    convert_rounding_mode_FPA(expr.rounding_mode());
+    out << ' ';
+    convert_expr(expr.op());
+    out << ")";
+  }
+  else
+    UNEXPECTEDCASE("TODO floatbv_round_to_integral without FPA");
+}
+
 void smt2_convt::convert_struct(const struct_exprt &expr)
 {
   const struct_typet &struct_type =
diff --git a/src/solvers/smt2/smt2_conv.h b/src/solvers/smt2/smt2_conv.h
index fa6d79aec05..a2a8168c8aa 100644
--- a/src/solvers/smt2/smt2_conv.h
+++ b/src/solvers/smt2/smt2_conv.h
@@ -32,6 +32,7 @@ Author: Daniel Kroening, kroening@kroening.com
 
 class floatbv_typecast_exprt;
 class ieee_float_op_exprt;
+class floatbv_round_to_integral_exprt;
 class union_typet;
 class update_bit_exprt;
 class update_bits_exprt;
@@ -145,6 +146,8 @@ class smt2_convt : public stack_decision_proceduret
   void convert_floatbv_div(const ieee_float_op_exprt &expr);
   void convert_floatbv_mult(const ieee_float_op_exprt &expr);
   void convert_floatbv_rem(const binary_exprt &expr);
+  void
+  convert_floatbv_round_to_integral(const floatbv_round_to_integral_exprt &);
   void convert_mod(const mod_exprt &expr);
   void convert_euclidean_mod(const euclidean_mod_exprt &expr);
   void convert_index(const index_exprt &expr);
diff --git a/src/solvers/smt2/smt2_parser.cpp b/src/solvers/smt2/smt2_parser.cpp
index add14229d3f..f8bd2b07099 100644
--- a/src/solvers/smt2/smt2_parser.cpp
+++ b/src/solvers/smt2/smt2_parser.cpp
@@ -1382,6 +1382,19 @@ void smt2_parsert::setup_expressions()
     return binary_exprt(op[0], ID_floatbv_rem, op[1]);
   };
 
+  expressions["fp.roundToIntegral"] = [this]
+  {
+    auto op = operands();
+
+    if(op.size() != 2)
+      throw error() << "fp.roundToIntegral takes two operands";
+
+    if(op[0].type().id() != ID_floatbv)
+      throw error() << "fp.roundToIntegral takes a FloatingPoint operand";
+
+    return floatbv_round_to_integral_exprt(op[0], op[1]);
+  };
+
   expressions["fp.eq"] = [this] {
     return function_application_ieee_float_eq(operands());
   };
diff --git a/src/util/floatbv_expr.h b/src/util/floatbv_expr.h
index c34f1c15af5..f0de4051c0e 100644
--- a/src/util/floatbv_expr.h
+++ b/src/util/floatbv_expr.h
@@ -83,6 +83,70 @@ inline floatbv_typecast_exprt &to_floatbv_typecast_expr(exprt &expr)
   return ret;
 }
 
+/// \brief Round a floating-point number to an integral value
+/// considering the given rounding mode
+class floatbv_round_to_integral_exprt : public binary_exprt
+{
+public:
+  floatbv_round_to_integral_exprt(exprt op, exprt rounding)
+    : binary_exprt(
+        op,
+        ID_floatbv_round_to_integral,
+        std::move(rounding),
+        op.type())
+  {
+  }
+
+  exprt &op()
+  {
+    return op0();
+  }
+
+  const exprt &op() const
+  {
+    return op0();
+  }
+
+  exprt &rounding_mode()
+  {
+    return op1();
+  }
+
+  const exprt &rounding_mode() const
+  {
+    return op1();
+  }
+};
+
+template <>
+inline bool can_cast_expr<floatbv_round_to_integral_exprt>(const exprt &base)
+{
+  return base.id() == ID_floatbv_round_to_integral;
+}
+
+/// \brief Cast an exprt to a \ref floatbv_round_to_integral_exprt
+///
+/// \a expr must be known to be \ref floatbv_round_to_integral_exprt.
+///
+/// \param expr: Source expression
+/// \return Object of type \ref floatbv_round_to_integral_exprt
+inline const floatbv_round_to_integral_exprt &
+to_floatbv_round_to_integral_expr(const exprt &expr)
+{
+  PRECONDITION(expr.id() == ID_floatbv_round_to_integral);
+  floatbv_round_to_integral_exprt::check(expr);
+  return static_cast<const floatbv_round_to_integral_exprt &>(expr);
+}
+
+/// \copydoc to_floatbv_round_to_integral_expr(const exprt &)
+inline floatbv_round_to_integral_exprt &
+to_floatbv_round_to_integral_expr(exprt &expr)
+{
+  PRECONDITION(expr.id() == ID_floatbv_round_to_integral);
+  floatbv_round_to_integral_exprt::check(expr);
+  return static_cast<floatbv_round_to_integral_exprt &>(expr);
+}
+
 /// \brief Evaluates to true if the operand is NaN
 class isnan_exprt : public unary_predicate_exprt
 {
diff --git a/src/util/irep_ids.def b/src/util/irep_ids.def
index 507043b0d7d..0fbb2a95263 100644
--- a/src/util/irep_ids.def
+++ b/src/util/irep_ids.def
@@ -560,6 +560,7 @@ IREP_ID_ONE(floatbv_div)
 IREP_ID_ONE(floatbv_mod)
 IREP_ID_ONE(floatbv_rem)
 IREP_ID_ONE(floatbv_typecast)
+IREP_ID_ONE(floatbv_round_to_integral)
 IREP_ID_ONE(compound_literal)
 IREP_ID_ONE(custom_bv)
 IREP_ID_ONE(custom_unsignedbv)