mp_real.h

Return whether |e| \le MP_REAL_EXP_MAX, respectively whether the exponent of x is safe.

Return the working precision n in limbs for a target relative accuracy \varepsilon = B^{-\mathit{prec\_limbs}} resp. 2^{-\mathit{prec\_bits}}: prec_limbs resp. \lceil \mathit{prec\_bits} / \mathtt{FLINT\_BITS} \rceil plus MP_REAL_PREC_GUARD.

Initialise x to exact zero, respectively free its memory.

Ensure that x has room for len mantissa limbs. The inline function calls _mp_real_grow when it has to reallocate; the contents of the mantissa are preserved.

Swap x and y efficiently.

Set the ball to exact zero, to a copy of x, or to the exact integer c resp. f.

Set res to x rounded to its top n \ge 1 limbs (a copy if x has at most n limbs). res may alias x.

Set x exactly to (p, \mathit{len})\, 2^e for the unsigned integer (p, \mathit{len}), which may have zero top limbs. The bit part of e costs one limb shift.

Set res to -x; in place this flips the sign bit.

Set res to x \cdot 2^e exactly (an exponent change and at most one limb shift of the mantissa).

Return whether x is exactly zero (zero midpoint and zero radius).

Return an e with |t| < 2^e for every t in the ball x.

Return an e such that the radius of x is below 2^e times the absolute value of its midpoint: -WORD_MAX / 2 for an exact x and WORD_MAX / 2 for a zero midpoint with a nonzero radius.

Add 2^e to the radius of x, respectively enlarge x to contain x(1 + t) for every |t| \le 2^e.

Add v ulps of the mantissa's last limb, respectively v B^{\mathit{anchor}}, to the radius of x, for v \ge 0.

Set res to a + b resp. a - b.

Set res to x + yc resp. x - yc for a word c.

Set res to ab resp. ac.

Set r_r + i r_i to (a_r + i a_i)(b_r + i b_i), to about n limbs relative to the larger part of the result. The outputs may alias the inputs but not each other.

Set res to a - bc, given that |a - bc| < B^E. Only the window of the product that reaches the result is computed, and the result is wrong if the bound does not hold.

Set res to a / b resp. a / c. The divisor ball must have a relative radius below 2^{-30}; mp_real_div throws for a wider or zero divisor, and mp_real_div_ui for c = 0.

Set res to \sqrt{x} resp. 1/\sqrt{x} for a positive ball x of relative radius below 2^{-30}. A wider ball, or a ball with a negative or zero midpoint, throws, except that mp_real_sqrt of an exact zero gives zero. res may alias x.

Set res to 1/\sqrt{c} for a word c \ge 1, by _mp_real_rsqrt_ui_newton; throws for c = 0.

Set res to x^{1/k} resp. x^{-1/k} for x > 0 and 1 \le k < MP_REAL_ROOT_K_MAX. They throw for k = 0, for x not positive (the third function for k < 2 as well). The third function takes the order r of the iteration (2 \le r \le 15, or 0 for the default) and recip = 1 for the reciprocal root.

Set res to the arithmetic-geometric mean \operatorname{agm}(x, y) of x, y \ge 0 to about n limbs; throws for a negative midpoint. The second function takes the order m of the final series (at most 16, 10 on 32-bit machines; 0 for the default).

Set res to the ball x without loss: the midpoint exactly and the radius rounded up.

If the ball x contains exactly one integer, set res to that integer (exactly) and return 1; otherwise return 0 and leave res unchanged. res may alias x, in which case no reallocation takes place.

Convert from and to the double expansions of dfloat: arrays of n \le 8 doubles (the midpoint components) and a double radius. mp_real_set_dfloat sets res to the ball \sum_i x_i \pm \mathit{rad} and throws for nonfinite components or radius. mp_real_get_dfloat rounds x greedily to n doubles and sets rad (if not NULL) to a bound for the error and the radius of x.

Print the fields size, exp and err of x and its value to about twenty digits.

Write a ball known to lie in [0, 1) as the fixed-point number (y, n), truncating, and set err (if not NULL) to a bound for the error in ulps B^{-n}, rounded up; UWORD_MAX means that no bound below B - 1 is known. A negative midpoint within the radius is clamped to zero and absorbed into the bound.

