nfloat.h

Initialises ctx to the ring of floating-point numbers with precision prec bits rounded up to a whole number of words (for example, prec = 53 gives 64-bit precision), with the given flags (a bitwise or of the flags below, or 0). Returns GR_UNABLE without initialising ctx if prec is not positive or exceeds the maximum NFLOAT_MAX_LIMBS * FLINT_BITS, and GR_SUCCESS otherwise. The context needs no memory and is cleared with gr_ctx_clear.

Set and get the function precision of ctx: the precision in bits to which the elementary functions (exponentials, logarithms, trigonometric and hyperbolic functions, their inverses and powers) are computed. It defaults to the full precision of the context; a value larger than that is clamped to it, and a value less than 1 makes nfloat_ctx_set_func_prec return GR_UNABLE. A smaller function precision makes these functions cheaper (for example, at 64 bits with a function precision of 50 bits, most functions evaluate their kernels with one word instead of two); the arithmetic is unaffected.

Writes a description of ctx to out, such as Floating-point numbers with prec = 128 (nfloat).

The precision of a context cannot be changed: nfloat_ctx_set_real_prec does nothing and returns GR_UNABLE. nfloat_ctx_get_real_prec sets res to the precision in bits and returns GR_SUCCESS.

Initialise res to zero, and clear it. Since elements own no memory, clearing does nothing.

Swap two elements word by word, and copy x to res.

Set res to 0, 1 or -1.

Set res to +\infty, -\infty or NaN if the context allows it. Otherwise the infinities return GR_DOMAIN and NaN returns GR_UNABLE, leaving res unchanged.

Sets res to \pi rounded to the precision of the context.

