quantiles: force step update with stable input values

The algorithm was designed for estimating quantiles in streams of
integer values. When the estimate is equal to the input value, the
step state variable does not change. This causes problems for the
floating-point adaptation used for measurents of delay in chrony.

One problem is numerical instability due to the strict comparison of
the input value and the current estimate.

Another problem is with signals that are so stable that the nanosecond
resolution of the system functions becomes the limitation. There is a
large difference in the value of the step state variable, which
determines how quickly the estimate will adapt to a new distribution,
between signals that are constant in the nanosecond resolution and
signals that can move in two nanoseconds.

Change the estimate update to never consider the input value equal to
the current estimate and don't set the estimate exactly to the input
value. Keep it off by a quarter of the minimum step to force jumping
around the input value if it's constant and decreasing the step variable
to negative values. Also fix the initial adjustment to step at least by
the minimum step (the original algorithm is described with ceil(), not
fabs()).

quantiles.c

@@ -138,22 +138,26 @@ static void
 update_estimate(struct Quantile *quantile, double value, double p, double rand,
                 double min_step)
 {
-  if (value > quantile->est && rand > (1.0 - p)) {
+  if (value >= quantile->est) {
+    if (rand < (1.0 - p))
+        return;
     quantile->step += quantile->sign > 0 ? min_step : -min_step;
-    quantile->est += quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    quantile->est += quantile->step > min_step ? quantile->step : min_step;
     if (quantile->est > value) {
       quantile->step += value - quantile->est;
-      quantile->est = value;
+      quantile->est = value + min_step / 4.0;
     }
     if (quantile->sign < 0 && quantile->step > min_step)
       quantile->step = min_step;
     quantile->sign = 1;
-  } else if (value < quantile->est && rand > p) {
+  } else {
+    if (rand < p)
+      return;
     quantile->step += quantile->sign < 0 ? min_step : -min_step;
-    quantile->est -= quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    quantile->est -= quantile->step > min_step ? quantile->step : min_step;
     if (quantile->est < value) {
       quantile->step += quantile->est - value;
-      quantile->est = value;
+      quantile->est = value - min_step / 4.0;
     }
     if (quantile->sign > 0 && quantile->step > min_step)
       quantile->step = min_step;

test/unit/quantiles.c

@@ -28,7 +28,7 @@ test_unit(void)
 {
   int i, j, k, min_k, max_k, q, r, in_order, out_order;
   QNT_Instance inst;
-  double x;
+  double x, x2;
 
   in_order = out_order = 0;
 
@@ -63,6 +63,32 @@ test_unit(void)
     QNT_Reset(inst);
     TEST_CHECK(inst->n_set == 0);
 
+    for (j = 0; j < 3; j++) {
+      QNT_Reset(inst);
+      x = (i + 950) / 1e9;
+      if (random() % 10)
+        x2 = (i + 951) / 1e9;
+      else
+        x2 = x;
+      for (k = 0; k < 1000; k++)
+        QNT_Accumulate(inst, random() % 2 ? x : x2);
+      for (k = 0; k < inst->n_quants; k++) {
+        TEST_CHECK(inst->quants[k].est > x - 0.4e-9);
+        TEST_CHECK(inst->quants[k].est < x2 + 0.4e-9);
+        TEST_CHECK(inst->quants[k].step < -15e-9);
+        TEST_CHECK(inst->quants[k].step > -1000e-9);
+        if (min_k * 2 == q && k < inst->repeat) {
+          if (x == x2) {
+            TEST_CHECK(inst->quants[k].step < -750e-9);
+            TEST_CHECK(inst->quants[k].step > -1000e-9);
+          } else {
+            TEST_CHECK(inst->quants[k].step < -350e-9);
+            TEST_CHECK(inst->quants[k].step > -600e-9);
+          }
+        }
+      }
+    }
+
     QNT_DestroyInstance(inst);
   }