How does this logic for the modulo operation implementation work?

Accelerated Computing CUDA CUDA Programming and Performance
1

Suppose we have an SM80 core with an unsigned integer x in SASS register R0, and an unsigned integer modulus q in register R7. Now, to obtain x % q in R0, we run:

 I2F.U32.RP R6, R7 
 IMAD.MOV R9, RZ, RZ, -R7 
 ISETP.NE.U32.AND P1, PT, R7, RZ, PT 
 MUFU.RCP R6, R6 
 IADD3 R4, R6, 0xffffffe, RZ 
 F2I.FTZ.U32.TRUNC.NTZ R5, R4 
 HFMA2.MMA R4, -RZ, RZ, 0, 0 
 IMAD R9, R9, R5, RZ 
 IMAD.HI.U32 R5, R5, R9, R4 
 IMAD.HI.U32 R5, R5, R0, RZ 
 IMAD.MOV R5, RZ, RZ, -R5 
 IMAD R0, R7, R5, R0 
 ISETP.GE.U32.AND P0, PT, R0, R7, PT 
 @P0 IADD3 R0, -R7, R0, RZ 
 ISETP.GE.U32.AND P0, PT, R0, R7, PT 
 @P0 IMAD.IADD R0, R0, 0x1, -R7 
 @!P1 LOP3.LUT R0, RZ, R7, RZ, 0x33, !PT

(I got this from compiling a modulo operation and assignment statement in GodBolt, see here).

I have a few questions about this code…

  • What are the @!Pn lines for n = 0, 1, 2 ? Are those predicated execution? Something else?
  • I don’t quite understand the use of the floating-point work we see here, followed be a conversion back to int: it seems to be taking 1.0f / q, and then turning off the least-significant-bit… huh? I kind of lost the logic there.
2

Yes, they indicate predicated execution.

3

MUFU.RCP provides a fast way to generate a reciprocal. Generating an accurate reciprocal approximation would be slower when using only integer operations and / or require big(-ish) tables.

When operating on integers, various variants of Newton-Raphson iteration require that the error is strictly to one side throughout to deliver correct results: either always an overestimate or always an underestimate. If I recall correctly, MUFU.RCP delivers results within ±1 ulp of the mathematical result, so subtracting 1 ulp would ensure the starting guess is always an underestimate. In uint32_t computation, addition of 0xFFFFFFFE is the same as subtracting 2.

While I did work on multiple versions of integer division emulation for CUDA while at NVIDIA, I do not recall all the details, nor could I either confirm or deny that the code shown here matches any of these versions. The HFMA2.MMA R4, -RZ, RZ, 0, 0 looks highly unusual, I cannot figure out what it does here. Looks like some sort of compiler magic, maybe some clever way of selecting either 0x0 vs 0xffffffff.

How did this question arise? Have you found a more efficient way of emulating 32-bit integer division?

4

How did this question arise? Have you found a more efficient way of emulating 32-bit integer division?

Oh, you give me too much credit… I’ve just been thinking about operations on very large integers, and specifically implementing an x % q operation with small q and large x, made up of multiple 32-bit values. So, I started by looking at what the small-q mod operation looks like, for some inspiration.

5

The implementation of % is typically / followed by multiply-subtract. If you search these forums, you should be able to locate various code versions of uint64_t division using a Halley iteration with cubic convergence that I posted here maybe five years ago or so.

Now if you are looking at very large integers, one can also formulate iteration versions with quartic and quintic, or generally n-th order, convergence:

http://numbers.computation.free.fr/Constants/Algorithms/inverse

As a general reference for arbitrary-precision arithmetic, I would recommend: R. P. Brent and P. Zimmermann, Modern Computer Arithmetic, Cambridge University Press 2010.

Ultimately the speed of division and modulo will depend on the speed of multiplication. My knowledge may be a bit dated, but in terms of increasing size of the operands the typical sequence of algorithms used to be: (1) Karatsuba (2) Toom-Cook, (3) FFT. Parallelizing that is a challenge. It may be more practical to search for open-source projects that already work on that.

This publication may be of interest:

A. P. Dieguez, M. Amor, R. Doallo, A. Nukada, S. Matsuoka: “Efficient high-precision integer multiplication on the GPU.” The International Journal of High Performance Computing Applications, Vol. 36, No. 3, May 2022, pp. 356-369.

6

I’m actually not interested in parallelizing over the length of the number; my interest is in that dimension being serial, and there being enough parallelism from having to compute % for lots of numbers. But thanks for the reference!

nor could I either confirm or deny that the code shown here matches any of these versions

I got this from GodBolt, you don’t need to confirm or deny anything… some of the secret sauce has to leaks out at the end :-P

7

If there are only a few different values of q, you could also use the division-by-constant => multiplication with pre-computed reciprocal (plus possibly a few shifts and adds) approach. The canonical reference is:

Torbjörn Granlund and Peter L. Montgomery, “Division by invariant integers using multiplication,” ACM SIGPLAN Notices, Vol. 29, June 1994, pp. 61-72

8

Certainly, the rule of thumb is: The best way to compute something is to avoid having to compute it (e.g. by precomputation).

9

huh ?

10

@CU_Steve Off-by-one error :-) 0xfffffffe = -2, of course. Thanks for spotting it. I should probably stop responding in these forums in the middle of the night …