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

[yk2ru]

С произведением разностей функция выдает те же результаты как и с суммами, никаких аккуратностей не потребовалось. Или тут квадраты сказываются?
Проверка:

p:=10000000;
for n:= 1 to 100000 do
 a:=random(p); b:=random(p);
 if km(a,b)<>kd(a,b) then
  writeln("km<>kd ",a," ",b);
 end;
end.

Никакого вывода не последовало, функции выдают одинаковый результат.

[feedinglight]

Проверил(проще чем думать), действительно unsigned int автоматом работает как нужно .

(проверка)

    unsigned int a1, a2, b1, b2, res;
    a1 = 10; a2 = 20;
    b1 = 100;b2 = 200;
    res = (a1 - a2)*(b1 - b2);
    printf("res = %u", res);  // res = 1000

[Munin]

Угу, баг, про который известно. У тега syntax есть своё "сворачивание-разворачивание", так что использовать off не требуется.

[feedinglight]

Ок, понятно. Убрал кракозябры.

[yk2ru]

Продолжение опытов. Используется данное в первом сообщении разложение в квадрат.

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

Далее нужно посчитать три произведения $A_{1}^2, A_{1}A_{2}, A_{2}^2$ . Для этого определяется, какое из чисел меньше. Пусть меньше первое число. Положим положительную разность этих чисел как $D$.
Для произведений трех чисел получаются такие формулы:

$$A_{1}A_{2} = A_{1}^2 + DA_{1} , A_{2}^2 = A_{1}A_{2} + DA_{1} + D^2$$

Если $A_{1}, A_{2}$ различаются значительно, то меньшее число можно домножить на нужную степень двойки и тогда разность станет маленькой и его квадрат и умножение с ним будут вычисляться быстрее. Последние формулы будут содержать поправки, множители со степенью два. Далее можно продолжить такую же рекурсию для нахождения уже произведений $A_{1}^2, DA_{1}, D^2$ . Сочинил следующий код:

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

function ab2(a:integer; b:integer; var a2:integer; var ab:integer; var b2:integer);
var ca,cb,la,lb,d,m,d2:integer;
begin
 ca:=bit_count(a);
 cb:=bit_count(b);
 if ca<2 then
  a2:=a*a;ab:=a*b;b2:=sq(b);
 elsif cb<2 then
  a2:=sq(a);ab:=a*b;b2:=b*b;
 end;
 ca:=0;cb:=0;
 while ca<la do
  if bit_test(a,ca)=1 then
   break;
  end;
  inc(ca);
 end;
 dec(la,ca);
 if la<2 then
   a2:=a*a;ab:=a*b;b2:=sq(b);
  return;
 end;
 while cb<lb do
  if bit_test(b,cb)=1 then
   break;
  end;
  inc(cb);
 end;
 dec(lb,cb);
 if lb<2 then
   a2:=sq(a);
   ab:=a*b;
   b2:=b*b;
  return;
 end;
 if a=b then
  a2:=sq(a);ab:=a2;b2:=a2;
  return;
 end;
 a:=bit_shift(a,-ca);
 b:=bit_shift(b,-cb);
 if a<b then
  la:=lb-la-1;
  if la>0 then
   d:=b-bit_shift(a,la)
  else
   d:=b-a;
  end;
  ab2(a,d,a2,m,d2);
  if la>0 then
   ab:=bit_shift(a2,la)+m;
   b2:=bit_shift(ab+m,la)+d2
  else
   ab:=a2+m;
   b2:=ab+m+d2;
  end;
 else
  lb:=la-lb-1;
  if lb>0 then
   d:=a-bit_shift(b,lb)
  else
   d:=a-b;
  end;
  ab2(b,d,b2,m,d2);
  if lb>0 then
   ab:=bit_shift(b2,lb)+m;
   a2:=bit_shift(ab+m,lb)+d2
  else
   ab:=b2+m;
   a2:=ab+m+d2;
  end;
 end;
 if ca>0 then
  a2:=bit_shift(a,bit_shift(ca,1));
 end;
 if cb>0 then
  b2:=bit_shift(a,bit_shift(cb,1));
 end;
 m:=ca+cb;
 if m>0 then
  ab:=bit_shift(ab,m);
 end;
end.

