ec.hpp

// Copyright © 2023 ZeroPass <zeropass@pm.me>
// Author: Crt Vavros
#pragma once
#include <cstdint>
#include <span>
#include <type_traits>

#include <ack/bigint.hpp>
#include <ack/fp.hpp>
#include <ack/types.hpp>
#include <ack/utils.hpp>

namespace ack {
    template<std::size_t N>
    using ec_fixed_bigint = fixed_bigint<bitsize_to_wordsize(N) * word_bit_size * 2>; // 2x size required for multiplication

    // Affine coordinates representation of an elliptic curve point
    template<typename PointT, typename CurveT>
    struct ec_point_base
    {
        using int_type = typename CurveT::int_type;

        constexpr ec_point_base() : // represents point at infinity
            curve_( nullptr )
        {}

        constexpr ec_point_base(const PointT& p) = delete;
        constexpr ec_point_base(PointT&& p) noexcept = delete;
        constexpr PointT& operator=(const PointT& p) = delete;
        constexpr PointT& operator=(PointT&& p) noexcept = delete;

        /**
         * Returns the curve this point belongs to.
         * @warning If this point is the identity element of the curve, the curve can be nullptr.
         *          Make sure to check if this point is not identity before calling this method.
         * @return the curve this point belongs to.
        */
        const CurveT& curve() const
        {
            check( curve_ != nullptr, "curve is null" );
            return *curve_;
        }

        /**
         * Checks if this point is the identity element of the curve, i.e. point at infinity.
         * @return true if this point is the identity element of the curve, false otherwise
        */
        inline constexpr bool is_identity() const
        {
            return underlying().is_identity();
        }

        /**
         * Checks if this point is on the curve.
         * @return true if this point is on the curve, false otherwise
        */
        [[nodiscard]] inline bool is_on_curve() const
        {
            return underlying().is_on_curve();
        }

        /**
         * Checks if this point is valid point generated by the generator point.
         * @note Should do SEC 1 section 3.2.2.1 like verification.
         * @return true if this point is valid, false otherwise
        */
        [[nodiscard]] inline bool is_valid() const
        {
            return underlying().is_valid();
        }

        /**
         * Returns the inverse of this point.
         * R = -this
         *
         * @return the inverse of this point
        */
        [[nodiscard]] inline constexpr PointT inverted() const
        {
            return underlying().inverted();
        }

        /**
         * Adds the given point to this point.
         * R = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        [[nodiscard]] inline PointT add(const PointT& a) const
        {
            return underlying().add( a );
        }

        /**
         * Adds the given point to this point.
         * R = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        [[nodiscard]] inline PointT add(const ec_point_base& a) const
        {
            return add( a.underlying() );
        }

        /**
         * Returns the double of this point.
         * R = 2 * this
         *
         * @return the double of this point
        */
        [[nodiscard]] inline PointT doubled() const
        {
            return underlying().doubled();
        }

        /**
         * Subtracts the given point from this point.
         * R = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        [[nodiscard]] inline PointT sub(const PointT& a) const
        {
            return underlying().sub( a );
        }

        /**
         * Subtracts the given point from this point.
         * R = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        [[nodiscard]] inline PointT sub(const ec_point_base& a) const
        {
            return sub( a.underlying() );
        }

        /**
         * Multiplies this point by the given scalar.
         * R = this * scalar
         *
         * @param scalar - the scalar to multiply this point by
         * @return the resulting point
        */
        [[nodiscard]] inline PointT mul(const int_type& scalar) const
        {
            return underlying().mul( scalar );
        }

        /**
         * Adds the given point to this point.
         * this = this + a
         *
         * @param a - the point to add other point
         * @param b - the point to add to the other point
         * @return reference to this point
        */
        [[nodiscard]] friend inline PointT operator + (const PointT& a, const PointT& b)
        {
            return a.add( b );
        }

        [[nodiscard]] friend inline PointT operator + (const ec_point_base& a, const ec_point_base& b)
        {
            return a.add( b );
        }

        /**
         * Adds the given point to this point.
         * this = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        inline PointT& operator += (const PointT& a)
        {
            return underlying() = add( a );
        }

        /**
         * Adds the given point to this point.
         * this = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        inline PointT& operator += (const ec_point_base& a)
        {
            return this-> operator += ( a.underlying() );
        }

        /**
         * Subtracts the given point from this point.
         * R = P - Q
         *
         * @param a - the point to subtract other point from
         * @param b - the point to subtract from the other point
         * @return the result point of the subtraction
        */
        [[nodiscard]] friend inline PointT operator - (const PointT& a, const PointT& b)
        {
            return a.sub( b );
        }

        /**
         * Subtracts the given point from this point.
         * R = P - Q
         *
         * @param a - the point to subtract other point from
         * @param b - the point to subtract from the other point
         * @return the result point of the subtraction
        */
        [[nodiscard]] friend inline PointT operator - (const ec_point_base& a, const ec_point_base& b)
        {
            return a.sub( b );
        }

        /**
         * Subtracts the given point from this point.
         * this = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        inline PointT& operator -= (const PointT& a)
        {
            return underlying() = sub( a );
        }

        /**
         * Subtracts the given point from this point.
         * this = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        inline PointT& operator -= (const ec_point_base& a)
        {
            return this->operator -= ( a.underlying() );
        }

        /**
         * Multiplies given point by the given scalar.
         * R = P * s
         *
         * @param p - the point to multiply
         * @param s - the scalar to multiply the point by
         * @return the result point of the multiplication
        */
        [[nodiscard]] friend inline PointT operator * (const PointT& p, const int_type& s)
        {
            return p.mul( s );
        }

        /**
         * Multiplies given point by the given scalar.
         * R = P * s
         *
         * @param p - the point to multiply
         * @param s - the scalar to multiply the point by
         * @return the result point of the multiplication
        */
        [[nodiscard]] friend inline PointT operator * (const ec_point_base& p, const int_type& s)
        {
            return p.mul( s );
        }

        /**
         * Multiplies the given scalar by this point.
         * R = s * P
         *
         * @param s - the scalar to multiply this point by
         * @param p - the point to multiply
         * @return the result point of the multiplication
        */
        [[nodiscard]] friend inline PointT operator * (const int_type& s, const PointT& p)
        {
            return p.mul( s );
        }

        /**
         * Multiplies the given scalar by this point.
         * R = s * P
         *
         * @param s - the scalar to multiply this point by
         * @param p - the point to multiply
         * @return the result point of the multiplication
        */
        [[nodiscard]] friend inline PointT operator * (const int_type& s, const ec_point_base& p)
        {
            return p.mul( s );
        }

        /**
         * Multiplies this point by the given scalar.
         * this = this * s
         *
         * @param s - the scalar to multiply this point by.
         * @return reference to this point
        */
        inline PointT& operator *= (const int_type& s)
        {
            return underlying() = mul( s );
        }

        /**
         * Returns the inverse of this point.
         * R = -this
         *
         * @return the inverse of this point
        */
        [[nodiscard]] inline constexpr PointT operator - () const
        {
            return inverted();
        }

        private:
            constexpr ec_point_base( const CurveT& curve ) :
                curve_( &curve )
            {}

            inline constexpr PointT& underlying()
            {
                return static_cast<PointT&>( *this );
            }

