Самый точный способ сделать операцию умножения и деления в 64-битной?

Question

Как наиболее точно можно выполнить операцию умножения и деления для 64-разрядных целых чисел, которая работает как в 32-разрядных, так и в 64-разрядных программах (в Visual C ++)?(В случае переполнения мне нужен мод результата 2 ⁶⁴ .)

(я ищу что-то вроде MulDiv64 , за исключением того, что этот использует встроенную сборку, который работает только в 32-битных программах.)

Очевидно, что приведение к double и обратно возможно, но мне интересно, есть ли более точный способ, который не слишком сложен.(т.е. я не ищу здесь арифметические библиотеки произвольной точности!)

Answer 1

Поскольку это помечено Visual C ++, я дам решение, которое использует MSVC-специфические свойства.

Этот пример довольно сложный. Это очень упрощенная версия того же алгоритма, который используется GMP и java.math.BigInteger для большого деления.

Хотя я имею в виду более простой алгоритм, он, вероятно, примерно в 30 раз медленнее.

Это решение имеет следующие ограничения / поведение:

• Требуется х64. Он не скомпилируется на x86.
• Коэффициент не равен нулю.
• Коэффициент насыщается, если он переполняется 64-битным.

Обратите внимание, что это для целого числа без знака. Тривиально создать оболочку для этого, чтобы она работала и для подписанных дел. Этот пример также должен давать правильно усеченные результаты.

Этот код не полностью протестирован. Тем не менее, он прошел все тестовые примеры, которые я к нему применил.
(Даже случаи, которые я намеренно создал, пытаясь сломать алгоритм.)

#include <intrin.h>

uint64_t muldiv2(uint64_t a, uint64_t b, uint64_t c){
    //  Normalize divisor
    unsigned long shift;
    _BitScanReverse64(&shift,c);
    shift = 63 - shift;

    c <<= shift;

    //  Multiply
    a = _umul128(a,b,&b);
    if (((b << shift) >> shift) != b){
        cout << "Overflow" << endl;
        return 0xffffffffffffffff;
    }
    b = __shiftleft128(a,b,shift);
    a <<= shift;


    uint32_t div;
    uint32_t q0,q1;
    uint64_t t0,t1;

    //  1st Reduction
    div = (uint32_t)(c >> 32);
    t0 = b / div;
    if (t0 > 0xffffffff)
        t0 = 0xffffffff;
    q1 = (uint32_t)t0;
    while (1){
        t0 = _umul128(c,(uint64_t)q1 << 32,&t1);
        if (t1 < b || (t1 == b && t0 <= a))
            break;
        q1--;
//        cout << "correction 0" << endl;
    }
    b -= t1;
    if (t0 > a) b--;
    a -= t0;

    if (b > 0xffffffff){
        cout << "Overflow" << endl;
        return 0xffffffffffffffff;
    }

    //  2nd reduction
    t0 = ((b << 32) | (a >> 32)) / div;
    if (t0 > 0xffffffff)
        t0 = 0xffffffff;
    q0 = (uint32_t)t0;

    while (1){
        t0 = _umul128(c,q0,&t1);
        if (t1 < b || (t1 == b && t0 <= a))
            break;
        q0--;
//        cout << "correction 1" << endl;
    }

//    //  (a - t0) gives the modulus.
//    a -= t0;

    return ((uint64_t)q1 << 32) | q0;
}

Обратите внимание, что если вам не нужен идеально усеченный результат, вы можете полностью удалить последний цикл. Если вы сделаете это, ответ будет не более чем на 2 больше, чем правильный коэффициент.

Контрольные примеры:

cout << muldiv2(4984198405165151231,6132198419878046132,9156498145135109843) << endl;
cout << muldiv2(11540173641653250113, 10150593219136339683, 13592284235543989460) << endl;
cout << muldiv2(449033535071450778, 3155170653582908051, 4945421831474875872) << endl;
cout << muldiv2(303601908757, 829267376026, 659820219978) << endl;
cout << muldiv2(449033535071450778, 829267376026, 659820219978) << endl;
cout << muldiv2(1234568, 829267376026, 1) << endl;
cout << muldiv2(6991754535226557229, 7798003721120799096, 4923601287520449332) << endl;
cout << muldiv2(9223372036854775808, 2147483648, 18446744073709551615) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 9223372036854775807) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 12) << endl;
cout << muldiv2(18446744073709551615, 18446744073709551615, 9223372036854775808) << endl;

Выход:

3337967539561099935
8618095846487663363
286482625873293138
381569328444
564348969767547451
1023786965885666768
11073546515850664288
1073741824
9223372032559808512
Overflow
18446744073709551615
Overflow
18446744073709551615

Answer 2

Вам просто нужно 64-битные целые числа. Есть несколько избыточных операций, но это позволяет использовать 10 в качестве базы и шага в отладчике.

uint64_t const base = 1ULL<<32;
uint64_t const maxdiv = (base-1)*base + (base-1);

