Умножение и сложение UInt64 с переполнением.

← →
DayGaykin © (2015-11-19 14:20) [0]

Delphi XE5. Проект 64 бита.

Есть большое число представленное в виде массива A: array[0..3] of UInt64.
Есть число M: UInt64;
Необходимо умножить A на M (и положить в A).
Подскажите как?

Как я понимаю это умножение в системе счисления 2^64.
Умножаем A[3] на M, результат кладем в A[3], а все что выше (в школе это называлось "в уме") прибавляем к A[2]. Если переполнилось, то прибавляем к A[1] и так далее.

Весь вопрос в том, как получить это значение "в уме", которое нужно прибавить к следующему разряду?

Причем как при сложении, так и при умножении.

← →
brother © (2015-11-19 14:27) [1]

Даже не считая: как учитывать переполнение после умножения?

← →
DayGaykin © (2015-11-19 14:29) [2]

Немного исследовав вопрос я нашел что при умножении в 32х битах командой mul результат кладется в EAX, а то что "в уме" в EDX.

Вопрос, справедливо ли это и для 64 бит и самое главное, как это использовать в delphi?

← →
DayGaykin © (2015-11-19 15:20) [3]

Оказывается на delphi под 64 можно писать функции на асме:
Написал такую функцию:

function Mul(A: UInt64; M: UInt64; out X: UInt64): UInt64; assembler; asm mov rax, rcx mul rdx mov [r8], rdx end;

Возвращает A*M, а то что осталось "в уме" в X.

Правильно ли?

← →
Rouse_ © (2015-11-19 17:10) [4]

Самое простое, (конечно не самое быстрое) это побитовое сложение, которое можно нарастить хоть до 1024-битных чисел.
Но так то вообще-то да, RDX - это переполнение.

mul rdx ; Multiply rax by rdx ; rax=low bits, rdx overflow

← →
DayGaykin © (2015-11-19 17:47) [5]

Спасибо.
Сань, глянь, пожалуйста еще:
var A: array[0..3] of UInt64; ... function Mul(AN: PUInt64; M: UInt64; A0: PUInt64): boolean; assembler; asm // rcx - указатель на текущий разряд. Вначале - последний слева разряд. // r8 - указатель на первый слева разряд. mov r9, rdx // Множитель xor r10, r10 // То что на уме xor rbx, rbx // Для определения переполнения @start: mov rax, [rcx] mul r9 // Умножаем текущий разряд на множитель. rax - результат, rdx - на ум. add rax, r10 // Прибавляем то что раннее на ум шло. mov [rcx], rax // Сохраняем результат. seto bl // А если при сложении было переполнение? add rdx, rbx // Добавялем к тому что на ум пойдет. mov r10, rdx // Записываем то что пошло на ум. sub rcx, 8 cmp rcx, r8 // Переходим к следующему разряду. jnz @start // Возвращаем True - если результат не влез в число. cmp r10, 0 setnz al end;
Умножаем там:
Mul(@A[High(A), M, @A[0]]) // M - некий множитель

Вопрос: с точки зрения асма правильно? или надо какие-то регистры в стеке сохранять и восстанавливать?
Если не лень, алгоритм тоже посмотри.

← →
DayGaykin © (2015-11-19 17:48) [6]

Mul(@A[High(A)], M, @A[0]) // M - некий множитель
Опечатался - тороплюсь!:)

← →
DayGaykin © (2015-11-19 18:16) [7]

Если можно как-то уменьшить тело цикла - буду очень благодарен!
Даже удаление 1-2 команд сократит мои расчеты на несколько часов.

← →
Rouse_ © (2015-11-19 18:31) [8]

На вскидку вроде правильно, но что-то мне не нравится. Не могу понять за что глаз зацепился.

← →
Rouse_ © (2015-11-19 18:33) [9]

Удалено модератором
Примечание: пардон, вру...

← →
Rouse_ © (2015-11-19 18:38) [10]

Нет, похоже не вру.
R8 регистр не контролируется, на входе. Если он не кратен 8 - то вечный цикл.

← →
DayGaykin © (2015-11-19 18:50) [11]

> Rouse_ © (19.11.15 18:38) [10]
> Нет, похоже не вру.
> R8 регистр не контролируется, на входе. Если он не кратен
> 8 - то вечный цикл.

Это не может случиться:
Mul(@A[High(A)], M, @A[0]) -- Разница между @A[High(A)] и @A[0] всегда кратна 8, т.к. элементы UInt64.

