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

[yk2ru]

Такой вариант умножения нормальный или возможно ускорить его? Наверное нужны сдвиги, а не деление нацело и (по модулю значение)

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

def mult(x,y):
    xl = (x.bit_length()+7) // 8
    yl = (y.bit_length()+7) // 8
    if xl < 2 or yl < 2:
        return x*y
    n = max(xl, yl)
    m = n // 2
    m256 = 256**m
    x_L = x % m256
    x_H  = x // m256
    y_L = y % m256
    y_H = y // m256
    a = mult(x_H,y_H)
    d = mult(x_L,y_L)
    e = mult(x_H + x_L, y_H + y_L) - a - d
    return a*(m256**2) + e*m256 + d

[arseniiv]

Так вы всё-таки хотите перебороть встроенный алгоритм умножения?

[slavav]

Да, нужны сдвиги вместо умножений и работа с байтами вместо делений.
Вам нужно извлечь байты: x_bytes = n.to_bytes(xl, byteorder='little').
Деление на $256^m$, взятие остатка от деления на $256^m$ и умножение на $256^m$ тогда делаются как операции над массивами байт.
Останется сложение и вычитание. Их предлагаю делать не на массивах а на самих числах. Массив байт в число переделывается так:
x = int.from_bytes(x_bytes, byteorder='little').
Преобразования число->байты и байты->число делаются за логарифм числа. Сложение/вычитание тоже. Это то что вам нужно.

[slavav]

arseniiv в сообщении #1433811 писал(а):

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

В поддержку ТС: Python делает преобразования целого в строку и наоборот за квадрат. Не боги горшки обжигают.

[arseniiv]

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

[slavav]

Умножение, возведение в квадрат и преобразование в другую систему счисления довольно близки.
Я решал задачу Перельмана про цифры $9^{9^9}$. Оказалось что Питону требуется несколько минут чтобы вычислить результат и месяц чтобы его напечатать. Это пример, который демонстрирует что не следует слепо доверять реализации очевидных вещей в очень популярном языке программирования.

[arseniiv]

(Оффтоп)

Просто преобразование чисел в строку и назад, насколько мне казалось, обычно не входит в используемые языками библиотеки длинной арифметики типа GMP (забыл, что используется питоном) (потому что это не часть собственно арифметики, и клиенты могут по-разному работать со строками, так что проще этого в библиотеку не включать), так что полагаться на арифметику и полагаться на преобразование в строку — разные вещи. Да, хорошо бы чтобы язык для любого случая давал достаточно оптимальный код, но понятно, что ни один этого по крайней мере пока сделать не может, потому что нужны человекочасы.

А вообще смысла нет спорить. Да, проверять полезно, и т. д. и т. п., но для этого нужна квалификация, которой у некоторых людей по-моему нету, раз они взялись писать умножение длинных чисел на питоне и не поминали ни PyPy, не Cython, ничего, на чистом питоне. Пожелаем им удачи за неимением других слов.

[B@R5uk]

Наткнулся я на эту тему, и показалась она мне интересной. В смысле, про метод Карацубы я давно слышал, это классический пример подхода "разделяй и властвуй", но никогда самому не доводилось его реализовывать. Решено: буду писать код, благо формула одна и очень простая. Оказалось быстро сказка сказывается, да долго дело делается. Провозился с отладкой и вылавливанием багов целую неделю. Плюс долго соображал как же лучше реализовать разные мелочи. Например, что взять за базовый разряд (я решил свою длинную арифметику делать в системе счисления с основанием $10^9$, это почти int в Java, и удвоенное значение не превышает величины $2^{31}$), как выделять память, чтобы не тратить лишку и всегда хватало на переполнение (я разрешил иметь ведущий ноль своим длинным числам), когда останавливаться расчленять множители по методу Карацубы и переходить на умножение в столбик (я выбрал пороговую величину 25). Я не стал заморачиваться с отрицательными числами, за то пришлось реализовать кучу вспомогательных функций типа печати числа, превращение строки в число, сравнение двух чисел (для проверки правильности счёта) и генерацию случайного числа. Это всё было необходимо для отладки.