            inline constexpr const PointT& underlying() const
            {
                return static_cast<const PointT&>( *this );
            }

        private:
            friend PointT;
            friend CurveT;
            const CurveT* curve_;
    };

    // Forward declaration of ec_curve_fp
    template<typename IntT, typename CurveTag>
    struct ec_curve_fp;

    /**
     * Type trait to check if the given curve type is derived from ec_curve_fp
     * @tparam CurveT - the curve type to check.
    */
    template<typename CurveT>
    constexpr bool is_ec_curve_fp = std::is_same_v<CurveT, ec_curve_fp<typename CurveT::int_type, typename CurveT::curve_tag>>;

    /**
     * Struct represents a affine point on an elliptic curve
     * over a prime finite field GF(p) with short Weierstrass equation:
     * y^2 = x^3 + ax + b
     *
     * The implementation follows the algorithms described in:
     *  - SECG standards SEC 1: Elliptic Curve Cryptography, Version 2.0
     *    https://www.secg.org/sec1-v2.pdf
     *  - RFC-6090: https://www.rfc-editor.org/rfc/rfc6090
     *
     * @warning The point's curve is stored as a pointer to the curve object.
     *          The curve object must outlive the point object.
     *
     * @tparam CurveT - the curve type. Required to be an instance of ec_curve_fp.
    */
    template<typename CurveT>
    struct ec_point_fp : ec_point_base<ec_point_fp<CurveT>, CurveT>
    {
        static_assert( is_ec_curve_fp<std::remove_cv_t<CurveT>> );

        using base_type          = ec_point_base<ec_point_fp<CurveT>, CurveT>;
        using int_type           = typename CurveT::int_type;
        using field_element_type = typename CurveT::field_element_type;
        using base_type::base_type;

        field_element_type x;
        field_element_type y;

        /**
         * Constructs a point at infinity
        */
        constexpr ec_point_fp() :
            base_type(),
            x( field_element_type::zero() ),
            y( field_element_type::zero() )
        {}

        /**
         * Checks if this point is the identity element of the curve, i.e. point at infinity.
         * @return true if this point is the identity element of the curve, false otherwise
        */
        inline constexpr bool is_identity() const
        {
            return this->curve_ == nullptr || ( x.is_zero() && y.is_zero() );
        }

        /**
         * Checks if this point is on the curve by calculating
         * the left and right hand side of the equation:  y^2 = x^3 + ax + b
         *
         * @return true if this point is on the curve, false otherwise
        */
        [[nodiscard]] bool is_on_curve() const
        {
            if ( is_identity() ) {
                return true;
            }
            return ( y.sqr() - ( ( x.sqr() + this->curve().a ) * x + this->curve().b )) == 0;
        }

        /**
         * Performs SEC 1 section 3.2.2.1 point check,
         * i.e. if point was generated using curve generator g and is not identity element.
         * Checks:
         *    - point is not the identity element
         *    - point is on the curve
         *    - point has order n
         *
         * @note Slow operation.
         *
         * @return true if this point is valid, false otherwise
        */
        [[nodiscard]] bool is_valid() const
        {
            if ( is_identity() ) {
                return false;
            }

            if ( !is_on_curve() ) {
                return false;
            }

            if ( !(this->curve().n * *this).is_identity() ) {
                return false;
            }

            return true;
        }

        /**
         * Returns the inverse of this point.
         * R = -this
         *
         * @return the inverse of this point
        */
        [[nodiscard]] constexpr ec_point_fp inverted() const
        {
            if ( is_identity() ) {
                return *this;
            }
            return ec_point_fp( this->curve(), x, -y );
        }

        /**
         * Adds the given point to this point.
         * R = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        [[nodiscard]] ec_point_fp add(const ec_point_fp& a) const
        {
            if ( is_identity() ) {
                return a;
            }

            if ( a.is_identity() ) {
                return *this;
            }

            // TODO: before this sub operation caused error in wasm:  memcpy with overlapping memory memcpy can only accept non-aliasing pointers
            auto s = a.x - x;
            if ( s.is_zero() ) {
                if ( y == a.y ) { // double point
                    return this->doubled();
                }
                return ec_point_fp(); // point at infinity
            }

            // Calculate tangent slope
            s = ( a.y - y ) / s;

            // Calculate new x and y
            auto x3 = s.sqr() - x - a.x;
            auto y3 = s * ( x - x3 ) - y;
            return ec_point_fp( this->curve(), std::move( x3 ), std::move( y3 ) );
        }

        /**
         * Returns the double of this point.
         * R = 2 * this
         *
         * @return the double of this point
        */
        [[nodiscard]] ec_point_fp doubled() const
        {
            if ( is_identity() || y.is_zero() ) { // check for y == 0 handles division by zero issue
                return ec_point_fp();
            }

            // Calculate tangent slope
            const auto s = ( 3 * x.sqr() + this->curve().a ) / ( 2 * y ) ;

            // Calculate new x and y
            auto x2 = s.sqr() - 2 * x;
            auto y2 = s * ( x - x2 ) - y;
            return ec_point_fp( this->curve(), std::move( x2 ), std::move( y2 ) );
        }

        /**
         * Subtracts the given point from this point.
         * R = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        [[nodiscard]] inline ec_point_fp sub(const ec_point_fp& a) const
        {
            return *this + (-a);
        }

        /**
         * Multiplies this point by the given scalar using double and add algorithm.
         * R = this * scalar
         *
         * @param scalar - the scalar to multiply this point by
         * @return the resulting point
        */
        [[nodiscard]] ec_point_fp mul(const int_type& scalar) const
        {
            if ( scalar.is_one() || is_identity() ) {
                return *this;
            }

            if ( scalar.is_zero() ) {
                return ec_point_fp();
            }

            if ( scalar < 0 ) {
                return inverted() * -scalar;
            }

            auto r   = ec_point_fp();
            auto tmp = *this;
            auto s   = scalar;
            while ( s != 0 ) {
                if ( s.is_odd() ) {
                    r += tmp;
                }
                tmp = tmp.doubled();
                s >>= 1;
            }
            return r;
        }

        /**
         * Compares 2 points for equality.
         * Points are equal if they have the same x and y coordinates.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are equal, false otherwise
        */
        constexpr friend bool operator == (const ec_point_fp& a, const ec_point_fp& b)
        {
            if ( a.x == b.x && a.y == b.y ) {
                return ( a.curve_ == b.curve_  ) || a.x == 0; // a.x == 0 means point at infinity
            }
            return false;
        }

        /**
         * Compares 2 points for inequality.
         * @note See operator == for details.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are not equal, false otherwise
        */
        constexpr friend inline bool operator != (const ec_point_fp& a, const ec_point_fp& b)
        {
            return !(a == b);
        }

        private:
            friend CurveT;
            template<typename>
            friend struct ec_point_fp_proj;
            template<typename>
            friend struct ec_point_fp_jacobi;

            constexpr ec_point_fp( const CurveT& curve, field_element_type x, field_element_type y ) :
                base_type( curve ),
                x( std::move(x) ),
                y( std::move(y) )
            {}
    };

