Быстрое умножение целых чисел с использованием таблиц

Хочу рассказать читателям о программистском трюке, с которым я познакомился в какой-то переводной книжке, содержащей подборку таких трюков, в те далёкие времена, когда ещё не изобрели не то что байт, а страшно сказать — стек, а великий Дейкстра ещё не проклял оператор GOTO (sic, in uppercase).

Трюк настолько мне понравился простотой и изяществом, что уже в этом тысячелетии я с удовольствием рассказывал о нём студентам в виде следующей задачи.

Представьте себе, что в процессе возвращения в 2030 году с Луны вы вдруг обнаружили, что ваш бортовой компьютер неправильно выполняет операции целочисленного умножения, и это непременно приведёт к аварии при посадке.

В таком сюжете нет ничего особо фантастического. Вспомним, например, какие проблемы случались когда-то с процессорами Pentium , а к моменту отправки на Луну вы ещё не достигли полного импортозамещения. И вообще надо проверить, а не были ли процессоры просверлены специально.

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

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

Вот только если цифры выбрать короткими (например 0 и 1), таблица умножения будет короткой, а столбики длинными, и их вычисление будет занимать много времени.

Если, наоборот, взять цифры длинные (например от 0 до 65535), то для 16-битной арифметики
результат берётся прямо из таблицы, а столбики отсутствуют. Однако размер классической таблицы Пифагора при этом оказывается около 17GB (4*65536*65536), если учесть симметрию относительно диагонали, то половина — около 8.5GB.

Может оказаться многовато.

Напрягаемся и вспоминаем алгебру.

$(a + b)^2 = a^2 + b^2 + 2ab$ (1)

$(a - b)^2 = a^2 + b^2 - 2ab $ (2)

Вычитаем (2) из (1)

$(a + b)^2 - (a - b)^2 = 4ab$

и далее

$a b = ((a + b)^2 - (a - b)^2)/4$

Таким образом, имея таблицу квадратов в массиве 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-битных процессорах, не имеющих аппаратного умножения.

Буду благодарен всем читателям этой заметки за исправления и замечания.
Поделиться публикацией

