Jak rozdělit 128-bit dividendy 64-bitový dělitel, kde dividendy jsou bity jsou 1, a tam, kde potřebuji pouze 64 LSBs kvocientu?

Question 1

Potřebuji vypočítat ((2^128) - 1) / x. Dělitel, xje unsigned 64-bitové číslo. Dividenda se skládá ze dvou unsigned 64-bit čísla (vysoké a nízké), kde obě čísla jsou UINT64_MAX. Mohu použít pouze 64-bit aritmetiky a potřebujete, aby byl přenosný (ne použití GNU __int128, MSCV je _udiv128, montáž, nebo něco takového). Nepotřebuji vysokou část kvocient, potřebuji pouze dolních 64 bitů.

Jak mohu udělat tuto operaci?

Také: x >= 3, x není mocninou 2.

Edit: vytvořil jsem si své vlastní řešení (odpověď níže). Ale vítám jakékoli jiné řešení, které provádí lépe :)

Question 2

Nejsem si vědom nějaké optimalizace, které se vztahují na celočíselné dělení s konstantní dividendy. Dvakrát zkontrolujte, zkoušel jsem testovací případ s all-ty dividendy s Compiler Explorer. Pomocí gcc, icc, a clang, s nejvyšší stupeň optimalizace uvedeno, vygenerovaný kód ukázala, žádné optimalizace se aplikuje do divize.

To je jistě možné vytvořit vysoce-výkon 128-bit rozdělení rutiny, ale z osobní zkušenosti vím, že je poměrně náchylný k chybám, a velmi sofistikované testování je potřeba k dosažení dobrého pokrytí testů včetně rohových případy, jako vyčerpávající test není možné v této operand velikost. Snaha o design a test snadno přesahuje to, co se zdá rozumné pro odpověď na Stackoverflow o dva desítkové řády.

Snadný způsob, jak provádět celočíselné dělení je použít algoritmus, který jsme se všichni naučili ve škole, pouze v binární. To je rozhodnutí o dalším kvocient trochu příliš snadné: je 1 když je aktuální dílčí zbytek je větší než, nebo rovno, dělitel, a 0 jinak. Pomocí ručně binární dělení, pouze celočíselné operace, které potřebujeme, jsou sčítání a odčítání.

Můžeme stavět přenosné primitiva pro provádění těchto na operandy z libovolné bitové délky napodobováním způsobu, jakým procesor počítače instrukce jsou použity pro efekt operace na multi-slovo celá čísla: PŘIDAT s carry-out, PŘIDEJTE se carry-in, ADD carry-in a carry-out; analogicky pro SUB. V kódu níže jsem pomocí jednoduché C makra pro to; určitě sofistikovanější přístupy jsou možné.

Protože systém jsem pracoval na teď není podpora pro 128-bitová celá čísla, já navrhnuté a testované tento přístup pro 64-bitové celá čísla. Na 128-bit verzi, pak bylo cvičení v jednoduché mechanické přejmenování. Na moderní 64-bitový procesor bych očekávat, že tato 128-bit divize funkci vykonávat ve zhruba 3000 cyklů.

#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 <[email protected]>
  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;
}

Rychlejší přístup, než je bit-wise výpočtu je výpočet vzájemné dělitel, vynásobte dividendy, což má za následek předběžné kvocient, pak vypočítat zbývající přesně nastavit kvocient. Celý výpočet lze provést ve fixed-point aritmetiky. Nicméně, na moderní procesory s rychle floating-point jednotky je výhodnější vytvořit počáteční aproximace pro reciproční s dvojitou přesností divize. Jeden Halley iterace s kubickou konvergenci pak výsledky ve full-přesnost reciproční.

Na Halley iterace pro reciproční je velmi celočíselné násobení intenzivní, s 64x64-bit násobení s 128-bit výsledek (umul64wide() v kódu níže) je stavební blok rozhodující pro výkon. Na moderních 64-bitových architekturách to je obvykle jeden stroj vykonávající instrukce v několika cyklech, nicméně to není přístupné přenosný kód. Přenosný kód emulovat instrukce vyžaduje asi 15 až 20 pokyny v závislosti na stylu a kompilátor.

Celou 128-bit divize by měla trvat zhruba 300 cyklů, nebo deset krát tak rychle jako jednoduché bit-wise výpočtu. Protože kód je poměrně složité, to vyžaduje značné množství testů, aby byla zajištěna správná funkce. V rámci níže jsem pomocí pattern-based a náhodné testy pro středně intenzivní testování, pomocí jednoduché bit-wise provádění jako referenční.

Provádění udiv128() níže předpokládá, že programovací prostředí používá IEEE-754 kompatibilní plovoucí-bod, aritmetika, že double typ je mapován na IEEE-754 je binary64 typu, a to rozdělení double operandy je správně zaoblené.

#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 <[email protected]>
  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;
}

Question 3

To je to, co jsem skončil kódování. Jsem si jistý, že tam jsou mnohem rychlejší alternativy, ale toto je alespoň funkční.

Na základě: https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(nepodepsaný)_with_remainder. Přizpůsobené pro tento konkrétní případ.

// 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;
}

njuffa · Answer 1 · 2021-11-22T01:55:43

