<repo_name>Ivan-Veselov/coding-hw<file_sep>/numbers.cpp
#include "numbers.h"

bool isPrime(Natural number) {
    if (number < 2) {
        return false;
    }

    for (Natural divisor = 2; Natural64(divisor) * Natural64(divisor) <= number; ++divisor) {
        if (number % divisor == 0) {
            return false;
        }
    }

    return true;
}
<file_sep>/main.cpp
#include <iostream>

#include "FinitePrimeFieldElement.h"
#include "Polynomial.h"

void findAllPrimePolynomials() {
    auto polynomials = Polynomial<GF<2>>::allPolynomials(3);

    for (const auto &polynomial : polynomials) {
        if (polynomial == Polynomial<GF<2>>::ZERO || polynomial[polynomial.deg()] != GF<2>::ONE) {
            continue;
        }

        bool isDivisible = false;
        for (const auto &other : polynomials) {
            if (other.deg() < 1 || other.deg() >= polynomial.deg()) {
                continue;
            }

            if (polynomial % other == Polynomial<GF<2>>::ZERO) {
                std::cout << polynomial.formulaString() << " = (" << other.formulaString() << ") * (" << (polynomial / other).formulaString() << ")" << std::endl;
                isDivisible = true;
                break;
            }
        }

        if (!isDivisible) {
            std::cout << polynomial.formulaString() << " is a prime polynomial" << std::endl;
        }
    }
}

void solve() {
    // x^3 + x + 1 for 2**3
    // x^2 + x + 2 for 3**2

    auto polynomials = Polynomial<GF<3>>::allPolynomials(1);

    auto x = Polynomial<GF<3>>::X;
    auto prime = x(2) + x(1) + 2;

    Natural power = 3 * 3 - 2;

    for (const auto &elem : polynomials) {
        std::cout << "The inverse of " << elem.formulaString() << " is " << std::endl;
        std::cout << "(" << elem.formulaString() << ")**" << power << " = " << std::endl;

        auto inverse = elem;
        for (int i = 0; i < power - 1; ++i) {
            inverse *= elem;
        }

        inverse %= prime;

        std::cout << inverse.formulaString() << std::endl;
        std::cout << "The product is " << ((elem * inverse) % prime).formulaString() << "\n" << std::endl;
    }
}

int main() {
    solve();

    /*auto x = Polynomial<GF<2>>::X;
    auto prime = x(3) + x(1) + 1;

    auto pol = x(1) + 1;
    std::cout << ((pol * pol * pol * pol * pol * pol) % prime).prettyString() << std::endl;
    std::cout << ((fastPower(pol, 6)) % prime).prettyString() << std::endl;*/

    return 0;
}<file_sep>/CMakeLists.txt
cmake_minimum_required(VERSION 3.13)
project(project)

set(CMAKE_CXX_STANDARD 17)

add_executable(project main.cpp FinitePrimeFieldElement.h numbers.h numbers.cpp Polynomial.h)<file_sep>/numbers.h
#pragma once

#include <cstdint>

typedef uint32_t Natural32;
typedef uint64_t Natural64;
typedef Natural32 Natural;

typedef int32_t Whole32;
typedef int64_t Whole64;
typedef Whole32 Whole;

bool isPrime(Natural number);

template<typename Number>
Number multiplicationNeutralElement() {
    return Number(1);
}

template <typename Number>
Number fastPower(Number value, Natural power) {
    if (power == 0) {
        return multiplicationNeutralElement<Number>();
    }

    if (power == 1) {
        return value;
    }

    Number result = fastPower(value, power / 2);
    result *= result;

    if (power % 2 == 1) {
        result *= value;
    }

    return result;
}
<file_sep>/Polynomial.h
#pragma once

#include <vector>

#include "numbers.h"

template <typename Coefficient>
class Polynomial;

template <typename Coefficient>
bool operator == (const Polynomial<Coefficient>&, const Polynomial<Coefficient>&);

template <typename Coefficient>
bool operator != (const Polynomial<Coefficient>&, const Polynomial<Coefficient>&);

template <typename Coefficient>
class Polynomial {
    friend bool operator == <>(const Polynomial<Coefficient>&, const Polynomial<Coefficient>&);
    friend bool operator != <>(const Polynomial<Coefficient>&, const Polynomial<Coefficient>&);

    private:
        std::vector<Coefficient> coefficients;
        // coefficients[i] is a coefficient for x^i

    public:
        class Variable {
            public:
                Polynomial<Coefficient> operator() (const Natural &power) const {
                    return Polynomial<Coefficient>(power);
                }
        };