function sqm1(a:integer):integer; # модификация №1 нахождения квадратов
var al,ar,a2,b2,ab,n:integer;
begin
 n:=bit_count(a);
 if n<2 then return a*a end;
 n:=bit_shift(bit_length(a),-1);
 al:=bit_shift(a,-n);
 ar:=a-bit_shift(al,n);
 ab2(al,ar,a2,ab,b2);
 return b2+bit_shift(ab,n+1)+bit_shift(a2,bit_shift(n,1));
end.

Также имеются рекурсии, как и в методе карацубы. Но чем длиннее числа, тем метод Карацубы все больше отстает.

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

==> p:=10**10;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 km(m,m);
end;
writeln(timer()-t).
27886
-: 1

==> p:=10**10;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 sqm1(m);
end;
writeln(timer()-t).
1190
-: 1

==> 27886/1190.
-: 23.4336134

==> p:=10**20;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 km(m,m);
end;
writeln(timer()-t).
82109
-: 1

==> p:=10**20;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 sqm1(m);
end;
writeln(timer()-t).
1462
-: 1

==> 82109/1462.
-: 56.1621067

==> p:=10**30;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 km(m,m);
end;
writeln(timer()-t).
153042
-: 1

==> p:=10**30;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 sqm1(m);
end;
writeln(timer()-t).
1675
-: 1

==> 153042/1675.
-: 91.3683582

==> p:=10**40;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 km(m,m);
end;
writeln(timer()-t).
243611
-: 1

==> p:=10**40;t:=timer();
for n:=1 to 100000 do
 m:=random(p)+p;
 sqm1(m);
end;
writeln(timer()-t).
1818
-: 1

==> 243611/1818.
-: 133.999450

[yk2ru]

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

[yk2ru]

В реализацию функции ab2 была ошибка, не были инициализированы до первого использования переменные la,lb. Поэтому числа там перемножались просто обычным способом без всяких рекурсий. С поправкой возведение в квадрат работает всё равно побыстрее, но не в 100 с лишним раз, а в примерно в 20 с лишним раз на числах порядка $10^{40}$ и в 10 раз на числах порядка $10^{10}$.

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

function ab2(a:integer; b:integer; var a2:integer; var ab:integer; var b2:integer);
var ca,cb,la,lb,d,m,d2:integer;
begin
 ca:=bit_count(a);
 cb:=bit_count(b);
 la:=bit_length(a);
 lb:=bit_length(b);
 if ca<2 then
  if a>0 then
   dec(la);
   a2:=bit_shift(1,bit_shift(la,1));
   ab:=bit_shift(b,la);b2:=sq(b);
  else
   a2:=0;ab:=0;b2:=sq(b);
  end;
  return;
 elsif cb<2 then
  if b>0 then
   dec(lb);
   a2:=sq(a);ab:=bit_shift(a,lb);
   b2:=bit_shift(1,bit_shift(lb,1))
  else
   a2:=sq(a);ab:=0;b2:=0;
  end;
  return;
 end;
 ca:=0;cb:=0;
 while ca<la do
  if bit_test(a,ca)=1 then
   break;
  end;
  inc(ca);
 end;
 dec(la,ca);
 if la<2 then
   a2:=a*a;ab:=a*b;b2:=sq(b);
  return;
 end;
 while cb<lb do
  if bit_test(b,cb)=1 then
   break;
  end;
  inc(cb);
 end;
 dec(lb,cb);
 if lb<2 then
   a2:=sq(a);
   ab:=a*b;
   b2:=b*b;
  return;
 end;
 a:=bit_shift(a,-ca);
 b:=bit_shift(b,-cb);
 if a=b then
  a2:=sq(a);ab:=a2;b2:=a2;
 elsif a<b then
  la:=lb-la-1;
  if la>0 then
   d:=b-bit_shift(a,la)
  else
   d:=b-a;
  end;
  ab2(d,a,d2,m,a2);
  if la>0 then
   ab:=bit_shift(a2,la)+m;
   b2:=bit_shift(ab+m,la)+d2
  else
   ab:=a2+m;
   b2:=ab+m+d2;
  end;
 else
  lb:=la-lb-1;
  if lb>0 then
   d:=a-bit_shift(b,lb)
  else
   d:=a-b;
  end;
  ab2(d,b,d2,m,b2);
  if lb>0 then
   ab:=bit_shift(b2,lb)+m;
   a2:=bit_shift(ab+m,lb)+d2
  else
   ab:=b2+m;
   a2:=ab+m+d2;
  end;
 end;
 if ca>0 then
  a2:=bit_shift(a2,bit_shift(ca,1));
 end;
 if cb>0 then
  b2:=bit_shift(b2,bit_shift(cb,1));
 end;
 m:=ca+cb;
 if m>0 then
  ab:=bit_shift(ab,m);
 end;