Дело в том, что в примере я указал массив A: array[0..3] of UInt64;, а на деле в массиве 6 млн элементов. Такое больше число я ожидаю на выходе. Перемножить мне нужно порядка 16 млн чисел (UInt64). Поэтому мне нужно экономить на каждом такте. Переписывание функции на ASM позволило увеличить скорость почти в 4 раза.

Спасибо за проверку! Ты очень помог!

← →
Rouse_ © (2015-11-19 18:57) [12]

> DayGaykin © (19.11.15 18:50) [11]
> Это не может случиться:

При таком вызове конечно не может, но ты сам попросил посмотреть :)
А ускорить, можно, но существенного прироста в данном частном случае не даст.

← →
DayGaykin © (2015-11-19 19:00) [13]

> А ускорить, можно, но существенного прироста в данном частном
> случае не даст.

Тогда оставлю как есть. В последствии мне нужно будет доказать, что этот алгоритм выполняет именно умножение. Чем проще - тем лучше.

← →
Rouse_ © (2015-11-19 19:42) [14]

О как, ты бы лучше тогда Sha подождал, он явно попридирчивей меня будет к коду

← →
Rouse_ © (2015-11-19 19:49) [15]

Зы, rbx не обязательно, но желательно на выходе иметь таким же что и на входе (push/pop ему сделай в эпилоге и прологе)

← →
Pavia © (2015-11-19 20:25) [16]

> Дело в том, что в примере я указал массив A: array[0..3]
> of UInt64;, а на деле в массиве 6 млн элементов. Такое больше
> число я ожидаю на выходе. Перемножить мне нужно порядка
> 16 млн чисел (UInt64).

Не от туда вы решаете задачу. Преждевременная оптимизация это ЗЛО.
Какую задачу вы решаете? Откуда взялись длинные числа? Почему используете школьный алгоритм умножения? А не например алгоритм Карацубы?

> А ускорить, можно, но существенного прироста в данном частном
> случае не даст.

Даст. Умножение легко паралелится. А векторность процессора в данном алгоритме не учтена ни каким образом.

← →
Rouse_ © (2015-11-19 22:08) [17]

> Pavia © (19.11.15 20:25) [16]
> Даст. Умножение легко паралелится. А векторность процессора
> в данном алгоритме не учтена ни каким образом.

Я знаю, но в именно данном случае сильного прироста не даст

← →
Sha © (2015-11-20 00:52) [18]

> DayGaykin © (19.11.15 17:47) [5]

Что сразу бросилось в глаза

1. Можно обойтись меньшим числом регистров, если переполнения обрабатывать командами adc или даже sbb.

2. Если развернуть цикл в 2 или в 4 раза, можно сэкономить на проверках и изменениях индекса.

P.S. Умножение немного странное - переполнение точно требуется выполнять от старших индексов к младшим?

P.P.S. Как-то привычнее передавать адрес первого и количество.

← →
DayGaykin © (2015-11-20 03:14) [19]

> Sha © (20.11.15 00:52) [18]
>

Спасибо за ответ.
Развернул цикл, заменил манипуляции с rbx на adc. Скорость чуть выросла, а расчеты сократились на несколько часов с примерно 38 часов, до 31, что уже хорошо и приемлемо. Спасибо еще раз!

> Pavia © (19.11.15 20:25) [16]

Я буду дольше разбираться с алгоритмом, чем он будет считать. Нынешний позволит найти произведение всех чисел примерно за 31 час - что приемлемо.
К тому же мне нужен точный и понятный алгоритм, т.к. в последствии нужно будет доказать, что полученный результат является произведением заданных чисел. Доказательством будет представленный исходный код. Поэтому мне нужно соблюсти баланс между скоростью и понятностью.

Если хотите, могу выложить массив с множителями.

← →
DayGaykin © (2015-11-20 03:26) [20]

Интересно! Насколько я знаю, когда я делаю mov rax, [rcx] или mov [rcx], rax или прыжки операционная система проверяет, есть ли у меня доступ к этой памяти. Можно ли как-то эту проверку отключить?

← →
Pavia © (2015-11-20 06:45) [21]

> Интересно! Насколько я знаю, когда я делаю mov rax, [rcx]
> или mov [rcx], rax или прыжки операционная система проверяет,
> есть ли у меня доступ к этой памяти. Можно ли как-то эту
> проверку отключить?

Не беспокойтесь такая проверка происходит за 1 такт. И не отключается.