Комментарии 47

    0

    "При цифрах размером в 1 байт, размер таблицы всего лишь 1022 байт, что может позволить использовать этот трюк в 8-битных процессорах, не имеющих аппаратного умножения."
    Выигрываем только в скорости вычислений, т.к. программная эмуляция даже для 32-битных операндов гораздо меньше памяти занимает, а для 8-битных тем более.

      0
      Вы правы, речь как раз и идёт об увеличении скорости при разумных затратах памяти.
        0
        Этот способ было бы любопытно прикинуть для ПЛИС по ресурсам и по частотам, например в контексте нейронных сетей. С одной стороны, в ПЛИС обычно валом DSP-слайсов как раз для умножения и сложения. Но с другой, использование блочной памяти (которой там примерно столько же, сколько DSP-слайсов, и каждая объемом 4Кб, т.е. прекрасно влезет такая табличка, а учитывая два независимых порта, к обоим квадратам можно делать одновременное обращение) для умножений таким способом потенциально сможет увеличить количество умножителей вдвое.
        0
        Спасибо за статью.
        А можно узнать предпосылки к «вычитаем 2 из 1»?
          0
          Если я правильно понял вопрос,
          то имелась в виду рутинная алгебра — вычитаем формулу (2) из формулы (1).
            0
            Нет, извините за не подробный вопрос. Я имею в виду — откуда возникла мысль вычесть одно выражение из другого и собственно откуда и как появилась идея начать с рассмотрения выражений для квадратов суммы и разности. Предполагаю, существуют какие-то работы или исследования, результатами которых вы и воспользовались. Не могли бы поделиться источниками или ссылками?
              0
              Так я и написал, что прочитал об этом в популярной книжке несколько десятков лет назад.
              Книжка давно утрачена. Я пытался найти ссылки перед написанием заметки, так как опасался наткнуться на статью в Википедии на эту тему и попасть пальцем в небо с актуальностью материала, но в поиске не преуспел, может быть неудачно формировал поисковые запросы. Может быть у кого-то есть ссылки на эту тему, буду благодарен если поделитесь.
                0
                Начальные условия очень похожи на объяснение метода умножения Карацубы.
                  +1
                  В англоязычной википедии это называется «Quarter square multiplication».
                    0
                    Огромное спасибо!
                    Ещё бы ссылку на книжку найти — там было ещё несколько интересных трюков.
              0
              Поскольку возведение в квадрат это тоже операция умножения, то интересно сколько элементов из таблицы можно вычислить из той-же таблицы за одно (а может и несколько) «умножений». Например любой квадрат нечётного числа сведётся к квадратам чётных, может быть половину таблицы можно так выкинуть.
                +1

                Для вычисления таблицы не нужно ни одного умножения, достаточно только сложений.
                (x + 1)^2 = x^2 + 2x + 1
                Соответственно таблицу можно проредить в N раз ценой N-1 дополнительных операций при каждом вычислении квадрата, как пример.

                  0
                  Ну с таким же успехом можно и предыдущее значение найти
                  (x - 1)^2 = x^2 - 2x + 1
                  Это меньше операций, если идти от ближайшего :)
                +1

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

                  +2
                  Современное железо с неисправным умножением скорее всего не сможет даже загрузить ОС. Статья скорее относится к неким бракованным микроконтроллерам, где неисправную инструкцию можно попробовать «обойти», а показатели скорости памяти и вычислительных ядер там могут очень сильно отличаться.
                    +3
                    Причём тут брак и неисправная инструкция? Умножения может тупо не быть на микроконтроллере.
                  0
                  Хорошая статья. Надо напомнить, что, к примеру, очень популярный z80 не имеет аппаратного умножения.
                    0
                    Оптимизация так себе на самом деле. Смотрите, для 8-битных аргументов аппаратное умножение представляет собой 64 AND-операции с последующим сумированием получаемых 64 бит.

                    В вашем случае, перед обращением в таблицу нужно вычислить a+b и a-b, что соответствует суммированию 2*2*8=32 бит. После обращения к таблице нужно еще раз вычислить разницу уже для 16-битных значений, это еще 32 бита которые необходимо просуммировать.

                    Получается в обоих случаях нужно выполнить сумирование 64 бит, но в вашем подходе 64 параллельных AND заменяются на lookup таблицу и 2 обращения к ней – сомнительная оптимизация как по скорости, так и по площади или энергопотреблению.
                      0
                      Хотя я пожалуй соглашусь что такая реализация вычислений может быть эфективной в ситуации когда «работаем с тем что имеем»: есть процессор который умеет только сумирование и не жалко тратить память, которая при этом относительно быстрая. В таком случае действительно можно получить умножение за меньшее количество тактов, заменив кучу сдвигов на 2 обращения к памяти, правда энергопотребление будет скорее всего выше (обращение в память сильно дороже арифметики и работы с регистрами).
                        0
                        Там ещё было сравнение с перестановкой аргументов :)
                        перед обращением в таблицу нужно вычислить a+b и a-b, что соответствует суммированию 2*2*8=32 бит

                        Почему 32, а не 16?
                        a + b = 8 бит + флаг переполнения/переноса
                        a - b = 8 бит (тут наверное «бесплатное» сравнение)
                          0
                          Я считаю количество бит-аргументов, которые нужно сложить, а не количество результирующих бит, по-этому 32. Просто я смотрю на эту оптимизацию как на аппаратную, а не софтварную. :)

                          Кстати, обратил внимание, что даже после срезания 2 битов, оставшиеся младшие 3 бита повторяются с перодом цикла 16. Так что если есть задача уменьшить размер таблицы – можно смело отсекать еще 3 бита от каждого значения и добавлять еще одну таблицу с 16 значениями.

                          Думаю есть еще какие-то паттерны, которые позволят декомпозировать задачу и тем самым уменьшить размер таблицы. Например, 4ый бит повторяется с шагом 32.
                            0
                            А если + и — сделать одной операцией с двумя выходами :)
                            Хитрая операция ab+- где на выходе 16 бит a+b и a-b
                              0

                              Вряд ли можно сделать настолько хитрую операцию. Переносы-то пойдут совсем по-разному, так что максимум что удастся сэкономить — первый слой схемы распространения переноса, и то не факт.

                        +3

                        Для доказательства возможности сэкономить два бита можно даже не рассматривать чётность или нечётность слагаемых в скобках, а также их младшие биты. Достаточно того факта, что всё выражение в скобках обязано нацело делиться на 4, по построению.

                          +1
                          Боюсь, что Вы неправы. Посмотрите ссылку, которую любезно прислал kinhort.
                            0
                            Ну почему же?
                            Если разность 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].
                              +1
                              Подумав, должен признать, что WinPooh73 совершенно прав.
                              Просто психологически при отбрасывании бит возникает ощущение, что может произойти потеря информации, и желание убедить себя, что этого не происходит. Я это и сделал, рассматривая чётности, примерно как и вавилонские математики 4000 лет назад.
                              Аргумент WinPooh73 проще и лучше.
                              Спасибо.
                          –1
                          На языке Lua, потому что с Луны. Вообще, Lua прекрасный язык.
                            0
                            Вот реализация данного алгоритма на языке 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;
                            }
                              0
                              Добавлю ещё полезных функций. Везде используется беззнаковое умножение.
                              #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;
                              }
                              
                                0
                                Поясню свой код. Для того, чтобы перемножить числа 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;
                                  0
                                  Выражение s = (x<0? -1: 0); компилятор и сам умеет вычислять без ветвлений. При необходимости для 16-битного x это выражение можно заменить следующим кодом:
                                  TUInt32 tmp;
                                  tmp.UIntPart = 0x7FFF;
                                  tmp.UIntPart -= x;
                                  s = (__int32)(__int16)tmp.HighWord;
                                  

                                    0
                                    Я упоминал про сложение с возвратом переноса в функции umul16x16to32_emul. Если не использовать ассемблер, то в 16-битной арифметике его можно выполнить с помощью следующей функции:
                                    BYTE add(WORD &Dest, WORD Src)
                                    {
                                        BYTE Carry = (WORD)~Dest < Src;
                                        Dest += Src;
                                        return Carry;
                                    }

                                    Данная функция увеличивает переменную Dest на Src и возвращает 0, если не было переноса в старшие разряды, и 1, если был перенос.
                                    Работает это так: перенос произойдёт, если Dest + Src > 0xFFFF. Преобразуем это выражение:
                                    Carry = 0xFFFF + (-Dest) < Src;
                                    При этом всегда выполняется тождество:
                                    -x = (~x)+1
                                    Следовательно:
                                    Carry = 0xFFFF + (~Dest) + 1 < Src
                                    Результат сложения 0xFFFF + 1 = 0x10000, следовательно, его можно не учитывать, так как нам нужны лишь младшие 16 битов, а они здесь равны 0.
                                    (~Dest) мы приводим к типу WORD, так как компилятор может неявно преобразовывать арифметические выражения к типу int.
                                    Получаем, что перенос:
                                    Carry = (WORD)~Dest < Src;
                                    0
                                    1. а если перед вычислением определять максимальное и его назначать на u и минимальное на v? тогда таблица из отрицательных значений уходит.
                                    2. при равенстве просто выдавать табличный квадрат?
                                      0
                                      Я привёл быструю версию этого алгоритма. Проверки аргументов замедлят программное выполнение, так как добавят лишний код.
                                      Обычно при умножении произвольных чисел равенство множителей возникает редко, а в тех формулах, где есть квадрат числа — лучше функцию для вычисления квадрата и использовать: быстрее передача аргументов и оптимальнее код. Пример функции вычисления квадрата я привёл (usqr8to16_emul).

                                      Вот сравнительные размеры таблиц, чтобы было нагляднее, сколько можно сэкономить, выбросив отрицательные значения.
                                      1. Умножение 1-байтовых чисел, результат умножения — 2-байтовое число.
                                        • Количество записей в таблице со знакопеременной индексацией: 255 + 1 + 2 * 255 = 766
                                        • Количество записей в таблице с неотрицательной индексацией: 1 + 2 * 255 = 511
                                        • Размер таблицы со знакопеременной индексацией: 1 532 байт
                                        • Размер таблицы с неотрицательной индексацией: 1 022 байт
                                        • Разница размеров таблиц: 510 байт
                                      2. Умножение 2-байтовых чисел, результат умножения — 4-байтовое число.
                                        • Количество записей в таблице со знакопеременной индексацией: 65 535 + 1 + 2 * 65 535 = 196 606
                                        • Количество записей в таблице с неотрицательной индексацией: 1 + 2 * 65 535 = 131 071
                                        • Размер таблицы со знакопеременной индексацией: 786 424 байт
                                        • Размер таблицы с неотрицательной индексацией: 524 284 байт
                                        • Разница размеров таблиц: 262 140 байт


                                      Отсюда видно, что разница размеров таблиц при умножении 1-байтовых чисел незначительна, и нет смысла сокращать размер таблицы, замедляя тем самым программу. А для 2-байтовых сомножителей таблица в любом случае получается слишком большой.

                                      Для ясности, как увеличится размер кода при отказе от отрицательной индексации, приведу примеры:
                                      // Эмуляция умножения с использованием ветвлений.
                                      // Отрицательная индексация в таблице не используется.
                                      DWORD umul8x8to16_emul_if(DWORD u, DWORD v)
                                      {
                                          DWORD tmp;
                                          if(u >= v) tmp = u - v;
                                          else tmp = v - u;
                                          return SqrDiv4[u + v + BIAS] - SqrDiv4[tmp + BIAS];
                                      }
                                      
                                      // Эмуляция умножения с использованием модуля разности.
                                      // Отрицательная индексация в таблице не используется.
                                      DWORD umul8x8to16_emul_abs(DWORD u, DWORD v)
                                      {
                                          TUInt32 tmp;
                                          tmp.UIntPart = u - v;
                                          // Если u >= v, то tmp.HighByte = 0
                                          // Если u < v, то tmp.HighByte = 0xFF
                                      
                                          // Вычисляем модуль
                                          tmp.LowByte = (BYTE)((tmp.LowByte + tmp.HighByte) ^ tmp.HighByte);
                                      
                                          // Находим произведение
                                          return SqrDiv4[u + v + BIAS] - SqrDiv4[(DWORD)tmp.LowByte + BIAS];
                                      }
                                      

                                      Как я уже упоминал, сложение с константой BIAS ничего не стоит, так как ещё на этапе компиляции она складывается со смещением начала массива. Здесь сложение с ней приведено для единообразия с предыдущим кодом.
                                        0
                                        Для удобства реализации на ассемблере, распишу умножение 8-битных чисел более подробно:
                                        // Эмуляция умножения двух 8-битных чисел.
                                        // Действия расписаны подробно.
                                        // Предполагается 32-битная адресация.
                                        inline DWORD umul8x8to16_verbose(DWORD u, DWORD v)
                                        {
                                            // Код данной функции эквивалентен следующему коду:
                                            // return SqrDiv4[u + v + BIAS] - SqrDiv4[u - v + BIAS];
                                        
                                            DWORD result;
                                        
                                            u *= sizeof(WORD); // эквивалентно u <<= 1
                                            // Сумма (SqrDiv4 + BIAS * sizeof(WORD))
                                            // вычисляется ещё на этапе компиляции
                                            u += (size_t)SqrDiv4 + BIAS * sizeof(WORD); 
                                            v *= sizeof(WORD); // эквивалентно v <<= 1
                                            u += v;
                                            result = *(WORD *)u;
                                            // Используем следующее тождество: u - v = (u + v) - 2v
                                            v *= 2; // эквивалентно v <<= 1
                                            u -= v;
                                            result -= *(WORD *)u;
                                            return result;
                                        }

                                        Чтобы приблизить код к настоящему, я добавил к атрибутам функции inline, чтобы компилятор точно знал, что эта функция — встраиваемая: для таких функций генерируется более быстрый код.
                                        В моих предыдущих постах в названии функции суффикс "_emul" был сокращением от «emulation» — здесь я решил его не использовать, чтобы не увеличивать название.

                                        Некоторые процессоры позволяют вычислять одной инструкцией выражение вида a * 2 + b, поэтому в данном примере я использовал умножение на 2, а не сдвиг влево на 1. Но, вообще говоря, следующие выражения эквивалентны:
                                        a * 2 + b = (a << 1) + b
                                        Скобки здесь нужны из-за того, что в языке C приоритет оператора сложения выше приоритета оператора сдвига. В языке Pascal, кстати, приоритеты операторов больше соответствуют интуиции:
                                        a * 2 + b = a shl 1 + b
                                        Если же процессор не поддерживает сдвиг, то можно записать так:
                                        a * 2 + b = a + a + b

                                        Таблица четвертей квадратов обычно вычисляется заранее, и подключается к исходному коду в виде константного массива — поэтому при вычислении этой таблицы сдвиг вполне можно использовать, даже если целевой процессор его не поддерживает.
                                –1
                                Если сломалось умножение — то алгоритм проверки не должен содержать функции умножения.
                                А проверка банально выполняется делением.
                                  0
                                  Этот трюк реально применялся на 8080, z80 PDP11 без сопроцессора.
                                    –5
                                    «результат берётся прямо из таблицы, а столбики отсутствуют. Однако размер классической таблицы Пифагора при этом оказывается около 17GB »

                                    — я вытираю лужи слёз от смеха,

                                    «Далее опишу некоторое улучшение, которого в книжке не было. Его я придумал сам, когда уже писал эту заметку. ))))
                                    Заметим, что два младших бита квадрата чётного числа всегда 00, а квадрата нечётного — всегда 01. С другой стороны для любой пары чисел их сумма и разность имеют одинаковую чётность.»
                                    Так это и заметно, что ты сам придумал.))
                                    слушай, прекрати, озёра, реки, да что там, моря слёз от смеха сейчас пролили программисты 60-х, 70-х и особенно 80-х.
                                    проблема только в том, что когда вы это придумываете, (а всё сделано ДО вас), то вы это слишком сильно цените в себе, и не любите здесь на Хабре, когда вам напоминают, что вы вообще-то последний в очереди. А тут ещё есть уважаемые люди. ))

                                    Короче, как вы сами высокомерно выражаетесь по отношению к инженерам электроники, с подключением вас. Вы только что родились на Земле, человеки с Марса )))

                                    Например, когда я здесь писал что таблицы синуса для 32 битной арифметики (да, именно таблица точности 0,4 ppb) достаточно 180 килобайт и она работает быстрее, чем любой сопроцессор, а для практической работы достаточно точности в 0,4 ppm таблица будет размером в 3,6 кбайт, то местные недоразвитые идиотики обсмеяли и выпилили своим любимым минусованием.

                                    Теперь у вас тут все толерантные и слушаются американских ботов бесприкословно ))) Адью, алтернативные. ) Изобретайте дальше свои велосипеды, кастрюли и поварёшки.

                                    PS, У меня давно составлена формула в Экселе, если пожертвовать точностью сильно, т.е. не выпендриваться, а сделать типа Брадиса, что хватает даже для космических полётов (с корректировкой всё равно), то при точности 1 промилле, хватит и таблицы в 72 байта через каждые 5 градусов )) что можно запихнуть вообще в любой современный проц прямо в программную память или даже в регистры.
                                    И это будет намнооого быстрее сопроцессора с плавающей точкой. А ты про умножение. ) Тут трансцедентная функция.
                                      0
                                      если пожертвовать точностью сильно, т.е. не выпендриваться, а сделать типа Брадиса, что хватает даже для космических полётов

                                      А специалисты по МДТТ и не знают, им даже float нужной точности не всегда даёт, не говоря уж о "таблицах в 72 байта"...

                                        0
                                        Советую Вам познакомиться с rybvv. Вы с ним явно найдете общий язык.
                                          0

                                          Кажется, Вы упомянутому товарищу устроили хабраэффект :)

                                          0
                                          Челябинские радиоинженеры настолько суровы…
                                            –4
                                            размер классической таблицы Пифагора при этом оказывается около 17GB
                                            Хе-хе… Наверное тогда этот тип женоподобный, как типичный программист, не смогший бороться с системой. «Ой не смогла, напрудила 17 GB.»
                                            В этом и разница между профи в 72 байта и 17 ГБ — от «альтернативно одарённых» (кстати эти слова, я подхватил у вас, на Хабре, и применяю их исключительно про вас, программеров, кесарю — кесарево, а слесарю — слесарево ))) хе-хе…
                                            «17 Гб на таблицу Пифагора» — это мне кажется, что некоторые из вас, программистов, уже не просто альтернативные дурачки для нормальных людей, а уже откровенно идиоты, нуждающиеся в психиатрическом лечении.

                                            Это я ещё не говорю про статистику по национальностям, которую вы тут боитесь ), по программерам. А она есть у меня. )) Пруф такой неплоххой про псих -заболевания, медицинский факт вполне.

                                            2 Cerberuser — «специалисты по МДТТ»
                                            Вы откуда свалились вообще, Вы понимаете для чего таблицы делают в системах? Судя по вашему высказыванию, нет.
                                            Их делают для скорости вычислений в реальном времени, когда времени на сопроцессор и ожидание — нет совсем и из железа выжимают максимум. А ваши конструкторские вычисления подождут в тишине. Их делают намного заранее до боя или действия, «до экшона». Вы если не понимаете, то не встревайте в тему, пожалуйста.
                                            Туда же — каждая пичужка норовит уколоть и высказаться.
                                            А специалисты по МДТТ и не знают, им даже float нужной точности не всегда даёт, не говоря уж о «таблицах в 72 байта»...
                                            писец какой-то. у него есть чувство сарказма. На Хабре? )) неее, не может быть… показалось… ))
                                            90 градусов через 5 — это 18 значений по 32 бита. Ноль вычислять не нужно, если вы не знали )) Как малолетние дети, ей-богу, сами восстановить ничего не могут. Реверс-инжиниринг — это удел богов наверное )
                                              0

                                              Так, может быть, расскажете нам, как же правильно сделать таблицу 65536 на 65535 с 32х-разрядными значениями, чтобы ее размер был меньше 16 ГиБ? (17 Гб — невнимательность автора, там их ровно 16)


                                              Лично я вижу только 1 вариант: не делать эту таблицу. Собственно, критикуемый вами автор именно так и поступил.

                                                0
                                                Их делают для скорости вычислений в реальном времени, когда времени на сопроцессор и ожидание — нет совсем и из железа выжимают максимум

                                                Прошу прощения. Из Вашего комментария было, мягко говоря, не очевидно (по крайне мере, не-специалисту в Вашей области), что Вы рассматриваете только частный случай, а не общеприменимый метод.


                                                Остальное лучше проигнорировать, а то ещё закончите, как автор по ссылке из комментария выше. Не хотелось бы.

                                              0
                                              А тут ещё есть уважаемые люди. ))
                                              Есть. Но вы к ним явно не относитесть.

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

                                            Самое читаемое