    /**
     * Struct represents a point on an elliptic curve in standard projective coordinates (homogeneous coordinates)
     * over a prime finite field GF(p) with short Weierstrass equation: y^2 * z = x^3 + ax * z^2 + b * z^3
     * Due to the nature of projective coordinates the arithmetic operations are more efficient than in affine coordinates.
     * i.e. addition and doubling of points are faster due to the lack of division, which is deferred to conversion to affine coordinates.
     *
     * Implementation follows the RFC 6090: https://www.rfc-editor.org/rfc/rfc6090
     *
     * @warning The point's curve is stored as a pointer to the curve object.
     *          The curve object must outlive the point object.
     *
     * @tparam CurveT - the curve type. Required to be an instance of ec_curve_fp.
    */
    template<typename CurveT>
    struct ec_point_fp_proj : ec_point_base<ec_point_fp_proj<CurveT>, CurveT>
    {
        static_assert( is_ec_curve_fp<std::remove_cv_t<CurveT>> );

        using base_type          = ec_point_base<ec_point_fp_proj<CurveT>, CurveT>;
        using int_type           = typename CurveT::int_type;
        using field_element_type = typename CurveT::field_element_type;
        using affine_point_type  = ec_point_fp<CurveT>;
        using base_type::base_type;

        field_element_type x;
        field_element_type y;
        field_element_type z;

        /**
         * Constructs a point at infinity
        */
        constexpr ec_point_fp_proj() :
            base_type(),
            x( field_element_type::zero() ),
            y( field_element_type::one()  ),
            z( field_element_type::zero() )
        {}

        /**
         * Constructs this point from the given affine point.
         * @warning The point's curve is stored as a pointer to the curve object.
         *          The curve object must outlive the point object.
         *
         * @param p - the affine point to construct this point from.
        */
        explicit constexpr ec_point_fp_proj(affine_point_type p) :
            ec_point_fp_proj()
        {
            if ( !p.is_identity() ) {
                this->curve_ = &p.curve();
                x = std::move(p.x);
                y = std::move(p.y);
                z = field_element_type( 1, p.curve().p );
            }
        }

        /**
         * Normalizes this point.
         * This ensures Z coordinate is 1 therefore the x, y coordinates reflect
         * those of the equivalent to point in an affine coordinate system.
         *
         * @return this point
        */
        ec_point_fp_proj& normalize()
        {
            if ( !is_identity() && !z.is_one() ) {
                const auto z_inv = z.inv();
                x *= z_inv;
                y *= z_inv;
                z = 1;
            }
            return *this;
        }

        /**
         * Returns the normalized representation of this point.
         * This ensures Z coordinate is 1 therefore the x, y coordinates reflect
         * those of the equivalent to point in an affine coordinate system.
        */
        [[nodiscard]] ec_point_fp_proj normalized() const
        {
            auto r = *this;
            r.normalize();
            return r;
        }

        /**
         * Returns the affine representation of this point.
         * @note No point verification is performed.
         * @note Slow operation due to division.
         *
         * @return this point in affine form.
        */
        [[nodiscard]] const affine_point_type to_affine() const
        {
            if ( is_identity() ) {
                return affine_point_type();
            }
            if ( z.is_one() ) {
                return affine_point_type( this->curve(), x, y );
            }

            // It is assumed that the point is on the curve and therefore calculated x and y are valid.
            // Calling affine_point_type() constructor will skip verification check.
            const auto z_inv = z.inv();
            return affine_point_type( this->curve(), x * z_inv, y * z_inv );
        }

        /**
         * Checks if this point is the identity element of the curve, i.e. point at infinity.
         * @return true if this point is the identity element of the curve, false otherwise
        */
        [[nodiscard]] constexpr bool is_identity() const
        {
            return this->curve_ == nullptr || ( z.is_zero() );
        }

        /**
         * Checks if this point is on the curve by calculating
         * the left and right hand side of the equation:  y^2 * z = x^3 + ax * z^2 + b * z^3
         *
         * @return true if this point is on the curve, false otherwise
        */
        [[nodiscard]] bool is_on_curve() const
        {
            if ( is_identity() ) {
                return true;
            }
            const auto z2 = z.sqr();
            return ( y.sqr() * z - ( ( ( x.sqr() + this->curve().a * z2 ) * x  + this->curve().b * z * z2 ) ) ) == 0;
        }

        /**
         * Performs SEC 1 section 3.2.2.1 point check,
         * i.e. if point was generated using curve generator g and is not identity element.
         * Checks:
         *    - point is not the identity element
         *    - point is on the curve
         *    - point has order n
         *
         * @note Slow operation.
         *
         * @return true if this point is valid, false otherwise
        */
        [[nodiscard]] bool is_valid() const
        {
            if ( is_identity() ) {
                return false;
            }

            if ( !is_on_curve() ) {
                return false;
            }

            if ( !(this->curve().n * *this).is_identity() ) {
                return false;
            }

            return true;
        }

        /**
         * Returns the inverse of this point.
         * R = -this
         *
         * @return the inverse of this point
        */
        [[nodiscard]] constexpr ec_point_fp_proj inverted() const
        {
            if ( is_identity() ) {
                return *this;
            }
            return ec_point_fp_proj( this->curve(), x, -y, z );
        }

        /**
         * Adds the given point to this point.
         * R = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        [[nodiscard]] ec_point_fp_proj add(const ec_point_fp_proj& q) const
        {
            const auto& p = *this;
            if ( p.is_identity() ) {
                return q;
            }

            if ( q.is_identity() ) {
                return p;
            }

           const auto t0 = p.y * q.z;
           const auto t1 = q.y * p.z;
           const auto u0 = p.x * q.z;
           const auto u1 = q.x * p.z;
            if ( u0 == u1 ) {
                if ( t0 == t1 ) {
                    return doubled();
                } else {
                    return ec_point_fp_proj();
                }
            }

            // Note: Wrapping the following code in 3 lambdas
            //       can make slightly faster execution time (few 10s of us)
            const auto t  = t0 - t1;
            const auto u  = u0 - u1;
            const auto u2 = u.sqr();
            const auto u3 = u * u2;

            const auto v  = p.z * q.z;
            const auto w  = t.sqr() * v - u2 * ( u0 + u1 );

            auto rx = u * w;
            auto ry = t * ( u0 * u2 - w ) - t0 * u3;
            auto rz = u3 * v;
            return make_point( std::move(rx), std::move(ry), std::move(rz) );
        }

        /**
         * Returns the double of this point.
         * R = 2 * this
         *
         * @return the double of this point
        */
        [[nodiscard]] ec_point_fp_proj doubled() const
        {
            const auto& p = *this;
            if ( p.is_identity() || p.y == 0 ) {
                return ec_point_fp_proj(); // identity
            }

            const auto t = []( const ec_point_fp_proj& p) {
                const auto x2 = p.x.sqr();
                if ( p.curve().a_is_zero ) {
                    return 3 * x2;
                }
                if ( p.curve().a_is_minus_3 ) {
                    return 3 * ( x2 - p.z.sqr() );
                }
                return 3 * x2 + p.curve().a * p.z.sqr();
            }( p );

            const auto dy = 2 * p.y;
            const auto u  = dy * p.z;
            const auto v  = u * p.x * dy;
            const auto w  = t.sqr() - v * 2;

            auto rx = u * w;

            const auto u2 = u.sqr();
            auto ry = t * ( v - w ) - u2 * p.y.sqr() * 2;

            auto rz = u2 * u;
            return make_point( std::move(rx), std::move(ry), std::move(rz) );
        }

