Подтвердить что ты не робот

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

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

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

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

4b9b3361

ОТВЕТЫ

Ответ 1

Так как это помечено Visual С++, я дам решение, которое нарушает специфические для MSVC встроенные функции.

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

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

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

  • Для этого требуется x64. Он не будет компилироваться на 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

Ответ 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
}

Ответ 3

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

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

Комментарий GregS верен: в разделе "Искусство компьютерного программирования, второе издание, том 2/" Семинумерные алгоритмы "в конце раздела 4.3.1" Множество прецизионных арифметических/классических алгоритмов "(стр. 255 - 265) копия). Это нелегко прочитать, по крайней мере, не для кого-то вроде меня, который забыл большинство математик за пределами алгебры 7-го класса. Как раз перед, Кнут также охватывает сторону умножения вещей.

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

  • Джек Креншоу раскрывает алгоритмы деления Кнута более читаемым образом в серии статей из журнала Embedded System Programming 1997 (к сожалению, в моих заметках нет точных проблем). К сожалению, статьи из старых вопросов ESP нелегко найти в Интернете. Если у вас есть доступ к университетской библиотеке, возможно, вам понадобятся некоторые проблемы с обратной связью или копия библиотеки CD-ROM ESP.
  • Томас Родеффер из исследования Microsoft опубликовал статью о подразделении Software Integer: http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
  • Статья Карла Хассельстрема "Быстрое разделение больших целых чисел": http://www.treskal.com/kalle/exjobb/original-report.pdf
  • Randall Hyde "Искусство языка ассемблера" (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. Это тоже медленное - выполнение бит-по-бит...

Ответ 4

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

Ответ 5

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

Ответ 6

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

В Delphi 2007 (последнее, что я установил - 32 бита), я бы реализовал MulDiv64 следующим образом:

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;

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

В любом случае, 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));

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

Ответ 7

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

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

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

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

Ответ 8

Здесь метод аппроксимации, который вы можете использовать: (полная точность, если 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.

Ответ 9

В C нет 32x32→64 или 64x64→128 умножений/делений. Результат всегда имеет тот же размер, что и самый большой из множителей и множитель. Это означает, что только int32 x int32 → int32 и int64 x int64 → int64. В случае переполнения результатом являются младшие биты. То же самое для делений. Таким образом, на самом деле вам нужно объявить переменные как __int64 или передать один из множителей/делителей в __int64, и результат будет ниже 64 бит по мере необходимости.

__int64 a, b, c;
c = a * b;
__int32 d, e;
b = __int64(e) / d;

Ответ 10

Считайте, что вы хотите умножить 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;.