Everywhere: Hoist the Libraries folder to the top-level

Author: https://github.com/trflynn89
Commit: 93712b24bf
Pull-request: https://github.com/LadybirdBrowser/ladybird/pull/2256
Reviewed-by: https://github.com/sideshowbarker

Libraries/LibCrypto/BigInt/Algorithms/BitwiseOperations.cpp

@@ -0,0 +1,268 @@
+/*
+ * Copyright (c) 2020, Itamar S. <itamar8910@gmail.com>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+#include <AK/BigIntBase.h>
+#include <AK/BuiltinWrappers.h>
+#include <AK/NumericLimits.h>
+
+namespace Crypto {
+
+/**
+ * Complexity: O(N) where N is the number of words in the shorter value
+ * Method:
+ * Apply <op> word-wise until words in the shorter value are used up
+ * then copy the rest of the words verbatim from the longer value.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::bitwise_or_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& output)
+{
+    // If either of the BigInts are invalid, the output is just the other one.
+    if (left.is_invalid()) {
+        output.set_to(right);
+        return;
+    }
+    if (right.is_invalid()) {
+        output.set_to(left);
+        return;
+    }
+
+    UnsignedBigInteger const *shorter, *longer;
+    if (left.length() < right.length()) {
+        shorter = &left;
+        longer = &right;
+    } else {
+        shorter = &right;
+        longer = &left;
+    }
+
+    output.m_words.resize_and_keep_capacity(longer->length());
+
+    size_t longer_offset = longer->length() - shorter->length();
+    for (size_t i = 0; i < shorter->length(); ++i)
+        output.m_words[i] = longer->words()[i] | shorter->words()[i];
+
+    __builtin_memcpy(output.m_words.data() + shorter->length(), longer->words().data() + shorter->length(), sizeof(u32) * longer_offset);
+}
+
+/**
+ * Complexity: O(N) where N is the number of words in the shorter value
+ * Method:
+ * Apply 'and' word-wise until words in the shorter value are used up
+ * and zero the rest.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::bitwise_and_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& output)
+{
+    // If either of the BigInts are invalid, the output is just the other one.
+    if (left.is_invalid()) {
+        output.set_to(right);
+        return;
+    }
+    if (right.is_invalid()) {
+        output.set_to(left);
+        return;
+    }
+
+    UnsignedBigInteger const *shorter, *longer;
+    if (left.length() < right.length()) {
+        shorter = &left;
+        longer = &right;
+    } else {
+        shorter = &right;
+        longer = &left;
+    }
+
+    output.m_words.resize_and_keep_capacity(longer->length());
+
+    size_t longer_offset = longer->length() - shorter->length();
+    for (size_t i = 0; i < shorter->length(); ++i)
+        output.m_words[i] = longer->words()[i] & shorter->words()[i];
+
+    __builtin_memset(output.m_words.data() + shorter->length(), 0, sizeof(u32) * longer_offset);
+}
+
+/**
+ * Complexity: O(N) where N is the number of words in the shorter value
+ * Method:
+ * Apply 'xor' word-wise until words in the shorter value are used up
+ * and copy the rest.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::bitwise_xor_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& output)
+{
+    // If either of the BigInts are invalid, the output is just the other one.
+    if (left.is_invalid()) {
+        output.set_to(right);
+        return;
+    }
+    if (right.is_invalid()) {
+        output.set_to(left);
+        return;
+    }
+
+    UnsignedBigInteger const *shorter, *longer;
+    if (left.length() < right.length()) {
+        shorter = &left;
+        longer = &right;
+    } else {
+        shorter = &right;
+        longer = &left;
+    }
+
+    output.m_words.resize_and_keep_capacity(longer->length());
+
+    size_t longer_offset = longer->length() - shorter->length();
+    for (size_t i = 0; i < shorter->length(); ++i)
+        output.m_words[i] = longer->words()[i] ^ shorter->words()[i];
+
+    __builtin_memcpy(output.m_words.data() + shorter->length(), longer->words().data() + shorter->length(), sizeof(u32) * longer_offset);
+}
+
+/**
+ * Complexity: O(N) where N is the number of words
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::bitwise_not_fill_to_one_based_index_without_allocation(
+    UnsignedBigInteger const& right,
+    size_t index,
+    UnsignedBigInteger& output)
+{
+    // If the value is invalid, the output value is invalid as well.
+    if (right.is_invalid()) {
+        output.invalidate();
+        return;
+    }
+
+    if (index == 0) {
+        output.set_to_0();
+        return;
+    }
+    size_t size = (index + UnsignedBigInteger::BITS_IN_WORD - 1) / UnsignedBigInteger::BITS_IN_WORD;
+
+    output.m_words.resize_and_keep_capacity(size);
+    VERIFY(size > 0);
+    for (size_t i = 0; i < size - 1; ++i)
+        output.m_words[i] = ~(i < right.length() ? right.words()[i] : 0);
+
+    index -= (size - 1) * UnsignedBigInteger::BITS_IN_WORD;
+    auto last_word_index = size - 1;
+    auto last_word = last_word_index < right.length() ? right.words()[last_word_index] : 0;
+
+    output.m_words[last_word_index] = (NumericLimits<UnsignedBigInteger::Word>::max() >> (UnsignedBigInteger::BITS_IN_WORD - index)) & ~last_word;
+}
+
+/**
+ * Complexity : O(N + num_bits % 8) where N is the number of words in the number
+ * Shift method :
+ * Start by shifting by whole words in num_bits (by putting missing words at the start),
+ * then shift the number's words two by two by the remaining amount of bits.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::shift_left_without_allocation(
+    UnsignedBigInteger const& number,
+    size_t num_bits,
+    UnsignedBigInteger& temp_result,
+    UnsignedBigInteger& temp_plus,
+    UnsignedBigInteger& output)
+{
+    // We can only do shift operations on individual words
+    // where the shift amount is <= size of word (32).
+    // But we do know how to shift by a multiple of word size (e.g 64=32*2)
+    // So we first shift the result by how many whole words fit in 'num_bits'
+    shift_left_by_n_words(number, num_bits / UnsignedBigInteger::BITS_IN_WORD, temp_result);
+
+    output.set_to(temp_result);
+
+    // And now we shift by the leftover amount of bits
+    num_bits %= UnsignedBigInteger::BITS_IN_WORD;
+
+    if (num_bits == 0) {
+        return;
+    }
+
+    for (size_t i = 0; i < temp_result.length(); ++i) {
+        u32 current_word_of_temp_result = shift_left_get_one_word(temp_result, num_bits, i);
+        output.m_words[i] = current_word_of_temp_result;
+    }
+
+    // Shifting the last word can produce a carry
+    u32 carry_word = shift_left_get_one_word(temp_result, num_bits, temp_result.length());
+    if (carry_word != 0) {
+
+        // output += (carry_word << temp_result.length())
+        // FIXME : Using temp_plus this way to transform carry_word into a bigint is not
+        //         efficient nor pretty. Maybe we should have an "add_with_shift" method ?
+        temp_plus.set_to_0();
+        temp_plus.m_words.append(carry_word);
+        shift_left_by_n_words(temp_plus, temp_result.length(), temp_result);
+        add_into_accumulator_without_allocation(output, temp_result);
+    }
+}
+
+FLATTEN void UnsignedBigIntegerAlgorithms::shift_right_without_allocation(
+    UnsignedBigInteger const& number,
+    size_t num_bits,
+    UnsignedBigInteger& output)
+{
+    output.m_words.resize_and_keep_capacity(number.length() - (num_bits / UnsignedBigInteger::BITS_IN_WORD));
+    Ops::shift_right(number.words_span(), num_bits, output.words_span());
+}
+
+void UnsignedBigIntegerAlgorithms::shift_left_by_n_words(
+    UnsignedBigInteger const& number,
+    size_t number_of_words,
+    UnsignedBigInteger& output)
+{
+    // shifting left by N words means just inserting N zeroes to the beginning of the words vector
+    output.set_to_0();
+    output.m_words.resize_and_keep_capacity(number_of_words + number.length());
+
+    __builtin_memset(output.m_words.data(), 0, number_of_words * sizeof(unsigned));
+    __builtin_memcpy(&output.m_words.data()[number_of_words], number.m_words.data(), number.m_words.size() * sizeof(unsigned));
+}
+
+void UnsignedBigIntegerAlgorithms::shift_right_by_n_words(
+    UnsignedBigInteger const& number,
+    size_t number_of_words,
+    UnsignedBigInteger& output)
+{
+    // shifting right by N words means just not copying the first words
+    output.set_to_0();
+    output.m_words.resize_and_keep_capacity(number.length() - number_of_words);
+    __builtin_memcpy(output.m_words.data(), &number.m_words.data()[number_of_words], (number.m_words.size() - number_of_words) * sizeof(unsigned));
+}
+
+/**
+ * Returns the word at a requested index in the result of a shift operation
+ */
+ALWAYS_INLINE UnsignedBigInteger::Word UnsignedBigIntegerAlgorithms::shift_left_get_one_word(
+    UnsignedBigInteger const& number,
+    size_t num_bits,
+    size_t result_word_index)
+{
+    // "<= length()" (rather than length() - 1) is intentional,
+    // The result index of length() is used when calculating the carry word
+    VERIFY(result_word_index <= number.length());
+    VERIFY(num_bits <= UnsignedBigInteger::BITS_IN_WORD);
+    u32 result = 0;
+
+    // we need to check for "num_bits != 0" since shifting right by 32 is apparently undefined behavior!
+    if (result_word_index > 0 && num_bits != 0) {
+        result += number.m_words[result_word_index - 1] >> (UnsignedBigInteger::BITS_IN_WORD - num_bits);
+    }
+    if (result_word_index < number.length() && num_bits < 32) {
+        result += number.m_words[result_word_index] << num_bits;
+    }
+    return result;
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/Division.cpp

@@ -0,0 +1,123 @@
+/*
+ * Copyright (c) 2020, Itamar S. <itamar8910@gmail.com>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+#include <AK/BigIntBase.h>
+#include <AK/BuiltinWrappers.h>
+
+namespace Crypto {
+
+using AK::Detail::div_mod_words;
+using AK::Detail::dword;
+
+/**
+ * Complexity: O(N^2) where N is the number of words in the larger number
+ * Division method:
+ * Knuth's Algorithm D, see UFixedBigIntDivision.h for more details
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::divide_without_allocation(
+    UnsignedBigInteger const& numerator,
+    UnsignedBigInteger const& denominator,
+    UnsignedBigInteger& quotient,
+    UnsignedBigInteger& remainder)
+{
+    size_t dividend_len = numerator.trimmed_length();
+    size_t divisor_len = denominator.trimmed_length();
+
+    VERIFY(divisor_len != 0);
+
+    // Fast paths
+    // Division by 1
+    if (divisor_len == 1 && denominator.m_words[0] == 1) {
+        quotient.set_to(numerator);
+        remainder.set_to_0();
+        return;
+    }
+
+    if (dividend_len < divisor_len) {
+        quotient.set_to_0();
+        remainder.set_to(numerator);
+        return;
+    }
+
+    if (divisor_len == 1 && dividend_len == 1) {
+        quotient.set_to(numerator.m_words[0] / denominator.m_words[0]);
+        remainder.set_to(numerator.m_words[0] % denominator.m_words[0]);
+        return;
+    }
+    // Division by Word
+    if (divisor_len == 1) {
+        quotient.resize_with_leading_zeros(dividend_len);
+        remainder.resize_with_leading_zeros(1);
+
+        // FIXME: Use a "DoubleWord" to allow increasing the Word size of
+        //        BigInt in the future
+        static_assert(UnsignedBigInteger::BITS_IN_WORD == 32);
+        auto u = dword(numerator.m_words[dividend_len - 2], numerator.m_words[dividend_len - 1]);
+        auto divisor = denominator.m_words[0];
+
+        auto top = u / divisor;
+        quotient.m_words[dividend_len - 1] = top >> UnsignedBigInteger::BITS_IN_WORD;
+        quotient.m_words[dividend_len - 2] = static_cast<UnsignedBigInteger::Word>(top);
+
+        auto carry = static_cast<UnsignedBigInteger::Word>(u % divisor);
+        for (size_t i = dividend_len - 2; i-- != 0;)
+            quotient.m_words[i] = div_mod_words(numerator.m_words[i], carry, divisor, carry);
+        remainder.m_words[0] = carry;
+        return;
+    }
+
+    // Knuth's algorithm D
+    auto dividend = numerator;
+    dividend.resize_with_leading_zeros(dividend_len + 1);
+    auto divisor = denominator;
+
+    quotient.resize_with_leading_zeros(dividend_len - divisor_len + 1);
+    remainder.resize_with_leading_zeros(divisor_len);
+
+    Ops::div_mod_internal<true>(
+        dividend.words_span(), divisor.words_span(),
+        quotient.words_span(), remainder.words_span(),
+        dividend_len, divisor_len);
+}
+
+/**
+ * Complexity : O(N) where N is the number of digits in the numerator
+ * Division method :
+ * Starting from the most significant one, for each half-word of the numerator, combine it
+ * with the existing remainder if any, divide the combined number as a u32 operation and
+ * update the quotient / remainder as needed.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::divide_u16_without_allocation(
+    UnsignedBigInteger const& numerator,
+    UnsignedBigInteger::Word denominator,
+    UnsignedBigInteger& quotient,
+    UnsignedBigInteger& remainder)
+{
+    VERIFY(denominator < (1 << 16));
+    UnsignedBigInteger::Word remainder_word = 0;
+    auto numerator_length = numerator.trimmed_length();
+    quotient.set_to_0();
+    quotient.m_words.resize(numerator_length);
+    for (int word_index = numerator_length - 1; word_index >= 0; --word_index) {
+        auto word_high = numerator.m_words[word_index] >> 16;
+        auto word_low = numerator.m_words[word_index] & ((1 << 16) - 1);
+
+        auto number_to_divide_high = (remainder_word << 16) | word_high;
+        auto quotient_high = number_to_divide_high / denominator;
+        remainder_word = number_to_divide_high % denominator;
+
+        auto number_to_divide_low = remainder_word << 16 | word_low;
+        auto quotient_low = number_to_divide_low / denominator;
+        remainder_word = number_to_divide_low % denominator;
+
+        quotient.m_words[word_index] = (quotient_high << 16) | quotient_low;
+    }
+    remainder.set_to(remainder_word);
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/GCD.cpp

@@ -0,0 +1,39 @@
+/*
+ * Copyright (c) 2020, Ali Mohammad Pur <mpfard@serenityos.org>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+
+namespace Crypto {
+
+void UnsignedBigIntegerAlgorithms::destructive_GCD_without_allocation(
+    UnsignedBigInteger& temp_a,
+    UnsignedBigInteger& temp_b,
+    UnsignedBigInteger& temp_quotient,
+    UnsignedBigInteger& temp_remainder,
+    UnsignedBigInteger& output)
+{
+    for (;;) {
+        if (temp_a == 0) {
+            output.set_to(temp_b);
+            return;
+        }
+
+        // temp_b %= temp_a
+        divide_without_allocation(temp_b, temp_a, temp_quotient, temp_remainder);
+        temp_b.set_to(temp_remainder);
+        if (temp_b == 0) {
+            output.set_to(temp_a);
+            return;
+        }
+
+        // temp_a %= temp_b
+        divide_without_allocation(temp_a, temp_b, temp_quotient, temp_remainder);
+        temp_a.set_to(temp_remainder);
+    }
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/ModularInverse.cpp

@@ -0,0 +1,90 @@
+/*
+ * Copyright (c) 2020, Ali Mohammad Pur <mpfard@serenityos.org>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+
+namespace Crypto {
+
+void UnsignedBigIntegerAlgorithms::modular_inverse_without_allocation(
+    UnsignedBigInteger const& a,
+    UnsignedBigInteger const& b,
+    UnsignedBigInteger& temp_1,
+    UnsignedBigInteger& temp_minus,
+    UnsignedBigInteger& temp_quotient,
+    UnsignedBigInteger& temp_d,
+    UnsignedBigInteger& temp_u,
+    UnsignedBigInteger& temp_v,
+    UnsignedBigInteger& temp_x,
+    UnsignedBigInteger& result)
+{
+    UnsignedBigInteger one { 1 };
+
+    temp_u.set_to(a);
+    if (!a.is_odd()) {
+        // u += b
+        add_into_accumulator_without_allocation(temp_u, b);
+    }
+
+    temp_v.set_to(b);
+    temp_x.set_to(0);
+
+    // d = b - 1
+    subtract_without_allocation(b, one, temp_d);
+
+    while (!(temp_v == 1)) {
+        while (temp_v < temp_u) {
+            // u -= v
+            subtract_without_allocation(temp_u, temp_v, temp_minus);
+            temp_u.set_to(temp_minus);
+
+            // d += x
+            add_into_accumulator_without_allocation(temp_d, temp_x);
+
+            while (!temp_u.is_odd()) {
+                if (temp_d.is_odd()) {
+                    // d += b
+                    add_into_accumulator_without_allocation(temp_d, b);
+                }
+
+                // u /= 2
+                divide_u16_without_allocation(temp_u, 2, temp_quotient, temp_1);
+                temp_u.set_to(temp_quotient);
+
+                // d /= 2
+                divide_u16_without_allocation(temp_d, 2, temp_quotient, temp_1);
+                temp_d.set_to(temp_quotient);
+            }
+        }
+
+        // v -= u
+        subtract_without_allocation(temp_v, temp_u, temp_minus);
+        temp_v.set_to(temp_minus);
+
+        // x += d
+        add_into_accumulator_without_allocation(temp_x, temp_d);
+
+        while (!temp_v.is_odd()) {
+            if (temp_x.is_odd()) {
+                // x += b
+                add_into_accumulator_without_allocation(temp_x, b);
+            }
+
+            // v /= 2
+            divide_u16_without_allocation(temp_v, 2, temp_quotient, temp_1);
+            temp_v.set_to(temp_quotient);
+
+            // x /= 2
+            divide_u16_without_allocation(temp_x, 2, temp_quotient, temp_1);
+            temp_x.set_to(temp_quotient);
+        }
+    }
+
+    // return x % b
+    divide_without_allocation(temp_x, b, temp_quotient, result);
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/ModularPower.cpp

@@ -0,0 +1,281 @@
+/*
+ * Copyright (c) 2020, Ali Mohammad Pur <mpfard@serenityos.org>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+
+namespace Crypto {
+
+void UnsignedBigIntegerAlgorithms::destructive_modular_power_without_allocation(
+    UnsignedBigInteger& ep,
+    UnsignedBigInteger& base,
+    UnsignedBigInteger const& m,
+    UnsignedBigInteger& temp_1,
+    UnsignedBigInteger& temp_2,
+    UnsignedBigInteger& temp_3,
+    UnsignedBigInteger& temp_multiply,
+    UnsignedBigInteger& temp_quotient,
+    UnsignedBigInteger& temp_remainder,
+    UnsignedBigInteger& exp)
+{
+    exp.set_to(1);
+    while (!(ep < 1)) {
+        if (ep.words()[0] % 2 == 1) {
+            // exp = (exp * base) % m;
+            multiply_without_allocation(exp, base, temp_1, temp_2, temp_3, temp_multiply);
+            divide_without_allocation(temp_multiply, m, temp_quotient, temp_remainder);
+            exp.set_to(temp_remainder);
+        }
+
+        // ep = ep / 2;
+        ep.set_to(ep.shift_right(1));
+
+        // base = (base * base) % m;
+        multiply_without_allocation(base, base, temp_1, temp_2, temp_3, temp_multiply);
+        divide_without_allocation(temp_multiply, m, temp_quotient, temp_remainder);
+        base.set_to(temp_remainder);
+
+        // Note that not clamping here would cause future calculations (multiply, specifically) to allocate even more unused space
+        // which would then persist through the temp bigints, and significantly slow down later loops.
+        // To avoid that, we can clamp to a specific max size, or just clamp to the min needed amount of space.
+        ep.clamp_to_trimmed_length();
+        exp.clamp_to_trimmed_length();
+        base.clamp_to_trimmed_length();
+    }
+}
+
+/**
+ * Compute (1/value) % 2^32.
+ * This needs an odd input value
+ * Algorithm from: Dumas, J.G. "On Newton–Raphson Iteration for Multiplicative Inverses Modulo Prime Powers".
+ */
+ALWAYS_INLINE static u32 inverse_wrapped(u32 value)
+{
+    VERIFY(value & 1);
+
+    u64 b = static_cast<u64>(value);
+    u64 k0 = (2 - b);
+    u64 t = (b - 1);
+    size_t i = 1;
+    while (i < 32) {
+        t = t * t;
+        k0 = k0 * (t + 1);
+        i <<= 1;
+    }
+    return static_cast<u32>(-k0);
+}
+
+/**
+ * Computes z = x * y + c. z_carry contains the top bits, z contains the bottom bits.
+ */
+ALWAYS_INLINE static void linear_multiplication_with_carry(u32 x, u32 y, u32 c, u32& z_carry, u32& z)
+{
+    u64 result = static_cast<u64>(x) * static_cast<u64>(y) + static_cast<u64>(c);
+    z_carry = static_cast<u32>(result >> 32);
+    z = static_cast<u32>(result);
+}
+
+/**
+ * Computes z = a + b. z_carry contains the top bit (1 or 0), z contains the bottom bits.
+ */
+ALWAYS_INLINE static void addition_with_carry(u32 a, u32 b, u32& z_carry, u32& z)
+{
+    u64 result = static_cast<u64>(a) + static_cast<u64>(b);
+    z_carry = static_cast<u32>(result >> 32);
+    z = static_cast<u32>(result);
+}
+
+/**
+ * Computes a montgomery "fragment" for y_i. This computes "z[i] += x[i] * y_i" for all words while rippling the carry, and returns the carry.
+ * Algorithm from: Gueron, "Efficient Software Implementations of Modular Exponentiation". (https://eprint.iacr.org/2011/239.pdf)
+ */
+UnsignedBigInteger::Word UnsignedBigIntegerAlgorithms::montgomery_fragment(UnsignedBigInteger& z, size_t offset_in_z, UnsignedBigInteger const& x, UnsignedBigInteger::Word y_digit, size_t num_words)
+{
+    UnsignedBigInteger::Word carry { 0 };
+    for (size_t i = 0; i < num_words; ++i) {
+        UnsignedBigInteger::Word a_carry;
+        UnsignedBigInteger::Word a;
+        linear_multiplication_with_carry(x.m_words[i], y_digit, z.m_words[offset_in_z + i], a_carry, a);
+        UnsignedBigInteger::Word b_carry;
+        UnsignedBigInteger::Word b;
+        addition_with_carry(a, carry, b_carry, b);
+        z.m_words[offset_in_z + i] = b;
+        carry = a_carry + b_carry;
+    }
+    return carry;
+}
+
+/**
+ * Computes the "almost montgomery" product : x * y * 2 ^ (-num_words * BITS_IN_WORD) % modulo
+ * [Note : that means that the result z satisfies z * 2^(num_words * BITS_IN_WORD) % modulo = x * y % modulo]
+ * assuming :
+ *  - x, y and modulo are all already padded to num_words
+ *  - k = inverse_wrapped(modulo) (optimization to not recompute K each time)
+ * Algorithm from: Gueron, "Efficient Software Implementations of Modular Exponentiation". (https://eprint.iacr.org/2011/239.pdf)
+ */
+void UnsignedBigIntegerAlgorithms::almost_montgomery_multiplication_without_allocation(
+    UnsignedBigInteger const& x,
+    UnsignedBigInteger const& y,
+    UnsignedBigInteger const& modulo,
+    UnsignedBigInteger& z,
+    UnsignedBigInteger::Word k,
+    size_t num_words,
+    UnsignedBigInteger& result)
+{
+    VERIFY(x.length() >= num_words);
+    VERIFY(y.length() >= num_words);
+    VERIFY(modulo.length() >= num_words);
+
+    z.set_to(0);
+    z.resize_with_leading_zeros(num_words * 2);
+
+    UnsignedBigInteger::Word previous_double_carry { 0 };
+    for (size_t i = 0; i < num_words; ++i) {
+        // z[i->num_words+i] += x * y_i
+        UnsignedBigInteger::Word carry_1 = montgomery_fragment(z, i, x, y.m_words[i], num_words);
+        // z[i->num_words+i] += modulo * (z_i * k)
+        UnsignedBigInteger::Word t = z.m_words[i] * k;
+        UnsignedBigInteger::Word carry_2 = montgomery_fragment(z, i, modulo, t, num_words);
+
+        // Compute the carry by combining all of the carries of the previous computations
+        // Put it "right after" the range that we computed above
+        UnsignedBigInteger::Word temp_carry = previous_double_carry + carry_1;
+        UnsignedBigInteger::Word overall_carry = temp_carry + carry_2;
+        z.m_words[num_words + i] = overall_carry;
+
+        // Detect if there was a "double carry" for this word by checking if our carry results are smaller than their components
+        previous_double_carry = (temp_carry < carry_1 || overall_carry < carry_2) ? 1 : 0;
+    }
+
+    if (previous_double_carry == 0) {
+        // Return the top num_words bytes of Z, which contains our result.
+        shift_right_by_n_words(z, num_words, result);
+        result.resize_with_leading_zeros(num_words);
+        return;
+    }
+
+    // We have a carry, so we're "one bigger" than we need to be.
+    // Subtract the modulo from the result (the top half of z), and write it to the bottom half of Z since we have space.
+    // (With carry, of course.)
+    UnsignedBigInteger::Word c { 0 };
+    for (size_t i = 0; i < num_words; ++i) {
+        UnsignedBigInteger::Word z_digit = z.m_words[num_words + i];
+        UnsignedBigInteger::Word modulo_digit = modulo.m_words[i];
+        UnsignedBigInteger::Word new_z_digit = z_digit - modulo_digit - c;
+        z.m_words[i] = new_z_digit;
+        // Detect if the subtraction underflowed - from "Hacker's Delight"
+        c = ((modulo_digit & ~z_digit) | ((modulo_digit | ~z_digit) & new_z_digit)) >> (UnsignedBigInteger::BITS_IN_WORD - 1);
+    }
+
+    // Return the bottom num_words bytes of Z (with the carry bit handled)
+    z.m_words.resize(num_words);
+    result.set_to(z);
+    result.resize_with_leading_zeros(num_words);
+}
+
+/**
+ * Complexity: still O(N^3) with N the number of words in the largest word, but less complex than the classical mod power.
+ * Note: the montgomery multiplications requires an inverse modulo over 2^32, which is only defined for odd numbers.
+ */
+void UnsignedBigIntegerAlgorithms::montgomery_modular_power_with_minimal_allocations(
+    UnsignedBigInteger const& base,
+    UnsignedBigInteger const& exponent,
+    UnsignedBigInteger const& modulo,
+    UnsignedBigInteger& temp_z,
+    UnsignedBigInteger& rr,
+    UnsignedBigInteger& one,
+    UnsignedBigInteger& z,
+    UnsignedBigInteger& zz,
+    UnsignedBigInteger& x,
+    UnsignedBigInteger& temp_extra,
+    UnsignedBigInteger& result)
+{
+    VERIFY(modulo.is_odd());
+
+    // Note: While this is a constexpr variable for clarity and could be changed in theory,
+    // various optimized parts of the algorithm rely on this value being exactly 4.
+    constexpr size_t window_size = 4;
+
+    size_t num_words = modulo.trimmed_length();
+    UnsignedBigInteger::Word k = inverse_wrapped(modulo.m_words[0]);
+
+    one.set_to(1);
+
+    // rr = ( 2 ^ (2 * modulo.length() * BITS_IN_WORD) ) % modulo
+    shift_left_by_n_words(one, 2 * num_words, x);
+    divide_without_allocation(x, modulo, temp_extra, rr);
+    rr.resize_with_leading_zeros(num_words);
+
+    // x = base [% modulo, if x doesn't already fit in modulo's words]
+    x.set_to(base);
+    if (x.trimmed_length() > num_words)
+        divide_without_allocation(base, modulo, temp_extra, x);
+    x.resize_with_leading_zeros(num_words);
+
+    one.set_to(1);
+    one.resize_with_leading_zeros(num_words);
+
+    // Compute the montgomery powers from 0 to 2^window_size. powers[i] = x^i
+    UnsignedBigInteger powers[1 << window_size];
+    almost_montgomery_multiplication_without_allocation(one, rr, modulo, temp_z, k, num_words, powers[0]);
+    almost_montgomery_multiplication_without_allocation(x, rr, modulo, temp_z, k, num_words, powers[1]);
+    for (size_t i = 2; i < (1 << window_size); ++i)
+        almost_montgomery_multiplication_without_allocation(powers[i - 1], powers[1], modulo, temp_z, k, num_words, powers[i]);
+
+    z.set_to(powers[0]);
+    z.resize_with_leading_zeros(num_words);
+    zz.set_to(0);
+    zz.resize_with_leading_zeros(num_words);
+
+    ssize_t exponent_length = exponent.trimmed_length();
+    for (ssize_t word_in_exponent = exponent_length - 1; word_in_exponent >= 0; --word_in_exponent) {
+        UnsignedBigInteger::Word exponent_word = exponent.m_words[word_in_exponent];
+        size_t bit_in_word = 0;
+        while (bit_in_word < UnsignedBigInteger::BITS_IN_WORD) {
+            if (word_in_exponent != exponent_length - 1 || bit_in_word != 0) {
+                almost_montgomery_multiplication_without_allocation(z, z, modulo, temp_z, k, num_words, zz);
+                almost_montgomery_multiplication_without_allocation(zz, zz, modulo, temp_z, k, num_words, z);
+                almost_montgomery_multiplication_without_allocation(z, z, modulo, temp_z, k, num_words, zz);
+                almost_montgomery_multiplication_without_allocation(zz, zz, modulo, temp_z, k, num_words, z);
+            }
+            auto power_index = exponent_word >> (UnsignedBigInteger::BITS_IN_WORD - window_size);
+            auto& power = powers[power_index];
+            almost_montgomery_multiplication_without_allocation(z, power, modulo, temp_z, k, num_words, zz);
+
+            swap(z, zz);
+
+            // Move to the next window
+            exponent_word <<= window_size;
+            bit_in_word += window_size;
+        }
+    }
+
+    almost_montgomery_multiplication_without_allocation(z, one, modulo, temp_z, k, num_words, zz);
+
+    if (zz < modulo) {
+        result.set_to(zz);
+        result.clamp_to_trimmed_length();
+        return;
+    }
+
+    // Note : Since we were using "almost montgomery" multiplications, we aren't guaranteed to be under the modulo already.
+    // So, if we're here, we need to respect the modulo.
+    // We can, however, start by trying to subtract the modulo, just in case we're close.
+    subtract_without_allocation(zz, modulo, result);
+
+    if (modulo < zz) {
+        // Note: This branch shouldn't happen in theory (as noted in https://github.com/rust-num/num-bigint/blob/master/src/biguint/monty.rs#L210)
+        // Let's dbgln the values we used. That way, if we hit this branch, we can contribute these values for test cases.
+        dbgln("Encountered the modulo branch during a montgomery modular power. Params : {} - {} - {}", base, exponent, modulo);
+        // We just clobber all the other temporaries that we don't need for the division.
+        // This is wasteful, but we're on the edgiest of cases already.
+        divide_without_allocation(zz, modulo, temp_extra, result);
+    }
+
+    result.clamp_to_trimmed_length();
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/Multiplication.cpp

@@ -0,0 +1,46 @@
+/*
+ * Copyright (c) 2020, Itamar S. <itamar8910@gmail.com>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+
+namespace Crypto {
+
+/**
+ * Complexity: O(N^2) where N is the number of words in the larger number
+ * Multiplication method:
+ * An integer is equal to the sum of the powers of two
+ * according to the indices of its 'on' bits.
+ * So to multiple x*y, we go over each '1' bit in x (say the i'th bit),
+ * and add y<<i to the result.
+ */
+FLATTEN void UnsignedBigIntegerAlgorithms::multiply_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& temp_shift_result,
+    UnsignedBigInteger& temp_shift_plus,
+    UnsignedBigInteger& temp_shift,
+    UnsignedBigInteger& output)
+{
+    output.set_to_0();
+
+    // iterate all bits
+    for (size_t word_index = 0; word_index < left.length(); ++word_index) {
+        for (size_t bit_index = 0; bit_index < UnsignedBigInteger::BITS_IN_WORD; ++bit_index) {
+            // If the bit is off - skip over it
+            if (!(left.m_words[word_index] & (1 << bit_index)))
+                continue;
+
+            size_t shift_amount = word_index * UnsignedBigInteger::BITS_IN_WORD + bit_index;
+
+            // output += (right << shift_amount);
+            shift_left_without_allocation(right, shift_amount, temp_shift_result, temp_shift_plus, temp_shift);
+            add_into_accumulator_without_allocation(output, temp_shift);
+        }
+    }
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/SimpleOperations.cpp

@@ -0,0 +1,106 @@
+/*
+ * Copyright (c) 2020, Itamar S. <itamar8910@gmail.com>
+ * Copyright (c) 2020-2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include "UnsignedBigIntegerAlgorithms.h"
+
+namespace Crypto {
+
+/**
+ * Complexity: O(N) where N is the number of words in the larger number
+ */
+void UnsignedBigIntegerAlgorithms::add_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& output)
+{
+    UnsignedBigInteger const* const longer = (left.length() > right.length()) ? &left : &right;
+    UnsignedBigInteger const* const shorter = (longer == &right) ? &left : &right;
+
+    output.set_to(*longer);
+    add_into_accumulator_without_allocation(output, *shorter);
+}
+
+/**
+ * Complexity: O(N) where N is the number of words in the larger number
+ */
+void UnsignedBigIntegerAlgorithms::add_into_accumulator_without_allocation(UnsignedBigInteger& accumulator, UnsignedBigInteger const& value)
+{
+    auto value_length = value.trimmed_length();
+
+    // If needed, resize the accumulator so it can fit the value.
+    accumulator.resize_with_leading_zeros(value_length);
+    auto final_length = accumulator.length();
+
+    // Add the words of the value into the accumulator, rippling any carry as we go
+    UnsignedBigInteger::Word last_carry_for_word = 0;
+    for (size_t i = 0; i < value_length; ++i) {
+        UnsignedBigInteger::Word current_carry_for_word = 0;
+        if (Checked<UnsignedBigInteger::Word>::addition_would_overflow(value.m_words[i], accumulator.m_words[i])) {
+            current_carry_for_word = 1;
+        }
+        UnsignedBigInteger::Word word_addition_result = value.m_words[i] + accumulator.m_words[i];
+        if (Checked<UnsignedBigInteger::Word>::addition_would_overflow(word_addition_result, last_carry_for_word)) {
+            current_carry_for_word = 1;
+        }
+        word_addition_result += last_carry_for_word;
+        last_carry_for_word = current_carry_for_word;
+        accumulator.m_words[i] = word_addition_result;
+    }
+
+    // Ripple the carry over the remaining words in the accumulator until either there is no carry left or we run out of words
+    while (last_carry_for_word && final_length > value_length) {
+        UnsignedBigInteger::Word current_carry_for_word = 0;
+        if (Checked<UnsignedBigInteger::Word>::addition_would_overflow(accumulator.m_words[value_length], last_carry_for_word)) {
+            current_carry_for_word = 1;
+        }
+        accumulator.m_words[value_length] += last_carry_for_word;
+        last_carry_for_word = current_carry_for_word;
+        value_length++;
+    }
+
+    if (last_carry_for_word) {
+        // Note : The accumulator couldn't add the carry directly, so we reached its end
+        accumulator.m_words.append(last_carry_for_word);
+    }
+}
+
+/**
+ * Complexity: O(N) where N is the number of words in the larger number
+ */
+void UnsignedBigIntegerAlgorithms::subtract_without_allocation(
+    UnsignedBigInteger const& left,
+    UnsignedBigInteger const& right,
+    UnsignedBigInteger& output)
+{
+    if (left < right) {
+        output.invalidate();
+        return;
+    }
+
+    u8 borrow = 0;
+    auto own_length = left.length();
+    auto other_length = right.length();
+
+    output.set_to_0();
+    output.m_words.resize_and_keep_capacity(own_length);
+
+    for (size_t i = 0; i < own_length; ++i) {
+        u32 other_word = (i < other_length) ? right.m_words[i] : 0;
+        i64 temp = static_cast<i64>(left.m_words[i]) - static_cast<i64>(other_word) - static_cast<i64>(borrow);
+        // If temp < 0, we had an underflow
+        borrow = (temp >= 0) ? 0 : 1;
+        if (temp < 0) {
+            temp += (UINT32_MAX + 1);
+        }
+        output.m_words[i] = temp;
+    }
+
+    // This assertion should not fail, because we verified that *this>=other at the beginning of the function
+    VERIFY(borrow == 0);
+}
+
+}