> Я буду дольше разбираться с алгоритмом, чем он будет считать.
> Нынешний позволит найти произведение всех чисел примерно
> за 31 час - что приемлемо.
Дело ваша за секунды-минуты считать или за часы. Да и разбираться с ним проще. И доказывать проще так как алгоритм известный и для него всё доказано.

Хотя честно мне интересно было узнать для чего понадобилось делать умножение? Похоже это учебная задача, на освоение алгоритмов и я думаю вас заставят делать через Карацубу. Единственная задача которую я знаю где нужны длинные числа это шифрование. Но там уже всё реализовано по 100 раз. Поэтому проще взять готовое.
В астрономии и гис может не хватать точности 64 бит и там да выгодно использовать 128 бит. В остальных случаях случаях длинные числа бесполезны.

← →
Sha © (2015-11-20 09:34) [22]

> DayGaykin

Код покажи, что в итоге получилось.
Может, там еще где-нибудь на спичках экономия возможна.

← →
DayGaykin © (2015-11-20 12:20) [23]

> Код покажи, что в итоге получилось.
> Может, там еще где-нибудь на спичках экономия возможна.

function Mul(M{RCX}: UInt64): boolean; assembler; asm mov r9, rcx // Множитель xor r10, r10 // То что на уме mov r8, offset MyNumber[0] // Ссылка на первый элемент числа mov rcx, Digits64 * 8 // Длина lea rcx, rcx + r8 - 8 // rcx - последний разряд @@1: mov rax, [rcx] mul r9 // Умножаем текущий разряд на множитель. rax - результат, rdx - на ум. add rax, r10 // Прибавляем то что раннее на ум шло. mov [rcx], rax // Сохраняем результат. adc rdx,0 // Если при прибавлении мы получили переполнение, прибавляем к следующему "на уме" единицу mov r10, rdx // И помещаем для следующей итерации sub rcx, 8 mov rax, [rcx] mul r9 // Умножаем текущий разряд на множитель. rax - результат, rdx - на ум. add rax, r10 // Прибавляем то что раннее на ум шло. mov [rcx], rax // Сохраняем результат. adc rdx,0 // Если при прибавлении мы получили переполнение, прибавляем к следующему "на уме" единицу mov r10, rdx // И помещаем для следующей итерации sub rcx, 8 cmp rcx, r8 jnb @@1 // Возвращаем True - если результат не влез в число. cmp r10, 0 setnz al end;
Немного переделаю еще так, чтобы не нужно было умножать нули вначале расчетов и чтобы функция не "портила" число, если оно не влезет.

> Pavia © (20.11.15 06:45) [21]

Одноразовый математический расчет. Ищут математики первый миллиард цифр числа Pi, вот и я кое чего ищу.

← →
DayGaykin © (2015-11-20 19:18) [24]

Еще вопрос:
Вот получил я число (на самом деле еще считается), записанное в массиве:
N: array[0..6000000] of UInt64 (В памяти оно лежит так-же, как оно было бы записано в мысленном типе данных UInt384000064).
Насколько сложно получить логарифм этого числа (в принципе по любому основанию, для простоты можно взять основание 2^64)?
Целую часть логарифма получить, думаю, вообще не составляет труда, а как получить дробную?

← →
Sha © (2015-11-20 21:16) [25]

> DayGaykin

Приближенно вычисляешь по любому удобному основанию, например, 2. Потом приводишь к нужному.

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

← →
DayGaykin © (2015-11-20 22:36) [26]

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

6000000 итераций для каждого умножения. 16000000 умножений. Приличное количество итераций выходит.

145 умножений в секунду.

А данные откуда? Столько в памяти не поместится.

> А данные откуда? Столько в памяти не поместится.

Напишу позже, сейчас хочу оставить задачу в тайне, чтобы избежать ненужного флуда в ветке.
16млн умножений - это столько у меня заготовлено множителей. Конечно, не влезет в 6 млн разрядов в системе счисления 2^64. Задача позволяет мне не умножать все, а умножать только до тех пор пока старший разряд из этих 6 млн не примет какое-либо ненулевое значение.
Я думаю уже можно догадаться что за задача, но все-же прошу воздержаться от догадок.