        static const Variable X;

    public:
        static const Polynomial ZERO;
        static const Polynomial ONE;

    public:
        Polynomial()
        : coefficients({Coefficient::ZERO}) {
        }

        Polynomial(const Coefficient &scalar)
        : coefficients({scalar}) {
        }

        explicit Polynomial(const Natural &power)
        : coefficients(power + 1, Coefficient::ZERO) {
            coefficients[power] = Coefficient::ONE;
        }

        static std::vector<Polynomial> allPolynomials(Natural maxDeg) {
            auto allCoefficients = Coefficient::allElements();

            std::vector<Polynomial> polynomials;

            if (maxDeg == 0) {
                for (const auto &coefficient : allCoefficients) {
                    polynomials.push_back(Polynomial(coefficient));
                }

                return polynomials;
            }

            std::vector<Polynomial> lowerDegreePolynomials = allPolynomials(maxDeg - 1);
            polynomials.insert(polynomials.end(), lowerDegreePolynomials.begin(), lowerDegreePolynomials.end());

            for (const auto &coefficient : allCoefficients) {
                if (coefficient == Coefficient::ZERO) {
                    continue;
                }

                for (auto lowerPolynomial : lowerDegreePolynomials) {
                    lowerPolynomial += Polynomial(maxDeg) * coefficient;
                    polynomials.push_back(lowerPolynomial);
                }
            }

            return polynomials;
        }

        Natural deg() const {
            return coefficients.size() - 1;
        }

        const Coefficient& operator [] (const Natural &degree) const {
            return coefficients[degree];
        }

        Polynomial& operator += (const Polynomial &other) {
            const auto numOfConsideredCoefficients = other.coefficients.size();
            addZeroCoefficientsUpTo(numOfConsideredCoefficients);

            for (auto i = 0; i < numOfConsideredCoefficients; ++i) {
                coefficients[i] += other[i];
            }

            normalize();
            return *this;
        }

        Polynomial& operator -= (const Polynomial &other) {
            const auto numOfConsideredCoefficients = other.coefficients.size();
            addZeroCoefficientsUpTo(numOfConsideredCoefficients);

            for (auto i = 0; i < numOfConsideredCoefficients; ++i) {
                coefficients[i] -= other[i];
            }

            normalize();
            return *this;
        }

        Polynomial& operator *= (const Coefficient &scalar) {
            for (auto &c : coefficients) {
                c *= scalar;
            }

            normalize();
            return *this;
        }

        Polynomial& operator *= (const Polynomial &other) {
            const auto numOfConsideredCoefficients = deg() + other.deg() + 1;
            addZeroCoefficientsUpTo(numOfConsideredCoefficients);

            for (int i = numOfConsideredCoefficients - 1; i >= 0; --i) {
                // i is an index of resulting multiplication
                for (int j = std::min(i, int(deg())); j >= 0 && i - j <= other.deg(); --j) {
                    // j is an index of original polynomial

                    if (j != i) {
                        coefficients[i] += coefficients[j] * other[i - j];
                    } else {
                        coefficients[i] = coefficients[j] * other[i - j];
                    }
                }
            }

            normalize();
            return *this;
        }

        // this may be passed as a 'remainder' or 'quotient'
        void divide(const Polynomial &divisor, Polynomial &quotient, Polynomial &remainder) const {
            assert(divisor != ZERO);

            remainder = *this;

            if (deg() < divisor.deg()) {
                quotient = ZERO;
                return;
            }

            Natural divisorDeg = divisor.deg();
            Natural quotientPower = deg() - divisorDeg + 1;
            quotient.coefficients.resize(quotientPower);

            Coefficient divisorHighCoefficientInv = divisor.coefficients.back().inv();

            do {
                --quotientPower;

                Coefficient quotientCoefficient = remainder[quotientPower + divisorDeg] * divisorHighCoefficientInv;
                quotient[quotientPower] = quotientCoefficient;

                for (auto i = 0; i <= divisorDeg; ++i) {
                    remainder[quotientPower + i] -= divisor[i] * quotientCoefficient;
                }

            } while (quotientPower != 0);

            remainder.normalize();
        }

        Polynomial& operator /= (const Polynomial &other) {
            Polynomial remainder;
            divide(other, *this, remainder);

            return *this;
        }

        Polynomial& operator %= (const Polynomial &other) {
            Polynomial quotient;
            divide(other, quotient, *this);

            return *this;
        }

