diff --git a/src/cdi/CDI.cc b/src/cdi/CDI.cc
index cfd58c8fa..4e294b496 100644
--- a/src/cdi/CDI.cc
+++ b/src/cdi/CDI.cc
@@ -533,7 +533,7 @@ double CDI::collapseMultigroupOpacitiesPlanck(std::vector<double> const &groupBo
   for (size_t g = 1; g < groupBounds.size(); ++g) {
     Check(planckSpectrum[g - 1] >= 0.0);
     Check(opacity[g - 1] >= 0.0);
-    sig_planck_sum += planckSpectrum[g - 1] * opacity[g - 1];
+    sig_planck_sum = FMA(planckSpectrum[g - 1], opacity[g - 1], sig_planck_sum);
     // Also collect some CDF data.
     emission_group_cdf[g - 1] = sig_planck_sum;
   }
@@ -597,7 +597,7 @@ double CDI::collapseMultigroupOpacitiesPlanck(std::vector<double> const &groupBo
   for (size_t g = 1; g < groupBounds.size(); ++g) {
     Check(planckSpectrum[g - 1] >= 0.0);
     Check(opacity[g - 1] >= 0.0);
-    sig_planck_sum += planckSpectrum[g - 1] * opacity[g - 1];
+    sig_planck_sum = FMA(planckSpectrum[g - 1], opacity[g - 1], sig_planck_sum);
   }
 
   //                         int_{\nu_0}^{\nu_G}{d\nu sigma(\nu,T) * B(\nu,T)}
diff --git a/src/cdi/CDI.hh b/src/cdi/CDI.hh
index 9749471b3..79e7e687f 100644
--- a/src/cdi/CDI.hh
+++ b/src/cdi/CDI.hh
@@ -16,6 +16,7 @@
 #include "GrayOpacity.hh"
 #include "MultigroupOpacity.hh"
 #include "ds++/Constexpr_Functions.hh"
+#include "ds++/FMA.hh"
 #include "ds++/Soft_Equivalence.hh"
 #include <algorithm>
 #include <array>
@@ -80,36 +81,16 @@ static inline double taylor_series_planck(double x) {
 
   double const xsqrd = x * x;
 
-  double taylor(coeff_21 * xsqrd);
-
-  taylor += coeff_19;
-  taylor *= xsqrd;
-
-  taylor += coeff_17;
-  taylor *= xsqrd;
-
-  taylor += coeff_15;
-  taylor *= xsqrd;
-
-  taylor += coeff_13;
-  taylor *= xsqrd;
-
-  taylor += coeff_11;
-  taylor *= xsqrd;
-
-  taylor += coeff_9;
-  taylor *= xsqrd;
-
-  taylor += coeff_7;
-  taylor *= xsqrd;
-
-  taylor += coeff_5;
-  taylor *= x;
-
-  taylor += coeff_4;
-  taylor *= x;
-
-  taylor += coeff_3;
+  double taylor = FMA(xsqrd, coeff_21, coeff_19);
+  taylor = FMA(taylor, xsqrd, coeff_17);
+  taylor = FMA(taylor, xsqrd, coeff_15);
+  taylor = FMA(taylor, xsqrd, coeff_13);
+  taylor = FMA(taylor, xsqrd, coeff_11);
+  taylor = FMA(taylor, xsqrd, coeff_9);
+  taylor = FMA(taylor, xsqrd, coeff_7);
+  taylor = FMA(taylor, xsqrd, coeff_5);
+  taylor = FMA(taylor, x, coeff_4);
+  taylor = FMA(taylor, x, coeff_3);
   taylor *= x * xsqrd * coeff;
 
   Ensure(taylor >= 0.0);