Как разделить 128-битное делимое на 64-битный делитель, где все биты делимого равны 1, а мне нужны только 64 младших бита частного?
Мне нужно посчитать
((2^128) - 1) / x. Делитель, представляет собой 64-битное число без знака. Делимое состоит из двух беззнаковых 64-битных чисел (старшего и младшего), где оба числа
UINT64_MAX. Я могу использовать только 64-битную арифметику, и мне нужно, чтобы она была переносимой (без использования GNU
__int128, MSCV's
_udiv128, сборка или что-то в этом роде). Мне не нужна высокая часть частного, мне нужны только младшие 64 бита.
Как мне проделать эту операцию?
Также:
x >= 3,
x не является степенью двойки.
2 ответа
Мне неизвестны какие-либо оптимизации, применимые к целочисленному делению с постоянным делимым. Чтобы перепроверить, я попробовал тестовый пример с дивидендом «все единицы» с помощью Compiler Explorer. При использовании gcc, icc и clang с заданным наивысшим уровнем оптимизации сгенерированный код показал, что оптимизация не применяется к подразделению.
Конечно, можно построить высокопроизводительные 128-битные подпрограммы деления, но из личного опыта я знаю, что это довольно подвержено ошибкам, и для достижения хорошего тестового покрытия, включая крайние случаи, требуется очень сложное тестирование, поскольку исчерпывающее тестирование невозможно при этот размер операнда. Усилия по проектированию и тестированию легко превышают то, что кажется разумным для ответа на Stackoverflow, на два десятичных порядка.
Простой способ выполнить целочисленное деление — использовать алгоритм, который мы все выучили в начальной школе, только в двоичном формате. Это делает решение о следующем бите частного особенно простым: он равен 1, когда текущий неполный остаток больше или равен делителю, и 0 в противном случае. При использовании обычного двоичного деления единственные операции с целыми числами, которые нам нужны, — это сложение и вычитание.
Мы можем создавать переносимые примитивы для выполнения этих действий над операндами любой битовой длины, имитируя способ, которым машинные инструкции процессора используются для выполнения операций над целыми числами, состоящими из нескольких слов: ADD с переносом, ADD с переносом, ADD с переносом. и проведение; аналогично для СУБ. В приведенном ниже коде я использую для этого простые макросы C; конечно, возможны более изощренные подходы.
Поскольку система, над которой я сейчас работаю, не поддерживает 128-битные целые числа, я создал прототип и протестировал этот подход для 64-битных целых чисел. Тогда 128-битная версия была упражнением в простом механическом переименовании. Я ожидаю, что на современном 64-битном процессоре эта 128-битная функция деления будет выполняться примерно за 3000 циклов.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>
#define SUBCcc(a,b,cy,t0,t1,t2) \
(t0=(b)+cy, t1=(a), cy=t0<cy, t2=t1<t0, cy=cy+t2, t1-t0)
#define SUBcc(a,b,cy,t0,t1) \
(t0=(b), t1=(a), cy=t1<t0, t1-t0)
#define SUBC(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), t1-t0)
#define ADDCcc(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), cy=t0<cy, t0=t0+t1, t1=t0<t1, cy=cy+t1, t0=t0)
#define ADDcc(a,b,cy,t0,t1) \
(t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)
#define ADDC(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), t0+t1)
typedef struct {
uint64_t l;
uint64_t h;
} my_uint128;
my_uint128 bitwise_division_128 (my_uint128 dvnd, my_uint128 dvsr)
{
my_uint128 quot, rem, tmp;
uint64_t cy, t0, t1, t2;
int bits_left = CHAR_BIT * sizeof (my_uint128);
quot.h = dvnd.h;
quot.l = dvnd.l;
rem.h = 0;
rem.l = 0;
do {
quot.l = ADDcc (quot.l, quot.l, cy, t0, t1);
quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
rem.l = ADDCcc (rem.l, rem.l, cy, t0, t1);
rem.h = ADDC (rem.h, rem.h, cy, t0, t1);
tmp.l = SUBcc (rem.l, dvsr.l, cy, t0, t1);
tmp.h = SUBCcc (rem.h, dvsr.h, cy, t0, t1, t2);
if (!cy) { // remainder >= divisor
rem.l = tmp.l;
rem.h = tmp.h;
quot.l = quot.l | 1;
}
bits_left--;
} while (bits_left);
return quot;
}
typedef struct {
uint32_t l;
uint32_t h;
} my_uint64;
my_uint64 bitwise_division_64 (my_uint64 dvnd, my_uint64 dvsr)
{
my_uint64 quot, rem, tmp;
uint32_t cy, t0, t1, t2;
int bits_left = CHAR_BIT * sizeof (my_uint64);
quot.h = dvnd.h;
quot.l = dvnd.l;
rem.h = 0;
rem.l = 0;
do {
quot.l = ADDcc (quot.l, quot.l, cy, t0, t1);
quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
rem.l = ADDCcc (rem.l, rem.l, cy, t0, t1);
rem.h = ADDC (rem.h, rem.h, cy, t0, t1);
tmp.l = SUBcc (rem.l, dvsr.l, cy, t0, t1);
tmp.h = SUBCcc (rem.h, dvsr.h, cy, t0, t1, t2);
if (!cy) { // remainder >= divisor
rem.l = tmp.l;
rem.h = tmp.h;
quot.l = quot.l | 1;
}
bits_left--;
} while (bits_left);
return quot;
}
/*
https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
From: geo <gmars...@gmail.com>
Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
Subject: 64-bit KISS RNGs
Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)
This 64-bit KISS RNG has three components, each nearly
good enough to serve alone. The components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64 (kiss64_t = (kiss64_x << 58) + kiss64_c, \
kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64 (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
kiss64_y ^= (kiss64_y << 43))
#define CNG64 (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
int main (void)
{
uint64_t a, b, res, ref;
my_uint64 aa, bb, rr;
do {
a = KISS64;
b = KISS64;
ref = a / b;
aa.l = (uint32_t)a;
aa.h = (uint32_t)(a >> 32);
bb.l = (uint32_t)b;
bb.h = (uint32_t)(b >> 32);
rr = bitwise_division_64 (aa, bb);
res = (((uint64_t)rr.h) << 32) + rr.l;
if (ref != res) {
printf ("a=%016llx b=%016llx res=%016llx ref=%016llx\n", a, b, res, ref);
return EXIT_FAILURE;
}
} while (a);
return EXIT_SUCCESS;
}
Более быстрый подход, чем побитовое вычисление, заключается в вычислении обратной величины делителя, умножении на делимое, в результате чего получается предварительное частное, а затем вычисление остатка для точной настройки частного. Все вычисления могут быть выполнены в арифметике с фиксированной точкой. Однако на современных процессорах с быстрыми модулями с плавающей запятой удобнее генерировать начальное приближение для обратной величины с делением с двойной точностью. Затем одна итерация Галлея с кубической сходимостью приводит к обратной величине с полной точностью.
Итерация Галлея для обратного умножения очень интенсивно использует целочисленное умножение с 64x64-битным умножением со 128-битным результатом (
umul64wide()в приведенном ниже коде), являясь строительным блоком, имеющим решающее значение для производительности. В современных 64-битных архитектурах это обычно одна машинная инструкция, выполняемая за несколько циклов, однако она недоступна для переносимого кода. Переносимый код, эмулирующий инструкцию, требует от 15 до 20 инструкций в зависимости от архитектуры и компилятора.
Все 128-битное деление должно занять примерно 300 циклов, или в десять раз быстрее, чем простое побитовое вычисление. Поскольку код довольно сложен, для обеспечения правильной работы требуется значительное количество тестов. В приведенной ниже структуре я использую тесты на основе шаблонов и случайные тесты для умеренно интенсивного тестирования, используя простую побитовую реализацию в качестве эталона.
Реализация
udiv128()ниже предполагается, что среда программирования использует совместимую с IEEE-754 арифметику с плавающей запятой, что тип сопоставлен с IEEE-754
binary64типа, и это подразделение
doubleоперанды правильно округлены.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <limits.h>
typedef struct {
uint64_t l;
uint64_t h;
} my_uint128;
my_uint128 make_my_uint128 (uint64_t h, uint64_t l);
my_uint128 add128 (my_uint128 a, my_uint128 b);
my_uint128 sub128 (my_uint128 a, my_uint128 b);
my_uint128 lsl128 (my_uint128 a, int sh);
my_uint128 lsr128 (my_uint128 a, int sh);
my_uint128 not128 (my_uint128 a);
my_uint128 umul128lo (my_uint128 a, my_uint128 b);
my_uint128 umul128hi (my_uint128 a, my_uint128 b);
double my_uint128_to_double (my_uint128 a);
int lt128 (my_uint128 a, my_uint128 b);
int eq128 (my_uint128 a, my_uint128 b);
uint64_t double_as_uint64 (double a);
double uint64_as_double (uint64_t a);
#define FP64_EXPO_BIAS (1023)
#define FP64_MANT_BITS (53)
#define FP64_MANT_IBIT (0x0010000000000000ULL)
#define FP64_MANT_MASK (0x000fffffffffffffULL)
#define FP64_INC_EXP_128 (0x0800000000000000ULL)
#define FP64_MANT_ADJ (2) // adjustment to ensure underestimate
my_uint128 udiv128 (my_uint128 dividend, my_uint128 divisor)
{
const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
const my_uint128 one = make_my_uint128 (0ULL, 1ULL);
const my_uint128 two = make_my_uint128 (0ULL, 2ULL);
my_uint128 recip, temp, quo, rem;
my_uint128 neg_divisor = sub128 (zero, divisor);
double r;
/* compute initial approximation for reciprocal; must be underestimate! */
r = 1.0 / my_uint128_to_double (divisor);
uint64_t i = double_as_uint64 (r) - FP64_MANT_ADJ + FP64_INC_EXP_128;
temp = make_my_uint128 (0ULL, (i & FP64_MANT_MASK) | FP64_MANT_IBIT);
int sh = (i >> (FP64_MANT_BITS-1)) - FP64_EXPO_BIAS - (FP64_MANT_BITS-1);
recip = (sh < 0) ? lsr128 (temp, -sh) : lsl128 (temp, sh);
/* perform Halley iteration with cubic convergence to refine reciprocal */
temp = umul128lo (neg_divisor, recip);
temp = add128 (umul128hi (temp, temp), temp);
recip = add128 (umul128hi (recip, temp), recip);
/* compute preliminary quotient and remainder */
quo = umul128hi (dividend, recip);
rem = sub128 (dividend, umul128lo (divisor, quo));
/* adjust quotient if too small; quotient off by 2 at most */
if (! lt128 (rem, divisor)) {
quo = add128 (quo, lt128 (sub128 (rem, divisor), divisor) ? one : two);
}
/* handle division by zero */
if (eq128 (divisor, zero)) quo = not128 (zero);
return quo;
}
#define SUBCcc(a,b,cy,t0,t1,t2) \
(t0=(b)+cy, t1=(a), cy=t0<cy, t2=t1<t0, cy=cy+t2, t1-t0)
#define SUBcc(a,b,cy,t0,t1) \
(t0=(b), t1=(a), cy=t1<t0, t1-t0)
#define SUBC(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), t1-t0)
#define ADDCcc(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), cy=t0<cy, t0=t0+t1, t1=t0<t1, cy=cy+t1, t0=t0)
#define ADDcc(a,b,cy,t0,t1) \
(t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)
#define ADDC(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), t0+t1)
uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}
double uint64_as_double (uint64_t a)
{
double r;
memcpy (&r, &a, sizeof r);
return r;
}
my_uint128 add128 (my_uint128 a, my_uint128 b)
{
uint64_t cy, t0, t1;
a.l = ADDcc (a.l, b.l, cy, t0, t1);
a.h = ADDC (a.h, b.h, cy, t0, t1);
return a;
}
my_uint128 sub128 (my_uint128 a, my_uint128 b)
{
uint64_t cy, t0, t1;
a.l = SUBcc (a.l, b.l, cy, t0, t1);
a.h = SUBC (a.h, b.h, cy, t0, t1);
return a;
}
my_uint128 lsl128 (my_uint128 a, int sh)
{
if (sh >= 64) {
a.h = a.l << (sh - 64);
a.l = 0ULL;
} else if (sh) {
a.h = (a.h << sh) + (a.l >> (64 - sh));
a.l = a.l << sh;
}
return a;
}
my_uint128 lsr128 (my_uint128 a, int sh)
{
if (sh >= 64) {
a.l = a.h >> (sh - 64);
a.h = 0ULL;
} else if (sh) {
a.l = (a.l >> sh) + (a.h << (64 - sh));
a.h = a.h >> sh;
}
return a;
}
my_uint128 not128 (my_uint128 a)
{
a.l = ~a.l;
a.h = ~a.h;
return a;
}
int lt128 (my_uint128 a, my_uint128 b)
{
uint64_t cy, t0, t1, t2;
a.l = SUBcc (a.l, b.l, cy, t0, t1);
a.h = SUBCcc (a.h, b.h, cy, t0, t1, t2);
return cy;
}
int eq128 (my_uint128 a, my_uint128 b)
{
return (a.l == b.l) && (a.h == b.h);
}
// derived from Hacker's Delight 2nd ed. figure 8-2
my_uint128 umul64wide (uint64_t u, uint64_t v)
{
my_uint128 r;
uint64_t u0, v0, u1, v1, w0, w1, w2, t;
u0 = (uint32_t)u; u1 = u >> 32;
v0 = (uint32_t)v; v1 = v >> 32;
w0 = u0 * v0;
t = u1 * v0 + (w0 >> 32);
w1 = (uint32_t)t;
w2 = t >> 32;
w1 = u0 * v1 + w1;
r.h = u1 * v1 + w2 + (w1 >> 32);
r.l = (w1 << 32) + (uint32_t)w0;
return r;
}
my_uint128 make_my_uint128 (uint64_t h, uint64_t l)
{
my_uint128 r;
r.h = h;
r.l = l;
return r;
}
my_uint128 umul128lo (my_uint128 a, my_uint128 b)
{
my_uint128 r;
r = umul64wide (a.l, b.l);
r.h = r.h + a.l * b.h + a.h * b.l;
return r;
}
my_uint128 umul128hi (my_uint128 a, my_uint128 b)
{
my_uint128 t0, t1, t2, t3;
t0 = umul64wide (a.l, b.l);
t3 = add128 (umul64wide (a.h, b.l), make_my_uint128 (0ULL, t0.h));
t1 = make_my_uint128 (0ULL, t3.l);
t2 = make_my_uint128 (0ULL, t3.h);
t1 = add128 (umul64wide (a.l, b.h), t1);
return add128 (add128 (umul64wide (a.h, b.h), t2), make_my_uint128 (0ULL, t1.h));
}
double my_uint128_to_double (my_uint128 a)
{
const int intbits = sizeof (a) * CHAR_BIT;
const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
my_uint128 rnd, i = a;
uint64_t j;
int sh = 0;
double r;
// normalize integer so MSB is set
if (lt128 (i, make_my_uint128(0x0000000000000001ULL, 0))) {i = lsl128 (i,64); sh += 64; }
if (lt128 (i, make_my_uint128(0x0000000100000000ULL, 0))) {i = lsl128 (i,32); sh += 32; }
if (lt128 (i, make_my_uint128(0x0001000000000000ULL, 0))) {i = lsl128 (i,16); sh += 16; }
if (lt128 (i, make_my_uint128(0x0100000000000000ULL, 0))) {i = lsl128 (i, 8); sh += 8; }
if (lt128 (i, make_my_uint128(0x1000000000000000ULL, 0))) {i = lsl128 (i, 4); sh += 4; }
if (lt128 (i, make_my_uint128(0x4000000000000000ULL, 0))) {i = lsl128 (i, 2); sh += 2; }
if (lt128 (i, make_my_uint128(0x8000000000000000ULL, 0))) {i = lsl128 (i, 1); sh += 1; }
// form mantissa with explicit integer bit
rnd = lsl128 (i, FP64_MANT_BITS);
i = lsr128 (i, intbits - FP64_MANT_BITS);
j = i.l;
// add in exponent, taking into account integer bit of mantissa
if (! eq128 (a, zero)) {
j += (uint64_t)(FP64_EXPO_BIAS + (intbits-1) - 1 - sh) << (FP64_MANT_BITS-1);
}
// round to nearest or even
rnd.h = rnd.h | (rnd.l != 0);
if ((rnd.h > 0x8000000000000000ULL) ||
((rnd.h == 0x8000000000000000ULL) && (j & 1))) j++;
// reinterpret bit pattern as IEEE-754 'binary64'
r = uint64_as_double (j);
return r;
}
my_uint128 bitwise_division_128 (my_uint128 dvnd, my_uint128 dvsr)
{
my_uint128 quot, rem, tmp;
uint64_t cy, t0, t1, t2;
int bits_left = CHAR_BIT * sizeof (dvsr);
quot.h = dvnd.h;
quot.l = dvnd.l;
rem.h = 0;
rem.l = 0;
do {
quot.l = ADDcc (quot.l, quot.l, cy, t0, t1);
quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
rem.l = ADDCcc (rem.l, rem.l, cy, t0, t1);
rem.h = ADDC (rem.h, rem.h, cy, t0, t1);
tmp.l = SUBcc (rem.l, dvsr.l, cy, t0, t1);
tmp.h = SUBCcc (rem.h, dvsr.h, cy, t0, t1, t2);
if (!cy) { // remainder >= divisor
rem.l = tmp.l;
rem.h = tmp.h;
quot.l = quot.l | 1;
}
bits_left--;
} while (bits_left);
return quot;
}
/*
https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
From: geo <gmars...@gmail.com>
Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
Subject: 64-bit KISS RNGs
Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)
This 64-bit KISS RNG has three components, each nearly
good enough to serve alone. The components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64 (kiss64_t = (kiss64_x << 58) + kiss64_c, \
kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64 (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
kiss64_y ^= (kiss64_y << 43))
#define CNG64 (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
my_uint128 v[100000]; /* FIXME: size appropriately */
int main (void)
{
const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
const my_uint128 one = make_my_uint128 (0ULL, 1ULL);
my_uint128 dividend, divisor, quot, ref;
int i, j, patterns, idx = 0, nbrBits = sizeof (v[0]) * CHAR_BIT;
int patterns_done = 0;
/* pattern class 1: 2**i */
for (i = 0; i < nbrBits; i++) {
v [idx] = lsl128 (one, i);
idx++;
}
/* pattern class 2: 2**i-1 */
for (i = 0; i < nbrBits; i++) {
v [idx] = sub128 (lsl128 (one, i), one);
idx++;
}
/* pattern class 3: 2**i+1 */
for (i = 0; i < nbrBits; i++) {
v [idx] = add128 (lsl128 (one, i), one);
idx++;
}
/* pattern class 4: 2**i + 2**j */
for (i = 0; i < nbrBits; i++) {
for (j = 0; j < nbrBits; j++) {
v [idx] = add128 (lsl128 (one, i), lsl128 (one, j));
idx++;
}
}
/* pattern class 5: 2**i - 2**j */
for (i = 0; i < nbrBits; i++) {
for (j = 0; j < nbrBits; j++) {
v [idx] = sub128 (lsl128 (one, i), lsl128 (one, j));
idx++;
}
}
patterns = idx;
/* pattern class 6: one's complement of pattern classes 1 through 5 */
for (i = 0; i < patterns; i++) {
v [idx] = not128 (v [i]);
idx++;
}
/* pattern class 7: two's complement of pattern classes 1 through 5 */
for (i = 0; i < patterns; i++) {
v [idx] = sub128 (zero, v[i]);
idx++;
}
patterns = idx;
printf ("Starting pattern-based tests. Number of patterns: %d\n", patterns);
for (long long int k = 0; k < 100000000000LL; k++) {
if (k < patterns * patterns) {
dividend = v [k / patterns];
divisor = v [k % patterns];
} else {
if (!patterns_done) {
printf ("Starting random tests\n");
patterns_done = 1;
}
dividend.l = KISS64;
dividend.h = KISS64;
divisor.h = KISS64;
divisor.l = KISS64;
}
/* exclude cases with undefined results: division by zero */
if (! eq128 (divisor, zero)) {
quot = udiv128 (dividend, divisor);
ref = bitwise_division_128 (dividend, divisor);
if (! eq128 (quot, ref)) {
printf ("@ (%016llx_%016llx, %016llx_%016llx): quot = %016llx_%016llx ref=%016llx_%016llx\n",
dividend.h, dividend.l, divisor.h, divisor.l,
quot.h, quot.l, ref.h, ref.l);
return EXIT_FAILURE;
}
}
}
printf ("unsigned 128-bit division: tests passed\n");
return EXIT_SUCCESS;
}
Это то, что я закончил кодировать. Я уверен, что есть гораздо более быстрые альтернативы, но, по крайней мере, это функционально.
На основе: https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(unsigned)_with_remainder . Адаптировано для этого конкретного варианта использования.
// q = (2^128 - 1) / d, where q is the 64 LSBs of the quotient
uint64_t two_pow_128_minus_1_div_d(uint64_t d) {
uint64_t q = 0, r_hi = 0, r_lo = 0;
for (int i = 127; i >= 0; --i) {
r_hi = (r_hi << 1) | (r_lo >> 63);
r_lo <<= 1;
r_lo |= 1UL;
if (r_hi || r_lo >= d) {
const uint64_t borrow = d > r_lo;
r_lo -= d;
r_hi -= borrow;
if (i < 64)
q |= 1UL << i;
}
}
return q;
}