        std::string formulaString() const {
            std::vector<std::string> terms;
            for (auto power = coefficients.size() - 1; power > 0; --power) {
                auto c = coefficients[power];
                if (c == Coefficient::ZERO) {
                    continue;
                }

                std::string termStr;
                if (c != Coefficient::ONE) {
                    termStr = c.formulaString() + " * ";
                }

                termStr += "x";

                if (power > 1) {
                    termStr += "^" + std::to_string(power);
                }

                terms.push_back(termStr);
            }

            auto bias = coefficients[0];
            if (bias != Coefficient::ZERO || terms.empty()) {
                terms.push_back(bias.formulaString());
            }

            std::string result;
            for (auto i = 0; i < terms.size(); ++i) {
                if (i > 0) {
                    result += " + ";
                }

                result += terms[i];
            }

            return result;
        }

        std::string prettyString() const {
            return formulaString() + ", a_i is from " + Coefficient::setPrettyString();
        }

    private:
        Coefficient& operator [] (const Natural &degree) {
            return coefficients[degree];
        }

        void normalize() {
            while (coefficients.back() == Coefficient::ZERO && coefficients.size() > 1) {
                coefficients.pop_back();
            }
        }

        void addZeroCoefficientsUpTo(const Natural numberOfCoefficients) {
            if (numberOfCoefficients > coefficients.size()) {
                coefficients.resize(numberOfCoefficients, Coefficient::ZERO);
            }
        }
};

template <typename Coefficient>
const Polynomial<Coefficient> Polynomial<Coefficient>::ZERO = Coefficient::ZERO;

template <typename Coefficient>
const Polynomial<Coefficient> Polynomial<Coefficient>::ONE  = Coefficient::ONE;

template <typename Coefficient>
Polynomial<Coefficient> operator + (Polynomial<Coefficient> element1, const Polynomial<Coefficient> &element2) {
    return element1 += element2;
}