Write y = \lfloor x B^n \rfloor for a ball known to lie in [0, 1) and return 1 if the radius determines that floor; otherwise return 0, so that the caller can retry at a higher precision.

Set res to a ball containing c = \pi/4, \log 2, Euler's constant \gamma, e, \log 10, Catalan's constant G, \zeta(3), \zeta(5), \Gamma(1/3), \Gamma(1/4), 2/\pi resp. 1/\log 2 - 1, accurate to about n limbs. With cache = 0 the value is computed from scratch; with cache = 1 it is read from (and if necessary added to) a process-wide cache of floors, and the ball is the cached floor with a radius of one ulp.

Set res to exactly \lfloor c B^n \rfloor for the constants c of mp_real_const_pi4 and its siblings, and err (if not NULL) to 1. The output has n limbs for the constants below 1 (\pi/4, \log 2, \gamma, G, 2/\pi, 1/\log 2 - 1) and n + 1 limbs, with a units limb, for the others. cache is as for the ball versions.

Free the shared caches of the constants. No other thread may be using them. flint_cleanup_master also frees them.

Set res to the constant computed by f with cache = 0, as an arb ball of prec bits. This is the entry point of the arb constants (arb_const_pi, arb_const_log2, arb_const_euler, ...), which keep their own caches.

Euler's constant with a forced choice of the logarithms used in the Brent–McMillan formula: set = 0, 1, 2, 3 selects the primes \{2\}, \{2, 3\}, \{2, 3, 5\}, \{2, 3, 5, 7\}, and -1 the automatic choice. For tuning.

Set res1 and res2 to balls containing \sin y and \cos y for every y in x. Either output may be NULL. The sine keeps its relative accuracy for small |x|; near the other zeros the loss of significance shows in the radius. Arguments with |x| \ge 2^{\max(65536, 4\,\mathit{prec})} give [0 \pm 1].

Set res to a ball containing \tan y for every y in x and return 1, or return 0 with res set to zero when the ball may contain a pole (or |x| \ge 2^{\max(65536, 4\,\mathit{prec})}). The result keeps its relative accuracy near the zeros and poles of the tangent.

The same as mp_real_sin_cos_bits and mp_real_tan_bits for \sin(\pi y), \cos(\pi y) and \tan(\pi y). The argument reduction is exact, so the outputs keep their relative accuracy near every zero (and the tangent near every pole) at any magnitude of x, and the values at multiples of 1/2 (and of the tangent at odd multiples of 1/4) are exact.

Set res1 and res2 to balls containing \sin(\pi p/q) and \cos(\pi p/q), respectively res to \tan(\pi p/q), for words p \ge 0 and q \ge 1, to a relative accuracy of about B^{-n}. Either output of the first function may be NULL; the two must be distinct. The tangent returns 0 with res set to zero at a pole, else 1. The rational values 0, \pm 1, \pm 1/2 are exact. Both throw for q = 0.

Set res to a ball containing \exp(y), \log(y) resp. \operatorname{atan}(y) for every y in x. The results keep their relative accuracy near \exp(y) = 1, \log(y) = 0 and \operatorname{atan}(y) = 0. mp_real_exp_bits throws for |y| \ge 2^{\mathtt{FLINT\_BITS} - 5}, where the exponent of the result would leave the safe range. mp_real_log_bits returns 0 with res set to zero unless the ball is strictly positive, and 1 otherwise.

Set res (rs to \sinh and rc to \cosh, either may be NULL) to a ball containing \log(1 + y), e^y - 1, \sinh y, \cosh y, \tanh y, \operatorname{asin} y, \operatorname{acos} y, \operatorname{asinh} y, \operatorname{acosh} y resp. \operatorname{atanh} y for every y in x. The results keep their relative accuracy at the zeros of the functions, and tiny and huge arguments are handled without squaring them. The functions returning int return 0 with res set to zero when the ball is not certified to lie in the domain, and 1 otherwise.

Set res to a ball containing \exp(y) for every y in x, to a relative accuracy of about B^{-n}, without the precomputed tables of mp_real_exp_bits. These are meant for a single exponential at huge precision, where building the tables does not pay. The ball handling and the range |y| < 2^{\mathtt{FLINT\_BITS} - 5} are those of mp_real_exp_bits.

