Hi everyone,

I know I was told what the algorithm was already, but I can't find it in my email history, so I'm forced to ask again... is there a reference to the algorithm that's used to calculate the reciprocal constants used in the integer division optimisations for x86 and AArch64? It's more to gain a personal understanding of the algorithm so I can possibly learn to write assembly language routines that both do it faster and produce the same results (my original implementation, such required double-width registers, produced magic constants that were bit-shifted to calc_magic_div_unsigned but which nonetheless produced correct results when the actual divisions were applied).

Gareth aka. Kit

Something tells me I should purchase that book - I sense it could reveal some interesting insights.

Note that while I understand the concept of turning integer division into multiplication (indeed, I implemented the first version into x86 before it was improved with the "calc_divconst_magic_unsigned" routine, and then implemented it for AArch64), the algorithm that is used in "calc_divconst_magic_unsigned" I don't quite get if just for the lack of comments and references, although I am studying it more closely.

Still, I figure I'll put down Hacker's Delight as a future purchase so I have a reputable source rather than just an online one (I know there's one somewhere that isn't behind a paywall).

Gareth aka. Kit

Well, you'll be pleased to know that I am now back from my camping and hiking trip in the Lake District and have actually purchased "Hacker's Delight", and I have to say it's extremely insightful. I would love to play around and implement some of the suggested optimisations.

Some of them will require some clever insight and manipulation at the node level though (e.g. converting "(x mod constant) = 0" into something that produces fast code unique to checking to see if the division is exact)... the main problem is seeing through all the typecasts and the different loads.

Gareth aka. Kit

Hey Marģers,

So I've been experimenting with your suggestion, and it looks like a resounding success!  I added some new tests to the "bdiv" bench test to see how it performs.  16-bit multiplications don't get improved as well as they could be on x86_64 because the intermediate values are all extended to 32-bit, but you can still see a massive improvement on all the unsigned divisions by 7.

I'm currently running the test suite for i386-win32 and x86_64-win64, then will implement the same system for AArch64. Thanks so much for spotting the improvement.

Gareth aka. Kit

P.S. I still need to work on the signed modulus at some point.

P.P.S. Timings:

Trunk:

Division compilation and timing test (using constants from System and Sysutils)
-------------------------------------------------------------------------------               Unsigned 16-bit division by 1 - Pass - average iteration duration: 0.349 ns
               Unsigned 16-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 16-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 16-bit modulus by 2 - Pass - average iteration duration: 0.349 ns
              Unsigned 16-bit division by 3 - Pass - average iteration duration: 0.640 ns
               Unsigned 16-bit modulus by 3 - Pass - average iteration duration: 0.757 ns
              Unsigned 16-bit division by 7 - Pass - average iteration duration: 0.873 ns
               Unsigned 16-bit modulus by 7 - Pass - average iteration duration: 0.990 ns
             Unsigned 16-bit division by 10 - Pass - average iteration duration: 0.582 ns
              Unsigned 16-bit modulus by 10 - Pass - average iteration duration: 0.815 ns
            Unsigned 16-bit division by 100 - Pass - average iteration duration: 0.582 ns
             Unsigned 16-bit modulus by 100 - Pass - average iteration duration: 0.815 ns
          Unsigned 16-bit division by 1,000 - Pass - average iteration duration: 0.582 ns
           Unsigned 16-bit modulus by 1,000 - Pass - average iteration duration: 0.757 ns
              Unsigned 32-bit division by 1 - Pass - average iteration duration: 0.291 ns
               Unsigned 32-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 32-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 32-bit modulus by 2 - Pass - average iteration duration: 0.291 ns
              Unsigned 32-bit division by 3 - Pass - average iteration duration: 0.582 ns
               Unsigned 32-bit modulus by 3 - Pass - average iteration duration: 0.698 ns
              Unsigned 32-bit division by 7 - Pass - average iteration duration: 0.815 ns
               Unsigned 32-bit modulus by 7 - Pass - average iteration duration: 0.990 ns
             Unsigned 32-bit division by 10 - Pass - average iteration duration: 0.582 ns
              Unsigned 32-bit modulus by 10 - Pass - average iteration duration: 0.757 ns
            Unsigned 32-bit division by 100 - Pass - average iteration duration: 0.582 ns
             Unsigned 32-bit modulus by 100 - Pass - average iteration duration: 0.698 ns
          Unsigned 32-bit division by 1,000 - Pass - average iteration duration: 0.524 ns
           Unsigned 32-bit modulus by 1,000 - Pass - average iteration duration: 0.757 ns
         Unsigned 32-bit division by 60,000 - Pass - average iteration duration: 0.524 ns
          Unsigned 32-bit modulus by 60,000 - Pass - average iteration duration: 0.698 ns
        Unsigned 32-bit division by 146,097 - Pass - average iteration duration: 0.582 ns
         Unsigned 32-bit modulus by 146,097 - Pass - average iteration duration: 0.698 ns
      Unsigned 32-bit division by 3,600,000 - Pass - average iteration duration: 0.582 ns
       Unsigned 32-bit modulus by 3,600,000 - Pass - average iteration duration: 0.698 ns
              Unsigned 64-bit division by 1 - Pass - average iteration duration: 0.349 ns
               Unsigned 64-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 64-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 64-bit modulus by 2 - Pass - average iteration duration: 0.349 ns
              Unsigned 64-bit division by 3 - Pass - average iteration duration: 0.582 ns
               Unsigned 64-bit modulus by 3 - Pass - average iteration duration: 0.698 ns
              Unsigned 64-bit division by 7 - Pass - average iteration duration: 0.873 ns
               Unsigned 64-bit modulus by 7 - Pass - average iteration duration: 0.990 ns
             Unsigned 64-bit division by 10 - Pass - average iteration duration: 0.524 ns
              Unsigned 64-bit modulus by 10 - Pass - average iteration duration: 0.698 ns
            Unsigned 64-bit division by 100 - Pass - average iteration duration: 0.815 ns
             Unsigned 64-bit modulus by 100 - Pass - average iteration duration: 0.990 ns
  Unsigned 64-bit division by 1,000,000,000 - Pass - average iteration duration: 0.815 ns
   Unsigned 64-bit modulus by 1,000,000,000 - Pass - average iteration duration: 0.990 ns
                Signed 32-bit division by 1 - Pass - average iteration duration: 0.349 ns
                 Signed 32-bit modulus by 1 - Pass - average iteration duration: 0.291 ns
              Signed 32-bit division by 100 - Pass - average iteration duration: 0.698 ns
               Signed 32-bit modulus by 100 - Pass - average iteration duration: 6.403 ns
                Signed 64-bit division by 1 - Pass - average iteration duration: 0.349 ns
                 Signed 64-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
               Signed 64-bit division by 10 - Pass - average iteration duration: 0.582 ns
                Signed 64-bit modulus by 10 - Pass - average iteration duration: 6.577 ns
               Signed 64-bit division by 18 - Pass - average iteration duration: 0.582 ns
                Signed 64-bit modulus by 18 - Pass - average iteration duration: 6.461 ns
               Signed 64-bit division by 24 - Pass - average iteration duration: 0.640 ns
                Signed 64-bit modulus by 24 - Pass - average iteration duration: 6.519 ns
              Signed 64-bit division by 100 - Pass - average iteration duration: 0.757 ns
               Signed 64-bit modulus by 100 - Pass - average iteration duration: 6.461 ns
              Signed 64-bit division by 153 - Pass - average iteration duration: 0.640 ns
               Signed 64-bit modulus by 153 - Pass - average iteration duration: 6.519 ns
            Signed 64-bit division by 1,461 - Pass - average iteration duration: 0.698 ns
             Signed 64-bit modulus by 1,461 - Pass - average iteration duration: 6.577 ns
Signed 64-bit division by 10,000 (Currency) - Pass - average iteration duration: 0.640 ns
 Signed 64-bit modulus by 10,000 (Currency) - Pass - average iteration duration: 6.694 ns
       Signed 64-bit division by 86,400,000 - Pass - average iteration duration: 0.640 ns
        Signed 64-bit modulus by 86,400,000 - Pass - average iteration duration: 6.636 ns

ok
- Sum of average durations: 96.275 ns
- Overall average duration: 1.375 ns

Extended magic-div:

Division compilation and timing test (using constants from System and Sysutils)
-------------------------------------------------------------------------------               Unsigned 16-bit division by 1 - Pass - average iteration duration: 0.349 ns
               Unsigned 16-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 16-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 16-bit modulus by 2 - Pass - average iteration duration: 0.291 ns
              Unsigned 16-bit division by 3 - Pass - average iteration duration: 0.524 ns
               Unsigned 16-bit modulus by 3 - Pass - average iteration duration: 0.640 ns
              Unsigned 16-bit division by 7 - Pass - average iteration duration: 0.466 ns
               Unsigned 16-bit modulus by 7 - Pass - average iteration duration: 0.640 ns
             Unsigned 16-bit division by 10 - Pass - average iteration duration: 0.582 ns
              Unsigned 16-bit modulus by 10 - Pass - average iteration duration: 0.698 ns
            Unsigned 16-bit division by 100 - Pass - average iteration duration: 0.524 ns
             Unsigned 16-bit modulus by 100 - Pass - average iteration duration: 0.582 ns
          Unsigned 16-bit division by 1,000 - Pass - average iteration duration: 0.524 ns
           Unsigned 16-bit modulus by 1,000 - Pass - average iteration duration: 0.640 ns
              Unsigned 32-bit division by 1 - Pass - average iteration duration: 0.349 ns
               Unsigned 32-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 32-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 32-bit modulus by 2 - Pass - average iteration duration: 0.291 ns
              Unsigned 32-bit division by 3 - Pass - average iteration duration: 0.466 ns
               Unsigned 32-bit modulus by 3 - Pass - average iteration duration: 0.582 ns
              Unsigned 32-bit division by 7 - Pass - average iteration duration: 0.524 ns
               Unsigned 32-bit modulus by 7 - Pass - average iteration duration: 0.582 ns
             Unsigned 32-bit division by 10 - Pass - average iteration duration: 0.524 ns
              Unsigned 32-bit modulus by 10 - Pass - average iteration duration: 0.640 ns
            Unsigned 32-bit division by 100 - Pass - average iteration duration: 0.524 ns
             Unsigned 32-bit modulus by 100 - Pass - average iteration duration: 0.582 ns
          Unsigned 32-bit division by 1,000 - Pass - average iteration duration: 0.466 ns
           Unsigned 32-bit modulus by 1,000 - Pass - average iteration duration: 0.582 ns
         Unsigned 32-bit division by 60,000 - Pass - average iteration duration: 0.524 ns
          Unsigned 32-bit modulus by 60,000 - Pass - average iteration duration: 0.582 ns
        Unsigned 32-bit division by 146,097 - Pass - average iteration duration: 0.524 ns
         Unsigned 32-bit modulus by 146,097 - Pass - average iteration duration: 0.524 ns
      Unsigned 32-bit division by 3,600,000 - Pass - average iteration duration: 0.524 ns
       Unsigned 32-bit modulus by 3,600,000 - Pass - average iteration duration: 0.524 ns
              Unsigned 64-bit division by 1 - Pass - average iteration duration: 0.291 ns
               Unsigned 64-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Unsigned 64-bit division by 2 - Pass - average iteration duration: 0.349 ns
               Unsigned 64-bit modulus by 2 - Pass - average iteration duration: 0.407 ns
              Unsigned 64-bit division by 3 - Pass - average iteration duration: 0.582 ns
               Unsigned 64-bit modulus by 3 - Pass - average iteration duration: 0.698 ns
              Unsigned 64-bit division by 7 - Pass - average iteration duration: 0.815 ns
               Unsigned 64-bit modulus by 7 - Pass - average iteration duration: 1.048 ns
             Unsigned 64-bit division by 10 - Pass - average iteration duration: 0.582 ns
              Unsigned 64-bit modulus by 10 - Pass - average iteration duration: 0.640 ns
            Unsigned 64-bit division by 100 - Pass - average iteration duration: 0.815 ns
             Unsigned 64-bit modulus by 100 - Pass - average iteration duration: 0.990 ns
  Unsigned 64-bit division by 1,000,000,000 - Pass - average iteration duration: 0.757 ns
   Unsigned 64-bit modulus by 1,000,000,000 - Pass - average iteration duration: 1.048 ns
                Signed 32-bit division by 1 - Pass - average iteration duration: 0.349 ns
                 Signed 32-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
              Signed 32-bit division by 100 - Pass - average iteration duration: 0.757 ns
               Signed 32-bit modulus by 100 - Pass - average iteration duration: 6.461 ns
                Signed 64-bit division by 1 - Pass - average iteration duration: 0.349 ns
                 Signed 64-bit modulus by 1 - Pass - average iteration duration: 0.349 ns
               Signed 64-bit division by 10 - Pass - average iteration duration: 0.640 ns
                Signed 64-bit modulus by 10 - Pass - average iteration duration: 6.694 ns
               Signed 64-bit division by 18 - Pass - average iteration duration: 0.698 ns
                Signed 64-bit modulus by 18 - Pass - average iteration duration: 6.636 ns
               Signed 64-bit division by 24 - Pass - average iteration duration: 0.698 ns
                Signed 64-bit modulus by 24 - Pass - average iteration duration: 6.461 ns
              Signed 64-bit division by 100 - Pass - average iteration duration: 0.757 ns
               Signed 64-bit modulus by 100 - Pass - average iteration duration: 6.519 ns
              Signed 64-bit division by 153 - Pass - average iteration duration: 0.640 ns
               Signed 64-bit modulus by 153 - Pass - average iteration duration: 6.577 ns
            Signed 64-bit division by 1,461 - Pass - average iteration duration: 0.698 ns
             Signed 64-bit modulus by 1,461 - Pass - average iteration duration: 6.694 ns
