completing x norm reconstruction

libs/math/docs/polynomial_recovery.md

@@ -2,13 +2,14 @@
 
 @tableofcontents
 
-Crypto3.Math provides three complementary recovery facilities:
+Crypto3.Math provides four complementary recovery facilities:
 
 * square testing and square roots in a finite polynomial quotient field;
-* bounded rational reconstruction from a residue modulo a polynomial; and
-* recovery of the fixed `X`-norm representation of an irreducible polynomial.
+* bounded rational reconstruction from a residue modulo a polynomial;
+* recovery of the fixed `X`-norm representation of an irreducible polynomial; and
+* factorization-aware recovery of the same representation for a general polynomial.
 
-Together they can recover representations of irreducible factors by the polynomial norm form
+Together they can recover representations by the polynomial norm form
 
     P(X)^2 - X * Q(X)^2.
 
@@ -209,10 +210,59 @@ This is the local representation of `g` by the norm from adjoining a square root
 
 All four polynomial products use the supplied arithmetic context, while multiplication by `X` is a coefficient shift.
 
-Consequently, a caller can factor a target polynomial, recover eligible irreducible factors independently, account for
-multiplicities and the scalar leading coefficient, and combine the local representations. Crypto3.Math deliberately
-keeps that application-level policy separate from the generic factorization, quotient-field square-root, bounded
-rational-reconstruction, and `X`-norm reconstruction facilities.
+### General polynomial recovery
+
+`recover_polynomial_x_norm_representation` applies the same norm construction to a canonical polynomial `H` without
+requiring the caller to decompose it first:
+
+```cpp
+// arithmetic_context and coefficient_generator are caller-owned.
+auto representation = math::recover_polynomial_x_norm_representation<backend_type>(
+    H, arithmetic_context, coefficient_generator);
+
+if (representation) {
+    polynomial_type exact_norm = math::evaluate_polynomial_x_norm<backend_type>(
+        *representation, arithmetic_context);
+    // exact_norm == H
+}
+```
+
+Zero is represented by `(0, 0)`. A constant `c` is represented by `(sqrt(c), 0)` when `c` is square and otherwise has
+no result. These cases do not invoke factorization or consume the coefficient generator.
+
+For a nonconstant input, the constant coefficient must be square. Its leading coefficient must be square when the
+degree is even, while the negated leading coefficient must be square when the degree is odd. These necessary tests
+reject impossible inputs before factorization, but passing them does not guarantee recovery.
+
+Complete factorization writes the remaining input as
+
+    H = c * product(g^e),
+
+where `c` is the leading coefficient and each `g` is monic and irreducible. An even factor power `g^(2r)` has the
+immediate representation `(g^r, 0)`. For an odd power `g^(2r+1)`, irreducible recovery first obtains `(P_g, Q_g)` for
+the unpaired copy of `g`; scaling both components by `g^r` then represents the complete factor power. Factor recovery
+stops as soon as an odd-multiplicity factor cannot be represented.
+
+The factor-power representations are combined with `multiply_polynomial_x_norm_representations` in a balanced product
+tree. Finally, both components are scaled by `sqrt(c)` to incorporate the factorization's leading coefficient. A
+nonsquare required scalar produces no result. Every successful path evaluates the completed norm and compares it
+exactly with `H` before returning.
+
+The coefficient generator is shared by complete factorization and odd-factor recovery. It must meet the factorization
+generator contract and must allow each square-root context to find a quotient-field nonsquare. All polynomial products
+and squares use the caller's compile-time-selected arithmetic context; no polynomial backend or random engine is
+constructed internally.
+
+| Outcome | Contract |
+|---|---|
+| Representation returned | `evaluate_polynomial_x_norm(result, context) == H` exactly. |
+| No value returned | A special-case or necessary square test fails, an odd factor is unrepresentable, bounded recovery fails, or scalar normalization is impossible. |
+| `std::invalid_argument` | `H` is empty or noncanonical, or a composed API contract is violated. |
+| `std::logic_error` | Completed factorization, recovery, combination, normalization, or final verification is internally inconsistent. |
+
+This factorization and factor-combination behavior is a generic polynomial operation: it depends only on finite-field
+polynomial arithmetic and the fixed norm map `P^2 - X * Q^2`. The API returns the two coefficient polynomials and does
+not impose a representation or policy beyond that algebraic identity.
 
 ### Worked example over F7
 

libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp

@@ -31,11 +31,13 @@
 #include <optional>
 #include <stdexcept>
 #include <utility>
+#include <vector>
 
 #include <nil/crypto3/algebra/fields/field_algorithms.hpp>
 
 #include <nil/crypto3/math/polynomial/operations/basic_operations.hpp>
 #include <nil/crypto3/math/polynomial/operations/shift.hpp>
+#include <nil/crypto3/math/polynomial/factorization/complete_factorization.hpp>
 #include <nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp>
 #include <nil/crypto3/math/polynomial/reconstruction/polynomial_rational_reconstruction.hpp>
 
@@ -244,6 +246,244 @@ namespace nil::crypto3::math {
         return std::move(representation);
     }
 
+    namespace detail {
+
+        /** Raise a canonical polynomial to a nonnegative integer power using the supplied arithmetic context. */
+        template<polynomial_arithmetic::PolynomialBackend Backend>
+        typename Backend::polynomial_type
+            polynomial_x_norm_power(const typename Backend::polynomial_type &base, std::size_t exponent,
+                                    polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context) {
+            using polynomial_type = typename Backend::polynomial_type;
+            using value_type = typename polynomial_type::value_type;
+
+            polynomial_type result = {value_type::one()};
+            polynomial_type current_power(base);
+            while (exponent != 0) {
+                if ((exponent & 1) != 0) {
+                    polynomial_type product;
+                    arithmetic_context.multiply(product, result, current_power);
+                    result = std::move(product);
+                }
+                exponent >>= 1;
+                if (exponent != 0) {
+                    polynomial_type square;
+                    arithmetic_context.square(square, current_power);
+                    current_power = std::move(square);
+                }
+            }
+            return result;
+        }
+
+        /**
+         * Construct an X-norm representation of one irreducible factor raised to its multiplicity. Complete
+         * factorization expresses the input as a leading scalar times a product of powers g^e, so each such power
+         * needs a representation before the factor representations can be combined.
+         *
+         * If e = 2r, then g^e is already a square. The pair (g^r, 0) represents it because
+         *
+         *     (g^r)^2 - X * 0^2 = g^(2r).
+         *
+         * If e = 2r + 1, recover (P_g, Q_g) for the one unpaired copy of g, where
+         *
+         *     P_g^2 - X * Q_g^2 = g.
+         *
+         * Scaling both components by g^r gives a representation of the complete factor power:
+         *
+         *     (g^r * P_g)^2 - X * (g^r * Q_g)^2 = g^(2r) * g = g^e.
+         *
+         * For example, if H contains g^3 * h^2 and (P_g, Q_g) represents g, then (g * P_g, g * Q_g)
+         * represents g^3, while (h, 0) represents h^2. Combining those two representations produces the
+         * representation of g^3 * h^2. Thus only the odd-multiplicity factor g requires irreducible recovery.
+         */
+        template<SupportsDivrem Backend, typename Generator>
+        std::optional<polynomial_x_norm_representation<typename Backend::polynomial_type>>
+            recover_polynomial_x_norm_factor_power(
+                const polynomial_factor<typename Backend::polynomial_type> &factor,
+                polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context,
+                Generator &coefficient_generator) {
+            using polynomial_type = typename Backend::polynomial_type;
+            using value_type = typename polynomial_type::value_type;
+            using representation_type = polynomial_x_norm_representation<polynomial_type>;
+
+            if (factor.multiplicity == 0) {
+                throw std::logic_error("complete factorization produced a factor with zero multiplicity");
+            }
+
+            const std::size_t half_multiplicity = factor.multiplicity / 2;
+            polynomial_type half_power =
+                polynomial_x_norm_power<Backend>(factor.polynomial, half_multiplicity, arithmetic_context);
+            if ((factor.multiplicity & 1) == 0) {
+                return representation_type {std::move(half_power), polynomial_type {value_type::zero()}};
+            }
+
+            auto odd_representation = recover_irreducible_polynomial_x_norm_representation<Backend>(
+                factor.polynomial, arithmetic_context, coefficient_generator);
+            if (!odd_representation) {
+                return std::nullopt;
+            }
+            if (half_multiplicity == 0) {
+                return odd_representation;
+            }
+
+            polynomial_type lifted_p;
+            polynomial_type lifted_q;
+            arithmetic_context.multiply(lifted_p, half_power, odd_representation->p);
+            arithmetic_context.multiply(lifted_q, half_power, odd_representation->q);
+            return representation_type {std::move(lifted_p), std::move(lifted_q)};
+        }
+
+        /**
+         * Combine X-norm representations in balanced levels so no left-deep product chain is formed. For example,
+         * five factor representations are combined as
+         *
+         *     [A, B, C, D, E]
+         *     [A * B, C * D, E]
+         *     [(A * B) * (C * D), E]
+         *     [(A * B) * (C * D) * E].
+         *
+         * Each product uses the X-norm product identity. If a level has an odd number of representations, its final
+         * representation is carried unchanged to the next level. An empty input represents the empty product and
+         * therefore returns the multiplicative identity (1, 0).
+         */
+        template<polynomial_arithmetic::PolynomialBackend Backend>
+        polynomial_x_norm_representation<typename Backend::polynomial_type>
+            combine_polynomial_x_norm_representations_balanced(
+                std::vector<polynomial_x_norm_representation<typename Backend::polynomial_type>>
+                    representations,
+                polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context) {
+            using polynomial_type = typename Backend::polynomial_type;
+            using value_type = typename polynomial_type::value_type;
+            using representation_type = polynomial_x_norm_representation<polynomial_type>;
+
+            if (representations.empty()) {
+                return representation_type {polynomial_type {value_type::one()}, polynomial_type {value_type::zero()}};
+            }
+
+            while (representations.size() > 1) {
+                std::vector<representation_type> next_level;
+                next_level.reserve((representations.size() + 1) / 2);
+                std::size_t index = 0;
+                for (; index + 1 < representations.size(); index += 2) {
+                    next_level.push_back(multiply_polynomial_x_norm_representations<Backend>(
+                        representations[index], representations[index + 1], arithmetic_context));
+                }
+                if (index < representations.size()) {
+                    next_level.push_back(std::move(representations[index]));
+                }
+                representations = std::move(next_level);
+            }
+            return std::move(representations.front());
+        }
+
+    }    // namespace detail
+
+    /**
+     * Recover P and Q satisfying
+     *
+     *     P^2 - X * Q^2 = h
+     *
+     * for a canonical polynomial h. Zero and square constants are handled directly. A nonconstant input is factored
+     * into monic irreducible factors. Even factor multiplicities are represented as polynomial squares; odd
+     * multiplicities use recover_irreducible_polynomial_x_norm_representation. Factor representations are combined in
+     * a balanced product tree, then scaled by the square root of the factorization's leading coefficient.
+     *
+     * The coefficient generator remains caller-owned and is shared by complete factorization and irreducible-factor
+     * recovery. It must satisfy the documented requirements of both operations.
+     *
+     * @return a representation whose evaluated norm is exactly h; no value if a necessary coefficient square test,
+     *         odd-factor recovery, or leading-scalar normalization fails.
+     * @throws std::invalid_argument if h is empty or noncanonical, or a composed factorization or recovery contract is
+     *         violated.
+     * @throws std::logic_error if completed internal operations produce an inconsistent identity.
+     */
+    template<detail::SupportsDivrem Backend, typename Generator>
+        requires algebra::FieldValue<typename Backend::polynomial_type::value_type> &&
+                 std::constructible_from<typename Backend::polynomial_type, std::size_t> &&
+                 requires(typename Backend::polynomial_type &polynomial, Generator &generator,
+                          const typename Backend::polynomial_type::value_type &value) {
+                     polynomial[0] = generator();
+                     { value.is_square() } -> std::convertible_to<bool>;
+                 }
+    std::optional<polynomial_x_norm_representation<typename Backend::polynomial_type>>
+        recover_polynomial_x_norm_representation(const typename Backend::polynomial_type &h,
+                                                 polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context,
+                                                 Generator &coefficient_generator) {
+        using polynomial_type = typename Backend::polynomial_type;
+        using value_type = typename polynomial_type::value_type;
+        using representation_type = polynomial_x_norm_representation<polynomial_type>;
+
+        if (h.empty() || (h.size() > 1 && h[h.size() - 1] == value_type::zero())) {
+            throw std::invalid_argument("polynomial X-norm recovery requires a canonical nonempty polynomial");
+        }
+
+        const auto verify_exact = [&](representation_type representation) -> std::optional<representation_type> {
+            if (evaluate_polynomial_x_norm<Backend>(representation, arithmetic_context) != h) {
+                throw std::logic_error("polynomial X-norm recovery failed exact verification");
+            }
+            return std::move(representation);
+        };
+
+        if (is_zero(h)) {
+            return verify_exact(
+                representation_type {polynomial_type {value_type::zero()}, polynomial_type {value_type::zero()}});
+        }
+        if (h.size() == 1) {
+            if (!h[0].is_square()) {
+                return std::nullopt;
+            }
+            return verify_exact(representation_type {polynomial_type {algebra::fields::sqrt_known_square(h[0])},
+                                                     polynomial_type {value_type::zero()}});
+        }
+
+        // square filters
+        const std::size_t degree = h.size() - 1;
+        if (!h[0].is_square()) {
+            return std::nullopt;
+        }
+        value_type signed_leading_coefficient = h[h.size() - 1];
+        if ((degree & 1) != 0) {
+            signed_leading_coefficient = value_type::zero() - signed_leading_coefficient;
+        }
+        if (!signed_leading_coefficient.is_square()) {
+            return std::nullopt;
+        }
+
+        std::vector<representation_type> factor_representations;
+        bool factor_recovery_failed = false;
+        const auto factorization = complete_factorization<Backend>(
+            h, arithmetic_context, coefficient_generator, [&](const polynomial_factor<polynomial_type> &factor) {
+                auto representation = detail::recover_polynomial_x_norm_factor_power<Backend>(
+                    factor, arithmetic_context, coefficient_generator);
+                if (!representation) {
+                    factor_recovery_failed = true;
+                    return factorization_control::stop_factorization;
+                }
+                factor_representations.push_back(std::move(*representation));
+                return factorization_control::continue_factorization;
+            });
+
+        if (factor_recovery_failed) {
+            return std::nullopt;
+        }
+        if (!factorization.complete || factor_representations.empty()) {
+            throw std::logic_error("complete factorization did not produce all nonconstant factors");
+        }
+
+        representation_type result = detail::combine_polynomial_x_norm_representations_balanced<Backend>(
+            std::move(factor_representations), arithmetic_context);
+        const value_type leading_coefficient = factorization.leading_coefficient;
+        if (leading_coefficient.is_zero()) {
+            throw std::logic_error("complete factorization produced a zero leading coefficient");
+        }
+        if (!leading_coefficient.is_square()) {
+            return std::nullopt;
+        }
+        const value_type scalar = algebra::fields::sqrt_known_square(leading_coefficient);
+        scalar_multiplication(result.p, result.p, scalar);
+        scalar_multiplication(result.q, result.q, scalar);
+        return verify_exact(std::move(result));
+    }
+
 }    // namespace nil::crypto3::math
 
 #endif    // CRYPTO3_MATH_POLYNOMIAL_X_NORM_RECONSTRUCTION_HPP