// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// Float.hpp
// Arbitrary-precision floating-point class
//
// This file provides a class for high-precision floating-point arithmetic
// that uses an arbitrary-precision integer (Int) as the mantissa.
//
// Main features:
// - Arbitrary-precision floating-point numbers
// - Standard arithmetic operations and mathematical functions
// - IEEE 754-compliant rounding modes
// - Advanced numerical operations

#ifndef SANGI_FLOAT_HPP
#define SANGI_FLOAT_HPP

#include <math/core/mp/Int.hpp>
#include <math/concepts/algebraic_concepts.hpp>
#include <math/core/common.hpp>
#include <math/core/traits.hpp>

#include <string>
#include <string_view>
#include <iosfwd>
#include <compare>  // C++20 three-way comparison operator
#include <span>
#include <climits>      // INT_MAX
#include <type_traits>  // std::is_integral_v
#include <utility>      // std::pair
#include <cstddef>      // std::size_t
#include <stdexcept>    // std::domain_error (toSizet)

namespace sangi {

    // Precision propagation policy selection
    // MIN_PROPAGATION: min — stepwise precision, saves computation
    // MAX_PROPAGATION: max/merge — preserves precision, prioritizes safety
    enum class PrecisionPolicy { MIN_PROPAGATION, MAX_PROPAGATION };

#ifdef SANGI_PRECISION_POLICY_MIN
    inline constexpr PrecisionPolicy PRECISION_POLICY = PrecisionPolicy::MIN_PROPAGATION;
#else
    inline constexpr PrecisionPolicy PRECISION_POLICY = PrecisionPolicy::MAX_PROPAGATION;
#endif

    // Forward declaration of the Float class
    class Float;

    // Forward declarations of mathematical functions
    [[nodiscard]] Float exp(const Float& x, int precision);
    [[nodiscard]] Float exp(Float&& x, int precision);
    [[nodiscard]] Float log(const Float& x, int precision);
    [[nodiscard]] Float log(Float&& x, int precision);
    [[nodiscard]] Float logUi(unsigned long long n, int precision);
    [[nodiscard]] Float sin(const Float& x, int precision);
    [[nodiscard]] Float sin(Float&& x, int precision);
    [[nodiscard]] Float cos(const Float& x, int precision);
    [[nodiscard]] Float cos(Float&& x, int precision);
    [[nodiscard]] Float tan(const Float& x, int precision);
    [[nodiscard]] Float tan(Float&& x, int precision);
    [[nodiscard]] Float sqr(const Float& x, int precision);
    [[nodiscard]] Float sqr(Float&& x, int precision);
    [[nodiscard]] Float sqrt(const Float& x, int precision);
    [[nodiscard]] Float sqrt(Float&& x, int precision);
    [[nodiscard]] Float sinh(const Float& x, int precision);
    [[nodiscard]] Float sinh(Float&& x, int precision);
    [[nodiscard]] Float cosh(const Float& x, int precision);
    [[nodiscard]] Float cosh(Float&& x, int precision);
    [[nodiscard]] Float tanh(const Float& x, int precision);
    [[nodiscard]] Float tanh(Float&& x, int precision);
    [[nodiscard]] Float pow(const Float& x, const Float& y, int precision);
    [[nodiscard]] Float pow(Float&& x, Float&& y, int precision);
    [[nodiscard]] Float pow(const Float& x, int n, int precision);
    [[nodiscard]] Float pow(Float&& x, int n, int precision);
    [[nodiscard]] Float abs(const Float& value);
    [[nodiscard]] Float abs(Float&& value);
    [[nodiscard]] Float ldexp(const Float& value, int exp);
    [[nodiscard]] Float ldexp(Float&& value, int exp);
    [[nodiscard]] Float frexp(const Float& value, int* exp);
    [[nodiscard]] Float modf(const Float& value, Float* iptr);

    // Rounding modes (IEEE 754-compliant)
    enum class RoundingMode {
        ToNearest,      // Round to nearest (default)
        TowardZero,     // Round toward zero (truncation)
        TowardPositive, // Round toward positive infinity (ceiling)
        TowardNegative, // Round toward negative infinity (floor)
        AwayFromZero    // Round away from zero
    };

    // IEEE 754 exception flags
    enum FloatException : unsigned {
        FE_NONE      = 0,
        FE_INEXACT   = 1u << 0,   // Inexact (rounding occurred)
        FE_UNDERFLOW = 1u << 1,   // Underflow
        FE_OVERFLOW  = 1u << 2,   // Overflow
        FE_DIVBYZERO = 1u << 3,   // Division by zero
        FE_INVALID   = 1u << 4,   // Invalid operation (NaN generated)
        FE_ALL       = 0x1Fu
    };

    /**
     * @brief Arbitrary-precision floating-point class
     *
     * This implementation uses an arbitrary-precision integer (Int) as the mantissa
     * and stores the exponent separately, enabling high-precision floating-point arithmetic.
     * It conforms to the IEEE 754 concepts of precision and rounding modes.
     */
    // Note: Float has thread_local static members, so SANGI_API (dllexport) cannot
    // be applied to the entire class (MSVC C2492). WINDOWS_EXPORT_ALL_SYMBOLS is used instead.
    // Future: once thread_local is moved to function-scope static, SANGI_API can be applied.
    class Float {
        // Friend class declarations
        friend struct numeric_traits<Float>;
        friend class FloatOps;

        // Friend operator declarations
        friend Float operator+(const Float& value);  // Unary plus operator
        friend Float operator-(const Float& value);  // Unary minus operator
        friend Float operator+(const Float& lhs, const Float& rhs);
        friend Float operator+(Float&& lhs, const Float& rhs);
        friend Float operator+(const Float& lhs, Float&& rhs);
        friend Float operator+(Float&& lhs, Float&& rhs);
        friend Float operator-(const Float& lhs, const Float& rhs);
        friend Float operator-(Float&& lhs, const Float& rhs);
        friend Float operator-(const Float& lhs, Float&& rhs);
        friend Float operator-(Float&& lhs, Float&& rhs);
        friend Float operator*(const Float& lhs, const Float& rhs);
        friend Float operator*(Float&& lhs, const Float& rhs);
        friend Float operator*(const Float& lhs, Float&& rhs);
        friend Float operator*(Float&& lhs, Float&& rhs);
        friend Float operator/(const Float& lhs, const Float& rhs);
        friend Float operator/(Float&& lhs, const Float& rhs);
        friend Float mulScalarF(const Float& lhs, uint64_t rhs);
        friend Float mulScalarF(Float&& lhs, uint64_t rhs);
        friend Float mulScalarF(const Float& lhs, int64_t rhs);
        friend Float mulScalarF(Float&& lhs, int64_t rhs);
        friend Float divScalarF(const Float& lhs, uint64_t rhs);
        friend Float divScalarF(Float&& lhs, uint64_t rhs);
        friend Float divScalarF(const Float& lhs, int64_t rhs);
        friend Float divScalarF(Float&& lhs, int64_t rhs);
        friend std::partial_ordering operator<=>(const Float& lhs, const Float& rhs);
        friend bool operator==(const Float& lhs, const Float& rhs);

        // Friend declarations for mathematical functions
        friend Float abs(const Float& value);
        friend Float abs(Float&& value);
        friend Float exp(const Float& x, int precision);
        friend Float log(const Float& x, int precision);
        friend Float sin(const Float& x, int precision);
        friend Float cos(const Float& x, int precision);
        friend Float tan(const Float& x, int precision);
        friend Float sqr(const Float& x, int precision);
        friend Float sqrt(const Float& x, int precision);
        friend Float importBinaryFloat(std::span<const uint8_t> data);
        friend Float sinh(const Float& x, int precision);
        friend Float cosh(const Float& x, int precision);
        friend Float tanh(const Float& x, int precision);
        friend Float pow(const Float& x, const Float& y, int precision);
        friend Float pow(const Float& x, int n, int precision);
        friend Float ldexp(const Float& value, int exp);
        friend Float ldexp(Float&& value, int exp);
        friend Float frexp(const Float& value, int* exp);

        // Friend declarations for utility functions
        friend void swap(Float& a, Float& b) noexcept;
        friend Float copySign(const Float& x, const Float& y);
        friend Float copySign(Float&& x, const Float& y);

    private:
        Int mantissa_;       // Mantissa (always non-negative; the sign is managed by is_negative_)
        int64_t exponent_;   // Exponent
        bool is_negative_;   // Sign flag (independent of the Int's sign)
            // mantissa_ is always kept with m_sign >= 0. Reasons for storing the sign separately:
            // (1) mpn functions assume unsigned arrays, eliminating signed branches
            // (2) To represent IEEE 754 -0.0 (in Int, the sign of zero is always 0)

        // Special-value flags
        bool is_infinity_;   // Infinity flag
        bool is_nan_;        // NaN flag

        // Precision management fields
        // effective_bits_: number of reliable bits actually contained in the mantissa
        //   - multiplication/division: min(lhs, rhs)
        //   - addition: min(lhs, rhs) (no catastrophic cancellation)
        //   - subtraction: min(lhs, rhs) - lost_bits (accounts for cancellation)
        //   - integer-derived (exact): INT_MAX
        //   - double-derived: 53
        // requested_bits_: target precision required for this value
        //   - operation result: max(lhs, rhs)
        //   - can be set explicitly via setPrecision()
        //   - integer-derived (exact): INT_MAX
        int effective_bits_; // Effective bit count (reliable precision)
        int requested_bits_; // Requested bit count (target precision)

        // Precision and mode settings
        // inline: ODR-compliant header definition. The variable exists in the caller's TU even across DLL boundaries.
        inline static thread_local int default_precision_ = 53;  // Default matches double precision
        // Ambient working precision (bits) for the current call scope. 0 = unset → fall back to
        // defaultPrecision(). Set via PrecisionScope at a special-function entry so that
        // context-free operations (exact÷exact, sqrt(exact), pow(exact, fractional)) compute at the
        // caller's working precision instead of defaulting to 53 digits (the exact/exact poison fix).
        inline static thread_local int working_precision_bits_ = 0;
        // ── exact/exact poison audit hook (opt-in; off by default, so the normal cost is just one bool read) ──
        // When ExactDivAudit is applied, it counts how many times a context-free operation
        // (exact÷exact / a unary function with exact input) fell back to default precision while no scope was set (working_precision_bits_==0). A safety net
        // for mechanically detecting from tests the "special functions that forgot to apply PrecisionScope" when rolling out broadly.
        inline static thread_local bool exact_audit_enabled_ = false;
        inline static thread_local long exact_audit_unscoped_hits_ = 0;
        inline static thread_local RoundingMode rounding_mode_ = RoundingMode::ToNearest;
        static thread_local int64_t emin_;                 // Minimum exponent (IEEE 754 emin)
        static thread_local int64_t emax_;                 // Maximum exponent (IEEE 754 emax)
        static thread_local unsigned exception_flags_;     // IEEE 754 exception flags

        // Exponent overflow guard constants
        // Values outside this range are rounded to ±∞ or 0
        static constexpr int64_t EXPONENT_MAX = (1LL << 60);
        static constexpr int64_t EXPONENT_MIN = -(1LL << 60);

        // Normalization function (adjusts the mantissa and updates the exponent)
        void normalize();

        // Exponent range check (overflow → ±∞, underflow → 0)
        void checkExponentBounds();

        // Common implementation that parses a decimal string and stores it in the body.
        // If requested_precision >= 0, non-exact values are rounded to that decimal precision;
        // if < 0, rounded to max(significant digits of the string, defaultPrecision()).
        void initFromDecimalString(std::string_view str, int requested_precision);

