|
| 1 | +// Unit tests for divlu and divllu. |
| 2 | +// To build standalone: |
| 3 | +// cc -O1 -o test_divlu test/test_divlu.c doc/divlu.c |
| 4 | + |
| 5 | +#include <stdint.h> |
| 6 | +#include <stdio.h> |
| 7 | +#include <stdlib.h> |
| 8 | + |
| 9 | +uint32_t divlu(uint32_t numhi, uint32_t numlo, uint32_t den, uint32_t *r); |
| 10 | +uint64_t divllu(uint64_t numhi, uint64_t numlo, uint64_t den, uint64_t *r); |
| 11 | + |
| 12 | +typedef unsigned long long ullong; |
| 13 | + |
| 14 | +static void check32( |
| 15 | + uint32_t numhi, uint32_t numlo, uint32_t den, uint32_t expected_quot, uint32_t expected_rem) { |
| 16 | + uint32_t rem; |
| 17 | + uint32_t quot = divlu(numhi, numlo, den, &rem); |
| 18 | + if (quot != expected_quot || rem != expected_rem) { |
| 19 | + fprintf(stderr, "divlu(%u, %u, %u): got q=%u r=%u, expected q=%u r=%u\n", numhi, numlo, den, |
| 20 | + quot, rem, expected_quot, expected_rem); |
| 21 | + abort(); |
| 22 | + } |
| 23 | +} |
| 24 | + |
| 25 | +static void check64( |
| 26 | + uint64_t numhi, uint64_t numlo, uint64_t den, uint64_t expected_quot, uint64_t expected_rem) { |
| 27 | + uint64_t rem; |
| 28 | + uint64_t quot = divllu(numhi, numlo, den, &rem); |
| 29 | + if (quot != expected_quot || rem != expected_rem) { |
| 30 | + fprintf(stderr, "divllu(%llu, %llu, %llu): got q=%llu r=%llu, expected q=%llu r=%llu\n", |
| 31 | + (ullong)numhi, (ullong)numlo, (ullong)den, (ullong)quot, (ullong)rem, |
| 32 | + (ullong)expected_quot, (ullong)expected_rem); |
| 33 | + abort(); |
| 34 | + } |
| 35 | +} |
| 36 | + |
| 37 | +// Compute expected result from reference arithmetic and check. |
| 38 | +static void verify32(uint32_t numhi, uint32_t numlo, uint32_t den) { |
| 39 | + uint32_t expected_quot, expected_rem; |
| 40 | + if (numhi >= den) { |
| 41 | + expected_quot = UINT32_MAX; |
| 42 | + expected_rem = UINT32_MAX; |
| 43 | + } else { |
| 44 | + uint64_t num = ((uint64_t)numhi << 32) | numlo; |
| 45 | + expected_quot = (uint32_t)(num / den); |
| 46 | + expected_rem = (uint32_t)(num % den); |
| 47 | + } |
| 48 | + check32(numhi, numlo, den, expected_quot, expected_rem); |
| 49 | +} |
| 50 | + |
| 51 | +// Compute hi:lo = a * b as a 128-bit product. |
| 52 | +static void mul64(uint64_t a, uint64_t b, uint64_t *phi, uint64_t *plo) { |
| 53 | + uint64_t alo = (uint32_t)a, ahi = a >> 32; |
| 54 | + uint64_t blo = (uint32_t)b, bhi = b >> 32; |
| 55 | + uint64_t ll = alo * blo; |
| 56 | + uint64_t lh = alo * bhi; |
| 57 | + uint64_t hl = ahi * blo; |
| 58 | + uint64_t hh = ahi * bhi; |
| 59 | + uint64_t mid = (ll >> 32) + (uint32_t)lh + (uint32_t)hl; |
| 60 | + *phi = hh + (lh >> 32) + (hl >> 32) + (mid >> 32); |
| 61 | + *plo = (mid << 32) | (uint32_t)ll; |
| 62 | +} |
| 63 | + |
| 64 | +// Rather than computing the expected value with 128-bit arithmetic, verify |
| 65 | +// the invariant: rem < den and quot * den + rem == numhi:numlo. |
| 66 | +static void verify64(uint64_t numhi, uint64_t numlo, uint64_t den) { |
| 67 | + uint64_t rem; |
| 68 | + uint64_t quot = divllu(numhi, numlo, den, &rem); |
| 69 | + if (numhi >= den) { |
| 70 | + if (quot != UINT64_MAX || rem != UINT64_MAX) { |
| 71 | + fprintf(stderr, "divllu(%llu, %llu, %llu): expected overflow, got q=%llu r=%llu\n", |
| 72 | + (ullong)numhi, (ullong)numlo, (ullong)den, (ullong)quot, (ullong)rem); |
| 73 | + abort(); |
| 74 | + } |
| 75 | + return; |
| 76 | + } |
| 77 | + if (rem >= den) { |
| 78 | + fprintf(stderr, "divllu(%llu, %llu, %llu): rem=%llu out of range (den=%llu)\n", |
| 79 | + (ullong)numhi, (ullong)numlo, (ullong)den, (ullong)rem, (ullong)den); |
| 80 | + abort(); |
| 81 | + } |
| 82 | + // Check quot * den + rem == numhi:numlo. |
| 83 | + uint64_t phi, plo; |
| 84 | + mul64(quot, den, &phi, &plo); |
| 85 | + plo += rem; |
| 86 | + if (plo < rem) phi++; |
| 87 | + if (phi != numhi || plo != numlo) { |
| 88 | + fprintf(stderr, "divllu(%llu, %llu, %llu): q=%llu r=%llu fails reconstruction\n", |
| 89 | + (ullong)numhi, (ullong)numlo, (ullong)den, (ullong)quot, (ullong)rem); |
| 90 | + abort(); |
| 91 | + } |
| 92 | +} |
| 93 | + |
| 94 | +int main(void) { |
| 95 | + // Trivial |
| 96 | + check32(0, 0, 1, 0, 0); |
| 97 | + check32(0, 1, 1, 1, 0); |
| 98 | + check32(0, 7, 3, 2, 1); |
| 99 | + // numlo == denom |
| 100 | + check32(0, UINT32_MAX, UINT32_MAX, 1, 0); |
| 101 | + // Divide across the 32-bit boundary: (1<<32) / 2 == 1<<31 |
| 102 | + check32(1, 0, 2, (uint32_t)1 << 31, 0); |
| 103 | + check32(1, 1, 2, (uint32_t)1 << 31, 1); |
| 104 | + // Near-maximum inputs |
| 105 | + check32(UINT32_MAX - 1, 0, UINT32_MAX, UINT32_MAX - 1, UINT32_MAX - 1); |
| 106 | + check32(UINT32_MAX - 1, UINT32_MAX, UINT32_MAX, UINT32_MAX, UINT32_MAX - 1); |
| 107 | + // Overflow: numhi >= den |
| 108 | + check32(1, 0, 1, UINT32_MAX, UINT32_MAX); |
| 109 | + check32(UINT32_MAX, 0, UINT32_MAX, UINT32_MAX, UINT32_MAX); |
| 110 | + // Divide by zero: numhi (0) >= den (0) |
| 111 | + check32(0, 1, 0, UINT32_MAX, UINT32_MAX); |
| 112 | + |
| 113 | + // --- divllu hardcoded cases --- |
| 114 | + |
| 115 | + check64(0, 0, 1, 0, 0); |
| 116 | + check64(0, 1, 1, 1, 0); |
| 117 | + check64(0, 5, 3, 1, 2); |
| 118 | + check64(0, UINT64_MAX, UINT64_MAX, 1, 0); |
| 119 | + // Divide across the 64-bit boundary |
| 120 | + check64(1, 0, 2, (uint64_t)1 << 63, 0); |
| 121 | + check64(1, 1, 2, (uint64_t)1 << 63, 1); |
| 122 | + // Near-maximum inputs |
| 123 | + check64(UINT64_MAX - 1, 0, UINT64_MAX, UINT64_MAX - 1, UINT64_MAX - 1); |
| 124 | + check64(UINT64_MAX - 1, UINT64_MAX, UINT64_MAX, UINT64_MAX, UINT64_MAX - 1); |
| 125 | + // Overflow |
| 126 | + check64(1, 0, 1, UINT64_MAX, UINT64_MAX); |
| 127 | + check64(UINT64_MAX, 0, UINT64_MAX, UINT64_MAX, UINT64_MAX); |
| 128 | + // Divide by zero |
| 129 | + check64(0, 1, 0, UINT64_MAX, UINT64_MAX); |
| 130 | + |
| 131 | + // NULL remainder pointer: quotient must still be correct. |
| 132 | + if (divlu(0, 7, 3, NULL) != 2) abort(); |
| 133 | + if (divllu(0, 5, 3, NULL) != 1) abort(); |
| 134 | + |
| 135 | + // --- Cases that trigger the qhat -= 1 correction --- |
| 136 | + // |
| 137 | + // The algorithm estimates the quotient digit qhat as (numhi / den1) where |
| 138 | + // den1 is the high half of the denominator. When den0 (the low half) is |
| 139 | + // large this estimate can be off by one. The correction fires when |
| 140 | + // qhat*den0 > rhat*base + num1. |
| 141 | + // |
| 142 | + // divlu (base 2^16): den=0x8000FFFF gives den1=0x8000, den0=0xFFFF. |
| 143 | + // With numhi=0x7FFFFFFF: qhat=0xFFFF, rhat=0x7FFF, |
| 144 | + // qhat*den0 = 0xFFFE0001 > rhat*base = 0x7FFF0000 → correction fires. |
| 145 | + verify32(0x7FFFFFFFu, 0, 0x8000FFFFu); |
| 146 | + verify32(0x7FFFFFFFu, 0xFFFF0000u, 0x8000FFFFu); |
| 147 | + // |
| 148 | + // divllu (base 2^32): den=0x80000000FFFFFFFF gives den1=0x80000000, den0=0xFFFFFFFF. |
| 149 | + // With numhi=0x7FFFFFFFFFFFFFFF: qhat=0xFFFFFFFF, rhat=0x7FFFFFFF, |
| 150 | + // qhat*den0 = 0xFFFFFFFE00000001 > rhat*base = 0x7FFFFFFF00000000 → correction fires. |
| 151 | + verify64(UINT64_C(0x7FFFFFFFFFFFFFFF), 0, UINT64_C(0x80000000FFFFFFFF)); |
| 152 | + verify64( |
| 153 | + UINT64_C(0x7FFFFFFFFFFFFFFF), UINT64_C(0xFFFFFFFF00000000), UINT64_C(0x80000000FFFFFFFF)); |
| 154 | + |
| 155 | + // --- Systematic sweeps --- |
| 156 | + |
| 157 | + // divlu: all combinations with den, numhi, numlo < 256 (~8M calls). |
| 158 | + for (uint32_t den = 1; den < 256; den++) |
| 159 | + for (uint32_t numhi = 0; numhi < den; numhi++) |
| 160 | + for (uint32_t numlo = 0; numlo < 256; numlo++) verify32(numhi, numlo, den); |
| 161 | + |
| 162 | + // divlu: medium denominators, probing numhi at 0, midpoint, and max (~353M calls). |
| 163 | + for (uint32_t den = 256; den < 0x70000; den++) |
| 164 | + for (uint32_t numlo = 0; numlo < 256; numlo++) { |
| 165 | + verify32(0, numlo, den); |
| 166 | + verify32(den / 2, numlo, den); |
| 167 | + verify32(den - 1, numlo, den); |
| 168 | + } |
| 169 | + |
| 170 | + // divlu: near UINT32_MAX |
| 171 | + for (uint32_t d = 0; d < 256; d++) { |
| 172 | + uint32_t den = UINT32_MAX - d; |
| 173 | + for (uint32_t numlo = 0; numlo < 256; numlo++) { |
| 174 | + verify32(den - 1, numlo, den); |
| 175 | + verify32(den - 1, UINT32_MAX - numlo, den); |
| 176 | + } |
| 177 | + } |
| 178 | + |
| 179 | + // divllu: all combinations with den, numhi, numlo < 64 |
| 180 | + for (uint64_t den = 1; den < 64; den++) |
| 181 | + for (uint64_t numhi = 0; numhi < den; numhi++) |
| 182 | + for (uint64_t numlo = 0; numlo < 64; numlo++) verify64(numhi, numlo, den); |
| 183 | + |
| 184 | + // divllu: medium denominators |
| 185 | + for (uint64_t den = 64; den < 0x70000; den++) |
| 186 | + for (uint64_t numlo = 0; numlo < 256; numlo++) { |
| 187 | + verify64(0, numlo, den); |
| 188 | + verify64(den / 2, numlo, den); |
| 189 | + verify64(den - 1, numlo, den); |
| 190 | + } |
| 191 | + |
| 192 | + return 0; |
| 193 | +} |
0 commit comments