В общем, программа получилась такая:

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

import java.util.*;

class BigNumber {
        
        private final static int base = 1000000000;
        private static int karatsubaThreshold = 25;
        
        private int[] digits;
        
        private BigNumber() {
        }
        
        public BigNumber(int arg) {
                if (0 > arg) {
                        throw new IllegalArgumentException("Numbers must be positive");
                }
                digits = new int[2];
                digits[0] = arg % base;
                digits[1] = arg / base;
        }
        
        public BigNumber(String arg) {
                BigNumber temp, ten;
                ten = new BigNumber(10);
                temp = new BigNumber(0);
                for (int k = 0; arg.length() > k; ++k) {
                        char c = arg.charAt(k);
                        if (('0' <= c) && (c <= '9')) {
                                temp = temp.multiplySimple(ten).add(new BigNumber(c - '0'));
                        }
                }
                digits = temp.digits;
        }
        
        public BigNumber(BigNumber arg) {
                digits = new int[arg.digits.length];
                for (int k = 0; arg.digits.length > k; ++k) {
                        digits[k] = arg.digits[k];
                }
        }
        /*
        private BigNumber(BigNumber arg, int padSize) {
                if (0 > padSize) {
                        throw new IllegalArgumentException("Pad size must be positive");
                }
                digits = new int[arg.digits.length + padSize];
                for (int k = 0; arg.digits.length > k; ++k) {
                        digits[k] = arg.digits[k];
                }
        }
        //*/
        public BigNumber(Random rng, int size) {
                digits = new int[size];
                for (int k = 0; size > k; ++k) {
                        digits[k] = rng.nextInt(base);
                }
        }
        
        public BigNumber add(BigNumber arg) {
                int k, size1, size2, sum;
                BigNumber augend, addend, result;
                size1 = this.digits.length - 1;
                if (0 != this.digits[size1]) {
                        ++size1;
                }
                size2 = arg.digits.length - 1;
                if (0 != arg.digits[size2]) {
                        ++size2;
                }
                if (size1 < size2) {
                        augend = arg;
                        addend = this;
                        sum = size1;
                        size1 = size2;
                        size2 = sum;
                } else {
                        augend = this;
                        addend = arg;
                }
                result = new BigNumber();
                result.digits = new int[size1 + 1];
                for (k = 0; size1 > k; ++k) {
                        result.digits[k] = augend.digits[k];
                }
                sum = 0;
                for (k = 0; size2 > k; ++k) {
                        sum += result.digits[k];
                        sum += addend.digits[k];
                        if (base > sum) {
                                result.digits[k] = sum;
                                sum = 0;
                        } else {
                                result.digits[k] = sum - base;
                                sum = 1;
                        }
                }
                while (0 != sum) {
                        sum += result.digits[k];
                        if (base > sum) {
                                result.digits[k] = sum;
                                sum = 0;
                        } else {
                                result.digits[k] = sum - base;
                                sum = 1;
                        }
                        ++k;
                }
                return result;
        }
        