end.

[FomaNeverov]

А сколько по времени занимает возведение в квадрат числа порядка $10^{40}$ ?

[Dmitriy40]

PARI/GP, который ну очень и очень не быстрый, справляется за пару мкс (2с на миллион итераций):

? t=10^42\71;t2=0;for(i=1,10^6,t2+=(t+113*i^3)^2);
  ***   last result computed in 1,873 ms.
? t=10^402\71;t2=0;for(i=1,10^6,t2+=(t+113*i^3)^2);
  ***   last result computed in 2,825 ms.
? t=10^4002\71;t2=0;for(i=1,10^6,t2+=(t+113*i^3)^2);
  ***   last result computed in 20,405 ms.

На С/С++/C# должно быть на порядок быстрее.
Число слишком короткое, всего 130 битов или 3-5 слов, алгоритму негде разогнаться, вот для чисел порядка $10^{4000}$ рост времени уже заметен.

[yk2ru]

На С# и C++ для замера можно воспользоваться System.Numerics.BigInteger .
Длинная арифметика от Microsoft - https://habr.com/ru/post/207754/

[yk2ru]

Сравнение на питоне внутри функции по вычислению чисел Фибоначчи, попрактиковался на новом языке программирования заодно. В функции sq на входе два числа Фибоначчи и вычисляется квадрат их суммы. Метод, который связан с разностями чисел - post1397310.html#p1397310 . Легко программируется без рекурсий, поэтому наверное и побыстрее отрабатывал.

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

import math
import timeit

def sq(x,y):
    if y<1073741824:
        xy = x+y
        return xy*xy
    k = 0
    while y>1073741824:
        d = y - x
        y = x
        x = d
        k += 1
    x2 = x*x
    xy = x*y
    y2 = y*y
    while k>0:
        k -= 1
        dx = xy
        xy = y2 + dx
        d2 = y2
        y2 = xy + dx + x2
        x2 = d2
    r = x2+2*xy+y2
    return r
    
