diff --git a/rtengine/EdgePreservingDecomposition.cc b/rtengine/EdgePreservingDecomposition.cc
index eac42ffaa..1ddd17107 100644
--- a/rtengine/EdgePreservingDecomposition.cc
+++ b/rtengine/EdgePreservingDecomposition.cc
@@ -42,7 +42,8 @@ float *SparseConjugateGradient(void Ax(float *Product, float *x, void *Pass), fl
     }
 
     //s is preconditionment of r. Without, direct to r.
-    float *s = r, rs = 0.0f;
+    float *s = r;
+    double rs = 0.0; // use double precision for large summations
 
     if(Preconditioner != nullptr) {
         s = new float[n];
@@ -77,7 +78,7 @@ float *SparseConjugateGradient(void Ax(float *Product, float *x, void *Pass), fl
 
     for(iterate = 0; iterate < MaximumIterates; iterate++) {
         //Get step size alpha, store ax while at it.
-        float ab = 0.0f;
+        double ab = 0.0; // use double precision for large summations
         Ax(ax, d, Pass);
 #ifdef _OPENMP
         #pragma omp parallel for reduction(+:ab)
@@ -94,7 +95,7 @@ float *SparseConjugateGradient(void Ax(float *Product, float *x, void *Pass), fl
         ab = rs / ab;
 
         //Update x and r with this step size.
-        float rms = 0.0;
+        double rms = 0.0; // use double precision for large summations
 #ifdef _OPENMP
         #pragma omp parallel for reduction(+:rms)
 #endif
diff --git a/rtengine/PF_correct_RT.cc b/rtengine/PF_correct_RT.cc
index 1a937b409..70c334928 100644
--- a/rtengine/PF_correct_RT.cc
+++ b/rtengine/PF_correct_RT.cc
@@ -70,7 +70,7 @@ void ImProcFunctions::PF_correct_RT(LabImage * src, LabImage * dst, double radiu
         gaussianBlur (src->b, tmp1->b, src->W, src->H, radius);
     }
 
-    float chromave = 0.0f;
+    double chromave = 0.0; // use double precision for large summations
 
 #ifdef _OPENMP
     #pragma omp parallel
@@ -382,7 +382,7 @@ void ImProcFunctions::PF_correct_RTcam(CieImage * src, CieImage * dst, double ra
         gaussianBlur (srbb, tmbb, src->W, src->H, radius);
     }
 
-    float chromave = 0.0f;
+    double chromave = 0.0; // use double precision for large summations
 
 #ifdef __SSE2__
 
@@ -1042,7 +1042,7 @@ void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, double radius
 
 
 // begin chroma badpixels
-    float chrommed = 0.f;
+    double chrommed = 0.0; // use double precision for large summations
 #ifdef _OPENMP
     #pragma omp parallel for reduction(+:chrommed)
 #endif
@@ -1646,7 +1646,7 @@ void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, double radius
 
     if(mode == 3) {
 // begin chroma badpixels
-        float chrommed = 0.f;
+        double chrommed = 0.0; // use double precision for large summations
 #ifdef _OPENMP
         #pragma omp parallel for reduction(+:chrommed)
 #endif
diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc
index ecbd4fa37..d274a056a 100644
--- a/rtengine/improcfun.cc
+++ b/rtengine/improcfun.cc
@@ -2002,8 +2002,7 @@ void ImProcFunctions::ciecam_02float (CieImage* ncie, float adap, int pW, int pw
                 hist16Q.clear();
             }
 
-            float sum = 0.f;
-//            float sumQ = 0.f;
+            double sum = 0.0; // use double precision for large summations
 
 #ifdef _OPENMP
             const int numThreads = min (max (width * height / 65536, 1), omp_get_max_threads());
@@ -2023,7 +2022,6 @@ void ImProcFunctions::ciecam_02float (CieImage* ncie, float adap, int pW, int pw
                     hist16Qthr.clear();
                 }
 
-                //    #pragma omp for reduction(+:sum,sumQ)
                 #pragma omp for reduction(+:sum)
 
 
diff --git a/rtengine/ipwavelet.cc b/rtengine/ipwavelet.cc
index 4767a7fba..cd1bbc5f5 100644
--- a/rtengine/ipwavelet.cc
+++ b/rtengine/ipwavelet.cc
@@ -1271,7 +1271,7 @@ void ImProcFunctions::Aver( float *  RESTRICT DataList, int datalen, float &aver
 
     //find absolute mean
     int countP = 0, countN = 0;
-    float averaP = 0.f, averaN = 0.f;
+    double averaP = 0.0, averaN = 0.0; // use double precision for large summations
 
     float thres = 5.f;//different fom zero to take into account only data large enough
     max = 0.f;
@@ -1332,7 +1332,7 @@ void ImProcFunctions::Aver( float *  RESTRICT DataList, int datalen, float &aver
 void ImProcFunctions::Sigma( float *  RESTRICT DataList, int datalen, float averagePlus, float averageNeg, float &sigmaPlus, float &sigmaNeg)
 {
     int countP = 0, countN = 0;
-    float variP = 0.f, variN = 0.f;
+    double variP = 0.0, variN = 0.0; // use double precision for large summations
     float thres = 5.f;//different fom zero to take into account only data large enough
 
 #ifdef _OPENMP
@@ -1687,7 +1687,7 @@ void ImProcFunctions::WaveletcontAllL(LabImage * labco, float ** varhue, float *
     float contrast = cp.contrast;
     float multL = (float)contrast * (maxl - 1.f) / 100.f + 1.f;
     float multH = (float) contrast * (maxh - 1.f) / 100.f + 1.f;
-    double avedbl = 0.f; // use double precision for big summations
+    double avedbl = 0.0; // use double precision for large summations
     float max0 = 0.f;
     float min0 = FLT_MAX;
 
diff --git a/rtengine/tmo_fattal02.cc b/rtengine/tmo_fattal02.cc
index dfdeb56b5..0a88958d1 100644
--- a/rtengine/tmo_fattal02.cc
+++ b/rtengine/tmo_fattal02.cc
@@ -276,7 +276,7 @@ float calculateGradients (Array2Df* H, Array2Df* G, int k, bool multithread)
     const int width = H->getCols();
     const int height = H->getRows();
     const float divider = pow ( 2.0f, k + 1 );
-    float avgGrad = 0.0f;
+    double avgGrad = 0.0; // use double precision for large summations
 
     #pragma omp parallel for reduction(+:avgGrad) if(multithread)