diff --git a/src/math.F90 b/src/math.F90
index a386747a4..6edf2f1c3 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -336,51 +336,7 @@ contains
     integer(C_INT), intent(in) :: n          ! number of components to polynomial
     real(C_DOUBLE), intent(out):: factors(n) ! output leading coefficient
 
-    integer :: i
-
-    real(8) :: sqrtE               ! sqrt(energy)
-    real(8) :: beta                ! sqrt(atomic weight ratio * E / kT)
-    real(8) :: half_inv_dopp2      ! 0.5 / dopp**2
-    real(8) :: quarter_inv_dopp4   ! 0.25 / dopp**4
-    real(8) :: erf_beta            ! error function of beta
-    real(8) :: exp_m_beta2         ! exp(-beta**2)
-
-    sqrtE = sqrt(E)
-    beta = sqrtE * dopp
-    half_inv_dopp2 = HALF / dopp**2
-    quarter_inv_dopp4 = half_inv_dopp2**2
-
-    if (beta > 6.0_8) then
-      ! Save time, ERF(6) is 1 to machine precision.
-      ! beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
-      erf_beta = ONE
-      exp_m_beta2 = ZERO
-    else
-      erf_beta = erf(beta)
-      exp_m_beta2 = exp(-beta**2)
-    end if
-
-    ! Assume that, for sure, we'll use a second order (1/E, 1/V, const)
-    ! fit, and no less.
-
-    factors(1) = erf_beta / E
-    factors(2) = ONE / sqrtE
-    factors(3) = factors(1) * (half_inv_dopp2 + E) &
-         + exp_m_beta2 / (beta * SQRT_PI)
-
-    ! Perform recursive broadening of high order components
-    do i = 1, n-3
-      if (i /= 1) then
-        factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 &
-             + factors(i+1) * (E + (ONE + TWO * i) * half_inv_dopp2)
-      else
-        ! Although it's mathematically identical, factors(0) will contain
-        ! nothing, and we don't want to have to worry about memory.
-        factors(i+3) = factors(i+1)*(E + (ONE + TWO * i) * half_inv_dopp2)
-      end if
-    end do
-
-    ! call broaden_wmp_polynomials_cc(E, dopp, n, factors)
+    call broaden_wmp_polynomials_cc(E, dopp, n, factors)
 
   end subroutine broaden_wmp_polynomials
 
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 9f698f42a..7ab2a8450 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -827,7 +827,7 @@ double watt_spectrum_c(double a, double b) {
 // double complex w_derivative_c(double complex z, int order){
 //   double complex wv; // The resultant w(z) value
 
-//   const double complex twoi_sqrtpi = 2.0 / std::sqrt(PI) * I;
+//   const double complex twoi_sqrtpi = 2.0 / SQRT_PI * I;
 
 //   switch(order) {
 //     case 0:
@@ -859,6 +859,7 @@ void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
   double erf_beta;            // error function of beta
   double exp_m_beta2;         // exp(-beta**2)
   int i;
+  double ip1_dbl;
 
   sqrtE = std::sqrt(E);
   beta = sqrtE * dopp;
@@ -881,17 +882,20 @@ void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
   factors[0] = erf_beta / E;
   factors[1] = 1. / sqrtE;
   factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 /
-       (beta * std::sqrt(PI));
+       (beta * SQRT_PI);
 
   // Perform recursive broadening of high order components
   for (i = 0; i < n - 3; i++) {
+    ip1_dbl = static_cast<double>(i + 1);
     if (i != 0) {
-      factors[i + 3] = -factors[i - 1] * (i - 1.) * i * quarter_inv_dopp4 +
-           factors[i + 1] * (E + (1. + 2. * i) * half_inv_dopp2);
+      factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl *
+           quarter_inv_dopp4 + factors[i + 1] *
+           (E + (1. + 2. * ip1_dbl) * half_inv_dopp2);
     } else {
       // Although it's mathematically identical, factors[0] will contain
       // nothing, and we don't want to have to worry about memory.
-      factors[i + 3] = factors[i + 1]*(E + (1. + 2. * i) * half_inv_dopp2);
+      factors[i + 3] = factors[i + 1] *
+           (E + (1. + 2. * ip1_dbl) * half_inv_dopp2);
     }
   }
 }
diff --git a/src/math_functions.h b/src/math_functions.h
index 4e7b2e8c6..ab35d140a 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -16,6 +16,8 @@ namespace openmc {
 // modifying test results
 const double PI = 3.1415926535898;
 
+const double SQRT_PI = std::sqrt(PI);
+
 //==============================================================================
 // NORMAL_PERCENTILE calculates the percentile of the standard normal
 // distribution with a specified probability level
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 03309bafe..99e16712c 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -199,17 +199,16 @@ def test_broaden_wmp_polynomials():
     # First lets do beta > 6
     test_E = 0.5
     test_dopp = 100.  # approximately U235 at room temperature
-    n = 4
-    ref_val = [2., 1.41421356, 1.0001, 0.70731891]
+    n = 6
+
+    ref_val = [2., 1.41421356, 1.0001, 0.70731891, 0.50030001, 0.353907]
     test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n)
 
     assert np.allclose(ref_val, test_val)
 
     # now beta < 6
     test_dopp = 5.
-    ref_val = [1.99999885, 1.41421356, 1.04, 0.79195959]
+    ref_val = [1.99999885, 1.41421356, 1.04, 0.79195959, 0.6224, 0.50346003]
     test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n)
 
     assert np.allclose(ref_val, test_val)
-
-test_calc_zn()
\ No newline at end of file