def karatsuba(x, y):
    if x < 1073741824 or y < 1073741824:
        return x*y

    #making the two numbers strings here
    str_x = str(x)
    str_y = str(y)

    #finding the mid point for each number here
    n = max(len(str(x)), len(str(y)))
    n_2 = int(n / 2)

    #higher bits of each number
    x_h = int(str_x[:n_2])
    y_h = int(str_y[:n_2])

    #lower bits for each number here
    x_l = int(str_x[n_2:])
    y_l = int(str_y[n_2:])

    a = karatsuba(x_h, y_h)
    d = karatsuba(x_l, y_l)
    e = karatsuba(x_h + x_l, y_h + y_l) - a - d

    return a*10**len(str_x) + e*10**(len(str_x) // 2) + d

def fibonacci(n):
    n2 = 2**100
    k = n2
    l = 100
    a = 0
    b = 1
    if n < 0: 
        print("Incorrect input") 
    elif n == 0: 
        return a 
    elif n == 1: 
        return b 
    else: 
        for i in range(2,n): 
            c = a + b 
            if a>k:
                k = k * n2
                t1 = timeit.default_timer()
                r1 = sq(a,b)
                t2 = timeit.default_timer()
                r2 = karatsuba(c,c)
                t3 = timeit.default_timer()
                print('%s: %s %s' % (l, t2-t1, t3-t2))
                l = l + 100
            a = b 
            b = c
        return b 

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

38500: 0.9804285950012854 1.8227088320018083
38600: 0.9961100830005307 1.8862310869990324
38700: 1.0019839939996018 2.267042440998921
38800: 1.3236269150002045 2.446670182998787
38900: 1.5489434560004156 2.3210190220015647
39000: 1.037292652999895 1.9048197089978203
39100: 1.5514899580011843 2.8959486429994286
39200: 1.0976877459979733 2.788560734999919
39300: 1.0404055819999485 1.9331318080003257
39400: 1.286475517998042 2.1424040460005926
39500: 1.298370137999882 3.030312638999021
39600: 2.335647541000071 2.772786024001107
39700: 3.8512898500011943 5.150986110998929
39800: 4.845729540997127 6.39639591300147
39900: 2.698005246998946 4.039674425999692
40000: 1.2129209699996863 2.1693502519992762
40100: 1.0955858200031798 2.9719087069970556
40200: 1.539383032999467 3.0746981619995495
40300: 1.1490264400017622 1.90916350200132
40400: 1.1330688980015111 1.9204274809999333
40500: 1.386953158002143 2.1585519689979265
40600: 1.1082464960018115 1.9867081789998338
40700: 1.1234650330006843 1.8799448360005044
40800: 1.76749957499851 2.1290064200002234
40900: 1.3289487290021498 2.674667433999275
41000: 1.0970706850021088 1.89014380899971
41100: 1.1742691679974087 1.923786764000397
41200: 1.2261482760004583 1.9756421630008845
41300: 1.2154652760000317 2.4787701009990997
41400: 1.3651241400002618 2.261001570997905
41500: 1.3201638459977403 1.9573589890023868
41600: 1.1669875270017656 1.9735675519987126
41700: 1.17371635499876 2.2550408839997544
41800: 1.263370140997722 3.0805593370023416
41900: 1.3819119850013521 2.04396842999995
42000: 1.2335806010014494 1.9092572100016696
42100: 1.2023409360008372 1.8780054380004003
42200: 1.1839886919988203 1.9098548500005563
42300: 1.3153181720008433 2.1218224769982044
42400: 1.2137951189979503 1.8586900199989032
42500: 1.2235551120029413 2.448624757998914
42600: 1.2095067420013947 1.9823751609983447
42700: 1.2054786320004496 1.9435984000010649
42800: 1.2420458490014425 1.8862686590000521
42900: 1.2291398950001167 2.2318529859985574
43000: 1.4651235760029522 1.9512128829992434
43100: 1.519923595002183 1.9853905579984712
43200: 1.2577896709990455 2.0976376489998074
43300: 1.2616500900003302 1.959719459999178
43400: 1.2820497730026545 1.9293508889968507
43500: 1.2838360550013022 2.071842208999442

[slavav]

В Питоне вызов str(n) выполняется за квадрат длины числа (или длины результирующей строки или числа бит в числе). Это делает вашу реализацию метода Карацубы медленнее чем она должна быть.
Чтобы разложить число в массив байт вам нужно сделать такую штуку: n.to_bytes((n.bit_length() + 7) // 8, 'little').

[yk2ru]

Код метода брал из интернета, все используют строки. Например тут http://www.marinamele.com/third-grade-karatsuba-multiplication-algorithms
Не нашел нормальной реализации.

[slavav]

В статье, как я её понял, не ставится цель получить высокую производительность в целом. Подсчитывается число элементарных сложений и умножений, другие операции не учитываются.
Питон хранит числа в двоичном представлении. Вызов str преобразует его в десятичное представление. Это нетривиальная операция. Исходный код отсылает к Кнуту, а там приведена схема Горнера для разрядов чисел, которая требует квадратичного времени.
Схему Горнера можно улучшить, но есть путь короче. Числа в Питоне хранятся как массивы байт, получить такой массив из числа (фактически его копию) можно за линейное время. Умножение Карацубы применяется тогда к числам записанным по основанию 256. Операции ведутся над массивами байт.
Не уверен что вы найдёте такой код в Сети. Во-первых сам Питон так умножает числа только гораздо быстрее. Во-вторых умножение двоичных чисел плохо подходит для учебного примера.

[arseniiv]

slavav в сообщении #1433739 писал(а):

не ставится цель получить высокую производительность в целом

Кстати если у ТС ставится (не очень понимаю его цели, хотя в теме уже и две страницы), следует или подумать насчёт использования Cython или PyPy, или насчёт другого языка. А если хочется посчитать число операций по-честному, то надо вместо вычисления сложений и умножений увеличивать какой-нибудь счётчик операций, и сами числа притом вообще не нужны, а нужны только их длины.

-- Вт янв 07, 2020 16:28:26 --

slavav в сообщении #1433739 писал(а):

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

Это тоже надо для ТС подчеркнуть: пользуясь любой библиотекой работы с большими числами (или языком, использующим какую-то такую реализацию под капотом), не надо велосипедить код для вычисления операций, которые предоставляются библиотекой. Лучше скорее всего не выйдет, или выйдет, но это будет бессмысленной тратой сил, потому что пользы будет меньше, чем усилий на корректную реализацию и её возможную поддержку.