[druntime] core.int128: Extract nested udivmod128_64 as udivmod overload

This commit is contained in:
Martin Kinkelin 2025-03-24 01:20:28 +01:00 committed by Nicholas Wilson
parent ea4b1e9607
commit 3a661c1f8c

View file

@ -485,10 +485,98 @@ Cent udivmod(Cent c1, Cent c2, out Cent modulus)
// Based on "Unsigned Doubleword Division" in Hacker's Delight
import core.bitop;
// Divides a 128-bit dividend by a 64-bit divisor.
// The result must fit in 64 bits.
static U udivmod128_64(Cent c1, U c2, out U modulus)
// Special cases
if (!tst(c2))
{
// Divide by zero
modulus = Zero;
return com(modulus);
}
if (c1.hi == 0 && c2.hi == 0)
{
// Single precision divide
const Cent rem = { lo:c1.lo % c2.lo };
modulus = rem;
const Cent ret = { lo:c1.lo / c2.lo };
return ret;
}
if (c1.hi == 0)
{
// Numerator is smaller than the divisor
modulus = c1;
return Zero;
}
if (c2.hi == 0)
{
// Divisor is a 64-bit value, so we just need one 128/64 division.
// If c1 / c2 would overflow, break c1 up into two halves.
const q1 = (c1.hi < c2.lo) ? 0 : (c1.hi / c2.lo);
if (q1)
c1.hi = c1.hi % c2.lo;
Cent rem;
const q0 = udivmod(c1, c2.lo, rem.lo);
modulus = rem;
const Cent ret = { lo:q0, hi:q1 };
return ret;
}
// Full cent precision division.
// Here c2 >= 2^^64
// We know that c2.hi != 0, so count leading zeros is OK
// We have 0 <= shift <= 63
const shift = (Ubits - 1) - bsr(c2.hi);
// Normalize the divisor so its MSB is 1
// v1 = (c2 << shift) >> 64
U v1 = shl(c2, shift).hi;
// To ensure no overflow.
Cent u1 = shr1(c1);
// Get quotient from divide unsigned operation.
U rem_ignored;
const Cent q1 = { lo:udivmod(u1, v1, rem_ignored) };
// Undo normalization and division of c1 by 2.
Cent quotient = shr(shl(q1, shift), 63);
// Make quotient correct or too small by 1
if (tst(quotient))
quotient = dec(quotient);
// Now quotient is correct.
// Compute rem = c1 - (quotient * c2);
Cent rem = sub(c1, mul(quotient, c2));
// Check if remainder is larger than the divisor
if (uge(rem, c2))
{
// Increment quotient
quotient = inc(quotient);
// Subtract c2 from remainder
rem = sub(rem, c2);
}
modulus = rem;
//printf("quotient "); print(quotient);
//printf("modulus "); print(modulus);
return quotient;
}
/****************************
* Unsigned divide 128-bit c1 / 64-bit c2. The result must fit in 64 bits.
* The remainder after division is stored to modulus.
* Params:
* c1 = dividend
* c2 = divisor
* modulus = set to c1 % c2
* Returns:
* quotient c1 / c2
*/
pure
U udivmod(Cent c1, U c2, out U modulus)
{
import core.bitop;
// We work in base 2^^32
enum base = 1UL << 32;
enum divmask = (1UL << (Ubits / 2)) - 1;
@ -540,83 +628,6 @@ Cent udivmod(Cent c1, Cent c2, out Cent modulus)
modulus = (rem * base + num0 - q0 * c2) >> shift;
return (cast(U)q1 << divshift) | q0;
}
// Special cases
if (!tst(c2))
{
// Divide by zero
modulus = Zero;
return com(modulus);
}
if (c1.hi == 0 && c2.hi == 0)
{
// Single precision divide
const Cent rem = { lo:c1.lo % c2.lo };
modulus = rem;
const Cent ret = { lo:c1.lo / c2.lo };
return ret;
}
if (c1.hi == 0)
{
// Numerator is smaller than the divisor
modulus = c1;
return Zero;
}
if (c2.hi == 0)
{
// Divisor is a 64-bit value, so we just need one 128/64 division.
// If c1 / c2 would overflow, break c1 up into two halves.
const q1 = (c1.hi < c2.lo) ? 0 : (c1.hi / c2.lo);
if (q1)
c1.hi = c1.hi % c2.lo;
Cent rem;
const q0 = udivmod128_64(c1, c2.lo, rem.lo);
modulus = rem;
const Cent ret = { lo:q0, hi:q1 };
return ret;
}
// Full cent precision division.
// Here c2 >= 2^^64
// We know that c2.hi != 0, so count leading zeros is OK
// We have 0 <= shift <= 63
const shift = (Ubits - 1) - bsr(c2.hi);
// Normalize the divisor so its MSB is 1
// v1 = (c2 << shift) >> 64
U v1 = shl(c2, shift).hi;
// To ensure no overflow.
Cent u1 = shr1(c1);
// Get quotient from divide unsigned operation.
U rem_ignored;
const Cent q1 = { lo:udivmod128_64(u1, v1, rem_ignored) };
// Undo normalization and division of c1 by 2.
Cent quotient = shr(shl(q1, shift), 63);
// Make quotient correct or too small by 1
if (tst(quotient))
quotient = dec(quotient);
// Now quotient is correct.
// Compute rem = c1 - (quotient * c2);
Cent rem = sub(c1, mul(quotient, c2));
// Check if remainder is larger than the divisor
if (uge(rem, c2))
{
// Increment quotient
quotient = inc(quotient);
// Subtract c2 from remainder
rem = sub(rem, c2);
}
modulus = rem;
//printf("quotient "); print(quotient);
//printf("modulus "); print(modulus);
return quotient;
}