template <typename Coefficient>
Polynomial<Coefficient> operator - (Polynomial<Coefficient> element1, const Polynomial<Coefficient> &element2) {
    return element1 -= element2;
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator + (Polynomial<Coefficient> element, const Scalar &scalar) {
    return element += Coefficient(scalar);
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator + (const Scalar &scalar, Polynomial<Coefficient> element) {
    return element += Coefficient(scalar);
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator - (Polynomial<Coefficient> element, const Scalar &scalar) {
    return element -= Coefficient(scalar);
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator - (const Scalar &scalar, Polynomial<Coefficient> element) {
    return element -= Coefficient(scalar);
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator * (Polynomial<Coefficient> element, const Scalar &scalar) {
    return element *= scalar;
}

template <typename Coefficient, typename Scalar>
Polynomial<Coefficient> operator * (const Scalar &scalar, Polynomial<Coefficient> element) {
    return element *= scalar;
}

template <typename Coefficient>
Polynomial<Coefficient> operator * (Polynomial<Coefficient> element1, const Polynomial<Coefficient> &element2) {
    return element1 *= element2;
}

template <typename Coefficient>
Polynomial<Coefficient> operator / (Polynomial<Coefficient> element1, const Polynomial<Coefficient> &element2) {
    return element1 /= element2;
}

template <typename Coefficient>
Polynomial<Coefficient> operator % (Polynomial<Coefficient> element1, const Polynomial<Coefficient> &element2) {
    return element1 %= element2;
}

template <typename Coefficient>
bool operator == (const Polynomial<Coefficient> &element1, const Polynomial<Coefficient> &element2) {
    return element1.coefficients == element2.coefficients;
}

template <typename Coefficient>
bool operator != (const Polynomial<Coefficient> &element1, const Polynomial<Coefficient> &element2) {
    return element1.coefficients != element2.coefficients;
}

namespace Polynomials {
    template <typename Coefficient>
    Polynomial<Coefficient> X(const Natural &power) {
        return Polynomial<Coefficient>(power);
    }
}

/*
https://www.fluentcpp.com/2017/08/15/function-templates-partial-specialization-cpp/

template <typename Coefficient>
Polynomial<Coefficient> multiplicationNeutralElement() {
    return Polynomial<Coefficient>(Coefficient::ONE);
}
*/
<file_sep>/FinitePrimeFieldElement.h
#pragma once

#include <string>
#include <cassert>
#include <vector>

#include "numbers.h"

template <Natural P>
class FinitePrimeFieldElement;

template <Natural P>
using GF = FinitePrimeFieldElement<P>;

template <Natural P>
bool operator == (const FinitePrimeFieldElement<P>&, const FinitePrimeFieldElement<P>&);

template <Natural P>
bool operator != (const FinitePrimeFieldElement<P>&, const FinitePrimeFieldElement<P>&);

template <Natural P>
class FinitePrimeFieldElement {
    friend bool operator == <>(const FinitePrimeFieldElement<P>&, const FinitePrimeFieldElement<P>&);
    friend bool operator != <>(const FinitePrimeFieldElement<P>&, const FinitePrimeFieldElement<P>&);

    private:
        static bool isChecked;

        Whole remainder;

    public:
        static const FinitePrimeFieldElement ZERO;
        static const FinitePrimeFieldElement ONE;

    public:
        FinitePrimeFieldElement()
        : remainder(0) {
            if (!isChecked) {
                check();
            }
        }

        FinitePrimeFieldElement(Whole number)
        : remainder(toRemainderClassMainRepresentative(number)) {
            if (!isChecked) {
                check();
            }
        }

        static std::vector<FinitePrimeFieldElement> allElements() {
            std::vector<FinitePrimeFieldElement> elements;
            for (Whole remainder = 0; remainder < P; ++remainder) {
                elements.push_back(FinitePrimeFieldElement(remainder));
            }

            return elements;
        }

        FinitePrimeFieldElement inv() const {
            assert(remainder != 0);
            return fastPower(*this, P - 2);
        }

        FinitePrimeFieldElement& operator += (const FinitePrimeFieldElement &other) {
            remainder += other.remainder;
            if (remainder >= P) {
                remainder -= P;
            }

            return *this;
        }

        FinitePrimeFieldElement& operator -= (const FinitePrimeFieldElement &other) {
            remainder -= other.remainder;
            if (remainder < 0) {
                remainder += P;
            }

            return *this;
        }

        FinitePrimeFieldElement& operator *= (const FinitePrimeFieldElement &other) {
            Whole64 tmp = Whole64(remainder) * other.remainder;
            remainder = Whole(tmp % P);

            return *this;
        }

        FinitePrimeFieldElement& operator /= (const FinitePrimeFieldElement &other) {
            return *this *= other.inv();
        }

        std::string prettyString() const {
            return std::to_string(remainder) + " (mod " + std::to_string(P) + ")";
        }

        std::string formulaString() const {
            return std::to_string(remainder);
        }

        static std::string setPrettyString() {
            return "GF(" + std::to_string(P) + ")";
        }

    private:
        static void check() {
            assert(isPrime(P));
        }

        Whole toRemainderClassMainRepresentative(Whole aRemainder) {
            Whole representative = aRemainder % P;
            if (representative < 0) {
                representative += P;
            }

            return representative;
        }
};

template <Natural P>
bool FinitePrimeFieldElement<P>::isChecked = false;

template <Natural P>
const FinitePrimeFieldElement<P> FinitePrimeFieldElement<P>::ZERO = 0;

template <Natural P>
const FinitePrimeFieldElement<P> FinitePrimeFieldElement<P>::ONE  = 1;

template <Natural P>
FinitePrimeFieldElement<P> operator + (FinitePrimeFieldElement<P> element1, const FinitePrimeFieldElement<P> &element2) {
    return element1 += element2;
}

template <Natural P>
FinitePrimeFieldElement<P> operator - (FinitePrimeFieldElement<P> element1, const FinitePrimeFieldElement<P> &element2) {
    return element1 -= element2;
}

template <Natural P>
FinitePrimeFieldElement<P> operator * (FinitePrimeFieldElement<P> element1, const FinitePrimeFieldElement<P> &element2) {
    return element1 *= element2;
}

template <Natural P>
FinitePrimeFieldElement<P> operator / (FinitePrimeFieldElement<P> element1, const FinitePrimeFieldElement<P> &element2) {
    return element1 /= element2;
}

template <Natural P>
bool operator == (const FinitePrimeFieldElement<P> &element1, const FinitePrimeFieldElement<P> &element2) {
    return element1.remainder == element2.remainder;
}

template <Natural P>
bool operator != (const FinitePrimeFieldElement<P> &element1, const FinitePrimeFieldElement<P> &element2) {
    return element1.remainder != element2.remainder;
}

/*
https://www.fluentcpp.com/2017/08/15/function-templates-partial-specialization-cpp/

template <Natural P>
FinitePrimeFieldElement<P> multiplicationNeutralElement() {
    return FinitePrimeFieldElement<P>(1);
}
*/