Set res to an integer, a rational number or a double, rounded to the precision of the context (toward zero, or in the direction of the context's directed rounding).

Convert from and to arf_t; the conversion to arf_t is exact. These are declared only if arf.h is included before nfloat.h.

Sets res to x if x is an integer and returns GR_SUCCESS; returns GR_DOMAIN if x is not an integer and GR_UNABLE if it is infinite or NaN.

Sets mant and expo such that x \approx \mathit{mant} \cdot 2^{\mathit{expo}} with |\mathit{mant}| \in [1/2, 1) (or both zero for x = 0), from the leading word of the mantissa, so with a relative error of about 2^{-53}. This avoids the overflow of a direct conversion to double. Returns GR_UNABLE for infinities and NaN.

Set res to (-1)^{\mathit{xsgnbit}} \cdot f \cdot 2^{\mathit{exp}}, where f = \{x, \mathit{xn}\} / 2^{\mathit{xn} \cdot \mathtt{FLINT\_BITS}} is the fraction in [0, 1) given by the xn words of x, rounded to the precision of the context. The version with an underscore requires that the top word x[xn - 1] is nonzero; the other removes leading zero words and handles f = 0.

Inline versions of nfloat_set_mpn_2exp for a two-, three- or four-word fraction (most significant word first in the argument list) and a context of exactly one word (nfloat_1_) or two words (nfloat_2_). The fraction is normalised and truncated (not rounded) to the precision; zero, underflow and overflow are handled as by the other functions.

Parses a decimal number with arb_set_str and rounds it to nearest (this ignores directed rounding); other expressions are parsed with the generic string parser of the gr interface.

Sets res to the element x of the ring x_ctx, rounded. Elements of other nfloat contexts, fmpz, fmpq, arf_t and arb_t (whose midpoint is used) are converted directly; other rings are converted through arf_t with gr_set_other, whose status is returned (for example GR_DOMAIN for a non-real number).

Write x in decimal as the arf_t ring of the same precision does, and set res to a random element (from arf_randtest) for testing.

Returns T_TRUE if x and y have the same representation, which for normal values means equal values, and T_FALSE otherwise. Unlike IEEE comparison, a NaN equals a NaN.

Return whether x is 0, 1 or -1; T_UNKNOWN for NaN.

Set res to -1, 0 or 1 according to the comparison of x and y, or of |x| and |y|, and return GR_SUCCESS; return GR_UNABLE if an operand is NaN.

Set res to the sign of x (-1, 0 or 1, and NaN for NaN), and to the imaginary part of x, which is zero.

Set res to -x or |x|, exactly.

Sum, difference, product and square. Without directed rounding the error is less than 2 ulp for the sum and difference (the smaller operand is truncated when it is aligned) and at most 1 ulp for the products (truncation, except that the approximate high product used from three words on can rarely give one ulp less in magnitude); with directed rounding the results are correctly rounded.

Set res to \mathit{res} + xy or \mathit{res} - xy, with two roundings, so the result can be far from the exact value (measured in ulp of the result) when the terms cancel.

Sets res to x \cdot 2^y, exactly unless the result overflows or underflows.

Inverse and quotients. Without directed rounding, the result is the quotient truncated toward zero or, for divisors of three or more nonzero words, possibly one ulp larger in magnitude (an error of less than 1 ulp); div_ui and div_si are always truncated. With directed rounding, the result is the correctly rounded floor or ceiling of the exact quotient. Division by zero gives NaN, so it returns GR_UNABLE unless NaN is allowed.

Square root and reciprocal square root. The square root is correctly rounded in every mode (truncated toward zero by default); the reciprocal square root has an error of less than 1 + 2^{-37} ulp by default and is correctly rounded with directed rounding. Negative arguments give NaN.

Sets res to \sqrt{x^2 + y^2}, correctly rounded in every mode (truncated by default), without intermediate overflow or underflow. This also implements nfloat_complex_abs.

Round x to an integer toward -\infty, toward +\infty, toward zero, or to the nearest integer (ties to even). The results are exact.

Exponentials e^x, e^x - 1, 2^x and logarithms \log x, \log(1 + x), \log_2 x. The exponentials reduce |x| modulo \log 2 in fixed point (with one integral word); expm1 and log1p have relative accuracy near zero, and exp2 and log2 are exact for integers and powers of two respectively.

Sine, cosine and tangent of x and of \pi x; the sin_cos versions compute both. The argument is reduced modulo \pi/2 in fixed point, repeated with more bits when the reduced argument is close to a zero of the result, so that the results have relative accuracy; the functions of \pi x reduce x modulo 2 exactly.

The arctangent, and the two-argument arctangent \operatorname{atan2}(y, x) \in (-\pi, \pi]. For |x| > 1, atan uses \pi/2 - \operatorname{atan}(1/|x|); atan2 takes the arctangent of the ratio of the smaller and the larger magnitude and adds the octant. \operatorname{atan2}(0, 0) = 0.

Hyperbolic functions, from E = e^{|x|} and 1/E, with the difference E - 1/E computed in fixed point for |x| < 1.

Sets res to x^y = \exp(y \log|x|), with the sign (-1)^y for x < 0 and integer y (NaN for x < 0 and noninteger y). Simple cases (y = \pm 1, 2, \pm 1/2, |x| = 1, integer powers of powers of two) are handled directly.

The less common elementary functions. After the exact special values (zeros, poles, |x| = 1 for the inverse functions, multiples of 1/2 for the functions of \pi x), they are evaluated as short compositions of the functions above in an internal context with a guard word (none when the function precision is sufficiently below the precision), using formulas without cancellation, such as \operatorname{asin}(x) = \operatorname{atan2}(x, \sqrt{(1 - x)(1 + x)}) and \operatorname{asinh}(x) = \operatorname{log1p}(|x| + x^2/(1 + \sqrt{1 + x^2})), and a rigorous bound for the composed error, which is added to obtain valid bounds with directed rounding.

The gamma function and the Riemann zeta function, computed by converting to arf_t and calling the generic ball implementation (gr_gamma, gr_zeta over arb_t), then rounding.

Initialise a vector of len elements to zero, clear it (no-op), copy and zero.

Elementwise sum, difference and product, product with the scalar y, and \mathit{res}_i \pm x_i y, with the same results as the scalar functions.

Elementwise quotients and quotients by a scalar, with the same results as the scalar division, except that _nfloat_vec_div_scalar for three or more words and a divisor of two or more nonzero words may differ by up to 1 ulp without directed rounding.

Set res to s \pm \sum_{i} x_i y_i (or with y reversed), where s is initial (zero if NULL) and the sign is negative if subtract is nonzero. The sum is accumulated in fixed point with one guard word after a pass that finds the largest exponent, so the result has an error of a few ulps relative to the largest term rather than one rounding per term; with directed rounding it is a valid bound. See Dot products.

Sets C = AB, choosing between classical multiplication with the dot product, nfloat_mat_mul_fixed and nfloat_mat_mul_block according to the dimensions and the precision (tuned cutoffs). Contexts with infinities, NaN or directed rounding, and matrices with a dimension of at most 2, always use classical multiplication. The methods are described in Matrix multiplication.

Sets C = AB by converting A and B to fixed-point matrices with common exponents and multiplying them with _nfixed_mat_mul. The fixed-point precision is the target precision plus the spread of the exponents in A and in B plus a few guard bits; returns GR_UNABLE if this exceeds the target by max_extra_prec bits or more, if a matrix has zero entries (but is not zero), or if infinities or NaN are enabled.

Sets C = AB by splitting the inner dimension into blocks in which the entries of each row of A and each column of B have similar exponents, scaling the blocks to integer matrices and multiplying them exactly with fmpz_mat_mul (multimodular or FFT methods), as arb_mat_mul_block [Joh2019dot]. Blocks shorter than min_block_size are handled with dot products.

Triangular solving, LU factorisation and LQ factorisation, choosing between the classical (Gram–Schmidt for LQ) and the recursive versions of the generic algorithms (gr_mat_nonsingular_solve_tril_classical, gr_mat_lu_recursive, gr_mat_lq_recursive, ...) by precision-dependent cutoffs, so that large problems run on fast matrix multiplication.

Initialises ctx to the ring of complex numbers whose parts have precision prec bits rounded up to a whole number of words, with the flags of nfloat_ctx_init. Returns GR_UNABLE if the precision is out of range.

Memory management and constants, as for the real numbers; i sets res to the imaginary unit.

Conversions from and to acf_t (exact for get) and acb_t (the midpoint for set). They are declared only if acf.h or acb.h is included before nfloat.h.

Conversions from real values (with zero imaginary part) and from other rings, output, and random elements.

Predicates and comparisons. cmp compares real numbers and returns GR_DOMAIN if an imaginary part is nonzero; cmpabs compares absolute values.

Conjugate, real and imaginary parts (as complex numbers), negation, sum, difference, product, square, multiplication by 2^y and absolute value (a real number, correctly rounded as by nfloat_hypot).

Inverse, quotient, square root and reciprocal square root (principal branches), with errors of a few ulps relative to the magnitude of the result. Purely real or imaginary operands are handled with the real functions.

Complex elementary functions, as compositions of the real functions, with the branch cuts of acb_t. Arguments in the real domain of a function use the real function.

Vector operations on complex numbers, as for the real vectors; the dot products accumulate the real and imaginary parts in fixed point.

Complex matrix multiplication, as for the real matrices. The fixed-point and block versions use three real products (the Karatsuba formula (AC - BD) + ((A + B)(C + D) - AC - BD)i); nfloat_complex_mat_mul_reorder separates the real and imaginary parts into real matrices and multiplies them with nfloat_mat_mul (with fewer products when A or B is real). nfloat_complex_mat_mul uses the reordering for real factors and for parts that are not well scaled.

Triangular solving and LU factorisation, choosing between the classical and recursive generic algorithms by precision-dependent cutoffs.

Prints the fixed-point number x times 2^{\mathit{exp}} in decimal.

Elementwise sum and difference of vectors of len fixed-point numbers (exact).

Dot products of fixed-point vectors with 2 to 8 words, accumulating truncated high products in separate positive and negative sums. The strides are the offsets in words between consecutive entries and may be negative.

Set the m \times p matrix C to the product of the m \times n matrix A and the n \times p matrix B (dense, row-major). The classical version uses the dot products; the _precise version computes with one extra word and is intended for testing. Waksman's algorithm [Waksman1970] halves the number of multiplications at the cost of more additions, which pays off from about 9 words. Strassen's algorithm [Str1969] recurses down to cutoff (a negative value selects the default). _nfixed_mat_mul chooses: Strassen for dimensions above about 26 to 57 (depending on the precision and parity), otherwise Waksman or classical.

For an m \times n \times p product at precision nlimbs with entries of A and B bounded by A and B, set bound to a bound for all intermediate values of the corresponding algorithm, which must be less than 1 for the computation to be valid, and error to a bound for the output error in ulps. _nfixed_mat_mul_bound follows the algorithm choice of _nfixed_mat_mul. The bounds are nondecreasing in A and B. For complex matrices, the entries are bounded by A + Bi and C + Di.

Set res to the result of an underflow (zero, returning GR_SUCCESS only with NFLOAT_ALLOW_UNDERFLOW) or an overflow (a signed infinity with NFLOAT_ALLOW_INF, otherwise zero and GR_UNABLE).

Return -1, 0 or 1 according to the comparison of x and y, or of |x| and |y|; neither may be NaN.

The fixed-size kernels of addition and subtraction of magnitudes: set res to (-1)^{\mathit{xsgnbit}} (|x| \pm |y|), where x has mantissa x and exponent xexp, and y has mantissa y and exponent xexp - delta with \mathit{delta} \ge 0 (so |x| \ge |y| in exponent, and for subtraction |x| \ge |y|). The mantissas have 1, 2, 3, 4 or nlimbs words.

Set res1 and res2 to the real and imaginary parts of (a + bi)^2 or (a + bi)(c + di) with real nfloat parts: naive with real products and sums, standard with fused dot products of two terms, and karatsuba with three real products instead of four. _nfloat_complex_sqr uses the naive version when infinities or NaN are enabled, real squarings when a part is zero, and otherwise the standard version below 20 words and the Karatsuba version from 20 words.