Умножить два 64-битных числа и поделить на третье

#c #алгоритм #арифметика

Есть ли быстрый способ найти значения этих выражений, используя 64-битные целочисленные
переменные?

x * 2**64 / 10**19
(x * 2**64) % 10**19


где 0 <= x < 10**19, **— возведение в степень. Нельзя использовать вещественные числа
и переменные размером больше 64 бит.
    

Ответы

Ответ 1

264 / 1019 = 245 / 519,


245 = 9*518 + 3490*512 + 2938*56 + 10707,

245 / 519 = 9/5 + 3490/57 + 2938/513 + 10707/519.


x = 519r +q, 
где r = x / 519; q = x % 519, q<519.
В первом случае:


(x*264)/1019 = ((519r +q) * 245) / 519,

(x*264)/1019  = 245r + r0 + r1 + r2 +r3,

r0 = (9q) / 5; q0 = (9q) % 5;
r1 = (q056 + 3490q) / 57; q1 = (q056 + 3490q) % 57;
r2 = (q156 + 2938q) / 513; 
q2 = (q156 + 2938q) % 513;
r3 = (q256 + 10707q) / 519; 


Bo втором случае можно выносить общий множитель:


(x*264) % 1019 = 219((519r +q)*245) % 519) = 219((q*245) % 519), 

(x*264) % 1019 = 219(518((9q) % 5) + 512((3490q) % 57) + 56((2938q) % 513) + 10707q)
%519). 

Алгоритм линейный, проблемы с программой исключены.

Пусть x = 2 345 678 912 345 678 912 = 122 981 * 519 + 2 490 226 538 287,
r = 122 981, q = 2 490 226 538 287. 

Прямые вычисления дают
(1) = 4 327 013 857 513 811 925,
(2) = 9 969 087 398 626 721 792.

Вычисления по предложенному алгоритму дают:

245r = 4 327 009 263 856 648 192;
r0 = 4 482 407 768 916, q0 = 3;
r1 = 111 243 399 918, q1 = 74 755;
r2 = 5 993 502, q1 = 116 440 331;
r3 = 1 397,

(1): 4 327 013 857 513 811 925;

(2): 219((518*3 + 512*27 880 + 56*169 096 581 + 17 195 145 048 284) % 519) = 219(((3*15
625 + 27 880)*15 625 + 169 096 581) *15 625 + 17 195 145 048 284) % 519) = 9 969 087
398 626 721 792.


Ответ 2

Творчески переработав Уоррена, вот такая функция для вашего деления и получения остатка:

unsigned long long divlu(unsigned long long x, unsigned long long * rem)
{
    unsigned long long y = 0ll, z = 10000000000000000000ull, t;
    long long i;
    for(i = 1; i <= 64; ++i)
    {
        t = (long long)x >> 63;
        x = (x << 1)| (y >> 63);
        y <<= 1;
        if ((x|t) >= z)
        {
            x -= z;
            ++y;
        }
    }
    if (rem) *rem = x;
    return y;
}


Передаете в нее свой x, и получаете на выходе частное, а в rem - делитель...

Выборочная проверка с длинной арифметикой подтверждает корректность.

Update
Полный код для общего случая деления 128-битного xy (т.е. понятно, что в x - старшие
64 бита, в y - младшие) на 64-битное z:

unsigned long long divlu(unsigned long long x, 
                         unsigned long long y,
                         unsigned long long z,
                         unsigned long long * rem)
{
    unsigned long long t;
    long long i;
    for(i = 1; i <= 64; ++i)
    {
        t = (long long)x >> 63;
        x = (x << 1)| (y >> 63);
        y <<= 1;
        if ((x|t) >= z)
        {
            x -= z;
            ++y;
        }
    }
    if (rem) *rem = x;
    return y;
}