Set res to the exact number of partitions p(n) of n = n_{\mathit{hi}} B + n_{\mathit{lo}} (first function) or of any fmpz n (second function), as an exact mp_real_t. The size of the result is not limited by that of an fmpz_t.

Set res to a ball containing S with about \mathtt{FLINT\_BITS}(n - 1) accurate bits. The second function takes the coefficients as int64_t, as they appear in y-cruncher's formula files. Both throw if the series is not geometrically convergent, if Q(k) = 0 for some k \ge 1, if c_D = 0, or if the tail cannot be bounded.

Set res to \operatorname{atan}(p/q) resp. \operatorname{atanh}(p/q) for the integers 0 \le p < q given as limb arrays, to about n limbs. The quotient p/q must be bounded away from 1 (at most about 0.99); p \ge q throws.

Set or read the level of progress reports on stderr for huge computations: 0 is silent, 1 prints at most one line every two seconds from the long loops and the end of every phase that took that long, and 2 prints the start and end of every phase above a few thousand limbs with its size, thread budget, time and (on Linux) memory use.

Evaluate the Taylor series of the function at (x, n) with x < 2^{-32} by rectangular splitting. The outputs have n + 1 limbs for exp, sin, cos, sinh and cosh and n limbs for atan and atanh; the combined versions allow either output to be NULL. err is at most 10 for exp and 15 for the others, and the bound is one-sided (the result never exceeds the true value) for exp, sinh, cosh and atanh.

Evaluate \exp(x), \log(1 + x), \sin x and \cos x, \tan x resp. \operatorname{atan} x for any x \in [0, 1) by bitwise (BKM-type) argument reduction below 2^{-r} with cached tables of \log(1 + 2^{-i}) resp. \operatorname{atan}(2^{-i}), followed by a series and a reconstruction. The outputs have n + 1 limbs for exp, sin, cos and tan (either of ysin, ycos may be NULL) and n limbs for log1p and atan. err is at most 9r + 100 (exp), 3r + 64 (log1p), 6r + 128 (sin and cos), 8r + 256 (tan) resp. 4r + 64 (atan), for the r used. r = 0 selects the tuned default, at most 768; explicit values require r \ge 32. On 32-bit machines n \ge 2 is required.

Return the reduction parameter that r = 0 selects at n limbs: the depth of the per-size kernel where one exists, and the tuned value beyond (from 64 at a dozen limbs up to 768). The last one serves _mp_real_sin_cos_bitwise_rs and _mp_real_tan_bitwise_rs.

Evaluate \exp(t), respectively \sin t and g = 1 - \cos t, for a reduced argument (t, n) with t < 2^{-r}, r \ge 16 (r \ge 32 for alg = 1, 2). The exponential is written to n + 1 limbs; the sine and g are pure fractions written to n limbs, but both buffers need room for n + 1 limbs. err is the rigorous bound computed at run time, at most 96.

Evaluate \exp(x), respectively \sin x and \cos x, for any (x, n) in [0, 1) without precomputed tables: the argument is halved down to the depth r, the reduced function evaluated by _mp_real_exp_reduced resp. _mp_real_sin_cos_reduced, and the result squared resp. the angle doubled back. Outputs have n + 1 limbs; err is at most 128. The _r versions take the depth r explicitly, for tuning.

Evaluate \exp(x), \sin x and \cos x (either may be NULL), resp. \tan x for (x, n) in [0, 1) by multi-prime argument reduction [HJ2024]: x minus an integer combination of \log p_j (resp. of \pi/2 and the arguments of Gaussian primes), the reduced function by _mp_real_exp_reduced resp. _mp_real_sin_cos_reduced, and a product by a rational (resp. Gaussian rational) number. Outputs have n + 1 limbs; err is at most 3 for exp and 4 for the others. The _tune versions take the number of primes (2 to 64) and the weight budget of the reduction.

Set (y, n) to -\log x for (x, n) in [1/2, 1), respectively to \operatorname{atan} x for (x, n) in [0, 1), by one Newton–Taylor step over the forward function: -\log x = t - \log(1 + w) with w = x e^t - 1, resp. \operatorname{atan} x = t + \operatorname{atan} w with w = (x\cos t - \sin t)/(x \sin t + \cos t), from a starting value t correct to a fraction of the precision. err is computed at run time and at most 2.

