Basically, to (naïvely) implement modular reduction by
M = (1 << 128) + 0x87, you need to:
- take the upper 128 bits $x_H$ of the 256-bit (actually, 255-bit, since the high bit can never be set)
PCLMULQDQ result $x = x_H \| x_L$,
- shift them down by 128 bits,
- XOR-multiply them by
- XOR the result with the lower 128 bits $x_L$ of the 256-bit result, and
- repeat (see note below).
Step 3 could be implemented with
PCLMULQDQ, but since one of the multiplicands is a small constant, it's probably more efficient to just implement it with discrete shifts and XORs.
Note that the result of step 3 (and thus of step 4) can be longer than 128 bits — up to 6 bits longer, in this case — and may thus need to be reduced a second time. In this second reduction, $x_H$ is at most 6 bits long (and thus the result of step 3 is at most 13 bits long), so this pass is simpler to implement, and guaranteed not to overflow and require a third pass.
(Also, if you can be sure that one of the original multiplicands is at most 122 bits long, then $x_H$ is at most 121 bits long, and so step 3 cannot overflow 128 bits. This can simplify the implementation, and in particular, makes the second reduction pass unnecessary.)
What we're really doing here, in steps 2–4, is calculating
x' = x - M * x_H, where
x_H = x >> 128 and the multiplication and the subtraction are both carryless. Of course, carryless binary subtraction is just the same as carryless binary addition, i.e. XOR. Treating it as subtraction here generalizes better, though, and also more clearly explains why this method works: we're effectively just subtracting multiples of
M from the result so as to clear the upper half of it.
I have not looked at the Intel paper linked by fgrieu in detail, but I suspect they're doing something smarter and faster than this naïve reduction. I'll let him summarize the paper, if and when he has more time. :)