Nejsem si vědom nějaké optimalizace, které se vztahují na celočíselné dělení s konstantní dividendy. Dvakrát zkontrolujte, zkoušel jsem testovací případ s all-ty dividendy s Compiler Explorer. Pomocí gcc, icc, a clang, s nejvyšší stupeň optimalizace uvedeno, vygenerovaný kód ukázala, žádné optimalizace se aplikuje do divize.

To je jistě možné vytvořit vysoce-výkon 128-bit rozdělení rutiny, ale z osobní zkušenosti vím, že je poměrně náchylný k chybám, a velmi sofistikované testování je potřeba k dosažení dobrého pokrytí testů včetně rohových případy, jako vyčerpávající test není možné v této operand velikost. Snaha o design a test snadno přesahuje to, co se zdá rozumné pro odpověď na Stackoverflow o dva desítkové řády.

Snadný způsob, jak provádět celočíselné dělení je použít algoritmus, který jsme se všichni naučili ve škole, pouze v binární. To je rozhodnutí o dalším kvocient trochu příliš snadné: je 1 když je aktuální dílčí zbytek je větší než, nebo rovno, dělitel, a 0 jinak. Pomocí ručně binární dělení, pouze celočíselné operace, které potřebujeme, jsou sčítání a odčítání.

Můžeme stavět přenosné primitiva pro provádění těchto na operandy z libovolné bitové délky napodobováním způsobu, jakým procesor počítače instrukce jsou použity pro efekt operace na multi-slovo celá čísla: PŘIDAT s carry-out, PŘIDEJTE se carry-in, ADD carry-in a carry-out; analogicky pro SUB. V kódu níže jsem pomocí jednoduché C makra pro to; určitě sofistikovanější přístupy jsou možné.

Protože systém jsem pracoval na teď není podpora pro 128-bitová celá čísla, já navrhnuté a testované tento přístup pro 64-bitové celá čísla. Na 128-bit verzi, pak bylo cvičení v jednoduché mechanické přejmenování. Na moderní 64-bitový procesor bych očekávat, že tato 128-bit divize funkci vykonávat ve zhruba 3000 cyklů.

#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 <[email protected]>
  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;
}

Rychlejší přístup, než je bit-wise výpočtu je výpočet vzájemné dělitel, vynásobte dividendy, což má za následek předběžné kvocient, pak vypočítat zbývající přesně nastavit kvocient. Celý výpočet lze provést ve fixed-point aritmetiky. Nicméně, na moderní procesory s rychle floating-point jednotky je výhodnější vytvořit počáteční aproximace pro reciproční s dvojitou přesností divize. Jeden Halley iterace s kubickou konvergenci pak výsledky ve full-přesnost reciproční.

Na Halley iterace pro reciproční je velmi celočíselné násobení intenzivní, s 64x64-bit násobení s 128-bit výsledek (umul64wide() v kódu níže) je stavební blok rozhodující pro výkon. Na moderních 64-bitových architekturách to je obvykle jeden stroj vykonávající instrukce v několika cyklech, nicméně to není přístupné přenosný kód. Přenosný kód emulovat instrukce vyžaduje asi 15 až 20 pokyny v závislosti na stylu a kompilátor.

Celou 128-bit divize by měla trvat zhruba 300 cyklů, nebo deset krát tak rychle jako jednoduché bit-wise výpočtu. Protože kód je poměrně složité, to vyžaduje značné množství testů, aby byla zajištěna správná funkce. V rámci níže jsem pomocí pattern-based a náhodné testy pro středně intenzivní testování, pomocí jednoduché bit-wise provádění jako referenční.

Provádění udiv128() níže předpokládá, že programovací prostředí používá IEEE-754 kompatibilní plovoucí-bod, aritmetika, že double typ je mapován na IEEE-754 je binary64 typu, a to rozdělení double operandy je správně zaoblené.

#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 <[email protected]>
  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;
}

Úžasné odpovědi, díky! Dal bych ti víc hlasů, kdybych mohl :)
@tgonzalez89, Aby vám představu: odpověď jak je to nyní představuje asi 10 hodin práce. Design prostor pro celočíselné dělení je poměrně velký, dalo by se snadno věnovat celou diplomovou práci. Dva hlavní přístupy jsou číslice-číslice výpočty (z binárních číslic high-radix, jak je navrženo v komentáři, např. Knuth je Algoritmus D v TAOCP Obj. 2, nebo Papež-Stein) a metody založené na funkční iteraci, které navíc zahrnují problém generování počáteční aproximace. Tady většinou jsem uskakoval ta může červů pomocí double divize jako výchozí bod.

tgonzalez89 · Answer 2 · 2021-11-21T07:08:09

To je to, co jsem skončil kódování. Jsem si jistý, že tam jsou mnohem rychlejší alternativy, ale toto je alespoň funkční.

Na základě: https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(nepodepsaný)_with_remainder. Přizpůsobené pro tento konkrétní případ.

// 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;
}

Jak rozdělit 128-bit dividendy 64-bitový dělitel, kde dividendy jsou bity jsou 1, a tam, kde potřebuji pouze 64 LSBs kvocientu?

Otázka

Nejlepší odpověď

V jiných jazycích

Tato stránka je v jiných jazycích

Populární v této kategorii

Oblíbené položky v této kategorii