Set (y, n) to -\log x for (x, n) in [1/2, 1) by the Sasaki–Kanada AGM formula [SK1982]; err is computed at run time and at most 2. The second function takes the number N of terms of the theta series (0 for the default, 4).

Set (y_s, n + 1) and (y_c, n + 1) (either may be NULL) to \sin(\pi p/q) and \cos(\pi p/q) for q \ge 1 and 0 \le p \le q/2. err is at most 4 for the closed forms and the Chebyshev iteration and the kernel's bound plus 2 (at most 6r + 130, r \le 768) when the sine and cosine kernel is used.

For a \in [1/B, 1) with a_{a_n - 1} \ne 0, set (q, n + 2) to an approximation of 1/a \in (1, B] with absolute error at most 4 B^{-n}/a.

For b \in [0, 1) with b_n \ge 1 fraction limbs and a as for _mp_real_inv_newton, set (q, n + 2) to an approximation of b/a with absolute error at most 4 B^{-n}/a. The first function multiplies b by _mp_real_inv_newton; the second performs a Karp–Markstein iteration [KM1997].

Set (\mathit{res}, n) to the fraction limbs of an approximation of 1/\sqrt{a} for a word 2 \le a < B, with error at most 2B^{-n}. The main function uses the third-order steps of _mp_real_rsqrt_newton, with the residual a y^2 - 1 exact and short.

For a \in [B^{-2}, 1) with one of its top two limbs nonzero, set (q, n + 2) to an approximation of 1/\sqrt{a} \in (1, B] with absolute error at most 4 B^{-n}/\sqrt{a}.

For a as for _mp_real_rsqrt_newton, set (q, n + 2) to an approximation of \sqrt{a} \in [1/B, 1) (the computed value can round up to 1) with absolute error at most 4 B^{-n}/\sqrt{a}, which is proportional to 1/\sqrt{a}, not to the output. The first function multiplies by the reciprocal square root; the second performs a Karp–Markstein iteration.

Return the relation table for the first num primes or Gaussian primes, 2 \le \mathit{num} \le MP_REAL_REL_MAX, respectively whether a call would find it ready.

_mp_real_machin_table returns the best set for num values and _mp_real_machin_table_max the size of the largest set. The accessors return the argument x_i, the coefficient c_{ij}, or a row of coefficients, the last one as the cached limbs (climbs limbs per coefficient) and sign bytes.

Set res to the angles \pi/2 and 2\arg \pi_j, j < \mathit{num}, as arb balls of prec bits.

The table-driven kernels at one, two and three fraction limbs (64-bit only): for v \in [0, 1), \exp(v) (n + 1 limbs), \log(1 + v) (n limbs), \sin v and \cos v (n + 1 limbs each) resp. \operatorname{atan} v (n limbs), with err between 8 and 64.

For v \in [0, 1) at two resp. three fraction limbs (most significant limb first), set y = \exp(v) and z = \exp(-v), each with an integral limb, at about 1.5 times the cost of one exponential. 64-bit only.

For three-limb fractions N < D with D normalised (top bit set), set w = \lfloor N B^3/D \rfloor by three 3-by-2 divisions in registers. 64-bit only.

For fractions X \ge 1/2 and Y \le X at two resp. three limbs (most significant limb first), set the result to \operatorname{atan}(Y/X) at the same number of fraction limbs, with one division instead of a division followed by the arctangent kernel. 64-bit only.

The dense-table kernels at n = 4, \ldots, 10 fraction limbs (64-bit only): for x \in [0, 1), \exp(x) (n + 1 limbs), \log(1 + x) and \operatorname{atan}(x) (n limbs), and \sin x and \cos x (n + 1 limbs each), within the bounds of the bitwise kernels for the reduction depths MP_REAL_*_OPT_R of impl.h.

The dense-table kernels as generic loops, for 2 \le n \le MP_REAL_DENSE_RT_MAX = 32 fraction limbs (64-bit only), with the deeper table levels computed at run time as they are needed and shared by all threads. Outputs as for the per-size kernels. The dispatchers call them from eleven limbs on.

