Возведение в квадрат в компьютерах.

[yk2ru]

Возведение в квадрат

$$(A_{1} + 2^{k}A_{2})^2 = A_{1}^2 + 2^{k+1}A_{1}A_{2} + 2^{2k}A_{2}^2$$

сводится к трем умножениям чисел половинной длины.

Метод умножения Карацубы

$$A=A_{1}+2^{k}A_{2},\quad B=B_{1}+2^{k}B$$

$$AB=A_{1}B_{1}+2^{k}{\Big (}(A_{1}+A_{2})(B_{1}+B_{2})-(A_{1}B_{1}+A_{2}B_{2}){\Big )}+2^{2k}A_{2}B_{2}$$

состоит в том, что четыре умножения сводятся к трём.

Правильно ли я понимаю, что для процессора не программируют возведение в квадрат, а используют тот же метод, что для перемножения разных чисел? Для них же используют метод умножения Карацубы.

Используется синтаксис Pascal

function km(a:integer; b:integer):integer;#karatsuba method
var l1,l2,r,d,d2,al,ar,bl,br:integer;
begin
 l1:=bit_length(a);
 if l1<2 then return a*b; end;
 l2:=bit_length(b);
 if l2<2 then return a*b; end;
 d:=max(l1,l2);
 d2:=bit_shift(d,-1);
 al:=bit_shift(a,-d2);
 ar:=a-bit_shift(al,d2);
 bl:=bit_shift(b,-d2);
 br:=b-bit_shift(bl,d2);
 l1:=km(al, bl);
 l2:=km(ar, br);
 r:=km(al+ar, bl+br);
 r := bit_shift((r-(l1+l2)),d2)+bit_shift(l1,bit_shift(d2,1))+l2;
 return r;
end.

Подумал и запрограммировал некий способ возведения в квадрат в функцию sq(n). И получилось, что он работает заметно быстрее. Обе функции дают одинаковый результат в результате исполнения.

for n:=1 to 1000000 do
 if sq(n)<>km(n,n) then
  writeln(n," sq<>km");
 end;
end.

Вывода строки "sq<>km" тут не случилось.
Результаты:

==> t := timer();m:=0;
for n:=1 to 1000 do
 m:=2*m+1;
 s:=km(m,m);
end;
writeln(timer()-t).
80414
-: 1

==> t := timer();m:=0;
for n:=1 to 1000 do
 m:=2*m+1;
 s:=sq(m);
end;
writeln(timer()-t).
550
-: 1

==> 80414/550.
-: 146.207273

При перемножении чисел Мерсенна более чем в 100 раз быстрее.

==>t := timer();
for n:=1 to 1000000 do
 s:=km(n,n);
end;
writeln(timer()-t).
226362
-: 1

==> t := timer();
for n:=1 to 1000000 do
 s:=sq(n);
end;
writeln(timer()-t).
10344
-: 1

==> 226362/10344.
-: 21.8834107