        /**
         * Subtracts the given point from this point.
         * R = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        [[nodiscard]] inline ec_point_fp_proj sub(const ec_point_fp_proj& a) const
        {
            return *this + (-a);
        }

        /**
         * Multiplies this point by the given scalar using the double and add algorithm.
         * R = this * scalar
         *
         * @param scalar - the scalar to multiply this point by
         * @return the resulting point
        */
        [[nodiscard]] ec_point_fp_proj mul(const int_type& scalar) const
        {
            if ( scalar.is_one() || is_identity() ) {
                return *this;
            }

            if ( scalar.is_zero() ) {
                return ec_point_fp_proj();
            }

            if ( scalar < 0 ) {
                return inverted() * -scalar;
            }

            auto r   = ec_point_fp_proj();
            auto tmp = *this;
            auto s   = scalar;
            while ( s != 0 ) {
                if ( s.is_odd() ) {
                    r += tmp;
                }
                tmp = tmp.doubled();
                s >>= 1;
            }
            return r;
        }

        /**
         * Compares 2 points for equality.
         * Points are equal if they are both identity or if
         * a.x * b.z == b.x * a.z and a.y * b.z == b.y * a.z
         *
         * @note Calling this function can be slow.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are equal, false otherwise
        */
        constexpr friend bool operator == (const ec_point_fp_proj& a, const ec_point_fp_proj& b)
        {
            if ( a.is_identity() || b.is_identity() ) {
                return a.is_identity() && b.is_identity();
            }
            return ( a.curve_  == b.curve_  ) &&
                   ( a.x * b.z == b.x * a.z ) &&
                   ( a.y * b.z == b.y * a.z );
        }

        /**
         * Compares 2 points for inequality.
         * @note See operator == for details.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are not equal, false otherwise
        */
        constexpr friend inline bool operator != (const ec_point_fp_proj& a, const ec_point_fp_proj& b)
        {
            return !( a == b );
        }

        private:
            friend CurveT;
            constexpr ec_point_fp_proj(const CurveT& curve, field_element_type x, field_element_type y, field_element_type z) :
                base_type( curve ),
                x( std::move(x) ),
                y( std::move(y) ),
                z( std::move(z) )
            {}

            ec_point_fp_proj make_point(field_element_type x, field_element_type y, field_element_type z) const
            {
                return ec_point_fp_proj( this->curve(), std::move(x), std::move(y), std::move(z) );
            }
    };

    /**
     * Struct representing a point on an elliptic curve in Jacobian coordinates
     * over a prime finite field GF(p) for a short Weierstrass-form elliptic curve satisfying the equation:
     *   y^2 = x^3 + ax * z^4 + b * z^6
     *
     * The use of Jacobian coordinates results in more efficient arithmetic operations
     * compared to projective or affine coordinates.
     *
     *
     * @warning The point's curve is stored as a pointer to the curve object.
     *          The curve object must outlive the point object.
     *
     * @tparam CurveT - the curve type. Required to be an instance of ec_curve_fp.
    */
    template<typename CurveT>
    struct ec_point_fp_jacobi : ec_point_base<ec_point_fp_jacobi<CurveT>, CurveT>
    {
        static_assert( is_ec_curve_fp<std::remove_cv_t<CurveT>> );

        using base_type          = ec_point_base<ec_point_fp_jacobi<CurveT>, CurveT>;
        using int_type           = typename CurveT::int_type;
        using field_element_type = typename CurveT::field_element_type;
        using affine_point_type  = ec_point_fp<CurveT>;
        using base_type::base_type;

        field_element_type x;
        field_element_type y;
        field_element_type z;

        /**
         * Constructs a point at infinity
        */
        constexpr ec_point_fp_jacobi() :
            base_type(),
            x( field_element_type::zero() ),
            y( field_element_type::one()  ),
            z( field_element_type::zero() )
        {}

        /**
         * Constructs this point from the given affine point.
         * @warning The point's curve is stored as a pointer to the curve object.
         *          The curve object must outlive the point object.
         *
         * @param p - the affine point to construct this point from.
        */
        explicit constexpr ec_point_fp_jacobi(affine_point_type p) :
            ec_point_fp_jacobi()
        {
            if ( !p.is_identity() ) {
                this->curve_ = &p.curve();
                x = std::move(p.x);
                y = std::move(p.y);
                z = field_element_type( 1, p.curve().p );
            }
        }

        /**
         * Normalizes this point.
         * This ensures Z coordinate is 1 therefore the x, y coordinates reflect
         * those of the equivalent to point in an affine coordinate system.
         *
         * @return this point
        */
        ec_point_fp_jacobi& normalize()
        {
            if ( !is_identity() && !z.is_one() ) {
                const auto z_inv  = z.inv();
                const auto z_inv2 = z_inv * z_inv;
                x *= z_inv2;
                y *= z_inv2 * z_inv;
                z = 1;
            }
            return *this;
        }

        /**
         * Returns the normalized representation of this point.
         * This ensures Z coordinate is 1 therefore the x, y coordinates reflect
         * those of the equivalent to point in an affine coordinate system.
        */
        [[nodiscard]] ec_point_fp_jacobi normalized() const
        {
            auto r = *this;
            r.normalize();
            return r;
        }

        /**
         * Returns the affine representation of this point.
         * @note No point verification is performed.
         * @note Slow operation due to division.
         *
         * @return this point in affine form.
        */
        [[nodiscard]] const affine_point_type to_affine() const
        {
            if ( is_identity() ) {
                return affine_point_type();
            }
            if ( z.is_one() ) {
                return affine_point_type( this->curve(), x, y );
            }

            // It is assumed that the point is on the curve and therefore calculated x and y are valid.
            // Calling affine_point_type() constructor will skip verification check.
            const auto z_inv  = z.inv();
            const auto z_inv2 = z_inv.sqr();
            return affine_point_type( this->curve(), x * z_inv2, y * z_inv2 * z_inv );
        }

        /**
         * Checks if this point is the identity element of the curve, i.e. point at infinity.
         * @return true if this point is the identity element of the curve, false otherwise
        */
        [[nodiscard]] inline constexpr bool is_identity() const
        {
            return this->curve_ == nullptr || ( z.is_zero() );
        }

        /**
         * Checks if this point is on the curve by calculating
         * the left and right hand side of the equation: y^2 = x^3 + ax * z^4 + b * z^6
         *
         * @return true if this point is on the curve, false otherwise
         * @note Slow operation
        */
        [[nodiscard]] bool is_on_curve() const
        {
            if ( is_identity() ) {
                return true;
            }
            const auto z2 = z.sqr();
            const auto z4 = z2.sqr();
            const auto z6 = z2 * z4;
            return  y.sqr() == ( ( ( x.sqr() + this->curve().a * z4 ) * x + this->curve().b * z6 ) ) ;
        }

        /**
         * Performs SEC 1 section 3.2.2.1 point check,
         * i.e. if point was generated using curve generator g and is not identity element.
         * Checks:
         *    - point is not the identity element
         *    - point is on the curve
         *    - point has order n
         *
         * @note Slow operation.
         *
         * @return true if this point is valid, false otherwise
        */
        [[nodiscard]] bool is_valid() const
        {
            if ( is_identity() ) {
                return false;
            }

            if ( !is_on_curve() ) {
                return false;
            }

            if ( !(this->curve().n * *this).is_identity() ) {
                return false;
            }

            return true;
        }

