+
+static double
+round__ (double x, double mult, double fuzzbits, double adjustment)
+{
+ if (fuzzbits <= 0)
+ fuzzbits = settings_get_fuzzbits ();
+ adjustment += exp2 (fuzzbits - DBL_MANT_DIG);
+
+ x /= mult;
+ x = x >= 0. ? floor (x + adjustment) : -floor (-x + adjustment);
+ return x * mult;
+}
+
+double
+round_nearest (double x, double mult, double fuzzbits)
+{
+ return round__ (x, mult, fuzzbits, .5);
+}
+
+double
+round_zero (double x, double mult, double fuzzbits)
+{
+ return round__ (x, mult, fuzzbits, 0);
+}
+
+struct substring
+replace_string (struct expression *e,
+ struct substring haystack,
+ struct substring needle,
+ struct substring replacement,
+ double n)
+{
+ if (!needle.length
+ || haystack.length < needle.length
+ || n <= 0
+ || n == SYSMIS)
+ return haystack;
+
+ struct substring result = alloc_string (e, MAX_STRING);
+ result.length = 0;
+
+ size_t i = 0;
+ while (i <= haystack.length - needle.length)
+ if (!memcmp (&haystack.string[i], needle.string, needle.length))
+ {
+ size_t copy_len = MIN (replacement.length, MAX_STRING - result.length);
+ memcpy (&result.string[result.length], replacement.string, copy_len);
+ result.length += copy_len;
+ i += needle.length;
+
+ if (--n < 1)
+ break;
+ }
+ else
+ {
+ if (result.length < MAX_STRING)
+ result.string[result.length++] = haystack.string[i];
+ i++;
+ }
+ while (i < haystack.length && result.length < MAX_STRING)
+ result.string[result.length++] = haystack.string[i++];
+
+ return result;
+}
+
+static int
+compare_doubles (const void *a_, const void *b_)
+{
+ const double *ap = a_;
+ const double *bp = b_;
+ double a = *ap;
+ double b = *bp;
+
+ /* Sort SYSMIS to the end. */
+ return (a == b ? 0
+ : a == SYSMIS ? 1
+ : b == SYSMIS ? -1
+ : a > b ? 1 : -1);
+}
+
+double
+median (double *a, size_t n)
+{
+ /* Sort the array in-place, sorting SYSMIS to the end. */
+ qsort (a, n, sizeof *a, compare_doubles);
+
+ /* Drop SYSMIS. */
+ n = count_valid (a, n);
+
+ return (!n ? SYSMIS
+ : n % 2 ? a[n / 2]
+ : (a[n / 2 - 1] + a[n / 2]) / 2.0);
+}