diff --git a/reference.c b/reference.c
index 71737c8..d6eea69 100644
--- a/reference.c
+++ b/reference.c
@@ -325,23 +325,6 @@ REF_GetLeapMode(void)
   return leap_mode;
 }
 
-/* ================================================== */
-
-static double
-Sqr(double x)
-{
-  return x*x;
-}
-
-/* ================================================== */
-#if 0
-static double
-Cube(double x)
-{
-  return x*x*x;
-}
-#endif
-
 /* ================================================== */
 /* Update the drift coefficients to the file. */
 
@@ -1055,8 +1038,8 @@ REF_SetReference(int stratum,
     /* Set new frequency based on weighted average of old and new skew. With
        manual reference the old frequency has no weight. */
 
-    old_weight = leap != LEAP_Unsynchronised ? 1.0 / Sqr(previous_skew) : 0.0;
-    new_weight = 3.0 / Sqr(new_skew);
+    old_weight = leap != LEAP_Unsynchronised ? 1.0 / SQUARE(previous_skew) : 0.0;
+    new_weight = 3.0 / SQUARE(new_skew);
 
     sum_weight = old_weight + new_weight;
 
@@ -1065,7 +1048,7 @@ REF_SetReference(int stratum,
     delta_freq1 = previous_freq - our_frequency;
     delta_freq2 = new_freq - our_frequency;
 
-    skew1 = sqrt((Sqr(delta_freq1) * old_weight + Sqr(delta_freq2) * new_weight) / sum_weight);
+    skew1 = sqrt((SQUARE(delta_freq1) * old_weight + SQUARE(delta_freq2) * new_weight) / sum_weight);
     skew2 = (previous_skew * old_weight + new_skew * new_weight) / sum_weight;
     our_skew = skew1 + skew2;
 
@@ -1125,11 +1108,11 @@ REF_SetReference(int stratum,
 
   /* Update the moving average of squares of offset, quickly on start */
   if (avg2_moving) {
-    avg2_offset += 0.1 * (our_offset * our_offset - avg2_offset);
+    avg2_offset += 0.1 * (SQUARE(our_offset) - avg2_offset);
   } else {
-    if (avg2_offset > 0.0 && avg2_offset < our_offset * our_offset)
+    if (avg2_offset > 0.0 && avg2_offset < SQUARE(our_offset))
       avg2_moving = 1;
-    avg2_offset = our_offset * our_offset;
+    avg2_offset = SQUARE(our_offset);
   }
 }
 
diff --git a/samplefilt.c b/samplefilt.c
index 2c5f936..af0b34c 100644
--- a/samplefilt.c
+++ b/samplefilt.c
@@ -76,7 +76,7 @@ SPF_CreateInstance(int min_samples, int max_samples, double max_dispersion, doub
   /* Set the first estimate to the system precision */
   filter->avg_var_n = 0;
   filter->avg_var = LCL_GetSysPrecisionAsQuantum() * LCL_GetSysPrecisionAsQuantum();
-  filter->max_var = max_dispersion * max_dispersion;
+  filter->max_var = SQUARE(max_dispersion);
   filter->combine_ratio = combine_ratio;
   filter->samples = MallocArray(NTP_Sample, filter->max_samples);
   filter->selected = MallocArray(int, filter->max_samples);
diff --git a/smooth.c b/smooth.c
index a21dcd8..4125ea0 100644
--- a/smooth.c
+++ b/smooth.c
@@ -144,7 +144,7 @@ update_stages(void)
      is equal to the offset that should be smoothed out */
 
   s1 = smooth_offset / max_wander;
-  s2 = smooth_freq * smooth_freq / (2.0 * max_wander * max_wander);
+  s2 = SQUARE(smooth_freq) / (2.0 * SQUARE(max_wander));
   
   /* Calculate the lengths of the 1st and 3rd stage assuming there is no
      frequency limit.  The direction of the 1st stage is selected so that
diff --git a/sources.c b/sources.c
index 1822912..7df6279 100644
--- a/sources.c
+++ b/sources.c
@@ -550,7 +550,7 @@ combine_sources(int n_sel_sources, struct timespec *ref_time, double *offset,
     src_offset += elapsed * src_frequency;
     src_offset_sd += elapsed * src_frequency_sd;
     offset_weight = 1.0 / sources[index]->sel_info.root_distance;
-    frequency_weight = 1.0 / (src_frequency_sd * src_frequency_sd);
+    frequency_weight = 1.0 / SQUARE(src_frequency_sd);
 
     DEBUG_LOG("combining index=%d oweight=%e offset=%e osd=%e fweight=%e freq=%e fsd=%e skew=%e",
               index, offset_weight, src_offset, src_offset_sd,
@@ -558,13 +558,13 @@ combine_sources(int n_sel_sources, struct timespec *ref_time, double *offset,
 
     sum_offset_weight += offset_weight;
     sum_offset += offset_weight * src_offset;
-    sum2_offset_sd += offset_weight * (src_offset_sd * src_offset_sd +
+    sum2_offset_sd += offset_weight * (SQUARE(src_offset_sd) +
         (src_offset - *offset) * (src_offset - *offset));
 
     sum_frequency_weight += frequency_weight;
     sum_frequency += frequency_weight * src_frequency;
-    inv_sum2_frequency_sd += 1.0 / (src_frequency_sd * src_frequency_sd);
-    inv_sum2_skew += 1.0 / (src_skew * src_skew);
+    inv_sum2_frequency_sd += 1.0 / SQUARE(src_frequency_sd);
+    inv_sum2_skew += 1.0 / SQUARE(src_skew);
 
     combined++;
   }
diff --git a/sourcestats.c b/sourcestats.c
index e13012a..2d28aff 100644
--- a/sourcestats.c
+++ b/sourcestats.c
@@ -553,7 +553,7 @@ SST_DoNewRegression(SST_Stats inst)
       sd_weight = 1.0;
       if (peer_distances[i] > min_distance)
         sd_weight += (peer_distances[i] - min_distance) / sd;
-      weights[i] = sd_weight * sd_weight;
+      weights[i] = SQUARE(sd_weight);
     }
   }
 
diff --git a/util.h b/util.h
index 4c55651..e3d6767 100644
--- a/util.h
+++ b/util.h
@@ -200,4 +200,6 @@ extern void UTI_GetRandomBytes(void *buf, unsigned int len);
 /* Macro to clamp a value between two values */
 #define CLAMP(min, x, max) (MAX((min), MIN((x), (max))))
 
+#define SQUARE(x) ((x) * (x))
+
 #endif /* GOT_UTIL_H */