        /**
         * Returns the inverse of this point.
         * R = -this
         *
         * @return the inverse of this point
        */
        [[nodiscard]] constexpr ec_point_fp_jacobi inverted() const
        {
            if ( is_identity() ) {
                return *this;
            }
            return ec_point_fp_jacobi( this->curve(), x, -y, z );
        }

        /**
         * Adds the given point to this point.
         * R = this + a
         *
         * @param a - the point to add to this point
         * @return reference to this point
        */
        [[nodiscard]] ec_point_fp_jacobi add(const ec_point_fp_jacobi& q) const
        {
            const auto& p = *this;
            if ( p.is_identity() ) {
                return q;
            }

            if ( q.is_identity() ) {
                return p;
            }

            // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-1998-cmo-2
            // note: faster than https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl

            const bool z1_is_one = p.z.is_one();
            const bool z2_is_one = q.z.is_one();
            if ( z1_is_one && z2_is_one ) {
                return add_z_1( p, q );
            }
            else if ( z1_is_one ) {
                return add_z2_1( q, p );
            }
            else if ( z2_is_one ) {
                return add_z2_1( p, q );
            }
            return add_ne( p, q );
        }

        /**
         * Returns the double of this point.
         * R = 2 * this
         *
         * @return the double of this point
        */
        [[nodiscard]] ec_point_fp_jacobi doubled() const
        {
            const auto& p = *this;
            if ( p.is_identity() ) {
                return p;
            }

            if ( p.y == 0 ) {
                return ec_point_fp_jacobi(); // identity
            }

            // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-1986-cc
            // note: this algo was measured to be the most efficient of them all.

            const auto M = [](const auto& p) {
                const bool z_is_one = p.z.is_one();
                if ( p.curve().a_is_zero ) {
                    return 3 * p.x.sqr();
                }
                else if ( p.curve().a_is_minus_3 ) {
                    const auto z2 = z_is_one ? p.z : p.z.sqr();
                    return 3 * ( p.x - z2 ) * ( p.x + z2 );
                }
                else {
                    const auto z4 = z_is_one ? p.z : p.z.sqr().sqr();
                    return 3 * p.x.sqr() + p.curve().a * z4;
                }
            }( p );

            const auto y2 = p.y.sqr();
            const auto S  = 4 * p.x * y2;
            auto X3 = M.sqr() - 2 * S;
            auto Y3 = M * ( S - X3 ) - 8 * y2.sqr();
            auto Z3 = 2 * p.y * p.z;
            return make_point( std::move(X3), std::move(Y3), std::move(Z3) );
        }

        /**
         * Subtracts the given point from this point.
         * R = this - a
         *
         * @param a - the point to subtract from this point
         * @return reference to this point
        */
        [[nodiscard]] inline ec_point_fp_jacobi sub(const ec_point_fp_jacobi& a) const
        {
            return *this + ( -a );
        }

        /**
         * Multiplies this point by the given scalar.
         * R = this * scalar
         *
         * @param scalar - the scalar to multiply this point by
         * @return the resulting point
        */
        [[nodiscard]] ec_point_fp_jacobi mul(const int_type& scalar) const
        {
            if ( scalar.is_one() || is_identity() ) {
                return *this;
            }

            if ( scalar.is_zero() ) {
                return ec_point_fp_jacobi();
            }

            if ( scalar < 0 ) {
                return inverted() * -scalar;
            }

            // Point scalar multiplication by the reversed NAF method
            // optimized for cases where this.z == 1.
            //
            // The base algorithm is described in:
            // "Comparative analysis of the scalar point multiplication
            // algorithms in the NIST FIPS 186 elliptic curve cryptography"
            // (Section 3.2: Scalar multiplication by the NAF method)
            // https://ceur-ws.org/Vol-2913/paper2.pdf
            const auto& p   = *this;
            const auto pneg = -p;
            auto r = make_point(0, 0, 1);
            for ( const char i : scalar.to_rnaf() ) {
                r = r.doubled();
                if ( i > 0 ) {
                    r += p;
                }
                else if ( i < 0 ) {
                    r += pneg;
                }
            }
            return r;
        }

        /**
         * Compares 2 points for equality.
         * Points are equal if they are both identity or if
         * a.x * b.z^2 == b.x * a.z^2 and a.y * b.z^3 == b.y * a.z^3
         *
         * @note Calling this function can be slow.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are equal, false otherwise
        */
        constexpr friend bool operator == (const ec_point_fp_jacobi& a, const ec_point_fp_jacobi& b)
        {
            if ( std::addressof(a) == std::addressof(b) ) {
                return true;
            }

            if ( a.is_identity() || b.is_identity() ) {
                return a.is_identity() && b.is_identity();
            }

            // Using short circuit in case of inequality
            const auto az2 = a.z.sqr();
            const auto bz2 = b.z.sqr();
            return ( a.curve_  == b.curve_  ) &&
                   ( a.x * bz2 == b.x * az2 ) &&
                   ( a.y * bz2 * b.z == b.y * az2 * a.z );
        }

        /**
         * Compares 2 points for inequality.
         * @note See operator == for details.
         *
         * @param a - the first point to compare
         * @param b - the second point to compare
         * @return true if the points are not equal, false otherwise
        */
        constexpr friend inline bool operator != (const ec_point_fp_jacobi& a, const ec_point_fp_jacobi& b)
        {
            return !( a == b );
        }

        private:
            friend CurveT;
            constexpr ec_point_fp_jacobi(const CurveT& curve, field_element_type x, field_element_type y, field_element_type z) :
                base_type( curve ),
                x( std::move(x) ),
                y( std::move(y) ),
                z( std::move(z) )
            {}

            ec_point_fp_jacobi make_point(field_element_type x, field_element_type y, field_element_type z) const
            {
                return ec_point_fp_jacobi( this->curve(), std::move(x), std::move(y), std::move(z) );
            }

            ec_point_fp_jacobi make_point(int_type x, int_type y, int_type z) const
            {
                return ec_point_fp_jacobi(
                    this->curve(),
                    field_element_type( std::move(x), this->curve().p ),
                    field_element_type( std::move(y), this->curve().p ),
                    field_element_type( std::move(z), this->curve().p )
                );
            }

            [[nodiscard]]
            __attribute__((always_inline)) // note: forced inline produces slightly more efficient computation. [[clang::always_inline]] doesn't work.
            static ec_point_fp_jacobi addex(const ec_point_fp_jacobi& p, const ec_point_fp_jacobi& q,
                                            const field_element_type& U1, const field_element_type& U2,
                                            const field_element_type& S1, const field_element_type& S2)
            {
                // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-1998-cmo-2
                //  note: faster than https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl
                const auto H  = U2 - U1;
                const auto R  = S2 - S1;
                if ( H.is_zero() ) {
                    if ( R.is_zero() ) {
                        return p.doubled();
                    } else {
                        return ec_point_fp_jacobi();
                    }
                }
                const auto H2 = H.sqr();
                const auto H3 = H2 * H;
                const auto V  = U1 * H2;

                auto X3 = R.sqr() - H3 - 2 * V;
                auto Y3 = R * ( V - X3 ) - S1 * H3;
                auto Z3 = std::move( H );
                if ( !p.z.is_one() ) {
                    Z3 *= p.z;
                }
                if ( !q.z.is_one() ) {
                    Z3 *= q.z;
                }
                return ec_point_fp_jacobi( p.curve(), std::move(X3), std::move(Y3), std::move(Z3) );
            }