Signed 64-bit division by 10,000 (Currency) - Pass - average iteration duration: 0.698 ns
 Signed 64-bit modulus by 10,000 (Currency) - Pass - average iteration duration: 6.694 ns
       Signed 64-bit division by 86,400,000 - Pass - average iteration duration: 0.698 ns
        Signed 64-bit modulus by 86,400,000 - Pass - average iteration duration: 6.577 ns

ok
- Sum of average durations: 93.540 ns
- Overall average duration: 1.336 ns

Very interesting, Gareth !

Is the div-by-7 related to 2 to the 3rd ? If yes, is it possible to design a div-by-3 with similar magics ?

Cheers, Ched'

Le 03.09.21 à 15:35, J. Gareth Moreton via fpc-devel a écrit :

Div by 3 has a very elegant implementation.  In 32-bit, it's

movl    $0xAAAAAAAB,%eax
mull    Input ; answer is in %edx:%eax
shrl    $1,%edx
movl    %edx, Result

As Marģers discovered, you can reduce the cycle count and minimise pipeline stalls by extending to 64-bit and applying the shift to the reciprocal, since the smaller input domain permits a larger rounding error before it starts affecting the result:

movq    $0x5555555555555556, %rax
mulq    Input ; answer is in %rdx:%rax.  Since input is a 32-bit register, the upper 32 bits of its 64-bit counterpart are zero. movl    %edx, Result

Now, 7 (in fact, any value of the form 7 × 2^k, where k ≥ 0) is a bit different because there's no magic number and shift value that will work for the entire 32-bit range; you need an extra bit of precision.  To simulate this, the input is added to the high word of the result, as if there were a 1 in the 33rd bit of the reciprocal.

movl    $0x24924925,%eax
mull    input
addl    input,%edx
rcrl    $1,%edx
shrl    $2,%edx
movl    %edx, Result

If the add operation overflows, you need to use rcrl (rotate with carry) to cycle the carry bit back into the result as you shift right, then you can shift right normally to get the result you need.

In 64-bit, this problem vanishes because you can insert that extra bit of precision as you shift the magic number - this sequence is about twice as fast, not just due to the reduced number of operations, but the reduction in pipeline stalls from each operation waiting on the previous one to complete (i.e. addl input,%edx cannot execute until the multiplication is finished, rcrl $1,%edx has to wait on addl, shrl $2,%edx has to wait on rcrl, and movl %edx, Result has to wait on shrl):

movq    0x2492492492492493
mulq   input
movl    %edx, Result

Modern x86 processors do out-of-order execution and can delay "movl %edx, Result" until it's actually needed, executing any instructions after it in parallel while the multiplication completes.  Unfortunately I received some test failures in tmoddiv4 and tmoddiv6 so I'm not out of the woods yet, and I'll have to see what's going on with those before I can continue.

Gareth aka. Kit

Fixed the problem.  I was deallocating a register too soon, so it got overwritten in some cases while still in use.  Running tests again.

Gareth aka. Kit

This one is for Marģers especially:

https://gitlab.com/freepascal.org/fpc/source/-/issues/39355

Gareth aka. Kit

The problem is that Free Pascal automatically scales up inputs to the word size, usually 32-bits.  If you do a node dump of a division node, you'll see that the input will often have a typecast node cast upon it.  It is possible to remove it, but it would have to be done at a different compiler stage.  I might see if it's possible to implement this at a later stage, since I've done just that for my latest merge request for "(x mod n) = 0".

Gareth aka. Kit

This one for Marģers specifically,

You'll be pleased to know that your insight has been partially implemented!

https://gitlab.com/freepascal.org/fpc/source/-/merge_requests/51

This only expands 32-bit divisions to 64-bit, since the smaller sizes requires more work at the node level, but is certainly possible and will happily continue to experiment with seeing if it can be implemented.  This was fun to do!

The test "bench/bdiv.pp" (also "tests/test/cg/tmoddiv6.pp" which just includes the given file) is my benchmark test for division and modulus operations if anyone wants to see if they can make the fundamental routine faster on any platform.

Gareth aka. Kit

Hi,

in my big decimal unit, I need to divide by 1000000000 all the time, e.g.
https://github.com/benibela/bigdecimalmath/blob/master/bigdecimalmath.pas#L1324-L1325

would the magic div help there much?

Bye,
Benito

If the code breaks down to "x div 1000000000" (with 1000000000 as an explicit constant and x being a 64-bit integer), then yes it will.  If the divisor is a variable though, then it won't.  However, the "magic div" algorithm only works for numerators up to a given bit size; if the number exceeds it, you'll lose precision and risk getting an incorrect answer (a warning for if you want to try to implement the algorithm yourself).

Gareth aka. Kit