[XiaohongGong]

zhi, the loop-carried borrow in Knuth Algorithm D, was stored as an sdigit (32-bit int) and then widened to 64-biy stwodigits on the next iteration. That extra conversion sat on the critical chain which is strong latency bounded:

    ldr     w0, [x12, x2, lsl #2]
    ldr     w1, [x27, x2, lsl #2]
    add     w0, w3, w0              // critical chain
    sxtw    x0, w0                  // critical chain (32 -> 64)
    umsubl  x0, w1, w14, x0         // critical chain
    and     w1, w0, #0x3fffffff
    str     w1, [x12, x2, lsl #2]
    asr     x0, x0, #30             // critical chain
    mov     w3, w0                  // critical chain (64 -> 32)

Store zhi as stwodigits instead to save the type conversions. In addition, the compiler can then compute vk[i] - q*w0[i] at first
and only then add the carried borrow, so this computation no longer depends on zhi:

    ldr     w4, [x3, x1, lsl #2]
    ldr     w0, [x28, x1, lsl #2]
    umsubl  x0, w0, w13, x4
    add     x2, x0, x2              // critical chain
    and     w4, w2, #0x3fffffff
    str     w4, [x3, x1, lsl #2]
    asr     x2, x2, #30             // critical chain

The loop-carried chain shrinks to two instructions. Knuth invariants are unchanged, so quotients and remainders are the same.

Below is the performance data from a micro python benchmark and the pidigits on my AArch64 and X86 machines respecitively. I didn't observe obvious performance regressions for other pyperformance benchmarks with this change.

AArch64 (NVIDIA Grace，Neoverse V):

| Benchmark      | base               | with this change                 |
|----------------|:------------------:|:--------------------------------:|
| div 500/70     | 50.1 us +- 0.1 us  | 22.4 us +- 0.0 us: 2.23x faster  |
| div 200/60     | 14.1 us +- 0.0 us  | 6.49 us +- 0.03 us: 2.18x faster |
| div 200/10     | 4.23 us +- 0.01 us | 2.59 us +- 0.01 us: 1.63x faster |
| div 60/10      | 1.16 us +- 0.00 us | 727 ns +- 6 ns: 1.59x faster     |
| div 61/60      | 262 ns +- 1 ns     | 207 ns +- 5 ns: 1.27x faster     |
| div 11/10      | 107 ns +- 0 ns     | 90.3 ns +- 0.9 ns: 1.18x faster  |
| div 60/2       | 606 ns +- 2 ns     | 541 ns +- 3 ns: 1.12x faster     |
| div 10/2       | 137 ns +- 4 ns     | 129 ns +- 4 ns: 1.06x faster     |
| div 3/2        | 73.3 ns +- 1.2 ns  | 72.4 ns +- 1.1 ns: 1.01x faster  |
| Geometric mean | (ref)              | 1.42x faster                     |
|----------------|--------------------|----------------------------------|
| pidigits       | 172 ms +- 0 ms     | 153 ms +- 0 ms: 1.13x faster     |

X86_64 (Intel Xeon Gold 6444Y):

| Benchmark      | base               | with this change                 |
|----------------|:------------------:|:--------------------------------:|
| div 500/70     | 36.9 us +- 0.0 us  | 22.5 us +- 0.0 us: 1.64x faster  |
| div 200/60     | 10.7 us +- 0.0 us  | 6.61 us +- 0.01 us: 1.62x faster |
| div 60/10      | 939 ns +- 1 ns     | 744 ns +- 1 ns: 1.26x faster     |
| div 200/10     | 3.48 us +- 0.01 us | 2.79 us +- 0.01 us: 1.25x faster |
| div 61/60      | 216 ns +- 0 ns     | 188 ns +- 1 ns: 1.15x faster     |
| div 11/10      | 88.3 ns +- 0.3 ns  | 79.3 ns +- 0.3 ns: 1.11x faster  |
| div 10/2       | 119 ns +- 1 ns     | 117 ns +- 0 ns: 1.01x faster     |
| div 3/2        | 56.0 ns +- 0.4 ns  | 55.2 ns +- 0.7 ns: 1.01x faster  |
| div 60/2       | 563 ns +- 1 ns     | 557 ns +- 2 ns: 1.01x faster     |
| Geometric mean | (ref)              | 1.21x faster                     |
|----------------|--------------------|----------------------------------|
| pidigits       | 147 ms +- 1 ms     | 134 ms +- 0 ms: 1.09x faster     |

Micro-benchmark:

import pyperf

BITS = 30

runner = pyperf.Runner()
CASES = [
    ("div 3/2", 3, 2),
    ("div 10/2", 10, 2),
    ("div 11/10", 11, 10),
    ("div 60/2", 60, 2),
    ("div 60/10", 60, 10),
    ("div 61/60", 61, 60),
    ("div 200/10", 200, 10),
    ("div 200/60", 200, 60),
    ("div 500/70", 500, 70),
]
for name, na, nb in CASES:
    setup = (
        f"import random; r = random.Random(1234); "
        f"a = r.getrandbits({na * BITS - 1}) | (1 << {na * BITS - 2}); "
        f"b = r.getrandbits({nb * BITS - 1}) | (1 << {nb * BITS - 2})"
    )
    runner.timeit(name, stmt="divmod(a, b)", setup=setup)

Fixes gh-160872.

Co-authored-by: Kyrylo Tkachov ktkachov@nvidia.com

[bedevere-app]

Most changes to Python require a NEWS entry. Add one using the blurb_it web app or the blurb command-line tool.

If this change has little impact on Python users, wait for a maintainer to apply the skip news label instead.

[XiaohongGong]

Hi @eendebakpt , thanks a lot for your review!

About integrating this PR, do I need to waiting for others to review? If not, could anyone please help merge this PR?

[eendebakpt]

Nothing more required, just some patience. I will give other devs some time to have a look as well, if no further comments I will merge.

[XiaohongGong]

Nothing more required, just some patience. I will give other devs some time to have a look as well, if no further comments I will merge.

Good to know this. Thanks a lot for your inputs and help!