При перемножении последовательных чисел более чем в 20 раз быстрее.
Код и проверка сделаны на Aribas ( http://progopedia.ru/language/aribas/ ).
Возможно, что неправильно запрограммировал метод умножения Карацубы и поэтому разрыв в скорости.
Предлагаю желающим попробовать сочинить две функции и повторить эти опыты с возведением в квадрат в любой из систем. Удастся ли вам повторить эти результаты?

[slavav]

Возведение в квадрат и умножение должны отличаться по времени исполнения не более чем в константное число раз. Ваши примеры очень маленькие. Попробуйте числа вида $10^n$, где $n$ достаточно велико чтобы время вычисления было порядка секунды. Варьируя $n$ получите несколько времён, затем восстановите регрессией константу $C$ и степень $a$ в формуле $Cn^a$.

integer тип из ARIBAS может оказаться неудачным контейнером для чисел. Алгоритм Карацубы рассчитан на работу с массивами ограниченных целых. Вы используете integer, следовательно вам надо проверять что элементарные операции (умножение на степень двойки, сложение) у вас имеют ту же сложность что и операции над массивами бит/ограниченных целых.

[yk2ru]

На чем поэкспериментировать, может Питон, там вроде бы целый тип не ограничен? В Арибас есть ограничение по целым числам.

[slavav]

Питон подойдёт, там нормальные неограниченные целые.
В Арибас я прочитал что числа не ограничены. Ошибся?

[12d3]

В питоне, а точнее, в cpython, уже используется алгоритм Карацубы при умножении больших целых чисел. Реализуйте лучше сами класс BigInteger, и на нем сравнивайте умножения.

[mihaild]

Тут вопрос в том, какая у вас реализация sq. Для одиночных чисел, квадрат которых влезает в регистры (скорее всего до $2^31 - 1$ в вашей реализации) ничего быстрее прямого умножения (которое скомпилируется в одну инструкцию) сделать не удастся. Для многих чисел нужно умножать пачками (и придется руками перемножать отдельные 32 битные части). Никакие софтверные реализации тут в любом случае не помогут.

Внутри для умножения используются схемы логарифмической глубины (зато с большим числом элементов), так что "асимптотика" умножения двух $n$-битных чисел - $O(\log n)$. Разумеется это утверждение особого смысла не имеет, т.к. на практике у нас одна схема, а не набор схем под каждую разрядность.

Код имеет смысл писать для перемножения длинных чисел (не влезающих в регистры). И тут какие-то оптимизации для возведения в квадрат по сравнению с умножением Карацубы могут иметь смысл. Но т.к. обычное умножение векторизуется сильно лучше, чем умножение Карацубы, скорее всего есть смысл сделать несколько итераций умножения Карацубы, а потом обычное.
Для достаточно длинных чисел (сотни тысяч разрядов) вообще лучше брать метод Шёнхаге-Штрассена (к сожалению его реализовать на "верхних уровнях" с переключением на обычное умножение внизу нельзя).

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

[slavav]

mihaild, в Питоне не лучшая реализация длинной арифметики. Там Карацуба, а не Шёнхаге-Штрассен, а GNU MP они оба. Печатает Питон длинные числа вообще за квадрат.
Как я понял ТС, он будет использовать длинные целые только для хранения, сложения, умножения на степень двойки, умножения произвольных чисел с ограниченной длиной в битах. При должной аккуратности можно смоделировать Карацубу довольно точно. Идея ведь реализовать возведение в квадрат и умножение и сравнить их О-большие, не более.

[Munin]

yk2ru в сообщении #1394590 писал(а):

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

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

Обычное дело - инструкция извлечения квадратного корня. Возведения в целочисленную степень не бывает. Странно, но не делают и извлечения кубического корня. Вообще возведение в произвольную степень реализуется через трансцендентные функции $2^x$ и $\operatorname{lb}(x)=\log_2 x.$ На трансцендентных функциях экономят (по числу функций, а не по скорости выполнения, при том что они обычно медленные): часто встречается набор $\sin,\cos,\operatorname{atan}_2,$ а остальные - делайте из них сами. (Правда, этот выбор может быть обусловлен и тем, что это "самые хорошие" функции по области определения и всяким особенностям.)

[yk2ru]

Код sq не содержит рекурсий, возможно поэтому работал быстрее. На длинных числах будет в два раза быстрее обычного умножения столбиком, как я понимаю. Возможна оптимизация, не на каждом нуле устраивая сдвиги, а потом на единице скопом. И в конце только результат вычислять.

Используется синтаксис Pascal

function sq(a:integer):integer;
var al,b,n,r:integer;
begin
 if a<2 then
  return a*a;
 end;
 r:=1; al:=1;
 n:=bit_length(a)-1;
 while n>0 do
  dec(n);
  b:=bit_test(a,n);
  if b=0 then
   al:=bit_shift(al,1);
   r:=bit_shift(r,2);
  else
   r:=1+bit_shift(r+al,2);
   al:=1+bit_shift(al,1);
  end;
 end;
 return r;
end.

[yk2ru]

Ещё раз взглянув на формулу вначале, на которой основан метод Карацубы, возник вопрос: почему там используются суммы, а не разности? Можно же через произведение разностей записать ту формулу.

yk2ru в сообщении #1394590 писал(а):

$$A=A_{1}+2^{k}A_{2},\quad B=B_{1}+2^{k}B_{2}$$

$$AB=A_{1}B_{1}+2^{k}{\Big (}A_{1}B_{1}+A_{2}B_{2}-(A_{1}-A_{2})(B_{1}-B_{2}){\Big )}+2^{2k}A_{2}B_{2}$$

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

Используется синтаксис Pascal

function kd(a:integer; b:integer):integer;#karatsuba method (differences)
var l1,l2,r,d,d2,al,ar,bl,br:integer;
begin
 l1:=bit_length(a);
 if l1<2 then return a*b; end;
 l2:=bit_length(b);
 if l2<2 then return a*b; end;
 d:=max(l1,l2);
 d2:=bit_shift(d,-1);
 al:=bit_shift(a,-d2);
 ar:=a-bit_shift(al,d2);
 bl:=bit_shift(b,-d2);
 br:=b-bit_shift(bl,d2);
 l1:=km(al, bl);
 l2:=km(ar, br);
 r:=km(ar-al, br-bl);
 r := bit_shift(l1+l2-r,d2)+bit_shift(l1,bit_shift(d2,1))+l2;
 return r;
end.

Примеры с квадратами чисел.

t := timer();
for n:=1 to 1000000 do
s:=kd(n,n);
end;
writeln(timer()-t).
191593
-: 1

==> t := timer();
for n:=1 to 1000000 do
s:=km(n,n);
end;
writeln(timer()-t).
224869
-: 1

==> (224869-191593)/224869*100.
-: 14.7979490

==> 224869/191593.
-: 1.17368067

С разностями отрабатывается быстрее на 15-17 процентов.

[feedinglight]

беззнаковые сложнее вычитать

[yk2ru]

Юрий Матиясевич дает объяснение метода именно через произведение разностей.
Быстрая арифметика ( http://book.etudes.ru/toc/fast-arithmetic/ )

[Munin]

feedinglight в сообщении #1396495 писал(а):

беззнаковые сложнее вычитать

Да вроде нет, с чего бы?

[feedinglight]

Munin
Если $A_{1} - A_{2} < 0$, то плохо же будет, не? Ну или делать так чтобы минусов не было тогда, но это тяжеловасто было бы.

(Оффтоп)

честно говоря вообще не понимаю о чём эта тема

[Munin]

Ну получится число в дополнительном коде. Потом оно просто сработает как отрицательное. (Да, вы правы, тут нужна бо́льшая аккуратность и внимательность.)