    public:
        // Type definitions (compatible with standard containers)
        using value_type = Float;
        using size_type = std::size_t;
        using difference_type = std::ptrdiff_t;
        using reference = Float&;
        using const_reference = const Float&;

        // Constructors and destructor

        /**
         * @brief Default constructor (initializes to zero)
         */
        Float();

        /**
         * @brief Constructor from an integer
         * @param value Initial value
         */
        explicit Float(int value);

        /**
         * @brief Constructor from a 64-bit integer
         * @param value Initial value
         */
        explicit Float(int64_t value);

        /**
         * @brief Constructor from an unsigned 64-bit integer (covers size_t)
         * @param value Initial value
         * @note Kept explicit (consistent with the other integer ctors, which are
         *       explicit by design — Float mixes with integers via the heterogeneous
         *       operator overloads, not implicit conversion). This ctor disambiguates
         *       `static_cast<Float>(size_t)`, otherwise ambiguous between the
         *       int64_t/double ctors. Delegates to Int so the full range is exact.
         */
        explicit Float(uint64_t value) : Float(Int(value)) {}

        /**
         * @brief Constructor taking an integer mantissa and an exponent
         * @param mantissa Integer mantissa
         * @param exponent Exponent
         * @param is_negative Sign flag
         */
        explicit Float(int64_t mantissa, int64_t exponent, bool is_negative = false);

        /**
         * @brief Constructor from an arbitrary-precision integer
         * @param value Initial value
         */
        explicit Float(const Int& value);

        /**
         * @brief Constructor from a floating-point number
         * @param value Initial value
         */
        explicit Float(double value);

        /**
         * @brief Constructor from a string
         * @param str String representing the numeric value
         *
         * A non-exact decimal string (e.g. "0.1" = 1/10) is rounded and stored at the precision of
         * the current defaultPrecision() (or the significant digits of the string, whichever is larger).
         * When high precision is required, use the 2-argument version with a precision argument
         * (the precision can be specified per literal without depending on the global defaultPrecision).
         */
        explicit Float(std::string_view str);

        /**
         * @brief Constructor from a string with an explicit target precision
         * @param str String representing the numeric value
         * @param precision Number of decimal digits of precision to which non-exact values are rounded
         *
         * `Float("0.1", 200)` stores 1/10 at 200 digits of precision.
         * Integers and binary finite decimals ("0.5" etc.) are mathematically exact, so regardless of precision
         * they are treated as INT_MAX (full precision). This
         * avoids a ceiling caused by insufficient precision of the input literal in calls
         * such as `gamma(Float("0.1", 200), 200)`.
         */
        Float(std::string_view str, int precision);

        /**
         * @brief Constructor taking a mantissa and an exponent
         * @param mantissa Mantissa
         * @param exponent Exponent
         * @param is_negative Sign flag
         */
        Float(const Int& mantissa, int64_t exponent, bool is_negative = false);
        Float(Int&& mantissa, int64_t exponent, bool is_negative = false);

        /**
         * @brief Copy constructor
         */
        Float(const Float& other) = default;

        /**
         * @brief Move constructor
         */
        Float(Float&& other) noexcept = default;

        /**
         * @brief Copy assignment operator
         */
        Float& operator=(const Float& other) = default;

        /**
         * @brief Move assignment operator
         */
        Float& operator=(Float&& other) noexcept = default;

        /**
         * @brief Assignment from an integer value
         */
        Float& operator=(int64_t value);

        /**
         * @brief Assignment from a floating-point value
         */
        Float& operator=(double value);

        /**
         * @brief Addition-assignment operator
         */
        Float& operator+=(const Float& rhs);

        /**
         * @brief Subtraction-assignment operator
         */
        Float& operator-=(const Float& rhs);

        /**
         * @brief Multiplication-assignment operator
         */
        Float& operator*=(const Float& rhs);

        /**
         * @brief Division-assignment operator
         */
        Float& operator/=(const Float& rhs);

        // Heterogeneous compound assignment (integer / double).
        // Float mixes with integers and doubles via the binary operator overloads
        // because the Float constructors are explicit by design. The matching
        // compound forms below let Float act as a drop-in scalar T in generic numeric
        // code: `x /= 2`, `x *= n`, `x += 0.5` etc. They delegate to the member
        // operator OP=(const Float&) above (declared earlier, so no dependency on the
        // later free binary operators) and build the Float through int64_t/uint64_t
        // to avoid the ambiguity of Float(unsigned). Purely additive: no Float
        // compound assignment with a scalar existed before.
        template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
        Float& operator+=(T rhs) {
            if constexpr (std::is_signed_v<T>) return *this += Float(static_cast<int64_t>(rhs));
            else return *this += Float(static_cast<uint64_t>(rhs));
        }
        template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
        Float& operator-=(T rhs) {
            if constexpr (std::is_signed_v<T>) return *this -= Float(static_cast<int64_t>(rhs));
            else return *this -= Float(static_cast<uint64_t>(rhs));
        }
        template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
        Float& operator*=(T rhs) {
            if constexpr (std::is_signed_v<T>) return *this *= Float(static_cast<int64_t>(rhs));
            else return *this *= Float(static_cast<uint64_t>(rhs));
        }
        template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
        Float& operator/=(T rhs) {
            if constexpr (std::is_signed_v<T>) return *this /= Float(static_cast<int64_t>(rhs));
            else return *this /= Float(static_cast<uint64_t>(rhs));
        }
        Float& operator+=(double rhs) { return *this += Float(rhs); }
        Float& operator-=(double rhs) { return *this -= Float(rhs); }
        Float& operator*=(double rhs) { return *this *= Float(rhs); }
        Float& operator/=(double rhs) { return *this /= Float(rhs); }

        /**
         * @brief Destructor
         */
        ~Float() = default;

        // Basic accessors

        /**
         * @brief Get the mantissa
         * @return Reference to the mantissa
         */
        [[nodiscard]] const Int& mantissa() const { return mantissa_; }

        /**
         * @brief Get the exponent
         * @return Exponent value
         */
        [[nodiscard]] int64_t exponent() const { return exponent_; }

        /**
         * @brief Get the sign
         * @return true if the value is negative
         */
        [[nodiscard]] bool isNegative() const { return is_negative_; }

        /**
         * @brief Check whether the value is zero
         * @return true if zero
         */
        [[nodiscard]] bool isZero() const;

        /**
         * @brief Check whether the value is infinity
         * @return true if infinity
         */
        [[nodiscard]] bool isInfinity() const { return is_infinity_; }

        /**
         * @brief Check whether the value is NaN
         * @return true if NaN
         */
        [[nodiscard]] bool isNaN() const { return is_nan_; }

        /**
         * @brief Check whether the value is exact (integer-derived)
         * @return true if exact (adapts to the other operand's precision in operations)
         */
        [[nodiscard]] bool isExact() const { return effective_bits_ >= INT_MAX && requested_bits_ >= INT_MAX; }

        /**
         * @brief Get the effective bit count
         * @return Effective bit count
         */
        [[nodiscard]] int effectiveBits() const { return effective_bits_; }

        /// Directly set the effective bit count (used when mathematical functions reflect input precision)
        void setEffectiveBits(int bits) { effective_bits_ = bits; }

        /**
         * @brief Get the requested bit count
         * @return Requested bit count
         */
        [[nodiscard]] int requestedBits() const { return requested_bits_; }

        /**
         * @brief Get the decimal precision (number of significant digits)
         * @return Number of significant digits
         */
        [[nodiscard]] int precision() const;

        /**
         * @brief Set the decimal precision (number of significant digits)
         * @param precision New number of significant digits
         * @return Reference to this object
         */
        Float& setPrecision(int precision);

        /**
         * @brief Set the precision fields on a mathematical-function result
         * @param precision Computation precision (decimal digits)
         * @return Reference to this object
         *
         * After calling setPrecision(precision), sets effective_bits_ to
         * precisionToBits(precision). Used just before a mathematical function returns its result.
         */
        Float& setResultPrecision(int precision);

        /**
         * @brief Fast word-aligned approximate truncation (for intermediate computations)
         * @param precision Target precision (decimal digits)
         *
         * Unlike setPrecision, this truncates the low-order part on word boundaries without rounding.
         * Up to 63 extra bits may be retained, but they are absorbed by guard bits.
         * Used for precision control of intermediate values inside loops; use setPrecision for final results.
         */
        void truncateToApprox(int precision);

        // Utilities for converting between precision and bit counts
        [[nodiscard]] static int precisionToBits(int precision);
        [[nodiscard]] static int bitsToPrecision(int bits);

        /// Helper that determines an appropriate decimal precision from the input's requestedBits
        /// Falls back to defaultPrecision for special values (infinity, NaN, zero) and exact values
        [[nodiscard]] static int requestedPrecision(const Float& x) {
            int rb = x.requestedBits();
            // contextPrecision(): the scope precision if a PrecisionScope exists, otherwise defaultPrecision().
            // This way unary math functions with exact input (sqrt(Float(3)) etc.) are also computed at scope precision,
            // eliminating the "exact input → default fallback" that is isomorphic to exact/exact poison.
            int dp = contextPrecision();
            if (rb <= 0 || rb >= INT_MAX) { auditContextFreeFallback(); return dp; }
            return std::max(bitsToPrecision(rb), dp);
        }

        /**
         * @brief Get the default decimal precision
         * @return Default number of significant digits
         */
        [[nodiscard]] static int defaultPrecision() { return default_precision_; }

        /**
         * @brief Set the default decimal precision
         * @param precision New default number of significant digits
         */
        static void setDefaultPrecision(int precision) { default_precision_ = precision; }

        /**
         * @brief Working precision (bits) to use for context-free operations.
         *
         * Returns the ambient PrecisionScope precision if one is active, otherwise falls back to
         * defaultPrecision(). Operations whose result is non-terminating yet have no precision to
         * inherit from their operands (exact÷exact, sqrt of a non-square exact, pow with a
         * fractional exact exponent) consult this instead of blindly using defaultPrecision().
         */
        [[nodiscard]] static int contextBits() {
            return working_precision_bits_ > 0 ? working_precision_bits_
                                               : precisionToBits(default_precision_);
        }

        /// Decimal-precision companion of contextBits(): the active PrecisionScope precision if set,
        /// otherwise defaultPrecision(). Used by the inline math overloads as the exact-input fallback.
        [[nodiscard]] static int contextPrecision() {
            return working_precision_bits_ > 0 ? bitsToPrecision(working_precision_bits_)
                                               : default_precision_;
        }

        /// Records one "context-free op fell back to default precision with no PrecisionScope active".
        /// No-op unless an ExactDivAudit is in scope. Called from operator/ and the exact-input math
        /// fallbacks so a test can assert that a special function scopes all of its irrational ops.
        static void auditContextFreeFallback() {
            if (exact_audit_enabled_ && working_precision_bits_ == 0) ++exact_audit_unscoped_hits_;
        }

        /**
         * @brief RAII guard that sets the ambient working precision for the current scope.
         *
         * Place one at the top of a special-function entry point:
         *   Float airyBi(const Float& x, int precision) {
         *       Float::PrecisionScope _ps(precision);
         *       ...                       // exact÷exact / sqrt(exact) / pow(exact,frac) now use `precision`
         *   }
         * Nested scopes stack (innermost wins) and the previous value is restored on destruction, so
         * the result stays bit-reproducible (the precision is set explicitly from the call argument).
         * Each thread carries its own value, so OpenMP workers must open their own scope.
         */
        struct PrecisionScope {
            int saved_;
            explicit PrecisionScope(int precision)
                : saved_(working_precision_bits_) {
                working_precision_bits_ = precisionToBits(precision);
            }
            ~PrecisionScope() { working_precision_bits_ = saved_; }
            PrecisionScope(const PrecisionScope&) = delete;
            PrecisionScope& operator=(const PrecisionScope&) = delete;
        };