P.S.
Пока у меня ведутся расчеты возник вопрос:
Есть у нас массив A[0..N-1] of UInt64, где N достаточно больше число (порядка 10^6--10^8).
Массив A содержит число, где в каждом элементе находится значение разряда в 2^64 системе счисления. Порядок большого значения не имеет.
Есть B: UInt64. - т.е. одноразрядное число.
Задача: умножить A на B и результат положить в A.
Неужели это умножение можно произвести быстрее чем в столбик?

При умножении, каждый разряд числа A в любом случае нужно будет изменить, каким бы не был алгоритм. В столбик на каждый разряд приходится всего одно умножение и два сложения. Я считаю, что быстрее чем в столбик не получится.

Встречаются задачи, которые могут быть решены без умножения.
Просто мысль )

Рассчитал число, теперь приступаю к его анализу и сразу возникла проблема/вопрос:

Может есть готовая библиотека для расчетов с плавающей точкой повышенной точности?

Мне нужно умножать, делить (или находить 1/x), находить логарифм. Желательно чтобы точность регулировалась, но хотелось бы не менее 1 млн значащих цифр.

> При умножении, каждый разряд числа A в любом случае нужно
> будет изменить, каким бы не был алгоритм. В столбик на каждый
> разряд приходится всего одно умножение и два сложения.
> Я считаю, что быстрее чем в столбик не получится.

Можно. Начните распаралеливать, и вы увидите что второе сложение можно "сократить".

> Может есть готовая библиотека для расчетов с плавающей точкой
> повышенной точности?

Только на Си, на паскале нету. Хотел делать но понял что это безперспективно достаточно чисел с плавающей точкой 64 бита. А если кому не хватает, то удвоить точность.

> Pavia © (21.11.15 18:31) [32]
> второе сложение можно "сократить"

Собственно операция сложения выполняется гораздо быстрее операции сложения с учетом переноса. А от переноса никуда не уйти. К тому же сильно "растащить" сложения с переносом не особо получится: мешают другие арифметические операции, а неарифметических там всего две. Добавить пару переименований, чтобы уйти от сложения - сомнительная идея.

У меня (на 32 битах) основной цикл столбика выглядит так:
@3: mul edi xor ebx, ebx add ebp, eax mov eax, [esi] mov [esi-4], ebp adc ebx, edx @0: mul edi xor ebp, ebp add ebx, eax mov eax, [esi+4] mov [esi], ebx adc ebp, edx @1: mul edi xor ebx, ebx add ebp, eax mov eax, [esi+8] mov [esi+4], ebp adc ebx, edx @2: mul edi xor ebp, ebp add ebx, eax mov eax, [esi+12] mov [esi+8], ebx adc ebp, edx add esi, 16 sub ecx, 4 jg @3

Мне не совсем видно, как здесь можно было бы повысить скорость.

С другой стороны, просто запустив несколько потоков вычислений по одному для каждой группы множителей, можно кратно ее увеличить.

С третьей стороны ) на 64 битах доступных регистров больше. Там можно в рамках одного потока на каждом ядре организовать 2 виртуальных потока вычислений, по одному для каждой из двух половинок большого числа, а потом "срастить" их. Возможно, это будет немного быстрее.

> Там можно в рамках одного потока на каждом ядре организовать
> 2 виртуальных потока вычислений

Это как так?

> DayGaykin © (21.11.15 20:18) [34]

Грубо говоря идти по циклу не так: 0, 1, 2, ... N-2, N-1
А так: 0, N/2, 1, N/2+1, ... N/2-1, N-1

Всем спасибо за помощь.

Число рассчитал, но моя идея по укорачиванию его записи провалилась, способ, который меня казался интересным, теперь мне кажется не эффективным.
Идея была перемножить первые N (как оказалось примерно 22 миллиона) простых чисел, а затем представить в виде небольшого количества слагаемых вида a^b. Поиск a и b я думал сделать следующим образом: a перебираем подряд, b рассчитываем с помощью логарифмирования. Перебирая "a" выбираем такие значение, для которых ln(M)/ln(a) максимально близко к натуральному. Другими словами дробная часть максимально близка к 0 или к 1.

Число у меня записано в бинарном виде и занимает чуть более 48 млн байт. Меня удивил тот факт, что этот файл не поддается сжатию ни 7z ни Zip. Этот факт меня еще больше убедил в сложности решаемой мной задачи. Можно, конечно, разбираться дальше, но мой любительский интерес на этом закончен.

Большое спасибо Sha и Rouse_ за комментарии и помощь.

В предыдущем значении M - то самое мое число.

Умножение и сложение UInt64 с переполнением. Найти похожие ветки