// 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 #include #include #include #include #include #include #include // C++20 three-way comparison operator #include #include // INT_MAX #include // std::is_integral_v #include // std::pair #include // std::size_t #include // 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; 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 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(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 && !std::is_same_v, int> = 0> Float& operator+=(T rhs) { if constexpr (std::is_signed_v) return *this += Float(static_cast(rhs)); else return *this += Float(static_cast(rhs)); } template && !std::is_same_v, int> = 0> Float& operator-=(T rhs) { if constexpr (std::is_signed_v) return *this -= Float(static_cast(rhs)); else return *this -= Float(static_cast(rhs)); } template && !std::is_same_v, int> = 0> Float& operator*=(T rhs) { if constexpr (std::is_signed_v) return *this *= Float(static_cast(rhs)); else return *this *= Float(static_cast(rhs)); } template && !std::is_same_v, int> = 0> Float& operator/=(T rhs) { if constexpr (std::is_signed_v) return *this /= Float(static_cast(rhs)); else return *this /= Float(static_cast(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. */ [[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(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(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 { 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 && !std::is_same_v, int> = 0> [[nodiscard]] Float operator*(const Float& lhs, T rhs) { if constexpr (std::is_signed_v) return mulScalarF(lhs, static_cast(rhs)); else return mulScalarF(lhs, static_cast(rhs)); } template && !std::is_same_v, int> = 0> [[nodiscard]] Float operator*(Float&& lhs, T rhs) { if constexpr (std::is_signed_v) return mulScalarF(std::move(lhs), static_cast(rhs)); else return mulScalarF(std::move(lhs), static_cast(rhs)); } template && !std::is_same_v, int> = 0> [[nodiscard]] Float operator*(T lhs, const Float& rhs) { if constexpr (std::is_signed_v) return mulScalarF(rhs, static_cast(lhs)); else return mulScalarF(rhs, static_cast(lhs)); } template && !std::is_same_v, int> = 0> [[nodiscard]] Float operator/(const Float& lhs, T rhs) { if constexpr (std::is_signed_v) return divScalarF(lhs, static_cast(rhs)); else return divScalarF(lhs, static_cast(rhs)); } template && !std::is_same_v, int> = 0> [[nodiscard]] Float operator/(Float&& lhs, T rhs) { if constexpr (std::is_signed_v) return divScalarF(std::move(lhs), static_cast(rhs)); else return divScalarF(std::move(lhs), static_cast(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); // -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 remquo(const Float& x, const Float& y); [[nodiscard]] std::pair 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 values, int precision); [[nodiscard]] Float dot(std::span a, std::span 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 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 exportBinary(const Float& value); [[nodiscard]] Float importBinaryFloat(std::span 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. // ================================================================ [[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 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/RVar // 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 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 { 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(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 --- // 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() #include #if defined(__cpp_lib_format) && __cpp_lib_format >= 202106L template<> struct std::formatter { 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() #endif // SANGI_FLOAT_HPP