        /**
         * @brief RAII audit guard that counts un-scoped context-free fallbacks within its lifetime.
         *
         * Enables the audit hook and zeroes the counter on construction; restores the previous
         * audit state on destruction. Use in a test to verify that a special function has wrapped
         * all of its irrational constant work (exact÷exact, sqrt(exact), pow(exact,frac)) in a
         * PrecisionScope — a properly-scoped function reports hits()==0:
         *   { Float::ExactDivAudit a; (void)someSpecialFn(x, 200); EXPECT_EQ(a.hits(), 0); }
         * A non-zero count means "this code path computed an irrational at defaultPrecision with no
         * scope active" → add a PrecisionScope at the entry. Opt-in: zero cost when not constructed.
         */
        struct ExactDivAudit {
            bool saved_enabled_;
            long saved_hits_;
            ExactDivAudit()
                : saved_enabled_(exact_audit_enabled_), saved_hits_(exact_audit_unscoped_hits_) {
                exact_audit_enabled_ = true;
                exact_audit_unscoped_hits_ = 0;
            }
            ~ExactDivAudit() {
                exact_audit_enabled_ = saved_enabled_;
                exact_audit_unscoped_hits_ = saved_hits_;
            }
            [[nodiscard]] long hits() const { return exact_audit_unscoped_hits_; }
            ExactDivAudit(const ExactDivAudit&) = delete;
            ExactDivAudit& operator=(const ExactDivAudit&) = delete;
        };

        /**
         * @brief Get the rounding mode
         * @return Current rounding mode
         */
        [[nodiscard]] static RoundingMode roundingMode() { return rounding_mode_; }

        /**
         * @brief Set the rounding mode
         * @param mode New rounding mode
         */
        static void setRoundingMode(RoundingMode mode) { rounding_mode_ = mode; }

        // ── Exponent-range control (IEEE 754 emin/emax) ──

        /// Get the minimum exponent
        [[nodiscard]] static int64_t emin() { return emin_; }
        /// Get the maximum exponent
        [[nodiscard]] static int64_t emax() { return emax_; }
        /// Set the minimum exponent
        static void setEmin(int64_t e) { emin_ = e; }
        /// Set the maximum exponent
        static void setEmax(int64_t e) { emax_ = e; }
        /// Reset the exponent range to the default (±2^60)
        static void resetExponentRange() { emin_ = EXPONENT_MIN; emax_ = EXPONENT_MAX; }

        /// Subnormal emulation: right-shifts the mantissa when the exponent falls below emin
        void subnormalize(RoundingMode mode = rounding_mode_);

        // ── IEEE 754 exception-flag management ──

        /// Get the exception flags
        [[nodiscard]] static unsigned getExceptionFlags() { return exception_flags_; }
        /// Set the specified exception flags
        static void raiseException(unsigned flags) { exception_flags_ |= (flags & FE_ALL); }
        /// Clear the specified exception flags
        static void clearExceptionFlags(unsigned flags = FE_ALL) { exception_flags_ &= ~(flags & FE_ALL); }
        /// Test whether the specified exceptions have been raised
        [[nodiscard]] static bool testException(unsigned flags) { return (exception_flags_ & flags) != 0; }

        /**
         * @brief Convert to a string
         * @param precision Output precision (default: current precision)
         * @return String representation
         */
        [[nodiscard]] std::string toString(int precision = -1) const;

        /// Base-N representation (e.g. binary "11.0101", hexadecimal "FF.A5")
        /// @param base      Radix (2-36)
        /// @param fracDigits Number of fractional digits (-1 to determine automatically from all mantissa bits)
        [[nodiscard]] std::string toString(int base, int fracDigits) const;

        /**
         * @brief Convert to a decimal string
         * @param precision Output precision (default: current precision)
         * @return Decimal string
         */
        [[nodiscard]] std::string toDecimalString(int precision = -1) const;

        /**
         * @brief Convert to a scientific-notation string
         * @param precision Output precision (default: current precision)
         * @return Scientific-notation string
         */
        [[nodiscard]] std::string toScientificString(int precision = -1) const;

        /**
         * @brief Convert to double (precision may be lost)
         * @return double value
         */
        [[nodiscard]] double toDouble() const;

        // Special-value constructors

        /**
         * @brief Construct positive infinity
         * @return Positive infinity
         */
        [[nodiscard]] static Float positiveInfinity();

        /**
         * @brief Construct negative infinity
         * @return Negative infinity
         */
        [[nodiscard]] static Float negativeInfinity();

        /**
         * @brief Construct NaN
         * @return NaN
         */
        [[nodiscard]] static Float nan();

        /**
         * @brief Construct machine epsilon
         * @param precision Precision (default: current precision)
         * @return Machine epsilon
         */
        [[nodiscard]] static Float epsilon(int precision = default_precision_);

        /**
         * @brief Construct zero with the specified precision
         * @param precision Precision
         * @return Zero at the specified precision
         */
        [[nodiscard]] static Float zero(int precision = default_precision_);

        /**
         * @brief Construct one with the specified precision
         * @param precision Precision
         * @return One at the specified precision
         */
        [[nodiscard]] static Float one(int precision = default_precision_);

        /**
         * @brief Construct π
         * @param precision Precision (default: current precision)
         * @return π at the specified precision
         */
        [[nodiscard]] static Float pi(int precision = default_precision_);

        /**
         * @brief Construct the base of the natural logarithm, e
         * @param precision Precision (default: current precision)
         * @return e at the specified precision
         */
        [[nodiscard]] static Float e(int precision = default_precision_);

        /**
         * @brief Construct log(2)
         * @param precision Precision (default: current precision)
         * @return log(2) at the specified precision
         */
        [[nodiscard]] static Float log2(int precision = default_precision_);

        /**
         * @brief Construct log(10)
         * @param precision Precision (default: current precision)
         * @return log(10) at the specified precision
         */
        [[nodiscard]] static Float log10(int precision = default_precision_);

        /**
         * @brief Construct the Euler-Mascheroni constant γ
         * @param precision Precision (default: current precision)
         * @return γ ≈ 0.5772156649... at the specified precision
         */
        [[nodiscard]] static Float euler(int precision = default_precision_);

        /**
         * @brief Construct Catalan's constant G
         * @param precision Precision (default: current precision)
         * @return G ≈ 0.9159655941... at the specified precision
         */
        [[nodiscard]] static Float catalan(int precision = default_precision_);

        /**
         * @brief Construct √2
         * @param precision Precision (default: current precision)
         * @return √2 ≈ 1.4142135623... at the specified precision
         */
        [[nodiscard]] static Float sqrt2(int precision = default_precision_);

        /**
         * @brief Construct the lemniscate constant ω
         * @param precision Precision (default: current precision)
         * @return ω ≈ 2.6220575542... at the specified precision
         *
         * Uses AGM: ω = 4·(a+b) / (t·√2)
         * where a, b, t are the convergence values of AGM(1, 1/√2)
         */
        [[nodiscard]] static Float lemniscate(int precision = default_precision_);

        /**
         * @brief Construct Γ(1/4)
         * @param precision Precision (default: current precision)
         * @return Γ(1/4) ≈ 3.6256099082... at the specified precision
         *
         * Uses AGM: Γ(1/4) = √(ω·√(2π))
         * = √(lemniscate · √(2π))
         */
        [[nodiscard]] static Float gamma14(int precision = default_precision_);

        /**
         * @brief Construct Apéry's constant ζ(3)
         * @param precision Precision (default: current precision)
         * @return ζ(3) ≈ 1.2020569031... at the specified precision
         *
         * Amdeberhan-Zeilberger second formula (Binary Splitting)
         */
        [[nodiscard]] static Float zeta3(int precision = default_precision_);

        /**
         * @brief Construct ζ(5)
         * @param precision Precision (default: current precision)
         * @return ζ(5) ≈ 1.0369277551... at the specified precision
         *
         * Broadhurst BBP formula (Binary Splitting)
         */
        [[nodiscard]] static Float zeta5(int precision = default_precision_);

        // --- π-derived constants ---
        [[nodiscard]] static Float half_pi(int precision = default_precision_);       // π/2
        [[nodiscard]] static Float quarter_pi(int precision = default_precision_);    // π/4
        [[nodiscard]] static Float two_pi(int precision = default_precision_);        // 2π (τ)
        [[nodiscard]] static Float inv_pi(int precision = default_precision_);        // 1/π
        [[nodiscard]] static Float two_inv_pi(int precision = default_precision_);    // 2/π
        [[nodiscard]] static Float inv_sqrt_pi(int precision = default_precision_);   // 1/√π
        [[nodiscard]] static Float two_inv_sqrt_pi(int precision = default_precision_); // 2/√π

        // --- Square roots and higher roots ---
        [[nodiscard]] static Float sqrt3(int precision = default_precision_);         // √3
        [[nodiscard]] static Float sqrt5(int precision = default_precision_);         // √5
        [[nodiscard]] static Float inv_sqrt2(int precision = default_precision_);     // 1/√2 = √2/2
        [[nodiscard]] static Float cbrt2(int precision = default_precision_);         // ∛2

        // --- Logarithms ---
        [[nodiscard]] static Float ln3(int precision = default_precision_);           // log(3)
        [[nodiscard]] static Float ln5(int precision = default_precision_);           // log(5)
        [[nodiscard]] static Float log2e(int precision = default_precision_);         // log₂(e) = 1/ln2
        [[nodiscard]] static Float log10e(int precision = default_precision_);        // log₁₀(e) = 1/ln10

        // --- Golden ratio and basic constants ---
        [[nodiscard]] static Float phi(int precision = default_precision_);           // φ = (1+√5)/2
        [[nodiscard]] static Float sin1(int precision = default_precision_);          // sin(1)
        [[nodiscard]] static Float cos1(int precision = default_precision_);          // cos(1)
        [[nodiscard]] static Float degree(int precision = default_precision_);        // π/180
        [[nodiscard]] static Float egamma_exp(int precision = default_precision_);    // e^γ

        // --- Zeta / series constants ---
        [[nodiscard]] static Float zeta7(int precision = default_precision_);         // ζ(7)