            [[nodiscard]]
            __attribute__((always_inline))
            static ec_point_fp_jacobi add_z_1(const ec_point_fp_jacobi& p, const ec_point_fp_jacobi& q)
            {
                return addex( p, q, p.x, q.x, p.y, q.y );
            }

            [[nodiscard]]
            __attribute__((always_inline))
            static ec_point_fp_jacobi add_z2_1(const ec_point_fp_jacobi& p, const ec_point_fp_jacobi& q)
            {
                const auto pz2 = p.z.sqr();
                return addex( p, q, p.x, q.x * pz2, p.y, q.y * pz2 * p.z );
            }

            [[nodiscard]]
            __attribute__((always_inline))
            static ec_point_fp_jacobi add_ne(const ec_point_fp_jacobi& p, const ec_point_fp_jacobi& q)
            {
                // This extra function, although inlined, produces slightly more efficient code than
                // it would if put directly into the calling scope.
                const auto pz2 = p.z.sqr();
                const auto qz2 = q.z.sqr();
                return addex( p, q, p.x * qz2, q.x * pz2, p.y * qz2 * q.z, q.y * pz2 * p.z );
            }
    };

    /**
     * Base struct for representing elliptic curve.
    */
    template<typename CurveT, typename FieldElementT, typename PointT>
    struct ec_curve_base {
        using int_type           = typename FieldElementT::int_type;
        using field_element_type = FieldElementT;
        using point_type         = PointT;

        /**
         * Generates a point from base point g and given scalar x.
         *
         * @tparam PointU - point type to create.
         *
         * @param x - scalar to multiply base point g with.
         * @return PointU - point on the curve.
        */
        template<typename PointU = point_type, typename IntT = int_type>
        [[nodiscard]] inline PointU generate_point(const IntT &x) const
        {
            return static_cast<const CurveT&>( *this )
                .template generate_point<PointU>( x );
        }

        /**
         * Generates a point from base point g and given hex string scalar sx.
         *
         * @tparam PointU - point type to create.
         *
         * @param sx - hex string scalar to multiply base point g with.
         * @return PointU - point on the curve.
        */
        template<typename PointU = point_type>
        [[nodiscard]] inline PointU generate_point(const std::string_view sx) const
        {
            return static_cast<const CurveT&>( *this )
                .template generate_point<PointU>( sx );
        }

        /**
         * Generates a point from base point g and given hex string literal scalar sx.
         *
         * @tparam PointU - point type to create.
         *
         * @param sx - hex string literal scalar to multiply base point g with.
         * @return PointU - point on the curve.
        */
        template<typename PointU = point_type, std::size_t N>
        [[nodiscard]] inline PointU generate_point(const char (&sx)[N]) const
        {
            return static_cast<const CurveT&>( *this )
                .template generate_point<PointU>( sx );
        }

        /**
         * Returns curve field element from the given integer.
         * @param x the integer to convert
         * @return curve field element
        */
        template<typename IntT = int_type>
        [[nodiscard]] inline constexpr field_element_type make_field_element(IntT&& x) const
        {
            return static_cast<const CurveT&>( *this )
                .make_field_element( std::move(x) );
        }

        /**
         * Returns curve field element from the given integer.
         * @param x the integer to convert
         * @return curve field element
        */
        template<typename IntT = int_type>
        [[nodiscard]] inline constexpr field_element_type make_field_element(const IntT& x) const
        {
            return static_cast<const CurveT&>( *this )
                .make_field_element( x );
        }

        /**
         * Returns curve point from the given big numbers.
         * @param x the x coordinate of the point
         * @param y the y coordinate of the point
         * @param verify if true, the point is verified to be valid point on the curve created by the curve generator point g.
         * @return curve point
        */
        template<typename PointU = point_type, typename IntT = int_type>
        [[nodiscard]] inline constexpr PointU make_point(IntT x, IntT y, bool verify = false) const
        {
            return static_cast<const CurveT&>( *this )
                .template make_point<PointU>( std::move(x), std::move(y), verify );
        }

        /**
         * Verifies curve parameters.
         * @return true if curve parameters are valid.
        */
        [[nodiscard]] inline bool verify() const
        {
            return static_cast<const CurveT&>( *this ).verify();
        }
    };