uint64_t multdiv(uint64_t a, uint64_t b, uint64_t c)
{
    // First get the easy thing
    uint64_t res = (a/c) * b + (a%c) * (b/c);
    a %= c;
    b %= c;
    // Are we done?
    if (a == 0 || b == 0)
        return res;
    // Is it easy to compute what remain to be added?
    if (c < base)
        return res + (a*b/c);
    // Now 0 < a < c, 0 < b < c, c >= 1ULL
    // Normalize
    uint64_t norm = maxdiv/c;
    c *= norm;
    a *= norm;
    // split into 2 digits
    uint64_t ah = a / base, al = a % base;
    uint64_t bh = b / base, bl = b % base;
    uint64_t ch = c / base, cl = c % base;
    // compute the product
    uint64_t p0 = al*bl;
    uint64_t p1 = p0 / base + al*bh;
    p0 %= base;
    uint64_t p2 = p1 / base + ah*bh;
    p1 = (p1 % base) + ah * bl;
    p2 += p1 / base;
    p1 %= base;
    // p2 holds 2 digits, p1 and p0 one

    // first digit is easy, not null only in case of overflow
    uint64_t q2 = p2 / c;
    p2 = p2 % c;

    // second digit, estimate
    uint64_t q1 = p2 / ch;
    // and now adjust
    uint64_t rhat = p2 % ch;
    // the loop can be unrolled, it will be executed at most twice for
    // even bases -- three times for odd one -- due to the normalisation above
    while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
        q1--;
        rhat += ch;
    }
    // subtract 
    p1 = ((p2 % base) * base + p1) - q1 * cl;
    p2 = (p2 / base * base + p1 / base) - q1 * ch;
    p1 = p1 % base + (p2 % base) * base;

    // now p1 hold 2 digits, p0 one and p2 is to be ignored
    uint64_t q0 = p1 / ch;
    rhat = p1 % ch;
    while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
        q0--;
        rhat += ch;
    }
    // we don't need to do the subtraction (needed only to get the remainder,
    // in which case we have to divide it by norm)
    return res + q0 + q1 * base; // + q2 *base*base
}

Answer 3

Это ответ вики-сообщества, поскольку на самом деле это просто набор указателей на другие статьи / ссылки (я не могу опубликовать соответствующий код).

Умножение двух 64-битных целых на128-битный результат довольно прост, если использовать простую технику карандаша и бумаги, которую все учат в начальной школе.

Комментарий Грегса верен: Кнут рассказывает о разделении «Искусство компьютерного программирования», второе издание, том 2 / Полу числовойАлгоритмы »в конце Раздела 4.3.1 Арифметика с множественной точностью / Классические алгоритмы (страницы 255 - 265 в моем экземпляре).Это не легко читать, по крайней мере, не для тех, кто, как я, забыл большинство математики за 7-й класс алгебры.Как и ранее, Кнут охватывает также и умножение.

Некоторые другие варианты идей (эти заметки предназначены для алгоритмов деления, но большинство также обсуждают умножение):

• Джек Креншоуболее подробно читается об алгоритмах деления Кнута в серии статей из журнала Embedded System Programming 1997 (к сожалению, в моих заметках нет точных проблем).К сожалению, статьи из старых проблем ESP не так легко найти в Интернете.Если у вас есть доступ к университетской библиотеке, возможно, у вас есть какие-то проблемы или копия библиотеки ESP CD-ROM.
• У Томаса Родехеффера из Microsoft Research есть статья о Software Integer Division: http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
• Статья Карла Хассельстрёма "Быстрое деление больших целых чисел": http://www.treskal.com/kalle/exjobb/original-report.pdf
• "Язык ассемблера" Рэндалла Хайда (http://webster.cs.ucr.edu/AoA/Windows/HTML/AoATOC.html), в частности, том четвертый, раздел 4.2.5 (Расширенное прецизионное деление)): http://webster.cs.ucr.edu/AoA/Windows/HTML/AdvancedArithmetica2.html#998729 это в варианте Hyde на ассемблере x86, но есть также некоторый псевдокод и достаточное объяснение, чтобы перенести алгоритм на C. Это тоже медленно - выполнение побитового деления ...

Answer 4

Вам не нужна произвольная точность арифметики для этого.Вам нужна только 128-битная арифметика.Т.е. вам нужно умножение 64 * 64 = 128 и деление 128/64 = 64 (с правильным поведением переполнения).Это не так сложно реализовать вручную.

Answer 5

Есть ли у вас в распоряжении VC ++ тип COMP (64-разрядный целочисленный тип на основе x87)?Я использовал это время от времени в Delphi, когда мне была нужна 64-битная целочисленная математика.В течение многих лет это было намного быстрее, чем основанная на библиотеке 64-битная целочисленная математика - конечно, когда делалось это деление.

В Delphi 2007 (последняя установленная мной версия - 32-битная версия) я реализовал бы MulDiv64 какthis:

function MulDiv64(const a1, a2, a3: int64): int64;
var
  c1: comp absolute a1;
  c2: comp absolute a2;
  c3: comp absolute a3;
  r: comp absolute result;
begin
  r := c1*c2/c3;
end;

(Эти странные операторы absolute накладывают переменные comp поверх их 64-битных целочисленных счетчиков. Я бы использовал приведение простых типов, за исключением того, что компилятор Delphi запуталсяоб этом - вероятно, потому что язык Delphi (или как его там сейчас называют) не имеет четкого синтаксического различия между приведением типов (переинтерпретация) и преобразованием типов значений.)

В любом случае, Delphi 2007 отображает вышеизложенное следующим образом:

0046129C 55               push ebp
0046129D 8BEC             mov ebp,esp
0046129F 83C4F8           add esp,-$08

004612A2 DF6D18           fild qword ptr [ebp+$18]
004612A5 DF6D10           fild qword ptr [ebp+$10]
004612A8 DEC9             fmulp st(1)
004612AA DF6D08           fild qword ptr [ebp+$08]
004612AD DEF9             fdivp st(1)
004612AF DF7DF8           fistp qword ptr [ebp-$08]
004612B2 9B               wait 

004612B3 8B45F8           mov eax,[ebp-$08]
004612B6 8B55FC           mov edx,[ebp-$04]
004612B9 59               pop ecx
004612BA 59               pop ecx
004612BB 5D               pop ebp
004612BC C21800           ret $0018

Следующий оператор выдает 256204778801521550, что представляется правильным.

writeln(MulDiv64($aaaaaaaaaaaaaaa, $555555555555555, $1000000000000000));

Если вы хотите реализовать это как встроенную сборку VC ++, возможно, вам потребуется выполнить некоторыенастройка флагов округления по умолчанию для достижения той же цели, я не знаю - мне не нужно былопока нет - :)