Libraries/LibCrypto/BigInt/Algorithms/UnsignedBigIntegerAlgorithms.h

@@ -0,0 +1,43 @@
+/*
+ * Copyright (c) 2021, Dex♪ <dexes.ttp@gmail.com>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#pragma once
+
+#include <LibCrypto/BigInt/UnsignedBigInteger.h>
+
+namespace Crypto {
+
+class UnsignedBigIntegerAlgorithms {
+    using Ops = AK::StorageOperations<UnsignedBigInteger::Word>;
+
+public:
+    static void add_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& output);
+    static void add_into_accumulator_without_allocation(UnsignedBigInteger& accumulator, UnsignedBigInteger const& value);
+    static void subtract_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& output);
+    static void bitwise_or_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& output);
+    static void bitwise_and_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& output);
+    static void bitwise_xor_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& output);
+    static void bitwise_not_fill_to_one_based_index_without_allocation(UnsignedBigInteger const& left, size_t, UnsignedBigInteger& output);
+    static void shift_left_without_allocation(UnsignedBigInteger const& number, size_t bits_to_shift_by, UnsignedBigInteger& temp_result, UnsignedBigInteger& temp_plus, UnsignedBigInteger& output);
+    static void shift_right_without_allocation(UnsignedBigInteger const& number, size_t num_bits, UnsignedBigInteger& output);
+    static void multiply_without_allocation(UnsignedBigInteger const& left, UnsignedBigInteger const& right, UnsignedBigInteger& temp_shift_result, UnsignedBigInteger& temp_shift_plus, UnsignedBigInteger& temp_shift, UnsignedBigInteger& output);
+    static void divide_without_allocation(UnsignedBigInteger const& numerator, UnsignedBigInteger const& denominator, UnsignedBigInteger& quotient, UnsignedBigInteger& remainder);
+    static void divide_u16_without_allocation(UnsignedBigInteger const& numerator, UnsignedBigInteger::Word denominator, UnsignedBigInteger& quotient, UnsignedBigInteger& remainder);
+
+    static void destructive_GCD_without_allocation(UnsignedBigInteger& temp_a, UnsignedBigInteger& temp_b, UnsignedBigInteger& temp_quotient, UnsignedBigInteger& temp_remainder, UnsignedBigInteger& output);
+    static void modular_inverse_without_allocation(UnsignedBigInteger const& a_, UnsignedBigInteger const& b, UnsignedBigInteger& temp_1, UnsignedBigInteger& temp_minus, UnsignedBigInteger& temp_quotient, UnsignedBigInteger& temp_d, UnsignedBigInteger& temp_u, UnsignedBigInteger& temp_v, UnsignedBigInteger& temp_x, UnsignedBigInteger& result);
+    static void destructive_modular_power_without_allocation(UnsignedBigInteger& ep, UnsignedBigInteger& base, UnsignedBigInteger const& m, UnsignedBigInteger& temp_1, UnsignedBigInteger& temp_2, UnsignedBigInteger& temp_3, UnsignedBigInteger& temp_multiply, UnsignedBigInteger& temp_quotient, UnsignedBigInteger& temp_remainder, UnsignedBigInteger& result);
+    static void montgomery_modular_power_with_minimal_allocations(UnsignedBigInteger const& base, UnsignedBigInteger const& exponent, UnsignedBigInteger const& modulo, UnsignedBigInteger& temp_z0, UnsignedBigInteger& temp_rr, UnsignedBigInteger& temp_one, UnsignedBigInteger& temp_z, UnsignedBigInteger& temp_zz, UnsignedBigInteger& temp_x, UnsignedBigInteger& temp_extra, UnsignedBigInteger& result);
+
+private:
+    static UnsignedBigInteger::Word montgomery_fragment(UnsignedBigInteger& z, size_t offset_in_z, UnsignedBigInteger const& x, UnsignedBigInteger::Word y_digit, size_t num_words);
+    static void almost_montgomery_multiplication_without_allocation(UnsignedBigInteger const& x, UnsignedBigInteger const& y, UnsignedBigInteger const& modulo, UnsignedBigInteger& z, UnsignedBigInteger::Word k, size_t num_words, UnsignedBigInteger& result);
+    static void shift_left_by_n_words(UnsignedBigInteger const& number, size_t number_of_words, UnsignedBigInteger& output);
+    static void shift_right_by_n_words(UnsignedBigInteger const& number, size_t number_of_words, UnsignedBigInteger& output);
+    ALWAYS_INLINE static UnsignedBigInteger::Word shift_left_get_one_word(UnsignedBigInteger const& number, size_t num_bits, size_t result_word_index);
+};
+
+}