Трюк настолько мне понравился простотой и изяществом, что уже в этом тысячелетии я с удовольствием рассказывал о нём студентам в виде следующей задачи.
Представьте себе, что в процессе возвращения в 2030 году с Луны вы вдруг обнаружили, что ваш бортовой компьютер неправильно выполняет операции целочисленного умножения, и это непременно приведёт к аварии при посадке.
В таком сюжете нет ничего особо фантастического. Вспомним, например, какие проблемы случались когда-то с процессорами Pentium , а к моменту отправки на Луну вы ещё не достигли полного импортозамещения. И вообще надо проверить, а не были ли процессоры просверлены специально.
Но к делу. Вам надо срочно реализовать умножение программно, чтоб работало быстро в реальном времени и укладывалось в доступный ресурс.
Из школьного курса арифметики вспоминаем, что многозначные числа умножать можно столбиком, а результат умножения отдельных цифр брать из таблицы умножения.
Вот только если цифры выбрать короткими (например 0 и 1), таблица умножения будет короткой, а столбики длинными, и их вычисление будет занимать много времени.
Если, наоборот, взять цифры длинные (например от 0 до 65535), то для 16-битной арифметики
результат берётся прямо из таблицы, а столбики отсутствуют. Однако размер классической таблицы Пифагора при этом оказывается около 17GB (4*65536*65536), если учесть симметрию относительно диагонали, то половина — около 8.5GB.
Может оказаться многовато.
Напрягаемся и вспоминаем алгебру.
(1)
(2)
Вычитаем (2) из (1)
и далее
Таким образом, имея таблицу квадратов в массиве sqr, получаем
a * b = ( sqr[a + b] — sqr[a — b]) / 4 (*)
Размер таблицы 8*(65535+65535) около 8.4MB, что уже может оказаться приемлемым.
Размер элемента таблицы в 8 байт связан с тем, что при максимальных a и b квадрат их суммы в 4 байта не влезает — не хватает 2-х бит.
Далее опишу некоторое улучшение, которого в книжке не было. Его я придумал сам, когда уже писал эту заметку.
Заметим, что два младших бита квадрата чётного числа всегда 00, а квадрата нечётного — всегда 01. С другой стороны для любой пары чисел их сумма и разность имеют одинаковую чётность.
Поэтому в нашей формуле (*) процесс вычитания в скобках не может иметь переносов,
связанных с двумя младшими битами. Поэтому содержимое элементов таблицы квадратов
можно заранее сдвинуть на два бита вправо и тем самым сэкономить половину памяти.
Окончательно имеем
a * b = sqr4[a + b] — sqr4[a — b] (**)
где sqr4 — модифицированная таблица квадратов.
В нашем примере её размер около 4.2 MB.
Ниже для иллюстрации этого подхода включен текст программы на языке Lua.
function create_multiplier(N) -- N это число разных цифр
local sqr4 = {} -- создаём таблицу
for i = 1, 2 * N - 1 do
local temp = 0
for j = 1, i - 1 do -- для вычисления квадрата
temp = temp + i - 1 -- используем многократное сложение
end -- т.к. умножение "неисправно"
sqr4[i] = math.floor(temp / 4) -- экономим 2 бита
end
return -- возвращаем мультипликатор (функцию)
function (i, j)
if i < j then i, j = j, i end
return sqr4[1 + i + j] - sqr4[1 + i - j]
end
end
N = 256
mpy = create_multiplier(N)
for i = 0, N - 1 do
for j = 0, N - 1 do
if i * j ~= mpy(i,j) then
print("Ошибка", i, j, i * j, mpy(i,j))
end
end
end
print("Всё работает.")
Для современных процессоров представляется разумным иметь размер цифр кратным размеру байта для лёгкого доступа к ним. При цифрах размером в 1 байт, размер таблицы всего лишь 1022 байт, что может позволить использовать этот трюк в 8-битных процессорах, не имеющих аппаратного умножения.
Буду благодарен всем читателям этой заметки за исправления и замечания.
Комментарии (34)
NightShad0w
20.06.2019 14:14Спасибо за статью.
А можно узнать предпосылки к «вычитаем 2 из 1»?
ladle Автор
20.06.2019 14:37Если я правильно понял вопрос,
то имелась в виду рутинная алгебра — вычитаем формулу (2) из формулы (1).NightShad0w
20.06.2019 14:40Нет, извините за не подробный вопрос. Я имею в виду — откуда возникла мысль вычесть одно выражение из другого и собственно откуда и как появилась идея начать с рассмотрения выражений для квадратов суммы и разности. Предполагаю, существуют какие-то работы или исследования, результатами которых вы и воспользовались. Не могли бы поделиться источниками или ссылками?
ladle Автор
20.06.2019 15:12Так я и написал, что прочитал об этом в популярной книжке несколько десятков лет назад.
Книжка давно утрачена. Я пытался найти ссылки перед написанием заметки, так как опасался наткнуться на статью в Википедии на эту тему и попасть пальцем в небо с актуальностью материала, но в поиске не преуспел, может быть неудачно формировал поисковые запросы. Может быть у кого-то есть ссылки на эту тему, буду благодарен если поделитесь.
kinhort
20.06.2019 17:50В англоязычной википедии это называется «Quarter square multiplication».
ladle Автор
20.06.2019 17:56Огромное спасибо!
Ещё бы ссылку на книжку найти — там было ещё несколько интересных трюков.
GCU
20.06.2019 14:21Поскольку возведение в квадрат это тоже операция умножения, то интересно сколько элементов из таблицы можно вычислить из той-же таблицы за одно (а может и несколько) «умножений». Например любой квадрат нечётного числа сведётся к квадратам чётных, может быть половину таблицы можно так выкинуть.
Flux
20.06.2019 15:14Для вычисления таблицы не нужно ни одного умножения, достаточно только сложений.
(x + 1)^2 = x^2 + 2x + 1
Соответственно таблицу можно проредить вNраз ценойN-1дополнительных операций при каждом вычислении квадрата, как пример.GCU
20.06.2019 15:41Ну с таким же успехом можно и предыдущее значение найти
(x - 1)^2 = x^2 - 2x + 1
Это меньше операций, если идти от ближайшего :)
0xd34df00d
20.06.2019 16:27Осталось побенчмаркать и выяснить, насколько это оправдано на современном железе с медленной памятью и быстрыми вычислительными ядрами.
GCU
20.06.2019 18:16+1Современное железо с неисправным умножением скорее всего не сможет даже загрузить ОС. Статья скорее относится к неким бракованным микроконтроллерам, где неисправную инструкцию можно попробовать «обойти», а показатели скорости памяти и вычислительных ядер там могут очень сильно отличаться.
Zenitchik
20.06.2019 19:45+1Причём тут брак и неисправная инструкция? Умножения может тупо не быть на микроконтроллере.
andy_p
20.06.2019 18:25Хорошая статья. Надо напомнить, что, к примеру, очень популярный z80 не имеет аппаратного умножения.
defuz
20.06.2019 18:40Оптимизация так себе на самом деле. Смотрите, для 8-битных аргументов аппаратное умножение представляет собой 64 AND-операции с последующим сумированием получаемых 64 бит.
В вашем случае, перед обращением в таблицу нужно вычислить a+b и a-b, что соответствует суммированию 2*2*8=32 бит. После обращения к таблице нужно еще раз вычислить разницу уже для 16-битных значений, это еще 32 бита которые необходимо просуммировать.
Получается в обоих случаях нужно выполнить сумирование 64 бит, но в вашем подходе 64 параллельных AND заменяются на lookup таблицу и 2 обращения к ней – сомнительная оптимизация как по скорости, так и по площади или энергопотреблению.defuz
20.06.2019 18:49Хотя я пожалуй соглашусь что такая реализация вычислений может быть эфективной в ситуации когда «работаем с тем что имеем»: есть процессор который умеет только сумирование и не жалко тратить память, которая при этом относительно быстрая. В таком случае действительно можно получить умножение за меньшее количество тактов, заменив кучу сдвигов на 2 обращения к памяти, правда энергопотребление будет скорее всего выше (обращение в память сильно дороже арифметики и работы с регистрами).
GCU
20.06.2019 19:02Там ещё было сравнение с перестановкой аргументов :)
перед обращением в таблицу нужно вычислить a+b и a-b, что соответствует суммированию 2*2*8=32 бит
Почему 32, а не 16?
a + b= 8 бит + флаг переполнения/переноса
a - b= 8 бит (тут наверное «бесплатное» сравнение)defuz
20.06.2019 19:13Я считаю количество бит-аргументов, которые нужно сложить, а не количество результирующих бит, по-этому 32. Просто я смотрю на эту оптимизацию как на аппаратную, а не софтварную. :)
Кстати, обратил внимание, что даже после срезания 2 битов, оставшиеся младшие 3 бита повторяются с перодом цикла 16. Так что если есть задача уменьшить размер таблицы – можно смело отсекать еще 3 бита от каждого значения и добавлять еще одну таблицу с 16 значениями.
Думаю есть еще какие-то паттерны, которые позволят декомпозировать задачу и тем самым уменьшить размер таблицы. Например, 4ый бит повторяется с шагом 32.GCU
20.06.2019 19:25А если + и — сделать одной операцией с двумя выходами :)
Хитрая операция ab+- где на выходе 16 бит a+b и a-bmayorovp
21.06.2019 15:40Вряд ли можно сделать настолько хитрую операцию. Переносы-то пойдут совсем по-разному, так что максимум что удастся сэкономить — первый слой схемы распространения переноса, и то не факт.
WinPooh73
20.06.2019 19:00+1Для доказательства возможности сэкономить два бита можно даже не рассматривать чётность или нечётность слагаемых в скобках, а также их младшие биты. Достаточно того факта, что всё выражение в скобках обязано нацело делиться на 4, по построению.
ladle Автор
20.06.2019 20:52+1Zenitchik
20.06.2019 21:34Ну почему же?
Если разность a и b нацело делится на 4, следовательно дробные части {a/4} и {b/4} — равны.
А, поскольку, нам нужна разность:
a/4 — b/4 = ([a/4]+{a/4}) — ([b/4]+{b/4}) = ([a/4] — [b/4]) — ({a/4} — {b/4}) = [a/4] — [b/4].ladle Автор
21.06.2019 00:39Подумав, должен признать, что WinPooh73 совершенно прав.
Просто психологически при отбрасывании бит возникает ощущение, что может произойти потеря информации, и желание убедить себя, что этого не происходит. Я это и сделал, рассматривая чётности, примерно как и вавилонские математики 4000 лет назад.
Аргумент WinPooh73 проще и лучше.
Спасибо.
kinhort
21.06.2019 09:06Вот реализация данного алгоритма на языке C. Может, кому пригодится.
#include <stdio.h> typedef unsigned __int8 BYTE; typedef unsigned __int16 WORD; typedef unsigned __int32 DWORD; #define MAX_BYTE 0xFF #define BIAS MAX_BYTE /* Таблица округлённых вниз четвертей квадратов: SqrDiv4[x + BIAS] = [x? /4] Состоит из MAX_BYTE записей для отрицательных значений аргумента, 1 записи для нулевого аргумента, и 2 * MAX_BYTE записей для положительных значений аргумента. */ WORD SqrDiv4[MAX_BYTE + 1 + 2 * MAX_BYTE]; // Инициализация таблицы четвертей квадратов void InitSqrDiv4() { // Используем формулу: // (x + 1)? = x? + 2x + 1 DWORD x = 1; DWORD q = 1; // = x? DWORD Delta = 3; // = 2x + 1 SqrDiv4[BIAS] = 0; for(; x <= MAX_BYTE; x++) { WORD tmp = (WORD)(q >> 2); SqrDiv4[BIAS + x] = tmp; SqrDiv4[BIAS - x] = tmp; q += Delta; Delta += 2; } // Так как макcимальное значение x = 2 * MAX_BYTE, // то максимальное значение q = 4 * MAX_BYTE?, // что займёт 3 байта, и в WORD не поместится, // поэтому q имеет тип DWORD for(; x <= 2 * MAX_BYTE; x++) { SqrDiv4[BIAS + x] = (WORD)(q >> 2); q += Delta; Delta += 2; } } // Эмуляция умножения двух 8-битных чисел DWORD umul8x8to16_emul(DWORD u, DWORD v) { return SqrDiv4[u + v + BIAS] - SqrDiv4[u - v + BIAS]; } int main() { InitSqrDiv4(); // Вычисляем произведение 10 * 15 printf("10 * 15 = %u\n", umul8x8to16_emul(10, 15)); return 0; }kinhort
21.06.2019 11:08Добавлю ещё полезных функций. Везде используется беззнаковое умножение.
#pragma pack(push, 1) union TUInt32 { DWORD UIntPart; BYTE Bytes[4]; struct { BYTE LowByte; union { BYTE HighByte; struct { WORD MidWord; BYTE SignByte; }; }; }; struct { WORD LowWord; WORD HighWord; }; }; #pragma pack(pop) // Эмуляция умножения двух 16-битных чисел DWORD umul16x16to32_emul(DWORD u, DWORD v) { TUInt32 U, V; TUInt32 q; TUInt32 Dest; U.UIntPart = u; V.UIntPart = v; Dest.LowWord = (WORD)umul8x8to16_emul(U.LowByte, V.LowByte); Dest.HighWord = (WORD)umul8x8to16_emul(U.HighByte, V.HighByte); q.UIntPart = 0; q.MidWord = (WORD)umul8x8to16_emul(U.LowByte, V.HighByte); Dest.UIntPart += q.UIntPart; q.MidWord = (WORD)umul8x8to16_emul(U.HighByte, V.LowByte); return Dest.UIntPart + q.UIntPart; } // Вычисление квадрата 8-битного числа DWORD usqr8to16_emul(DWORD u) { return SqrDiv4[2 * u + BIAS]; } // Вычисление квадрата 16-битного числа DWORD usqr16to32_emul(DWORD u) { TUInt32 U; TUInt32 q; TUInt32 Dest; U.UIntPart = u; Dest.LowWord = (WORD)usqr8to16_emul(U.LowByte); Dest.HighWord = (WORD)usqr8to16_emul(U.HighByte); q.UIntPart = 0; q.MidWord = (WORD)umul8x8to16_emul(U.LowByte, U.HighByte); return Dest.UIntPart + q.UIntPart * 2; }kinhort
23.06.2019 15:56Поясню свой код. Для того, чтобы перемножить числа u и v, нам нужно найти сумму u+v и разность u-v. Для 1-байтовых аргументов сумма будет принимать значения в диапазоне [0, 2*255], а разность будет в диапазоне [-255, 255]. Объединив эти диапазоны мы получаем диапазон [-255, 2*255]. Так как язык C поддерживает массивы с индексацией лишь от нуля, то сдвигаем индексацию в этом массиве на 255 с помощью константы BIAS. Затем находим значения элементов массива.
В коде инициализации массива — два цикла. Сначала мы находим значения [(-x)?/4]=[x?/4], а затем, во втором цикле, значения, специфичные только для u+v.
В приведённой здесь функции умножения не нужно сравнивать операнды и обменивать их между собой, так как для отрицательных значений u-v в таблице тоже есть записи. Что касается сложения с константой BIAS, то оптимизатор ещё на этапе компиляции сложит её со смещением массива, и во время исполнения никакого явного сложения с этой константой не будет.
Код использует только сложение, вычитание и сдвиг. Аппаратное умножение не используется. Хотя в примере для вычисления квадрата есть умножение на 2, но оно эквивалентно сдвигу влево на 1.
Приведенные функции оптимизированы для 32-битной архитектуры. Хотя, с формальной точки зрения, функцию умножения можно было бы определить так:
WORD umul8x8to16_emul(BYTE u, BYTE v)
однако, если мы попытаемся использовать это определение на 32-битной архитектуре, то компилятор может начать делать неприятные вещи. Он всё равно будет использовать 32-битные регистры, однако при передаче параметров добавит лишний код вроде: u &= 0xFF, v &= 0xFF, а при возврате тоже сделает return (...) & 0xFFFF; Поэтому я предпочёл использовать DWORD, а тип параметров и результата указать в названии функции. На 16-битной архитектуре это следует скорректировать, соответственно, вместо DWORD написать WORD, кроме инициализации таблицы: там я в комментариях написал, почему.
Функцию умножения 8-битных чисел можно использовать как составную часть функции умножения 16-битных чисел. Я тоже привёл пример этого. Функция umul16x16to32_emul оптимизирована для 32-битной архитектуры. В случае 16-битной архитектуры строка
q.UIntPart = 0;
не нужна. Вместо этого мы складываем umul8x8to16_emul(U.LowByte, V.HighByte) с Dest.MidWord, а возможный перенос от сложения складываем с байтом Dest.SignByte. То же самое делаем и со вторым сложением с umul8x8to16_emul(U.HighByte, V.LowByte).
Такие странные длинные названия объясняются тем, что это — фрагмент кода из другой моей программы. Там ещё было умножение со знаком, которое называлось imul16x16to32.
Для умножения со знаком я первоначально определял знак сомножителей, при необходимости инвертировал знак, чтобы сомножители всегда были беззнаковыми, вызывал функцию беззнакового умножения, затем, при необходимости, инвертировал знак результата.
Это можно сделать без ветвлений, если использовать следующий трюк:
(x+(-1))^(-1)=-x
(x+0)^0=x
То есть, мы создаём переменную s, которая хранит -1, если сомножитель x отрицательный, и 0, если x не отрицательный:
s = (x<0? -1: 0);
Затем вычисляем модуль x:
x=(x+s)^s;
То же самое делаем со вторым сомножителем, знак которого сохраняем в переменную s2.
Затем делаем так:
s3 = s^s2;
После чего s3 содержит -1, если сомножители имеют противоположный знак, и 0, если знак у сомножителей одинаковый.
Выполняем беззнаковое умножение.
У результата умножения r знак корректируем с помощью выражения:
r = (r+s3)^s3;
return r;kinhort
23.06.2019 18:02Выражение s = (x<0? -1: 0); компилятор и сам умеет вычислять без ветвлений. При необходимости для 16-битного x это выражение можно заменить следующим кодом:
TUInt32 tmp; tmp.UIntPart = 0x7FFF; tmp.UIntPart -= x; s = (__int32)(__int16)tmp.HighWord;
AVI-crak
21.06.2019 11:31-1Если сломалось умножение — то алгоритм проверки не должен содержать функции умножения.
А проверка банально выполняется делением.
Vanellope
"При цифрах размером в 1 байт, размер таблицы всего лишь 1022 байт, что может позволить использовать этот трюк в 8-битных процессорах, не имеющих аппаратного умножения."
Выигрываем только в скорости вычислений, т.к. программная эмуляция даже для 32-битных операндов гораздо меньше памяти занимает, а для 8-битных тем более.
ladle Автор
Вы правы, речь как раз и идёт об увеличении скорости при разумных затратах памяти.
Brak0del
Этот способ было бы любопытно прикинуть для ПЛИС по ресурсам и по частотам, например в контексте нейронных сетей. С одной стороны, в ПЛИС обычно валом DSP-слайсов как раз для умножения и сложения. Но с другой, использование блочной памяти (которой там примерно столько же, сколько DSP-слайсов, и каждая объемом 4Кб, т.е. прекрасно влезет такая табличка, а учитывая два независимых порта, к обоим квадратам можно делать одновременное обращение) для умножений таким способом потенциально сможет увеличить количество умножителей вдвое.