Answer 6

Ну, вы можете разделить свои 64-битные операнды на 32-битные блоки (младшая и старшая части).Чем делать операции на них, что вы хотите.Все промежуточные результаты будут менее 64 бит и поэтому могут быть сохранены в ваших типах данных.

Answer 7

Если вам нужна только поддержка Windows 7 и новее, один хороший способ это:

#include <mfapi.h>
#include <assert.h>
#pragma comment( lib, "mfplat.lib" )

uint64_t mulDiv64( uint64_t a, uint64_t b, uint64_t c )
{
    assert( a <= LLONG_MAX && b <= LLONG_MAX && c <= LLONG_MAX );
    // https://docs.microsoft.com/en-us/windows/desktop/api/Mfapi/nf-mfapi-mfllmuldiv
    return (uint64_t)MFllMulDiv( (__int64)a, (__int64)b, (__int64)c, (__int64)c / 2 );
}

Этот метод намного проще, чем в других ответах, он округляет результат вместо усечения и работает на всех платформах Windows, включая ARM.

Answer 8

Учтите, что вы хотите умножить a на b, а затем разделить на d:

uint64_t LossyMulDiv64(uint64_t a, uint64_t b, uint64_t d)
{
    long double f = long double(b)/d;
    uint64_t highPart = uint64_t((a & ~0xffffffff) * f + 0.5);
    uint64_t lowPart = uint64_t((a & 0xffffffff) * f + 0.5);
    return highPart + lowPart;
}

Этот код разбивает значение a на более высокие и более низкие 32-битные частиЗатем 32-битные части умножаются отдельно на 52-битное точное соотношение от b до d, округляются частичные умножения и суммируются обратно до целого числа.Некоторая точность все еще теряется, но результат более точен, чем просто return a * double(b) / d;.

Answer 9

Вот метод аппроксимации, который вы можете использовать: (полная точность, если a> 0x7FFFFFFF или b> 0x7FFFFFFF и c больше, чем a или b)

constexpr int64_t muldiv(int64_t a, int64_t b, int64_t c, unsigned n = 0) {
  return (a < 0x7FFFFFFF && b < 0x7FFFFFFF) ? (a * b) / c : (n != 2) ? (c <= a) ? ((a / c) * b + muldiv(b, a % c, c, n + 1)) : muldiv(a, b, c / 2) / 2 : 0;
}

Модуль используется для определения потери точности, который затем подключается обратно в функцию.Это похоже на классический алгоритм деления.

2 было выбрано, потому что (x % x) % x = x % x.

Answer 10

Для кода 64-битного режима вы можете реализовать умножение 64 * 64 = 128 аналогично реализации 128/64 = 64: 64 деления здесь .

Для 32-битного кодаэто будет более сложным, потому что нет инструкции ЦП, которая бы выполняла умножение или деление таких длинных операндов в 32-битном режиме, и вам придется объединить несколько меньших умножений в большее и переопределить длинное деление.

Вы можете использовать код из этого ответа в качестве основы для построения длинного деления.

Конечно, если ваши делители всегда меньше 2 ³² (или лучшеТем не менее, 2 ¹⁶ ), вы можете сделать намного более быстрое деление в 32-битном режиме путем объединения нескольких делений (делить старшие 64 (или 32) бита делимого на 32-битный (или 16-битный)) делитель, затем объедините остаток от этого деления со следующими 64 (32) битами дивиденда, разделите его на делитель и продолжайте делать это, пока не израсходуете весь дивиденд).Кроме того, если делитель большой, но может быть разложен на достаточно малые числа, деление на его факторы с использованием этого цепочечного деления будет лучше, чем классическое решение с циклическим вычислением.