        // --- Special constants ---
        [[nodiscard]] static Float glaisher(int precision = default_precision_);      // Glaisher-Kinkelin A
        [[nodiscard]] static Float khinchin(int precision = default_precision_);      // Khinchin K
        [[nodiscard]] static Float omega(int precision = default_precision_);         // Ω (Lambert W(1))
        [[nodiscard]] static Float plastic(int precision = default_precision_);       // plastic number ρ
        [[nodiscard]] static Float twin_prime(int precision = default_precision_);    // twin prime C₂
        [[nodiscard]] static Float landau_ramanujan(int precision = default_precision_); // Landau-Ramanujan K
        [[nodiscard]] static Float meissel_mertens(int precision = default_precision_); // Meissel-Mertens M
        [[nodiscard]] static Float bernstein(int precision = default_precision_);     // Bernstein β
        [[nodiscard]] static Float gauss_kuzmin(int precision = default_precision_);  // Gauss-Kuzmin-Wirsing λ
        [[nodiscard]] static Float feigenbaum_delta(int precision = default_precision_); // Feigenbaum δ
        [[nodiscard]] static Float feigenbaum_alpha(int precision = default_precision_); // Feigenbaum α
        [[nodiscard]] static Float erdos_borwein(int precision = default_precision_); // Erdős-Borwein E
        [[nodiscard]] static Float laplace_limit(int precision = default_precision_); // Laplace limit λ*
        [[nodiscard]] static Float soldner(int precision = default_precision_);       // Ramanujan-Soldner μ
        [[nodiscard]] static Float backhouse(int precision = default_precision_);     // Backhouse B
        [[nodiscard]] static Float porter(int precision = default_precision_);        // Porter C
        [[nodiscard]] static Float lieb_square_ice(int precision = default_precision_); // Lieb (8√3/9)
        [[nodiscard]] static Float niven(int precision = default_precision_);         // Niven C
        [[nodiscard]] static Float reciprocal_fibonacci(int precision = default_precision_); // ψ = Σ 1/F(k)
        [[nodiscard]] static Float sierpinski(int precision = default_precision_);    // Sierpiński K
        [[nodiscard]] static Float mills(int precision = default_precision_);         // Mills θ
        [[nodiscard]] static Float dottie(int precision = default_precision_);        // Dottie d
        [[nodiscard]] static Float golomb_dickman(int precision = default_precision_); // Golomb-Dickman λ
        [[nodiscard]] static Float salem(int precision = default_precision_);         // smallest Salem (Lehmer) τ
        [[nodiscard]] static Float cahen(int precision = default_precision_);         // Cahen C
        [[nodiscard]] static Float levy(int precision = default_precision_);          // Lévy β
        [[nodiscard]] static Float copeland_erdos(int precision = default_precision_); // Copeland-Erdős C_CE
        [[nodiscard]] static Float champernowne(int precision = default_precision_);   // Champernowne C_10
        [[nodiscard]] static Float liouville(int precision = default_precision_);      // Liouville Σ10^{-k!}
        [[nodiscard]] static Float alladi_grinstead(int precision = default_precision_); // Alladi-Grinstead
        [[nodiscard]] static Float hafner_sarnak_mccurley(int precision = default_precision_); // HSM
        [[nodiscard]] static Float lengyel(int precision = default_precision_);        // Lengyel L
        [[nodiscard]] static Float prime_quadruplet(int precision = default_precision_); // Hardy-Littlewood prime quadruplet
        [[nodiscard]] static Float pi_squared_over_6(int precision = default_precision_); // π²/6 = ζ(2)
        [[nodiscard]] static Float pi_squared_over_12(int precision = default_precision_); // π²/12 = ζ(2)/2

        /**
         * @brief Convert to Int
         * @return Int value (fractional part is truncated)
         */
        [[nodiscard]] Int toInt() const;

        /**
         * @brief Convert to a machine int64_t (fractional part truncated toward zero)
         * @note For converting a Float coordinate/count to a machine integer
         *       (e.g. an array index) WITHOUT going through double — keeps the
         *       conversion exact for in-range values. Throws on NaN/Inf (via toInt()).
         */
        [[nodiscard]] int64_t toInt64() const { return toInt().toInt64(); }

        /**
         * @brief Convert to size_t (fractional part truncated toward zero)
         * @throws std::domain_error if the value is negative (size_t is unsigned)
         * @note For array-index use. Avoids the double round-trip of static_cast<double>.
         */
        [[nodiscard]] std::size_t toSizet() const {
            if (is_negative_ && !isZero()) {
                throw std::domain_error(
                    "Float::toSizet: negative value cannot be converted to size_t");
            }
            return static_cast<std::size_t>(toInt().toInt64());
        }

        /**
         * @brief Get the bit length
         * @return Bit length
         */
        [[nodiscard]] int bitLength() const;

        /**
         * @brief Get the minimum precision (in bits) needed to represent x exactly
         * @return Minimum precision in bits (0 for 0, NaN, Inf)
         */
        [[nodiscard]] int minPrec() const;

        /**
         * @brief Test whether the value is positive
         * @return true if positive (>= 0)
         */
        [[nodiscard]] bool isPositive() const { return !isNegative(); }

        /**
         * @brief Test whether the value is an integer
         * @return true if integer (NaN/∞ return false)
         */
        [[nodiscard]] bool isInteger() const;

        /**
         * @brief Test whether the value fits in int (fractional part truncated)
         * @return true if within int range
         */
        [[nodiscard]] bool fitsInt() const;

        /**
         * @brief Test whether the value fits in int64_t (fractional part truncated)
         * @return true if within int64_t range
         */
        [[nodiscard]] bool fitsInt64() const;

        /**
         * @brief Test whether the value can be approximately represented as a double (range check)
         * @return true if within double range
         */
        [[nodiscard]] bool fitsDouble() const;

        /**
         * @brief Left-shift assignment operator
         * @param shift Shift amount
         * @return Reference to this object
         */
        Float& operator<<=(int shift);

        /**
         * @brief Right-shift assignment operator
         * @param shift Shift amount
         * @return Reference to this object
         */
        Float& operator>>=(int shift);

        /**
         * @brief Left-shift operator
         * @param shift Shift amount
         * @return Shift result
         */
        friend Float operator<<(const Float& value, int shift);

        /**
         * @brief Right-shift operator
         * @param shift Shift amount
         * @return Shift result
         */
        friend Float operator>>(const Float& value, int shift);

    private:
        // Internal implementation helper functions

        /**
         * @brief Implementation of unsigned addition
         * @param rhs Value to add
         * @return Result
         */
        // Tag type used to construct a Float directly from an mpn computation result
        struct PreNormalizedTag {};

        // Precondition: mantissa is already MSB/LSB normalized (no need for normalize()/trimZeroWords())
        Float(Int&& mantissa, int64_t exponent, bool is_negative, PreNormalizedTag)
            : mantissa_(std::move(mantissa)), exponent_(exponent), is_negative_(is_negative),
              is_infinity_(false), is_nan_(false),
              effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
            checkExponentBounds();
        }

        Float addUnsigned(const Float& rhs) const &;
        Float addUnsigned(const Float& rhs) &&;

        /**
         * @brief Implementation of unsigned subtraction
         * @param rhs Value to subtract
         * @return Result
         */
        Float subtractUnsigned(const Float& rhs) const &;
        Float subtractUnsigned(const Float& rhs) &&;

        // In-place versions: overwrite this->mantissa_ buffer directly.
        // Called from operator+= / -=. effective_bits_ / requested_bits_ are set by the caller.
        // is_negative_ is not modified (same semantics as addUnsigned / subtractUnsigned).
        void addUnsignedInPlace(const Float& rhs);
        void subtractUnsignedInPlace(const Float& rhs);

        /**
         * @brief Unsigned comparison
         * @param rhs Value to compare against
         * @return -1: this is smaller, 0: equal, 1: this is larger
         */
        int compareUnsigned(const Float& rhs) const;

        // P1-1: helper for the small-precision fast path
        // Constructs a Float directly from a pre-normalized limb buffer (skipping normalize)
        // Precondition: data[0] != 0 && data[n-1] != 0 (no trailing/leading zero limbs)
        static Float fromRawLimbs(const uint64_t* data, size_t n,
                                  int64_t exponent, bool negative,
                                  int eff, int req) {
            auto mint = Int::fromRawWordsPreNormalized(
                std::span<const uint64_t>(data, n), 1);
            Float result(std::move(mint), exponent, negative, PreNormalizedTag{});
            result.effective_bits_ = eff;
            result.requested_bits_ = req;
            return result;
        }