        public BigNumber multiplyKaratsuba(BigNumber arg) {
                int k, l, size1, size2, offset, limit, sum;
                BigNumber result, a, b, a1, a2, b1, b2, axb1, axb2, apa, bpb, temp;
                size1 = digits.length;
                if (0 == digits[size1 - 1]) {
                        --size1;
                }
                size2 = arg.digits.length;
                if (0 == arg.digits[size2 - 1]) {
                        --size2;
                }
                if (karatsubaThreshold > size1 || karatsubaThreshold > size2) {
                        return this.multiplySimple(arg);
                }
                if (size2 < size1) {
                        a = this;
                        b = arg;
                        offset = (size1 + 1) / 2;
                } else {
                        a = arg;
                        b = this;
                        offset = size2;
                        size2 = size1;
                        size1 = offset;
                        offset = (size1 + 1) / 2;
                }
                //==============================
                result = new BigNumber();
                result.digits = new int[size1 + size2];
                a1 = new BigNumber();
                a2 = new BigNumber();
                a1.digits = new int[offset];
                a2.digits = new int[size1 - offset];
                for (k = 0; offset > k; ++k) {
                        a1.digits[k] = a.digits[k];
                }
                for (l = 0; size1 > k; ++k, ++l) {
                        a2.digits[l] = a.digits[k];
                }
                //==============================
                if (size2 > offset) {
                        b1 = new BigNumber();
                        b2 = new BigNumber();
                        b1.digits = new int[offset];
                        b2.digits = new int[size2 - offset];
                        for (k = 0; offset > k; ++k) {
                                b1.digits[k] = b.digits[k];
                        }
                        for (l = 0; size2 > k; ++k, ++l) {
                                b2.digits[l] = b.digits[k];
                        }
                        //==============================
                        axb1 = a1.multiplyKaratsuba(b1);
                        axb2 = a2.multiplyKaratsuba(b2);
                        apa = a1.add(a2);
                        bpb = b1.add(b2);
                        temp = apa.multiplyKaratsuba(bpb);
                        //==============================
                        limit = axb1.digits.length - 1;
                        if (0 != axb1.digits[limit]) {
                                ++limit;
                        }
                        sum = 0;
                        for (k = 0; limit > k; ++k) {
                                sum += temp.digits[k];
                                sum -= axb1.digits[k];
                                if (0 > sum) {
                                        temp.digits[k] = sum + base;
                                        sum = -1;
                                } else {
                                        temp.digits[k] = sum;
                                        sum = 0;
                                }
                        }
                        while (0 != sum) {
                                sum += temp.digits[k];
                                if (0 > sum) {
                                        temp.digits[k] = sum + base;
                                        sum = -1;
                                } else {
                                        temp.digits[k] = sum;
                                        sum = 0;
                                }
                                ++k;
                        }
                        //==============================
                        limit = axb2.digits.length - 1;
                        if (0 != axb2.digits[limit]) {
                                ++limit;
                        }
                        sum = 0;
                        for (k = 0; limit > k; ++k) {
                                sum += temp.digits[k];
                                sum -= axb2.digits[k];
                                if (0 > sum) {
                                        temp.digits[k] = sum + base;
                                        sum = -1;
                                } else {
                                        temp.digits[k] = sum;
                                        sum = 0;
                                }
                        }
                        while (0 != sum) {
                                sum += temp.digits[k];
                                if (0 > sum) {
                                        temp.digits[k] = sum + base;
                                        sum = -1;
                                } else {
                                        temp.digits[k] = sum;
                                        sum = 0;
                                }
                                ++k;
                        }
                        //==============================
                        limit = axb1.digits.length;
                        for (k = 0; limit > k; ++k) {
                                result.digits[k] = axb1.digits[k];
                        }
                        //==============================
                        limit = temp.digits.length - 1;
                        while (0 == temp.digits[limit]) {
                                --limit;
                        }
                        ++limit;
                        sum = 0;
                        for (k = 0, l = offset; limit > k; ++k, ++l) {
                                sum += result.digits[l];
                                sum += temp.digits[k];
                                if (base > sum) {
                                        result.digits[l] = sum;
                                        sum = 0;
                                } else {
                                        result.digits[l] = sum - base;
                                        sum = 1;
                                }
                        }
                        if (0 != sum) {
                                result.digits[l] = sum;
                        }
                        //==============================
                        limit = axb2.digits.length - 1;
                        if (0 != axb2.digits[limit]) {
                                ++limit;
                        }
                        sum = 0;
                        for (k = 0, l = 2 * offset; limit > k; ++k, ++l) {
                                sum += result.digits[l];
                                sum += axb2.digits[k];
                                if (base > sum) {
                                        result.digits[l] = sum;
                                        sum = 0;
                                } else {
                                        result.digits[l] = sum - base;
                                        sum = 1;
                                }
                        }
                        if (0 != sum) {
                                result.digits[l] = sum;
                        }
                } else {
                        axb1 = a1.multiplyKaratsuba(b);
                        axb2 = a2.multiplyKaratsuba(b);
                        //==============================
                        limit = axb1.digits.length;
                        for (k = 0; limit > k; ++k) {
                                result.digits[k] = axb1.digits[k];
                        }
                        //==============================
                        limit = axb2.digits.length - 1;
                        if (0 != axb2.digits[limit]) {
                                ++limit;
                        }
                        sum = 0;
                        for (k = 0, l = offset; limit > k; ++k, ++l) {
                                sum += result.digits[l];
                                sum += axb2.digits[k];
                                if (base > sum) {
                                        result.digits[l] = sum;
                                        sum = 0;
                                } else {
                                        result.digits[l] = sum - base;
                                        sum = 1;
                                }
                        }
                        if (0 != sum) {
                                result.digits[l] = sum;
                        }
                }
                return result;
        }
        