Fill every entry not yet computed of the run-time dense tables of a family (MP_REAL_DENSE_RT_EXP, MP_REAL_DENSE_RT_LOG1P or MP_REAL_DENSE_RT_ATAN, the last shared with sine and cosine) in the tier serving n limbs, and return the width of that tier: 21 for n \le 21, else MP_REAL_DENSE_RT_MAX. 64-bit only.

Set (y, n + 1) to \tan v for v \in [0, 1) and n \le 32 from the sine and cosine kernels and one division (64-bit only). With both kernel outputs within E ulps, err is \lfloor 19E/4 \rfloor + 2. The tangent kernel uses it up to 22 limbs.

Tapered rectangular splitting for a fraction 0 \le x < 2^{-r}, r \ge 8, of n limbs. The first function computes \exp(\pm x) (n + 1 limbs, within 5 ulps) or \sin x, \sinh x, \operatorname{atan} x, \operatorname{atanh} x, 1 - \cos x, \cosh x - 1 (n limbs, within 2 ulps) as selected by func. The second computes the sine and 1 - \cos (or \sinh and \cosh - 1 when hyperbolic is set) sharing the powers; either output may be NULL. Outputs may alias x.

Set (\mathit{res}, n) to \tan x, \operatorname{atan} x, \operatorname{atanh} x, \sin x or 1 - \cos x (selected by func) for 0 \le x < 2^{-r}, within 6 ulps, by Horner's rule in x^2 over tapered widths with static coefficient tables. The supported ranges are n \le 80, r \ge 32 for the tangent, n \le 13, r \ge 10 for atan and atanh, and n \le 16, r \ge 18 for the sine and cosine.

Set (\mathit{res}, n) to \tan x for 0 \le x < 2^{-r}, r \ge 32, within 3 ulps, for n \le 80 and beyond as long as _mp_real_series_rs_tan_ok returns nonzero. res may alias x.

The bitwise reductions of _mp_real_exp_bitwise_rs, _mp_real_log1p_bitwise_rs, _mp_real_atan_bitwise_rs and of the sine, cosine and tangent, with signed digits and the full-precision arithmetic in radix 2^{56} on vec4n lanes, for x \in [0, 1) and 32 \le r \le 64n - 16, within the bounds of the corresponding _bitwise_rs functions (64-bit only). Outputs: (\mathit{res}, n + 1) for exp, (\mathit{res}, n) for log1p and atan, and n + 1 limbs for each of \sin, \cos, \tan (each may be NULL).

Set (\mathit{res}, n) to the floor of \log(1 + 2^{-i}) resp. \operatorname{atan}(2^{-i}) scaled by B^n, or one less, for i \ge 1, by binary splitting in ball arithmetic. These compute entries of the bitwise tables.

Set res to the balls \log p_j, resp. \pi/2 followed by 2\arg \pi_j, j < \mathit{num}, to about n limbs, by the Machin-type sets of _mp_real_machin_table and followup series, the series running on the available threads.

Set res to \log(u/v) for integers u > v > 0, resp. to \operatorname{atan}(d/t) for integers 0 < d < t, by Zuniga's Ramanujan-type series [Zun2025] [Zun2025b], to about n limbs. See Zuniga's series.

The exact binary splittings of the bit-burst steps. The first computes T/(Q B^{QE}) = \sum_{k=1}^N x^k/(k!\,B^{Dk}) for the integer x given by (x_p, x_n), x_n \le D; the second the sine and 1 - \cos of x B^{-D} over one common denominator, 1 - \cos = A B^{ae}/(Q B^{QE}) and \sin = x B^{-D}(Q B^{QE} - B B^{be})/(Q B^{QE}), with A and B truncated to lmax limbs. The third returns the number of terms for a reduced argument below 2^{-r} and prec bits.

For a binary splitting of N terms with fixed midpoints (ranges of one or two terms being leaves) whose split steps are the len entries of xexp, set uses[i] to the number of merges using the step xexp[i] and rmax[i] to the longest right child of those merges.

Prepare, free, and use the transforms of the power-table entries xpow[i] of xlen[i] limbs for products by partners of up to ymax[i] limbs. _mp_real_powtab_pre_mul sets z to \mathit{xpow}[i] \cdot y (\mathit{xlen}[i] + y_n limbs, z may overlap y) through the transform of entry i and returns 1, or returns 0 without writing anything when the entry has no transform or the product is outside its range.