        /**
         * @brief Implementation of rounding
         * @param precision_bits Precision (in bits)
         * @param mode Rounding mode
         */
        void round(int precision_bits, RoundingMode mode = rounding_mode_);
    };

    // Declaration of numeric_traits specialization
    template<>
    struct numeric_traits<Float> {
        using value_type = Float;
        using category = floating_point_tag;

        static constexpr bool is_supported = true;
        static constexpr bool is_complex = false;
        static constexpr bool is_integer = false;
        static constexpr bool is_floating_point = true;

        static Float zero() { return Float::zero(); }
        static Float one() { return Float::one(); }
        static Float epsilon() { return Float::epsilon(); }

        static Float abs(const Float& value);
        static Float conj(const Float& value) { return value; }  // The conjugate of a real number is itself
        static Float norm(const Float& value);  // |x|^2
        static bool pivotBetter(const Float& a, const Float& b) {
            return abs(a) > abs(b);
        }
    };

    // Basic operators
    [[nodiscard]] Float operator+(const Float& lhs, const Float& rhs);
    [[nodiscard]] Float operator+(Float&& lhs, const Float& rhs);
    [[nodiscard]] Float operator+(const Float& lhs, Float&& rhs);
    [[nodiscard]] Float operator+(Float&& lhs, Float&& rhs);
    [[nodiscard]] Float operator-(const Float& lhs, const Float& rhs);
    [[nodiscard]] Float operator-(Float&& lhs, const Float& rhs);
    [[nodiscard]] Float operator-(const Float& lhs, Float&& rhs);
    [[nodiscard]] Float operator-(Float&& lhs, Float&& rhs);
    [[nodiscard]] Float operator-(const Float& value);  // Unary minus
    [[nodiscard]] Float operator*(const Float& lhs, const Float& rhs);
    [[nodiscard]] Float operator*(Float&& lhs, const Float& rhs);
    [[nodiscard]] Float operator*(const Float& lhs, Float&& rhs);
    [[nodiscard]] Float operator*(Float&& lhs, Float&& rhs);
    [[nodiscard]] Float operator/(const Float& lhs, const Float& rhs);
    [[nodiscard]] Float operator/(Float&& lhs, const Float& rhs);
    // Single-word multiplication/division — implementation functions (Float.cpp)
    [[nodiscard]] Float mulScalarF(const Float& lhs, uint64_t rhs);
    [[nodiscard]] Float mulScalarF(Float&& lhs, uint64_t rhs);
    [[nodiscard]] Float mulScalarF(const Float& lhs, int64_t rhs);
    [[nodiscard]] Float mulScalarF(Float&& lhs, int64_t rhs);
    [[nodiscard]] Float divScalarF(const Float& lhs, uint64_t rhs);
    [[nodiscard]] Float divScalarF(Float&& lhs, uint64_t rhs);
    [[nodiscard]] Float divScalarF(const Float& lhs, int64_t rhs);
    [[nodiscard]] Float divScalarF(Float&& lhs, int64_t rhs);

    // Template operators: Float * scalar, Float / scalar
    template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
    [[nodiscard]] Float operator*(const Float& lhs, T rhs) {
        if constexpr (std::is_signed_v<T>) return mulScalarF(lhs, static_cast<int64_t>(rhs));
        else return mulScalarF(lhs, static_cast<uint64_t>(rhs));
    }
    template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
    [[nodiscard]] Float operator*(Float&& lhs, T rhs) {
        if constexpr (std::is_signed_v<T>) return mulScalarF(std::move(lhs), static_cast<int64_t>(rhs));
        else return mulScalarF(std::move(lhs), static_cast<uint64_t>(rhs));
    }
    template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
    [[nodiscard]] Float operator*(T lhs, const Float& rhs) {
        if constexpr (std::is_signed_v<T>) return mulScalarF(rhs, static_cast<int64_t>(lhs));
        else return mulScalarF(rhs, static_cast<uint64_t>(lhs));
    }
    template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
    [[nodiscard]] Float operator/(const Float& lhs, T rhs) {
        if constexpr (std::is_signed_v<T>) return divScalarF(lhs, static_cast<int64_t>(rhs));
        else return divScalarF(lhs, static_cast<uint64_t>(rhs));
    }
    template<typename T, std::enable_if_t<std::is_integral_v<T> && !std::is_same_v<T, bool>, int> = 0>
    [[nodiscard]] Float operator/(Float&& lhs, T rhs) {
        if constexpr (std::is_signed_v<T>) return divScalarF(std::move(lhs), static_cast<int64_t>(rhs));
        else return divScalarF(std::move(lhs), static_cast<uint64_t>(rhs));
    }

    // Mixed operations with double — provided because Float(double) is explicit
    // No float version is needed: the float mantissa (23 bits) fits in one limb, so
    // Float(float) and Float(double) have the same internal representation.
    // C++'s implicit float→double conversion routes through these double overloads.
    [[nodiscard]] inline Float operator+(const Float& lhs, double rhs) { return lhs + Float(rhs); }
    [[nodiscard]] inline Float operator+(Float&& lhs, double rhs)      { return std::move(lhs) + Float(rhs); }
    [[nodiscard]] inline Float operator+(double lhs, const Float& rhs) { return Float(lhs) + rhs; }
    [[nodiscard]] inline Float operator-(const Float& lhs, double rhs) { return lhs - Float(rhs); }
    [[nodiscard]] inline Float operator-(Float&& lhs, double rhs)      { return std::move(lhs) - Float(rhs); }
    [[nodiscard]] inline Float operator-(double lhs, const Float& rhs) { return Float(lhs) - rhs; }
    [[nodiscard]] inline Float operator*(const Float& lhs, double rhs) { return lhs * Float(rhs); }
    [[nodiscard]] inline Float operator*(Float&& lhs, double rhs)      { return std::move(lhs) * Float(rhs); }
    [[nodiscard]] inline Float operator*(double lhs, const Float& rhs) { return Float(lhs) * rhs; }
    [[nodiscard]] inline Float operator/(const Float& lhs, double rhs) { return lhs / Float(rhs); }
    [[nodiscard]] inline Float operator/(Float&& lhs, double rhs)      { return std::move(lhs) / Float(rhs); }
    [[nodiscard]] inline Float operator/(double lhs, const Float& rhs) { return Float(lhs) / rhs; }
    [[nodiscard]] inline bool operator==(const Float& lhs, double rhs) { return lhs == Float(rhs); }
    [[nodiscard]] inline bool operator==(double lhs, const Float& rhs) { return Float(lhs) == rhs; }
    [[nodiscard]] inline std::partial_ordering operator<=>(const Float& lhs, double rhs) { return lhs <=> Float(rhs); }
    [[nodiscard]] inline std::partial_ordering operator<=>(double lhs, const Float& rhs) { return Float(lhs) <=> rhs; }

    [[nodiscard]] std::partial_ordering operator<=>(const Float& lhs, const Float& rhs);
    [[nodiscard]] bool operator==(const Float& lhs, const Float& rhs);

    // Stream I/O
    std::ostream& operator<<(std::ostream& os, const Float& value);
    std::istream& operator>>(std::istream& is, Float& value);

    // Floating-point utility functions
    [[nodiscard]] Float fmin(const Float& a, const Float& b);
    [[nodiscard]] Float fmin(Float&& a, Float&& b);
    [[nodiscard]] Float fmax(const Float& a, const Float& b);
    [[nodiscard]] Float fmax(Float&& a, Float&& b);
    [[nodiscard]] Float fdim(const Float& a, const Float& b);
    [[nodiscard]] Float fdim(Float&& a, Float&& b);
    [[nodiscard]] Float copySign(const Float& x, const Float& y);
    [[nodiscard]] Float copySign(Float&& x, const Float& y);
    [[nodiscard]] bool signBit(const Float& x);

    // Linear interpolation and midpoint
    [[nodiscard]] Float lerp(const Float& a, const Float& b, const Float& t, int precision);
    [[nodiscard]] Float midpoint(const Float& a, const Float& b);

    // <cmath>-compatible utilities
    [[nodiscard]] Float modf(const Float& x, Float& iptr);        // Splits into integer and fractional parts
    [[nodiscard]] int64_t ilogb(const Float& x);                   // floor(log2(|x|))
    [[nodiscard]] Float logb(const Float& x);                      // floor(log2(|x|)) as Float
    [[nodiscard]] Float scalbn(const Float& x, int n);             // x * 2^n (= ldexp)
    [[nodiscard]] Float scalbn(Float&& x, int n);
    [[nodiscard]] Float nearbyint(const Float& x);                 // Same as round
    [[nodiscard]] Float nearbyint(Float&& x);
    [[nodiscard]] Float rint(const Float& x);                      // Same as round
    [[nodiscard]] Float rint(Float&& x);

    // Rounding and integer-conversion functions
    [[nodiscard]] Float floor(const Float& x);
    [[nodiscard]] Float floor(Float&& x);
    [[nodiscard]] Float ceil(const Float& x);
    [[nodiscard]] Float ceil(Float&& x);
    [[nodiscard]] Float trunc(const Float& x);
    [[nodiscard]] Float trunc(Float&& x);
    [[nodiscard]] Float round(const Float& x);
    [[nodiscard]] Float round(Float&& x);
    [[nodiscard]] Float roundEven(const Float& x);               // Round half to even (banker's rounding)
    [[nodiscard]] Float roundEven(Float&& x);
    [[nodiscard]] Float frac(const Float& x);
    [[nodiscard]] Float frac(Float&& x);

    // Inverse trigonometric functions
    [[nodiscard]] Float asin(const Float& x, int precision);
    [[nodiscard]] Float asin(Float&& x, int precision);
    [[nodiscard]] Float acos(const Float& x, int precision);
    [[nodiscard]] Float acos(Float&& x, int precision);
    [[nodiscard]] Float atan(const Float& x, int precision);
    [[nodiscard]] Float atan(Float&& x, int precision);
    [[nodiscard]] Float atan2(const Float& y, const Float& x, int precision);
    [[nodiscard]] Float atan2(Float&& y, Float&& x, int precision);

    // Inverse trigonometric functions (result expressed in units of π)
    [[nodiscard]] Float asinPi(const Float& x, int precision);
    [[nodiscard]] Float asinPi(Float&& x, int precision);
    [[nodiscard]] Float acosPi(const Float& x, int precision);
    [[nodiscard]] Float acosPi(Float&& x, int precision);
    [[nodiscard]] Float atanPi(const Float& x, int precision);
    [[nodiscard]] Float atanPi(Float&& x, int precision);
    [[nodiscard]] Float atan2Pi(const Float& y, const Float& x, int precision);
    [[nodiscard]] Float atan2Pi(Float&& y, Float&& x, int precision);

    // Inverse hyperbolic functions
    [[nodiscard]] Float asinh(const Float& x, int precision);
    [[nodiscard]] Float asinh(Float&& x, int precision);
    [[nodiscard]] Float acosh(const Float& x, int precision);
    [[nodiscard]] Float acosh(Float&& x, int precision);
    [[nodiscard]] Float atanh(const Float& x, int precision);
    [[nodiscard]] Float atanh(Float&& x, int precision);

    // Logarithm/exponential variants
    [[nodiscard]] Float log2(const Float& x, int precision);
    [[nodiscard]] Float log2(Float&& x, int precision);
    [[nodiscard]] Float log10(const Float& x, int precision);
    [[nodiscard]] Float log10(Float&& x, int precision);
    [[nodiscard]] Float log1p(const Float& x, int precision);
    [[nodiscard]] Float log1p(Float&& x, int precision);
    [[nodiscard]] Float exp2(const Float& x, int precision);
    [[nodiscard]] Float exp2(Float&& x, int precision);
    [[nodiscard]] Float exp10(const Float& x, int precision);
    [[nodiscard]] Float exp10(Float&& x, int precision);
    [[nodiscard]] Float expm1(const Float& x, int precision);
    [[nodiscard]] Float expm1(Float&& x, int precision);
    [[nodiscard]] Float exp2m1(const Float& x, int precision);
    [[nodiscard]] Float exp2m1(Float&& x, int precision);
    [[nodiscard]] Float exp10m1(const Float& x, int precision);
    [[nodiscard]] Float exp10m1(Float&& x, int precision);
    [[nodiscard]] Float log2p1(const Float& x, int precision);
    [[nodiscard]] Float log2p1(Float&& x, int precision);
    [[nodiscard]] Float log10p1(const Float& x, int precision);
    [[nodiscard]] Float log10p1(Float&& x, int precision);
    [[nodiscard]] Float compound(const Float& x, int n, int precision);
    [[nodiscard]] Float compound(Float&& x, int n, int precision);

    // n-th root and reciprocal square root
    [[nodiscard]] Float cbrt(const Float& x, int precision);
    [[nodiscard]] Float cbrt(Float&& x, int precision);
    [[nodiscard]] Float nthRoot(const Float& x, int n, int precision);
    [[nodiscard]] Float nthRoot(Float&& x, int n, int precision);
    [[nodiscard]] Float recSqrt(const Float& x, int precision);
    [[nodiscard]] Float recSqrt(Float&& x, int precision);

    // Fused multiply-add
    [[nodiscard]] Float fma(const Float& a, const Float& b, const Float& c, int precision);
    [[nodiscard]] Float fma(Float&& a, Float&& b, Float&& c, int precision);
    [[nodiscard]] Float fms(const Float& a, const Float& b, const Float& c, int precision);
    [[nodiscard]] Float fms(Float&& a, Float&& b, Float&& c, int precision);

    // Sum/difference of two products (equivalent to MPFR mpfr_fmma/mpfr_fmms)
    [[nodiscard]] Float fmma(const Float& a, const Float& b, const Float& c, const Float& d, int precision);
    [[nodiscard]] Float fmms(const Float& a, const Float& b, const Float& c, const Float& d, int precision);

    // fmod / remainder / hypot
    [[nodiscard]] Float fmod(const Float& x, const Float& y);
    [[nodiscard]] Float fmod(Float&& x, Float&& y);
    [[nodiscard]] Float remainder(const Float& x, const Float& y);
    [[nodiscard]] Float remainder(Float&& x, Float&& y);
    [[nodiscard]] std::pair<Float, int> remquo(const Float& x, const Float& y);
    [[nodiscard]] std::pair<Float, int> remquo(Float&& x, Float&& y);
    [[nodiscard]] Float hypot(const Float& x, const Float& y, int precision);
    [[nodiscard]] Float hypot(Float&& x, Float&& y, int precision);

    // Joint computation
    void sinCos(const Float& x, Float& sin_result, Float& cos_result, int precision);
    void sinCos(Float&& x, Float& sin_result, Float& cos_result, int precision);
    void sinhCosh(const Float& x, Float& sinh_result, Float& cosh_result, int precision);
    void sinhCosh(Float&& x, Float& sinh_result, Float& cosh_result, int precision);

    // Reciprocal trigonometric and hyperbolic functions
    [[nodiscard]] Float sec(const Float& x, int precision);
    [[nodiscard]] Float sec(Float&& x, int precision);
    [[nodiscard]] Float csc(const Float& x, int precision);
    [[nodiscard]] Float csc(Float&& x, int precision);
    [[nodiscard]] Float cot(const Float& x, int precision);
    [[nodiscard]] Float cot(Float&& x, int precision);
    [[nodiscard]] Float sech(const Float& x, int precision);
    [[nodiscard]] Float sech(Float&& x, int precision);
    [[nodiscard]] Float csch(const Float& x, int precision);
    [[nodiscard]] Float csch(Float&& x, int precision);
    [[nodiscard]] Float coth(const Float& x, int precision);
    [[nodiscard]] Float coth(Float&& x, int precision);

    // Factorial
    [[nodiscard]] Float factorial(int n, int precision);

    // π-based trigonometric functions
    [[nodiscard]] Float sinPi(const Float& x, int precision);
    [[nodiscard]] Float sinPi(Float&& x, int precision);
    [[nodiscard]] Float cosPi(const Float& x, int precision);
    [[nodiscard]] Float cosPi(Float&& x, int precision);
    [[nodiscard]] Float tanPi(const Float& x, int precision);
    [[nodiscard]] Float tanPi(Float&& x, int precision);

    // Trigonometric functions with arbitrary angle units (IEEE 754-2019)
    // sinu(x, u) = sin(2πx/u), cosu(x, u) = cos(2πx/u), tanu(x, u) = tan(2πx/u)
    // u is the number of units in one full turn (360 = degrees, 400 = gradians)
    [[nodiscard]] Float sinu(const Float& x, int u, int precision);
    [[nodiscard]] Float sinu(Float&& x, int u, int precision);
    [[nodiscard]] Float cosu(const Float& x, int u, int precision);
    [[nodiscard]] Float cosu(Float&& x, int u, int precision);
    [[nodiscard]] Float tanu(const Float& x, int u, int precision);
    [[nodiscard]] Float tanu(Float&& x, int u, int precision);

    // nextAbove / nextBelow
    [[nodiscard]] Float nextAbove(const Float& x);
    [[nodiscard]] Float nextAbove(Float&& x);
    [[nodiscard]] Float nextBelow(const Float& x);
    [[nodiscard]] Float nextBelow(Float&& x);

    // Rounding-direction test (compatible with MPFR mpfr_can_round)
    // x: an approximation computed at err_bits of precision under rounding mode rnd1
    // target_prec: target precision (decimal digits)
    // Returns true if the rounding direction is uniquely determined when rounding to target_prec under rnd2
    [[nodiscard]] bool canRound(const Float& x, int err_bits,
                                RoundingMode rnd1, RoundingMode rnd2,
                                int target_prec);

    // AGM (arithmetic-geometric mean)
    [[nodiscard]] Float agm(const Float& a, const Float& b, int precision);
    [[nodiscard]] Float agm(Float&& a, Float&& b, int precision);

    // High-precision sum and dot product
    [[nodiscard]] Float sum(std::span<const Float> values, int precision);
    [[nodiscard]] Float dot(std::span<const Float> a, std::span<const Float> b, int precision);

    // Error functions
    [[nodiscard]] Float erf(const Float& x, int precision);
    [[nodiscard]] Float erf(Float&& x, int precision);
    [[nodiscard]] Float erfc(const Float& x, int precision);
    [[nodiscard]] Float erfc(Float&& x, int precision);
    [[nodiscard]] Float erfcx(const Float& x, int precision);    // exp(x²)·erfc(x)
    [[nodiscard]] Float erfcx(Float&& x, int precision);

    // Gamma function and related functions
    [[nodiscard]] Float lnGamma(const Float& x, int precision);
    [[nodiscard]] Float lnGamma(Float&& x, int precision);
    [[nodiscard]] Float gamma(const Float& x, int precision);
    [[nodiscard]] Float gamma(Float&& x, int precision);
    [[nodiscard]] Float beta(const Float& a, const Float& b, int precision);
    [[nodiscard]] Float beta(Float&& a, Float&& b, int precision);
    [[nodiscard]] Float digamma(const Float& x, int precision);
    [[nodiscard]] Float digamma(Float&& x, int precision);
    [[nodiscard]] Float trigamma(const Float& x, int precision);
    [[nodiscard]] Float trigamma(Float&& x, int precision);
    [[nodiscard]] Float polygamma(int n, const Float& x, int precision);
    [[nodiscard]] Float polygamma(int n, Float&& x, int precision);

    // Incomplete gamma functions
    [[nodiscard]] Float gammaP(const Float& a, const Float& x, int precision);     // P(a,x) = γ(a,x)/Γ(a) regularized lower
    [[nodiscard]] Float gammaP(Float&& a, Float&& x, int precision);
    [[nodiscard]] Float gammaQ(const Float& a, const Float& x, int precision);     // Q(a,x) = Γ(a,x)/Γ(a) regularized upper
    [[nodiscard]] Float gammaQ(Float&& a, Float&& x, int precision);
    [[nodiscard]] Float gammaLower(const Float& a, const Float& x, int precision); // γ(a,x) lower incomplete gamma
    [[nodiscard]] Float gammaLower(Float&& a, Float&& x, int precision);
    [[nodiscard]] Float gammaUpper(const Float& a, const Float& x, int precision); // Γ(a,x) upper incomplete gamma
    [[nodiscard]] Float gammaUpper(Float&& a, Float&& x, int precision);

    // Incomplete beta function
    [[nodiscard]] Float betaRegularized(const Float& x, const Float& a, const Float& b, int precision); // I_x(a,b)
    [[nodiscard]] Float betaRegularized(Float&& x, Float&& a, Float&& b, int precision);

    // Elliptic integrals — Carlson symmetric form
    [[nodiscard]] Float carlsonRF(const Float& x, const Float& y, const Float& z, int precision);
    [[nodiscard]] Float carlsonRF(Float&& x, Float&& y, Float&& z, int precision);
    [[nodiscard]] Float carlsonRC(const Float& x, const Float& y, int precision);
    [[nodiscard]] Float carlsonRC(Float&& x, Float&& y, int precision);
    [[nodiscard]] Float carlsonRD(const Float& x, const Float& y, const Float& z, int precision);
    [[nodiscard]] Float carlsonRD(Float&& x, Float&& y, Float&& z, int precision);
    [[nodiscard]] Float carlsonRJ(const Float& x, const Float& y, const Float& z, const Float& p, int precision);
    [[nodiscard]] Float carlsonRJ(Float&& x, Float&& y, Float&& z, Float&& p, int precision);

    // Elliptic integrals — Legendre form (complete)
    [[nodiscard]] Float ellipticK(const Float& k, int precision);                                      // K(k) first kind
    [[nodiscard]] Float ellipticK(Float&& k, int precision);
    [[nodiscard]] Float ellipticE(const Float& k, int precision);                                      // E(k) second kind
    [[nodiscard]] Float ellipticE(Float&& k, int precision);
    [[nodiscard]] Float ellipticPi(const Float& n, const Float& k, int precision);                     // Π(n,k) third kind
    [[nodiscard]] Float ellipticPi(Float&& n, Float&& k, int precision);

    // Elliptic integrals — Legendre form (incomplete)
    [[nodiscard]] Float ellipticF(const Float& phi, const Float& k, int precision);                    // F(φ,k) first kind
    [[nodiscard]] Float ellipticF(Float&& phi, Float&& k, int precision);
    [[nodiscard]] Float ellipticE(const Float& phi, const Float& k, int precision);                    // E(φ,k) second kind
    [[nodiscard]] Float ellipticE(Float&& phi, Float&& k, int precision);
    [[nodiscard]] Float ellipticPi(const Float& n, const Float& phi, const Float& k, int precision);   // Π(n,φ,k) third kind
    [[nodiscard]] Float ellipticPi(Float&& n, Float&& phi, Float&& k, int precision);

    // Jacobi elliptic functions
    [[nodiscard]] Float jacobiSn(const Float& u, const Float& k, int precision);
    [[nodiscard]] Float jacobiSn(Float&& u, Float&& k, int precision);
    [[nodiscard]] Float jacobiCn(const Float& u, const Float& k, int precision);
    [[nodiscard]] Float jacobiCn(Float&& u, Float&& k, int precision);
    [[nodiscard]] Float jacobiDn(const Float& u, const Float& k, int precision);
    [[nodiscard]] Float jacobiDn(Float&& u, Float&& k, int precision);

    // Bernoulli number sequence B_{2k} (Akiyama-Tanigawa, result[k]=B_{2k}, result[0]=1).
    // Implemented in FloatMath.cpp. Shared by asymptotic expansions of the gamma/zeta family and constant computations.
    [[nodiscard]] std::vector<Float> computeBernoulliNumbers(int max_k, int precision);

    // Zeta function
    [[nodiscard]] Float zeta(const Float& s, int precision);
    [[nodiscard]] Float zeta(Float&& s, int precision);
    [[nodiscard]] Float hurwitzZeta(const Float& s, const Float& a, int precision);  // ζ(s,a)
    [[nodiscard]] Float hurwitzZeta(Float&& s, Float&& a, int precision);
    [[nodiscard]] Float dirichletEta(const Float& s, int precision);                 // η(s)
    [[nodiscard]] Float dirichletEta(Float&& s, int precision);

    // Exponential integral and dilogarithm
    [[nodiscard]] Float expint(const Float& x, int precision);   // Ei(x)
    [[nodiscard]] Float expint(Float&& x, int precision);
    [[nodiscard]] Float li(const Float& x, int precision);        // li(x) = Ei(ln(x))
    [[nodiscard]] Float li(Float&& x, int precision);
    [[nodiscard]] Float dilog(const Float& x, int precision);     // Li₂(x)
    [[nodiscard]] Float dilog(Float&& x, int precision);

    // Bessel functions
    [[nodiscard]] Float besselJ(int n, const Float& x, int precision);  // J_n(x)
    [[nodiscard]] Float besselJ(int n, Float&& x, int precision);
    [[nodiscard]] Float besselY(int n, const Float& x, int precision);  // Y_n(x)
    [[nodiscard]] Float besselY(int n, Float&& x, int precision);
    [[nodiscard]] Float besselI(int n, const Float& x, int precision);  // I_n(x) modified first kind
    [[nodiscard]] Float besselI(int n, Float&& x, int precision);
    [[nodiscard]] Float besselK(int n, const Float& x, int precision);  // K_n(x) modified second kind
    [[nodiscard]] Float besselK(int n, Float&& x, int precision);
    [[nodiscard]] Float sphericalBesselJ(int n, const Float& x, int precision);  // j_n(x) spherical Bessel, first kind
    [[nodiscard]] Float sphericalBesselJ(int n, Float&& x, int precision);
    [[nodiscard]] Float sphericalBesselY(int n, const Float& x, int precision);  // y_n(x) spherical Bessel, second kind
    [[nodiscard]] Float sphericalBesselY(int n, Float&& x, int precision);

    // Airy functions
    [[nodiscard]] Float airyAi(const Float& x, int precision);       // Ai(x)
    [[nodiscard]] Float airyAi(Float&& x, int precision);
    [[nodiscard]] Float airyBi(const Float& x, int precision);       // Bi(x)
    [[nodiscard]] Float airyBi(Float&& x, int precision);
    [[nodiscard]] Float airyAiPrime(const Float& x, int precision);  // Ai'(x)
    [[nodiscard]] Float airyAiPrime(Float&& x, int precision);
    [[nodiscard]] Float airyBiPrime(const Float& x, int precision);  // Bi'(x)
    [[nodiscard]] Float airyBiPrime(Float&& x, int precision);

    // Hypergeometric functions
    [[nodiscard]] Float confHyperg(const Float& a, const Float& b, const Float& z, int precision);                         // ₁F₁(a; b; z)
    [[nodiscard]] Float confHyperg(Float&& a, Float&& b, Float&& z, int precision);
    [[nodiscard]] Float hyperg(const Float& a, const Float& b, const Float& c, const Float& z, int precision);             // ₂F₁(a, b; c; z)
    [[nodiscard]] Float hyperg(Float&& a, Float&& b, Float&& c, Float&& z, int precision);
    [[nodiscard]] Float hyperg0F1(const Float& b, const Float& z, int precision);                                          // ₀F₁(; b; z)
    [[nodiscard]] Float hyperg0F1(Float&& b, Float&& z, int precision);

    // Legendre polynomials
    [[nodiscard]] Float legendreP(int n, const Float& x, int precision);              // P_n(x)
    [[nodiscard]] Float legendreP(int n, Float&& x, int precision);
    [[nodiscard]] Float assocLegendreP(int n, int m, const Float& x, int precision);  // P_n^m(x)
    [[nodiscard]] Float assocLegendreP(int n, int m, Float&& x, int precision);

    // Hermite polynomials (physicists' convention: H_n(x))
    [[nodiscard]] Float hermite(int n, const Float& x, int precision);
    [[nodiscard]] Float hermite(int n, Float&& x, int precision);

    // Laguerre polynomials
    [[nodiscard]] Float laguerre(int n, const Float& x, int precision);               // L_n(x)
    [[nodiscard]] Float laguerre(int n, Float&& x, int precision);
    [[nodiscard]] Float assocLaguerre(int n, int m, const Float& x, int precision);   // L_n^m(x)
    [[nodiscard]] Float assocLaguerre(int n, int m, Float&& x, int precision);

    // Lambert W function
    [[nodiscard]] Float lambertW0(const Float& x, int precision);    // W₀(x)
    [[nodiscard]] Float lambertW0(Float&& x, int precision);
    [[nodiscard]] Float lambertWm1(const Float& x, int precision);   // W₋₁(x)
    [[nodiscard]] Float lambertWm1(Float&& x, int precision);

    // Exponential and trigonometric integrals
    [[nodiscard]] Float expintN(int n, const Float& x, int precision);        // E_n(x)
    [[nodiscard]] Float expintN(int n, Float&& x, int precision);
    [[nodiscard]] Float sinIntegral(const Float& x, int precision);           // Si(x)
    [[nodiscard]] Float sinIntegral(Float&& x, int precision);
    [[nodiscard]] Float cosIntegral(const Float& x, int precision);           // Ci(x)
    [[nodiscard]] Float cosIntegral(Float&& x, int precision);

    // Random-number generation
    [[nodiscard]] Float randomFloat(int precision);                                    // [0, 1) uniform distribution
    [[nodiscard]] Float randomFloat(const Float& lo, const Float& hi, int precision);  // [lo, hi) uniform distribution
    [[nodiscard]] Float randomFloat(Float&& lo, Float&& hi, int precision);
    [[nodiscard]] Float normalRandom(int precision);                                   // N(0,1) normal distribution
    [[nodiscard]] Float exponentialRandom(int precision);                              // Exp(1) exponential distribution

    // ================================================================
    // Binary serialization
    // Format: [tag:1B][sign:1B][exponent:8B LE][eff:4B LE][req:4B LE][mantissa...]
    //   tag: 0x00=zero, 0x01=normal, 0x02=infinity, 0x03=nan
    //   mantissa: output of Int::exportBinary (only when tag==0x01)
    // ================================================================
    [[nodiscard]] std::vector<uint8_t> exportBinary(const Float& value);
    [[nodiscard]] Float importBinaryFloat(std::span<const uint8_t> data);

    // ================================================================
    // Precision-elided overloads
    // The precision is derived from the input's requestedBits(). For exact (INT_MAX)
    // values it falls back to defaultPrecision().
    // Intended to be invoked via ADL from inside Complex<Float>.
    // ================================================================

    [[nodiscard]] inline Float exp(const Float& x)  { return exp(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float exp(Float&& x)       { int p = Float::requestedPrecision(x); return exp(std::move(x), p); }
    [[nodiscard]] inline Float log(const Float& x)  { return log(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log(Float&& x)       { int p = Float::requestedPrecision(x); return log(std::move(x), p); }
    [[nodiscard]] inline Float logUi(unsigned long long n) { return logUi(n, Float::defaultPrecision()); }
    [[nodiscard]] inline Float sin(const Float& x)  { return sin(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float sin(Float&& x)       { int p = Float::requestedPrecision(x); return sin(std::move(x), p); }
    [[nodiscard]] inline Float cos(const Float& x)  { return cos(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float cos(Float&& x)       { int p = Float::requestedPrecision(x); return cos(std::move(x), p); }
    [[nodiscard]] inline Float tan(const Float& x)  { return tan(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float tan(Float&& x)       { int p = Float::requestedPrecision(x); return tan(std::move(x), p); }
    [[nodiscard]] inline Float sqr(const Float& x)   { return sqr(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float sqr(Float&& x)       { int p = Float::requestedPrecision(x); return sqr(std::move(x), p); }
    [[nodiscard]] inline Float sqrt(const Float& x) { return sqrt(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float sqrt(Float&& x)      { int p = Float::requestedPrecision(x); return sqrt(std::move(x), p); }
    [[nodiscard]] inline Float sinh(const Float& x) { return sinh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float sinh(Float&& x)      { int p = Float::requestedPrecision(x); return sinh(std::move(x), p); }
    [[nodiscard]] inline Float cosh(const Float& x) { return cosh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float cosh(Float&& x)      { int p = Float::requestedPrecision(x); return cosh(std::move(x), p); }
    [[nodiscard]] inline Float tanh(const Float& x) { return tanh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float tanh(Float&& x)      { int p = Float::requestedPrecision(x); return tanh(std::move(x), p); }
    [[nodiscard]] inline Float asin(const Float& x) { return asin(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float asin(Float&& x)      { int p = Float::requestedPrecision(x); return asin(std::move(x), p); }
    [[nodiscard]] inline Float acos(const Float& x) { return acos(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float acos(Float&& x)      { int p = Float::requestedPrecision(x); return acos(std::move(x), p); }
    [[nodiscard]] inline Float atan(const Float& x) { return atan(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float atan(Float&& x)      { int p = Float::requestedPrecision(x); return atan(std::move(x), p); }
    [[nodiscard]] inline Float atan2(const Float& y, const Float& x) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return atan2(y, x, prec);
    }
    [[nodiscard]] inline Float atan2(Float&& y, Float&& x) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return atan2(std::move(y), std::move(x), prec);
    }
    [[nodiscard]] inline Float asinPi(const Float& x) { return asinPi(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float asinPi(Float&& x)     { int p = Float::requestedPrecision(x); return asinPi(std::move(x), p); }
    [[nodiscard]] inline Float acosPi(const Float& x) { return acosPi(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float acosPi(Float&& x)     { int p = Float::requestedPrecision(x); return acosPi(std::move(x), p); }
    [[nodiscard]] inline Float atanPi(const Float& x) { return atanPi(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float atanPi(Float&& x)     { int p = Float::requestedPrecision(x); return atanPi(std::move(x), p); }
    [[nodiscard]] inline Float atan2Pi(const Float& y, const Float& x) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return atan2Pi(y, x, prec);
    }
    [[nodiscard]] inline Float atan2Pi(Float&& y, Float&& x) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return atan2Pi(std::move(y), std::move(x), prec);
    }
    [[nodiscard]] inline Float asinh(const Float& x) { return asinh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float asinh(Float&& x)      { int p = Float::requestedPrecision(x); return asinh(std::move(x), p); }
    [[nodiscard]] inline Float acosh(const Float& x) { return acosh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float acosh(Float&& x)      { int p = Float::requestedPrecision(x); return acosh(std::move(x), p); }
    [[nodiscard]] inline Float atanh(const Float& x) { return atanh(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float atanh(Float&& x)      { int p = Float::requestedPrecision(x); return atanh(std::move(x), p); }
    [[nodiscard]] inline Float pow(const Float& x, const Float& y) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return pow(x, y, prec);
    }
    [[nodiscard]] inline Float pow(Float&& x, Float&& y) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return pow(std::move(x), std::move(y), prec);
    }
    [[nodiscard]] inline Float pow(const Float& x, int n) {
        return pow(x, n, Float::requestedPrecision(x));
    }
    [[nodiscard]] inline Float pow(Float&& x, int n) {
        int p = Float::requestedPrecision(x); return pow(std::move(x), n, p);
    }
    [[nodiscard]] inline Float log2(const Float& x)  { return log2(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log2(Float&& x)       { int p = Float::requestedPrecision(x); return log2(std::move(x), p); }
    [[nodiscard]] inline Float log10(const Float& x) { return log10(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log10(Float&& x)      { int p = Float::requestedPrecision(x); return log10(std::move(x), p); }
    [[nodiscard]] inline Float log1p(const Float& x) { return log1p(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log1p(Float&& x)      { int p = Float::requestedPrecision(x); return log1p(std::move(x), p); }
    [[nodiscard]] inline Float exp2(const Float& x)  { return exp2(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float exp2(Float&& x)       { int p = Float::requestedPrecision(x); return exp2(std::move(x), p); }
    [[nodiscard]] inline Float exp10(const Float& x) { return exp10(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float exp10(Float&& x)      { int p = Float::requestedPrecision(x); return exp10(std::move(x), p); }
    [[nodiscard]] inline Float expm1(const Float& x)  { return expm1(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float expm1(Float&& x)       { int p = Float::requestedPrecision(x); return expm1(std::move(x), p); }
    [[nodiscard]] inline Float exp2m1(const Float& x)  { return exp2m1(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float exp2m1(Float&& x)       { int p = Float::requestedPrecision(x); return exp2m1(std::move(x), p); }
    [[nodiscard]] inline Float exp10m1(const Float& x) { return exp10m1(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float exp10m1(Float&& x)      { int p = Float::requestedPrecision(x); return exp10m1(std::move(x), p); }
    [[nodiscard]] inline Float log2p1(const Float& x)  { return log2p1(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log2p1(Float&& x)       { int p = Float::requestedPrecision(x); return log2p1(std::move(x), p); }
    [[nodiscard]] inline Float log10p1(const Float& x) { return log10p1(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float log10p1(Float&& x)      { int p = Float::requestedPrecision(x); return log10p1(std::move(x), p); }
    [[nodiscard]] inline Float compound(const Float& x, int n) { return compound(x, n, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float compound(Float&& x, int n)      { int p = Float::requestedPrecision(x); return compound(std::move(x), n, p); }
    [[nodiscard]] inline Float cbrt(const Float& x)  { return cbrt(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float cbrt(Float&& x)       { int p = Float::requestedPrecision(x); return cbrt(std::move(x), p); }
    // Single-argument (precision-omitted) forms of special functions + std-compatible aliases.
    // tgamma = gamma (Γ), lgamma = lnGamma (ln|Γ|). So that generic_xxx (generic_math.hpp)
    // and autodiff's Dual<Float> can resolve them via ADL.
    [[nodiscard]] inline Float erf(const Float& x)    { return erf(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float erfc(const Float& x)   { return erfc(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float tgamma(const Float& x, int precision) { return gamma(x, precision); }
    [[nodiscard]] inline Float tgamma(const Float& x) { return gamma(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float lgamma(const Float& x, int precision) { return lnGamma(x, precision); }
    [[nodiscard]] inline Float lgamma(const Float& x) { return lnGamma(x, Float::requestedPrecision(x)); }
    // Single-argument form of ψ(x) = Γ'(x)/Γ(x) (digamma). Used so that autodiff's Dual<Float>/RVar<Float>
    // can resolve the exact derivative ψ(x)·Γ(x) of tgamma/lgamma via ADL.
    [[nodiscard]] inline Float digamma(const Float& x) { return digamma(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float nthRoot(const Float& x, int n) { return nthRoot(x, n, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float nthRoot(Float&& x, int n)      { int p = Float::requestedPrecision(x); return nthRoot(std::move(x), n, p); }
    [[nodiscard]] inline Float recSqrt(const Float& x) { return recSqrt(x, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float recSqrt(Float&& x)      { int p = Float::requestedPrecision(x); return recSqrt(std::move(x), p); }
    [[nodiscard]] inline Float hypot(const Float& x, const Float& y) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return hypot(x, y, prec);
    }
    [[nodiscard]] inline Float hypot(Float&& x, Float&& y) {
        int ry = y.requestedBits(), rx = x.requestedBits();
        int req = (ry >= INT_MAX && rx >= INT_MAX) ? INT_MAX
                : (ry >= INT_MAX) ? rx : (rx >= INT_MAX) ? ry : std::max(ry, rx);
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return hypot(std::move(x), std::move(y), prec);
    }
    [[nodiscard]] inline Float fma(const Float& a, const Float& b, const Float& c) {
        int ra = a.requestedBits(), rb = b.requestedBits(), rc = c.requestedBits();
        int req = INT_MAX;
        if (ra < INT_MAX) req = ra;
        if (rb < INT_MAX) req = (req < INT_MAX) ? std::max(req, rb) : rb;
        if (rc < INT_MAX) req = (req < INT_MAX) ? std::max(req, rc) : rc;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fma(a, b, c, prec);
    }
    [[nodiscard]] inline Float fma(Float&& a, Float&& b, Float&& c) {
        int ra = a.requestedBits(), rb = b.requestedBits(), rc = c.requestedBits();
        int req = INT_MAX;
        if (ra < INT_MAX) req = ra;
        if (rb < INT_MAX) req = (req < INT_MAX) ? std::max(req, rb) : rb;
        if (rc < INT_MAX) req = (req < INT_MAX) ? std::max(req, rc) : rc;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fma(std::move(a), std::move(b), std::move(c), prec);
    }
    [[nodiscard]] inline Float fms(const Float& a, const Float& b, const Float& c) {
        int ra = a.requestedBits(), rb = b.requestedBits(), rc = c.requestedBits();
        int req = INT_MAX;
        if (ra < INT_MAX) req = ra;
        if (rb < INT_MAX) req = (req < INT_MAX) ? std::max(req, rb) : rb;
        if (rc < INT_MAX) req = (req < INT_MAX) ? std::max(req, rc) : rc;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fms(a, b, c, prec);
    }
    [[nodiscard]] inline Float fms(Float&& a, Float&& b, Float&& c) {
        int ra = a.requestedBits(), rb = b.requestedBits(), rc = c.requestedBits();
        int req = INT_MAX;
        if (ra < INT_MAX) req = ra;
        if (rb < INT_MAX) req = (req < INT_MAX) ? std::max(req, rb) : rb;
        if (rc < INT_MAX) req = (req < INT_MAX) ? std::max(req, rc) : rc;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fms(std::move(a), std::move(b), std::move(c), prec);
    }
    [[nodiscard]] inline Float fmma(const Float& a, const Float& b, const Float& c, const Float& d) {
        int ra = a.requestedBits(), rb = b.requestedBits();
        int rc = c.requestedBits(), rd = d.requestedBits();
        int req = INT_MAX;
        for (int v : {ra, rb, rc, rd})
            if (v < INT_MAX) req = (req < INT_MAX) ? std::max(req, v) : v;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fmma(a, b, c, d, prec);
    }
    [[nodiscard]] inline Float fmms(const Float& a, const Float& b, const Float& c, const Float& d) {
        int ra = a.requestedBits(), rb = b.requestedBits();
        int rc = c.requestedBits(), rd = d.requestedBits();
        int req = INT_MAX;
        for (int v : {ra, rb, rc, rd})
            if (v < INT_MAX) req = (req < INT_MAX) ? std::max(req, v) : v;
        int prec = (req >= INT_MAX) ? Float::contextPrecision() : Float::bitsToPrecision(req);
        return fmms(a, b, c, d, prec);
    }
    inline void sinCos(const Float& x, Float& sin_result, Float& cos_result) {
        sinCos(x, sin_result, cos_result, Float::requestedPrecision(x));
    }
    inline void sinCos(Float&& x, Float& sin_result, Float& cos_result) {
        int p = Float::requestedPrecision(x); sinCos(std::move(x), sin_result, cos_result, p);
    }
    [[nodiscard]] inline Float sinu(const Float& x, int u) { return sinu(x, u, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float sinu(Float&& x, int u)      { int p = Float::requestedPrecision(x); return sinu(std::move(x), u, p); }
    [[nodiscard]] inline Float cosu(const Float& x, int u) { return cosu(x, u, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float cosu(Float&& x, int u)      { int p = Float::requestedPrecision(x); return cosu(std::move(x), u, p); }
    [[nodiscard]] inline Float tanu(const Float& x, int u) { return tanu(x, u, Float::requestedPrecision(x)); }
    [[nodiscard]] inline Float tanu(Float&& x, int u)      { int p = Float::requestedPrecision(x); return tanu(std::move(x), u, p); }

} // namespace sangi

// Add Float's mathematical functions and numeric_limits to the std namespace.
// This bridges generic algorithms in LinAlg / statistics / ODE / signal processing etc.
// that call `std::xxx(T)` directly so that Matrix<Float> works. Delegates to sangi::xxx.
namespace std {
    // Norm and absolute value
    inline sangi::Float abs(const sangi::Float& x) { return sangi::abs(x); }
    inline sangi::Float sqrt(const sangi::Float& x) { return sangi::sqrt(x); }
    inline sangi::Float cbrt(const sangi::Float& x) { return sangi::cbrt(x); }
    inline sangi::Float hypot(const sangi::Float& x, const sangi::Float& y) {
        return sangi::hypot(x, y);
    }

    // Exponential and logarithm
    inline sangi::Float exp(const sangi::Float& x)   { return sangi::exp(x); }
    inline sangi::Float log(const sangi::Float& x)   { return sangi::log(x); }
    inline sangi::Float log2(const sangi::Float& x)  { return sangi::log2(x); }
    inline sangi::Float log10(const sangi::Float& x) { return sangi::log10(x); }
    inline sangi::Float pow(const sangi::Float& x, const sangi::Float& y) {
        return sangi::pow(x, y);
    }
    inline sangi::Float pow(const sangi::Float& x, int n) {
        return sangi::pow(x, n);
    }
    // pow(Float, double): promote the double exponent to Float before calling
    inline sangi::Float pow(const sangi::Float& x, double y) {
        return sangi::pow(x, sangi::Float(y));
    }

    // Trigonometric functions
    inline sangi::Float sin(const sangi::Float& x)  { return sangi::sin(x); }
    inline sangi::Float cos(const sangi::Float& x)  { return sangi::cos(x); }
    inline sangi::Float tan(const sangi::Float& x)  { return sangi::tan(x); }
    inline sangi::Float asin(const sangi::Float& x) { return sangi::asin(x); }
    inline sangi::Float acos(const sangi::Float& x) { return sangi::acos(x); }
    inline sangi::Float atan(const sangi::Float& x) { return sangi::atan(x); }
    inline sangi::Float atan2(const sangi::Float& y, const sangi::Float& x) {
        return sangi::atan2(y, x);
    }

    // Hyperbolic functions
    inline sangi::Float sinh(const sangi::Float& x) { return sangi::sinh(x); }
    inline sangi::Float cosh(const sangi::Float& x) { return sangi::cosh(x); }
    inline sangi::Float tanh(const sangi::Float& x) { return sangi::tanh(x); }

    // Rounding
    inline sangi::Float floor(const sangi::Float& x)     { return sangi::floor(x); }
    inline sangi::Float ceil(const sangi::Float& x)      { return sangi::ceil(x); }
    inline sangi::Float round(const sangi::Float& x)     { return sangi::round(x); }
    inline sangi::Float trunc(const sangi::Float& x)     { return sangi::trunc(x); }
    inline sangi::Float nearbyint(const sangi::Float& x) { return sangi::nearbyint(x); }

    // Remainder
    inline sangi::Float fmod(const sangi::Float& x, const sangi::Float& y) {
        return sangi::fmod(x, y);
    }

    // Numeric-state predicates (member function → free-function bridge)
    inline bool isnan(const sangi::Float& x)    { return x.isNaN(); }
    inline bool isinf(const sangi::Float& x)    { return x.isInfinity(); }
    inline bool isfinite(const sangi::Float& x) { return !x.isInfinity() && !x.isNaN(); }

    // Sign manipulation: copysign(x, y) returns |x| with the sign of y
    inline sangi::Float copysign(const sangi::Float& x, const sangi::Float& y) {
        sangi::Float ax = sangi::abs(x);
        // Negative sign when y < 0; otherwise (including y == 0) positive
        return (y < sangi::Float(0)) ? -ax : ax;
    }
    inline bool signbit(const sangi::Float& x) { return x < sangi::Float(0); }

    // ldexp: x * 2^n
    inline sangi::Float ldexp(const sangi::Float& x, int n) { return sangi::ldexp(x, n); }

    template<>
    struct numeric_limits<sangi::Float> {
        static constexpr bool is_specialized = true;
        static sangi::Float epsilon() {
            // Machine epsilon based on the current default precision
            int bits = sangi::Float::defaultPrecision();
            // bits is a decimal digit count; convert to bits via bits * log2(10) ≈ bits * 3.3219
            int prec_bits = static_cast<int>(bits * 3.3219) + 1;
            return sangi::ldexp(sangi::Float(1), -prec_bits);
        }
        static sangi::Float min() { return sangi::ldexp(sangi::Float(1), -1000000); }
        static sangi::Float max() { return sangi::ldexp(sangi::Float(1), 1000000); }
        static sangi::Float lowest() { return -max(); }
        static constexpr bool is_signed = true;
        static constexpr bool is_integer = false;
        static constexpr bool is_exact = false;
        static constexpr bool has_infinity = false;
        static constexpr bool has_quiet_NaN = true;
        static sangi::Float quiet_NaN() { return sangi::Float::nan(); }
        static sangi::Float infinity() { return sangi::Float::positiveInfinity(); }
    };
}

// --- std::formatter specialization ---
// <format> is available on GCC 13+, Clang 17+, MSVC 19.29+.
// Use __cpp_lib_format (C++20 feature test macro) to check for availability.
#if __has_include(<format>)
#include <format>
#if defined(__cpp_lib_format) && __cpp_lib_format >= 202106L

template<>
struct std::formatter<sangi::Float> {
    int precision_ = -1;
    char type_ = 'g';  // 'g'=general, 'f'=fixed, 'e'=scientific

    constexpr auto parse(std::format_parse_context& ctx) {
        auto it = ctx.begin();
        auto end = ctx.end();
        if (it != end && *it == '.') {
            ++it;
            precision_ = 0;
            while (it != end && *it >= '0' && *it <= '9') {
                precision_ = precision_ * 10 + (*it - '0');
                ++it;
            }
        }
        if (it != end && *it != '}') {
            if (*it == 'f' || *it == 'e' || *it == 'E' || *it == 'g' || *it == 'G') {
                type_ = *it;
                ++it;
            } else {
                throw std::format_error("invalid format spec for Float");
            }
        }
        return it;
    }

    auto format(const sangi::Float& val, std::format_context& ctx) const {
        std::string result;
        switch (type_) {
            case 'f': result = val.toDecimalString(precision_); break;
            case 'e': case 'E': result = val.toScientificString(precision_); break;
            default: result = val.toString(precision_); break;
        }
        return std::format_to(ctx.out(), "{}", result);
    }
};

#endif // __cpp_lib_format
#endif // __has_include(<format>)

#endif // SANGI_FLOAT_HPP