        public BigNumber multiplySimple(BigNumber arg) {
                int k, l, size1, size2;
                long mult, sum;
                BigNumber result;
                size1 = digits.length;
                if (0 == digits[size1 - 1]) {
                        --size1;
                }
                size2 = arg.digits.length;
                if (0 == arg.digits[size2 - 1]) {
                        --size2;
                }
                result = new BigNumber();
                result.digits = new int[size1 + size2];
                for (k = 0; size1 > k; ++k) {
                        mult = digits[k];
                        sum = 0;
                        for (l = 0; size2 > l; ++l) {
                                sum += result.digits[k + l];
                                sum += mult * arg.digits[l];
                                result.digits[k + l] = (int) (sum % base);
                                sum /= base;
                        }
                        if (0 != sum) {
                                result.digits[k + l] = (int) sum;
                        }
                }
                return result;
        }
        
        public int compareTo(BigNumber arg) {
                int k, flag;
                BigNumber a, b;
                if (this.digits.length > arg.digits.length) {
                        a = this;
                        b = arg;
                        flag = 1;
                } else {
                        a = arg;
                        b = this;
                        flag = -1;
                }
                k = a.digits.length;
                while (b.digits.length > k) {
                        --k;
                        if (0 != a.digits[k]) {
                                return flag;
                        }
                }
                do {
                        --k;
                        if (b.digits[k] < a.digits[k]) {
                                return flag;
                        } else if (b.digits[k] > a.digits[k]) {
                                return -flag;
                        }
                } while (0 != k);
                return 0;
        }
        
        public String toString() {
                int k, l, value;
                l = 12 * digits.length;
                char[] buffer = new char[l];
                for (k = 0; digits.length > k; ++k) {
                        value = digits[k];
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = ' ';
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = ' ';
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value % 10);
                        value /= 10;
                        buffer[--l] = (char) ('0' + value);
                        buffer[--l] = ' ';
                }
                return new String(buffer);
        }
}

public class Study_Karatsuba_method {
        
        private final static Random rng = new Random();
        
        public static void main(String[] args) {
                int k, numberSize, tableSize;
                long time1, time2, time3;
                double temp, mean1, std1, mean2, std2;
                double[] table1, table2;
                BigNumber a, b, c, d;
                Locale.setDefault(new Locale("en", "US"));
                tableSize = 12;
                table1 = new double[tableSize];
                table2 = new double[tableSize];
                numberSize = 29000;
                while (30 < numberSize) {
                        mean1 = 0.0;
                        mean2 = 0.0;
                        for (k = 0; tableSize > k; ++k) {
                                a = new BigNumber(rng, numberSize);
                                b = new BigNumber(rng, numberSize);
                                time1 = System.nanoTime();
                                c = a.multiplySimple(b);
                                time2 = System.nanoTime();
                                d = a.multiplyKaratsuba(b);
                                time3 = System.nanoTime();
                                if (0 != d.compareTo(c)) {
                                        System.out.println("ERROR");
                                }
                                temp = 1e-9 * (time2 - time1);
                                mean1 += temp;
                                table1[k] = temp;
                                temp = 1e-9 * (time3 - time2);
                                mean2 += temp;
                                table2[k] = temp;
                        }
                        mean1 /= tableSize;
                        mean2 /= tableSize;
                        std1 = 0.0;
                        std2 = 0.0;
                        for (k = 0; tableSize > k; ++k) {
                                temp = table1[k] - mean1;
                                std1 += temp * temp;
                                temp = table2[k] - mean2;
                                std2 += temp * temp;
                        }
                        std1 = Math.sqrt(std1 / tableSize / (tableSize - 1));
                        std2 = Math.sqrt(std2 / tableSize / (tableSize - 1));
                        System.out.println(String.format("%6d %12.9f %12.9f %12.9f %12.9f", numberSize, mean1, std1, mean2, std2));
                        numberSize -= numberSize / 8;
                }
        }
}
 

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