    /**
     * Struct defines curve over a prime finite field GF(p)
     * with Weierstrass equation y^2 = x^3 + ax + b.
     *
     * @tparam IntT  - integer type
     * @tparam CurveTag - the curve tag
    */
    template<typename IntT, typename CurveTag>
    struct ec_curve_fp :
        ec_curve_base<
            ec_curve_fp<IntT, CurveTag>,
            fp_element<IntT, CurveTag>,
            ec_point_fp<ec_curve_fp<IntT, CurveTag>>
        >
    {
        using base_type = ec_curve_base<ec_curve_fp<IntT, CurveTag>,
                            fp_element<IntT, CurveTag>,
                            ec_point_fp<ec_curve_fp<IntT, CurveTag>>>;

        using curve_tag          = CurveTag;
        using int_type           = IntT;
        using field_element_type = typename base_type::field_element_type;
        using point_type         = typename base_type::point_type;

        const IntT       p;      // curve prime
        const IntT       a;      // curve coefficient
        const IntT       b;      // curve coefficient
        const point_type g;      // generator point
        const IntT       n;      // order of g
        const uint32_t   h;      // cofactor, i.e.: h = #E(Fp) / n
                                 //     #E(Fp) - number of points on the curve
        const bool a_is_minus_3; // cached constant a == p - 3
        const bool a_is_zero;    // cached constant a == 0
        const IntT p_minus_n;    // cached constant p - n; used for checking the maximum negative point coordinate

        /**
         * Creates a curve from the given parameters.
         * @param p - curve prime
         * @param a - curve coefficient
         * @param b - curve coefficient
         * @param g - generator point coordinates
         * @param n - order of g
         * @param h - cofactor, i.e.: h = p / n
        */
        inline constexpr ec_curve_fp(IntT p, IntT a, IntT b, std::pair<IntT, IntT> g, IntT n, uint32_t h):
            p( std::move(p) ),
            a( std::move(a) ),
            b( std::move(b) ),
            g( make_point( std::move(g.first), std::move(g.second) )),
            n( std::move(n) ),
            h( h ),
            a_is_minus_3( this->a == ( this->p - 3) ),
            a_is_zero( this->a.is_zero() ),
            p_minus_n( this->p - this->n )
        {}

        /**
         * Generates a point from base point g and given scalar x.
         *
         * @tparam PointT - EC point type to create.
         *
         * @param x - scalar to multiply base point g with.
         * @return PointT - point on the curve.
        */
        template<typename PointT = point_type, typename IntNumT = IntT>
        [[nodiscard]] PointT generate_point(const IntNumT& x) const
        {
            check( x > 0 && x < n, "x must be in range [1, n-1]" );
            if constexpr (std::is_same_v<PointT, point_type>) {
                return g * x;
            } else {
                return PointT( g ) * x;
            }
        }

        /**
         * Generates a point from base point g and given hex string scalar sx.
         *
         * @tparam PointU - EC point type to create.
         *
         * @param sx - hex string scalar to multiply base point g with.
         * @return PointU - point on the curve.
        */
        template<typename PointU = point_type>
        [[nodiscard]] PointU generate_point(const std::string_view sx) const
        {
            return generate_point<PointU>( int_type( sx ) );
        }

        /**
         * Generates a point from base point g and given hex string literal scalar sx.
         *
         * @tparam PointU - EC point type to create.
         *
         * @param sx - hex string literal scalar to multiply base point g with.
         * @return PointU - point on the curve.
        */
        template<typename PointU = point_type, std::size_t N>
        [[nodiscard]] inline PointU generate_point(const char (&sx)[N]) const
        {
            return generate_point<PointU>( std::string_view( sx, N ) );
        }

        /**
         * Creates a field element from a given integer.
         * @warning Returned field element stores pointer to curve prime.
         *          The curve must outlive the field element.
         *
         * @note Expects x > 0 and x < p
         *
         * @param x - Integer to convert to field element
         * @return Field element
         *
        */
        [[nodiscard]] constexpr field_element_type make_field_element(IntT&& x) const
        {
            return make_field_element( std::move(x), /*verify=*/ true );
        }

        /**
         * Creates a field element from a given integer.
         * @warning Returned field element stores pointer to curve prime.
         *          The curve must outlive the field element.
         *
         * @note Expects x > 0 and x < p
         *
         * @param x - Integer to convert to field element
         * @return Field element
         *
        */
        [[nodiscard]] inline constexpr field_element_type make_field_element(const IntT& x) const
        {
            return make_field_element( x, /*verify=*/ true );
        }

        /**
         * Creates a point from provided affine coordinates x & y.
         * @warning Returned point stores pointer to curve prime.
         *          The curve must outlive the point.
         * @note Expects x|y >= 0 and x|y < p
         * @note The returned point can be invalid, since the function allows
         *        creating points that were not generated with the generator point.
         *
         * @tparam PointU - EC point type to make. Default affine point_type.
         *
         * @param x - point x coordinate
         * @param y - point y coordinate
         * @param verify - If true, the point is verified to be valid point on the curve created by the curve generator point g.
         *                 Default is false. Slow operation, can be be performed also with call to point.is_valid() function.
         * @return Curve point
        */
        template<typename PointU = point_type>
        [[nodiscard]] constexpr PointU make_point(IntT&& x, IntT&& y, bool verify = false) const
        {
            check_integer( x, "Invalid point x coordinate" );
            check_integer( y, "Invalid point y coordinate" );
            auto point = make_point<PointU>(
                make_field_element( std::move(x) ),
                make_field_element( std::move(y) ),
                verify
            );
            return point;
        }

        /**
         * Creates a point from provided affine coordinates x & y.
         * @warning Returned point stores pointer to curve prime.
         *          The curve must outlive the point.
         * @note Expects x|y >= 0 and x|y < p
         * @note The returned point can be invalid, since the function allows
         *        creating points that were not generated with the generator point.
         *
         * @tparam PointU - EC point type to make. Default affine point_type.
         *
         * @param x - point x coordinate
         * @param y - point y coordinate
         * @param verify - If true, the point is verified to be valid point on the curve created by the curve generator point g.
         *                 Default is false. Slow operation, can be be performed also with call to point.is_valid() function.
         * @return Curve point
        */
        template<typename PointU = point_type>
        [[nodiscard]] constexpr PointU make_point(const IntT& x, const IntT& y, bool verify = false) const
        {
            check_integer( x, "Invalid point x coordinate" );
            check_integer( y, "Invalid point y coordinate" );
            auto point = make_point<PointU>(
                make_field_element( x, /*verify=*/ false ),
                make_field_element( y, /*verify=*/ false ),
                verify
            );
            return point;
        }

        /**
         * Solves the equation y = sqrt( x^3 + ax + b ) (mod p) and returns y.
         *
         * @param x   - field element x-coordinate.
         * @param odd - True if the calculated y-coordinate should be odd.
         * @return The y-coordinate, or 0 if no square root exists.
        */
        [[nodiscard]]
        inline field_element_type compute_y(const field_element_type& x, const bool odd) const {
            auto y = (( x.sqr() + a ) * x + b ).sqrt();
            if ( odd != y.value().test_bit( 0 )) {
                y = -y;
            }
            return y;
        }

        /**
         * Solves the equation y = sqrt( x^3 + ax + b ) (mod p) and returns y.
         *
         * @param x      - integer x-coordinate.
         * @param odd    - True if the calculated y-coordinate should be odd.
         * @param verify - True if x coordinate should be verified before calculating y. Default false.
         * @return The y-coordinate, or 0 if no square root exists.
        */
        [[nodiscard]]
        inline field_element_type compute_y(const IntT& x, const bool odd, bool verify = false) const {
            return compute_y( make_field_element( x, verify ), odd );
        }

        /**
         * Decompress point from provided x-coordinate.
         * @tparam PointT - returned EC point type.
         *
         * @param x    - field element x-coordinate.
         * @param yodd - True if the calculated y-coordinate should be odd.
         * @return The decompressed EC point. If the y-coordinate cannot be computed, the point at infinity is returned.
        */
        template<typename PointT = point_type>
        [[nodiscard]]
        inline PointT decompress_point(field_element_type x, const bool yodd) const {
            auto y = compute_y( x, yodd );
            if ( y.is_zero() ) {
                return PointT{};
            }
            return make_point<PointT>( std::move( x ), std::move( y ) );
        }

        /**
         * Decompress point from provided x-coordinate.
         * @tparam PointT - returned EC point type.
         *
         * @param x       - integer x-coordinate.
         * @param yodd    - True if the calculated y-coordinate should be odd.
         * @param verify  - True if x coordinate should be verified before calculating y. Default false.
         * @return The decompressed EC point. If the y-coordinate cannot be computed, the point at infinity is returned.
        */
        template<typename PointT = point_type>
        [[nodiscard]]
        inline PointT decompress_point(const IntT& x, const bool yodd, bool verify = false) const {
            return decompress_point<PointT>( make_field_element( x, verify ), yodd );
        }

        /**
         * Encodes an EC point into its byte-encoded form.
         *
         * This function follows the SEC1-v2 standard for EC point encoding.
         * The 'compress' parameter determines whether to use compressed or uncompressed
         * point encoding. Compressed encoding includes only the x-coordinate, while uncompressed encoding
         * includes both the x and y coordinates.
         *
         * Encoded Format: <type> | [x] | [y]
         *   type:
         *      0 - point at infinity, not coordinates is encoded
         *      2 - compressed form  with only x-coordinate (even y-coordinate is derived)
         *      3 - compressed form with only x-coordinate (odd y-coordinate is derived)
         *      4 - uncompressed form with both x and y coordinates
         *
         * @see SEC1-v2 section '2.3.3 Elliptic-Curve-Point-to-Octet-String Conversion'
         * https://www.secg.org/sec1-v2.pdf#page=16
         *
         * @param point    - The EC point to be encoded, represented as a `point_type` structure.
         * @param compress - If true the encoded point will be in compressed form; otherwise, encodes it in uncompressed form.
         * @return The byte-encoded form of the EC point.
         */
        [[nodiscard]] bytes encode_point(const point_type& point, bool compress) const
        {
            bytes epoint = { 0 };
            if ( point.is_identity() ) {
                return epoint;
            }

            if ( compress ) {
                epoint[0] = 2 + point.y.value().test_bit( 0 );
                epoint += point.x.to_bytes();
            }
            else {
                epoint[0] = 4;
                epoint += point.x.to_bytes() + point.y.to_bytes();
            }
            return epoint;
        }

        /**
         * Decodes an EC point from its byte-encoded form.
         *
         * This function follows the SEC1-v2 standard for EC point encoding. The encoded point
         * is expected to be in the format (type | x [| y]), where 'type' represents the encoding type or
         * parity of the y-coordinate (if compressed form), and 'x' & 'y' are coordinates of the point.
         *
         * @see SEC1-v2 section '2.3.4 Octet-String-to-Elliptic-Curve-Point Conversion'
         * https://www.secg.org/sec1-v2.pdf#page=17
         *
         * @see encode_point
         *
         * @note internally function halts if epoint is in invalid format or encoded coordinates are invalid.
         *
         * @tparam PointT - returned EC point type.
         *
         * @param epoint A bytes view of the encoded EC point (type | x [| y]).
         * @return The decoded EC point of type PointT.
         */
        template<typename PointT = point_type>
        [[nodiscard]] PointT decode_point(const bytes_view epoint) const
        {
            check( !epoint.empty(), "invalid encoded point" );

            const auto type = epoint[0];
            check( type <= 4 && type != 1, "invalid encoded point" );

            const std::size_t plen = p.byte_length();
            PointT P; // type 0 = point at infinity
            switch ( type ) {
                // 2.3.4.2
                case 2:
                case 3: {
                    check( epoint.size() == plen + 1, "invalid encoded point" );
                    P = decompress_point<PointT>( /*x=*/epoint.subspan( 1, plen ), /*yodd=*/type & 1 );
                } break;
                // 2.3.4.3
                case 4: {
                    check( epoint.size() == 2 * plen + 1, "invalid encoded point" );
                    P = make_point<PointT>(
                        epoint.subspan( 1, plen ),
                        epoint.subspan( 1 + plen, plen )
                    );
                } break;
            }
            return P;
        }

        /**
         * Verifies curve parameters.
         * @note Very basic verification is performed, no check for primality of p and n;
         *       Very slow operation.
         *
         * @return true if curve parameters are valid.
        */
        [[nodiscard]] bool verify() const
        {
            // Verify that p is prime
            // if ( !p.is_prime() ) {
            //     return false;
            // }

            // Verify p and n
            if ( p <= 3 || n <= 0 ) {
                return false;
            }

            // Verify that a and b are in the range [0, p-1]
            if ( a.is_negative() || a >= p ) {
                return false;
            }
            if ( b.is_negative() || b >= p ) {
                return false;
            }

            // Verify that a^3 + 27 * b^2 != 0
            auto afe = make_field_element( a, /*verify =*/ false );
            auto bfe = make_field_element( b, /*verify =*/ false );
            auto y2  = 4 * afe.sqr() * afe + 27 * bfe.sqr();

            // check that discriminant is nonzero. If zero, the curve is singular.
            if ( ( -16 * y2 ) == 0) {
                return false;
            }

            // Verify G point
            if ( !g.is_valid() ) {
                return false;
            }

            // Verify that n is prime
            // if ( !n.is_prime() ) {
            //     return false;
            // }

            // Verify that h is in the range [1, p-1]
            if ( h < 1 || h >= p ) {
                return false;
            }

            // Verify that n * h = p - 1
            if ( n * h == p - 1 ) {
                return false;
            }

            return true;
        }

        private:
            [[nodiscard]] inline constexpr bool is_valid_integer(const IntT& x) const
            {
                return !x.is_negative() && x < p;
            }

            inline constexpr void check_integer(const IntT& x, const char* error) const
            {
                check( is_valid_integer( x ), error );
            }

            [[nodiscard]] constexpr field_element_type make_field_element(IntT&& x, bool verify) const
            {
                if ( verify ) {
                    check_integer( x, "Invalid field element value" );
                }
                return field_element_type( std::move(x), p );
            }

            [[nodiscard]] constexpr field_element_type make_field_element(const IntT& x, bool verify) const
            {
                if ( verify ) {
                    check_integer( x, "Invalid field element value" );
                }
                return field_element_type( x, p );
            }

            template<typename PointU>
            [[nodiscard]] constexpr PointU make_point(field_element_type x, field_element_type y, bool verify = false) const
            {
                auto point = PointU{ point_type {
                    *this,
                    std::move( x ),
                    std::move( y )
                }};

                if ( verify ) {
                    check( point.is_valid(), "Invalid point" );
                }
                return point;
            }
    };

