diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 7b1a86d728ef4f181112cef1d8fe5f20379cade6..7b8c41518219728eb8c07e38d89d64cbf0394513 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -24,6 +24,7 @@
 #include <utility>
 #include <limits>
 #include <numeric>
+#include <numbers>
 
 #include "ExprNode.hh"
 #include "DataTree.hh"
@@ -6118,9 +6119,9 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
   switch (op_code)
     {
     case TrinaryOpcode::normcdf:
-      return (0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
+      return (0.5*(1+erf((v1-v2)/v3/numbers::sqrt2)));
     case TrinaryOpcode::normpdf:
-      return (1/(v3*sqrt(2*M_PI)*exp(pow((v1-v2)/v3, 2)/2)));
+      return (1/(v3*sqrt(2*numbers::pi)*exp(pow((v1-v2)/v3, 2)/2)));
     }
   // Suppress GCC warning
   exit(EXIT_FAILURE);
diff --git a/src/macro/Expressions.cc b/src/macro/Expressions.cc
index 54f747008cfc8e1c791c546276036cd6c3ff7bee..83fd7b3a55dce2bb45d3bc315002e97d8c75715d 100644
--- a/src/macro/Expressions.cc
+++ b/src/macro/Expressions.cc
@@ -18,6 +18,7 @@
  */
 
 #include <utility>
+#include <numbers>
 
 #include "Expressions.hh"
 
@@ -240,7 +241,7 @@ Real::normpdf(const BaseTypePtr &btp1, const BaseTypePtr &btp2) const
   auto btp22 = dynamic_pointer_cast<Real>(btp2);
   if (!btp12 || !btp22)
     throw StackTrace("Type mismatch for operands of `normpdf` operator");
-  return make_shared<Real>((1/(btp22->value*std::sqrt(2*M_PI)*std::exp(pow((value-btp12->value)/btp22->value, 2)/2))));
+  return make_shared<Real>((1/(btp22->value*std::sqrt(2*numbers::pi)*std::exp(pow((value-btp12->value)/btp22->value, 2)/2))));
 }
 
 RealPtr
@@ -250,7 +251,7 @@ Real::normcdf(const BaseTypePtr &btp1, const BaseTypePtr &btp2) const
   auto btp22 = dynamic_pointer_cast<Real>(btp2);
   if (!btp12 || !btp22)
     throw StackTrace("Type mismatch for operands of `normpdf` operator");
-  return make_shared<Real>((0.5*(1+std::erf((value-btp12->value)/btp22->value/M_SQRT2))));
+  return make_shared<Real>((0.5*(1+std::erf((value-btp12->value)/btp22->value/numbers::sqrt2))));
 }
 
 BaseTypePtr