int32_t stmesh_utils_round(float x, float eps, int32_t r)
  { 
    demand(eps > 0.0, "length unit {eps} must be positive");
    double div; /* 1 if round to any, 2 for round to even/odd. */
    double rem; /* 0 if round to any/even, 1 if round to odd. */
    if (r >= 0)
      { demand((r == 0) || (r == 1), "invalid remainder {r}");
        div = 2.0; rem = (double) r;
      }
    else
      { div = 1.0; rem = 0.0; }
    
    /* Compute {qd = (x/eps - rem)/div}, the amount that should be rounded to integer: */
    double qd = (((double)x)/((double)eps) - rem)/div;
    demand(isinf(qd) == 0, "overflow in rounding");
    
    /* Round {qd} to an integer {qn}, still represented as double: */
    double qlo = floor(qd + 0.5); /* Rounds to nearest; half-integers up. */
    double qhi = ceil(qd - 0.5); /* Rounds to nearest; half-integers down. */
    double qn;
    if (qlo == qhi)
      { qn = qlo; }
    else 
      { /* Tie in rounding: round {qn} to even. */
        if (fmod(qlo,2.0) != 0)
          { qn = qlo; }
        else
          { qn = qhi; }
      }
    assert(qn == floor(qn));
    
    /* Compute the result, still as double: */
    double qr = div*qn + rem;
    
    /* Convert to integer: */
    demand(qr >= ((double)INT32_MIN), "overflow (-) in rounding");
    demand(qr <= ((double)INT32_MAX), "overflow (+) in rounding");
    int32_t qk = (int32_t)qr;
    assert(qr == (double)qk);

    return qk;
  }