Re: [PATCH 2/5] lib/sort: Use more efficient bottom-up heapsort variant

[Rasmus Villemoes]

On 21/02/2019 09.21, George Spelvin wrote:
>  
> +/**
> + * parent - given the offset of the child, find the offset of the parent.
> + * @i: the offset of the heap element whose parent is sought.  Non-zero.
> + * @lsbit: a precomputed 1-bit mask, equal to "size & -size"
> + * @size: size of each element
> + *
> + * In terms of array indexes, the parent of element j = i/size is simply
> + * (j-1)/2.  But when working in byte offsets, we can't use implicit
> + * truncation of integer divides.
> + *
> + * Fortunately, we only need one bit of the quotient, not the full divide.
> + * size has a least significant bit.  That bit will be clear if i is
> + * an even multiple of size, and set if it's an odd multiple.
> + *
> + * Logically, we're doing "if (i & lsbit) i -= size;", but since the
> + * branch is unpredictable, it's done with a bit of clever branch-free
> + * code instead.
> + */
> +__attribute_const__ __always_inline
> +static size_t parent(size_t i, unsigned int lsbit, size_t size)
> +{
> +	i -= size;
> +	i -= size & -(i & lsbit);
> +	return i / 2;
> +}
> +

Really nice :) I had to work through this by hand, but it's solid.

>  /**
>   * sort - sort an array of elements
>   * @base: pointer to data to sort
> @@ -125,21 +151,26 @@ static void generic_swap(void *a, void *b, int size)
>   * @cmp_func: pointer to comparison function
>   * @swap_func: pointer to swap function or NULL
>   *
> - * This function does a heapsort on the given array. You may provide a
> - * swap_func function optimized to your element type.
> + * This function does a heapsort on the given array.  You may provide a
> + * swap_func function if you need to do something more than a memory copy
> + * (e.g. fix up pointers or auxiliary data), but the built-in swap isn't
> + * usually a bottleneck.
>   *
>   * Sorting time is O(n log n) both on average and worst-case. While
>   * qsort is about 20% faster on average, it suffers from exploitable
>   * O(n*n) worst-case behavior and extra memory requirements that make
>   * it less suitable for kernel use.
>   */
> -
>  void sort(void *base, size_t num, size_t size,
>  	  int (*cmp_func)(const void *, const void *),
>  	  void (*swap_func)(void *, void *, int size))
>  {
>  	/* pre-scale counters for performance */
> -	int i = (num/2 - 1) * size, n = num * size, c, r;
> +	size_t n = num * size, a = (num/2) * size;
> +	unsigned const lsbit = size & -size;	/* Used to find parent */
> +

Nit: qualifier before type, "const unsigned". And this sets ZF, so a
paranoid check for zero size (cf. the other mail) by doing "if (!lsbit)
return;" is practically free. Though it's probably a bit obscure doing
it that way...

> +	if (!n)
> +		return;

I'd make that n <= 1. Shouldn't be much more costly.

> -		}
> -	}
> +	/*
> +	 * Loop invariants:
> +	 * 1. elements [a,n) satisfy the heap property (compare greater than
> +	 *    all of their children),
> +	 * 2. elements [n,num*size) are sorted, and
> +	 * 3. a <= b <= c <= d <= n (whenever they are valid).
> +	 */
> +	for (;;) {
> +		size_t b, c, d;
>  

> +		if (a)			/* Building heap: sift down --a */
> +			a -= size;
> +		else if (n -= size)	/* Sorting: Extract root to --n */
> +			swap_func(base, base + n, size);
> +		else			/* Sort complete */
> +			break;
> +
> +		/*
> +		 * Sift element at "a" down into heap.  This is the
> +		 * "bottom-up" variant, which significantly reduces
> +		 * calls to cmp_func(): we find the sift-down path all
> +		 * the way to the leaves (one compare per level), then
> +		 * backtrack to find where to insert the target element.
> +		 *
> +		 * Because elements tend to sift down close to the leaves,
> +		 * this uses fewer compares than doing two per level
> +		 * on the way down.  (A bit more than half as many on
> +		 * average, 3/4 worst-case.)
> +		 */
> +		for (b = a; c = 2*b + size, (d = c + size) < n;)
> +			b = cmp_func(base + c, base + d) >= 0 ? c : d;
> +		if (d == n)	/* Special case last leaf with no sibling */
> +			b = c;
> +
> +		/* Now backtrack from "b" to the correct location for "a" */
> +		while (b != a && cmp_func(base + a, base + b) >= 0)
> +			b = parent(b, lsbit, size);
> +		c = b;			/* Where "a" belongs */
> +		while (b != a) {	/* Shift it into place */
> +			b = parent(b, lsbit, size);
> +			swap_func(base + b, base + c, size);
>  		}
>  	}
>  }
> 

Nice!

Rasmus