Результаты экспериментальной проверки времени работы на моём компе (с частотой процессора 2 ГГц) оказались такие:

[slavav]

B@R5uk в сообщении #1439184 писал(а):

Интересный вопрос, кстати, что быстрее: считать в двоичной, а потом как-то хитро конвертировать или же так, как сделал я?

Вы сделали правильно.
Смотрите секцию BASE CONVERSION в THE COMPLEXITY OF MULTIPLE-PRECISION
ARITHMETIC (https://pdfs.semanticscholar.org/e70f/be46bf4deb2918af73717b31e62ad1b49717.pdf?_ga=2.244133389.96897432.1581340337-1381867169.1581340337).

[Dmitriy40]

B@R5uk в сообщении #1439184 писал(а):

Интересный вопрос, кстати, что быстрее: считать в двоичной, а потом как-то хитро конвертировать, или же так, как сделал я?

По моему ответ на него очевиден: если надо сделать лишь одно возведение в квадрат (и плюс-минусы), то конечно лучше степень 10, если же надо делать ещё какие-то затратные действия, то однократный перевод из десятичной в двоичную и один раз обратно для вывода будет быстрее. Можете кстати сами легко проверить, установите base в степень двойки (например 29-ю) и сравните скорость, наверняка компилятор заменит деление на AND и сдвиги.

[slavav]

Преобразование системы счисления $10 \to 2$ - сложная операция если говорим о длинных числах и заботимся об асимптотике.
Обратное преобразование ещё сложнее так как использует длинную арифметику по основанию 10. Что и сделал ТС.

[Dmitriy40]

Преобразование $10\to 2$ сводится к умножениям (на 10, которое тоже легко сводится к двум сдвигам и сложению), вот обратное требует делений, что и плохо. Но и то и другое выполняется лишь для входных и выходных чисел, причём строго по одному разу для каждого, а не для всех промежуточных если таковые есть. Потому если промежуточных вычислений больше пары-тройки, выгоднее их ускорить за счёт замедления преобразований входных/выходных чисел, которых обычно меньше на порядки. Впрочем при большом желании всё это можно и точно подсчитать.

[B@R5uk]

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

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

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

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

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

На 0-ом уровне: возводится в квадрат 1 число
На 1-ом уровне: возводится в квадрат 2 числа, находится 1 произведение
На 2-ом уровне: возводится в квадрат 4 числа, находится 5 произведений
На 3-ем уровне: возводится в квадрат 8 чисел, находится 19 произведений
На 4-ом уровне: возводится в квадрат 16 чисел, находится 65 произведений
На 5-ом уровне: возводится в квадрат 32 числа, находится 211 произведение
...
На n-ом уровне: возводится в квадрат $q_n$ чисел, находится $p_n$ произведений
На (n+1)-ом уровне: возводится в квадрат $q_{n+1}=2q_n$ чисел, находится $p_{n+1}=3p_n+q_n$ произведений

Поскольку $p_n=3^n-2^n$ растёт быстрее, чем $q_n=2^n$, то число возведений в квадрат будет пренебрежимо мало по сравнению с числом умножений при достаточно большой глубине ветвления. Так что, если бы и был быстрый метод возведения в квадрат малых чисел, то он бы ничего не дал для ускорения счёта больших чисел.

[slavav]

При возведении в квадрат методом Карацубы произведения тоже выражаются через квадраты. Благодаря этому константа может быть лучше.
$4ab = (a + b)^2 - (a - b)^2$.

[B@R5uk]

slavav, не взлетит. $a^2$ и $b^2$ всё равно считать придётся.