    /**
     * Fast multiplication of points and addition of points, i.e.: a*P + b*Q
     * Function uses combination of non-adjacent form of number (NAF) and
     * Shamir's trick to calculate a*P + b*Q in one batch, thus the speedup shoule be ~2.5x.
     *
     * @tparam PointT - point type
     * @tparam CurveT - curve type
     *
     * @param a - First multiplier
     * @param p - First point
     * @param b - Second multiplier
     * @param q - Second point
     * @return PointT - Result of a*P + b*Q
    */
    template <typename PointT, typename CurveT, typename IntT = typename CurveT::int_type>
    [[nodiscard]] static PointT ec_mul_add_fast(const IntT &a, const ec_point_base<PointT, CurveT> &p,
                                                const IntT &b, const ec_point_base<PointT, CurveT> &q)
    {
        // Get the NAF representations of a and b.
        auto a_naf = a.to_rnaf();
        auto b_naf = b.to_rnaf();
        if ( a_naf.size() < b_naf.size() ) {
            a_naf.insert( a_naf.begin(), b_naf.size() - a_naf.size(), 0 );
        } else if ( a_naf.size() > b_naf.size() ) {
            b_naf.insert( b_naf.begin(), a_naf.size() - b_naf.size(), 0 );
        }

        // Precompute required points
        const auto pneg         = -p;
        const auto qneg         = -q;
        const auto pnegqneg     = pneg + qneg;
        const auto pqneg        = p + qneg;
        const auto pqneg_inv    = -pqneg;
        const auto pnegqneg_inv = -pnegqneg;

        // Iterate reversed NAF representations of a and b,
        // optimized for cases where p.z == 1 or q.z == 1.
        PointT r;
        for ( std::size_t i = 0; i < a_naf.size(); i++ ) {
            r = r.doubled();
            if ( a_naf[i] == 0 ) {
                if ( b_naf[i] == 0 ) {
                    continue;
                }
                else if ( b_naf[i] < 0 ) {
                    r += qneg;
                }
                else {
                    r += q;
                }
            }
            else if ( a_naf[i] < 0 ) {
                if ( b_naf[i] == 0 ) {
                    r += pneg;
                }
                else if ( b_naf[i] < 0 ) {
                    r += pnegqneg;
                }
                else {
                    r += pqneg_inv;
                }
            }
            else { // a_naf[i] > 0
                if ( b_naf[i] == 0 ) {
                    r += p;
                }
                else if ( b_naf[i] < 0 ) {
                    r += pqneg;
                }
                else {
                    r += pnegqneg_inv;
                }
            }
        }

        return r;
    }
}