radix.h

Initialises radix with digit radix b and limb radix B = b^e where e = exp. Requires b \ge 2, e < \mathtt{FLINT\_BITS} and b^e < \beta, and throws an exception otherwise. With e = 0, the largest e with b^e < \beta is chosen.

Frees the memory used by radix.

Initialise radix with a random digit radix b (a random prime b for the second function) and a random exponent, for testing.

Return b, B and e.

Return the number of digits of the limb c (the least k with c < b^k). The underscore version requires c > 0; the other returns 0 for c = 0.

Return the number of digits of (x, n). The underscore version requires n > 0 and x_{n-1} \ne 0; the other allows n = 0 and zero top limbs, which it skips, and returns 0 for zero.

Returns the number of trailing zero digits of the limb c > 0, that is, the largest k such that b^k divides c.

Set (\mathit{res}, n) to a uniformly random integer in [0, B^n - 1], respectively (\mathit{res}, \lceil n/e \rceil) to a uniformly random integer in [0, b^n - 1].

Like radix_rand_limbs and radix_rand_digits, but with a non-uniform distribution producing long runs of limbs equal to 0 or B - 1, for testing.

Set (\mathit{res}, \mathit{an}) to a + b, respectively a - b, modulo B^{an} and return the carry, respectively the borrow (0 or 1). Requires \mathit{an} \ge \mathit{bn} \ge 1.

Sets (\mathit{res}, \mathit{an}) to the B-complement (-a) \bmod B^{an} and returns the borrow (1 if a \ne 0, 0 otherwise).

Shift (a, n) by k = e digits, for 1 \le k less than the limb exponent: the left shift sets (\mathit{res}, n) to a b^k \bmod B^n and returns the k digits shifted out (as the low digits of the returned limb); the right shift sets (\mathit{res}, n) to \lfloor a / b^k \rfloor and returns the k digits shifted out as the high digits of the returned limb, that is, (a \bmod b^k) b^{e - k}. These are faster than radix_mul_1 and radix_divrem_1 with b^k. Shifts by whole limbs are pointer offsets.

Fused addition and subtraction with a shift of b by 1 \le \mathit{sh} < e digits: set (\mathit{res}, L) with L = \max(\mathit{an}, \mathit{bn} + 1) to a + b \cdot b^{\mathit{sh}}, a - b \cdot b^{\mathit{sh}} and b \cdot b^{\mathit{sh}} - a, and return the carry (0 or 1), respectively the borrow, which is 0 since the result must be nonnegative. Requires \mathit{bn} \ge 1. res may alias a but not b. The shifted operand is never formed in memory.

Returns the sign of a - b \cdot b^{\mathit{sh}} for normalised (a, \mathit{an}) and (b, \mathit{bn}) and 1 \le \mathit{sh} < e.

Sets (\mathit{res}, n) to ac \bmod B^n and returns the carry \lfloor ac / B^n \rfloor. Requires n \ge 1 and c < B.

Sets (\mathit{res}, \mathit{an}) to 2a \bmod B^{an} and returns the carry.

Set (\mathit{res}, \mathit{an} + \mathit{bn}) to ab, respectively (\mathit{res}, 2\mathit{an}) to a^2. Requires \mathit{an}, \mathit{bn} \ge 1; res must not overlap the inputs. Both are calls to radix_mulmid with the full window (there is no special squaring code yet).

Windowed product: for \mathit{an}, \mathit{bn} \ge 1 and 0 \le \mathit{lo} < \mathit{hi} \le \mathit{an} + \mathit{bn}, sets (\mathit{res}, \mathit{hi} - \mathit{lo}) to the limbs [\mathit{lo}, \mathit{hi}) of the product, computed without the carries from the limbs below \mathit{lo}. res must not overlap the inputs.

The algorithms behind radix_mulmid, with its interface: the schoolbook algorithm with delayed reduction of the column sums, efficient when an operand or the window is short; multiplication of the limb sequences as polynomials with the small-prime FFT of fft_small (only available if FLINT is built with it), reconstructing only the coefficients that reach the window and propagating carries in radix B; Kronecker substitution into one product in the binary radix \beta (each limb padded to two words); and an unoptimised reference implementation through fmpz_poly.

Returns the number of pairs (i, j) with 0 \le i < \mathit{an}, 0 \le j < \mathit{bn} and i + j < k, as a double (it may exceed the range of slong on 32-bit machines). It counts the limb products of the classical algorithm and is used for tuning decisions.

The length of the convolution that radix_mulmid_fft_small computes for the window [\mathit{lo}, \mathit{hi}), used by radix_mulmid to choose the algorithm: with \mathit{an}' and \mathit{bn}' the numbers of limbs of the operands that reach the window, it is \min(\mathit{an}', \mathit{bn}') + (\mathit{hi}' - \mathit{lo}') if the window (clipped to the product) is shorter than \max(\mathit{an}', \mathit{bn}'), and \mathit{an}' + \mathit{bn}' otherwise; 0 for an empty window.

The radix B analogue of _flint_mpn_mulhigh_known_low: sets (\mathit{out}, \mathit{khi} - \mathit{klo}) to the exact limbs [\mathit{klo}, \mathit{khi}) of xy, given that its low \mathit{klo} limbs are (\mathit{kl}, \mathit{kl\_len}) (zero above kl_len), using a windowed product with guard limbs. scratch must have room for \mathit{khi} limbs. Used by the Hensel division and p-adic square roots.

Set (\mathit{res}, \mathit{an}) to \lfloor a/d \rfloor and return the remainder, for 1 \le d < B; the same assuming that d divides a (the result is undefined otherwise); and division by 2.

Given (a, \mathit{an}) with a_{an-1} \ne 0, read as a fixed-point number a \in [1/B, 1) with \mathit{an} fraction limbs, set (q, n + 2) to an approximation of 1/a \in (1, B] with n fraction limbs and two integral limbs (the top one may be zero). The absolute error is less than 4 B^{-n}/a (relative error less than 4 B^{-n}). The general version uses Newton iteration from the basecase.

Given the divisor (A, \mathit{An}) with nonzero top limb, read as a fixed-point number in [1/B, 1) with \mathit{An} \ge 1 fraction limbs, and the dividend (B, \mathit{Bn}), read as a fixed-point number in [0, 1) with \mathit{Bn} \ge 1 fraction limbs, set (Q, n + 2) to an approximation of the quotient with n fraction limbs and two integral limbs; the absolute error is less than 4 B^{-n} divided by the divisor. Note the order of the arguments: the dividend comes first. The invmul version multiplies by an approximate reciprocal; the general version uses Karp–Markstein iteration [KM1997].

Sets (q, \mathit{an} - \mathit{bn} + 1) to \lfloor a/b \rfloor and, if r is not NULL, (r, \mathit{bn}) to the remainder. Requires \mathit{an} \ge \mathit{bn} \ge 1 and b_{bn-1} \ne 0.

The algorithms behind radix_divrem, with its interface: conversion to binary, division by a Newton reciprocal, and Karp–Markstein division (quotient from a half-precision reciprocal and one correction step).

radix_divrem with a precomputed inverse of b given as (\mathit{binv}, \mathit{binvn} + 2) with \mathit{binvn} fraction limbs, as computed by radix_inv_approx. Requires \mathit{binvn} \ge \mathit{an} - \mathit{bn} + 1; a few extra limbs can improve performance. r may be NULL.

If b divides a, sets (q, \mathit{an} - \mathit{bn} + 1) to the quotient and returns 1; otherwise returns 0, the contents of q being undefined.

Sets (q, \mathit{an} - \mathit{bn} + 1) to a/b, assuming that the division is exact; the result is undefined otherwise.

Returns -1, 0 or 1 according to whether (x, n) is less than, equal to or greater than \lfloor B^n/2 \rfloor. Requires n \ge 1. Used for rounding to nearest and for symmetric remainders.

Sets (s, \lceil \mathit{an}/2 \rceil) to \lfloor \sqrt a \rfloor and, if r is not NULL, (r, \lceil \mathit{an}/2 \rceil + 1) to the remainder a - s^2 \le 2s. Requires \mathit{an} \ge 1 and a_{an-1} \ne 0. Inputs shorter than a cutoff depending on the size of B (43 limbs for B = 10^{19}) are converted to binary and use flint_mpn_sqrtrem; longer ones use the Newton square root.

Sets (s, \lceil \mathit{an}/2 \rceil) to \lfloor \sqrt a \rfloor and returns 1 if a is a perfect square and 0 otherwise. Requires \mathit{an} \ge 1 and a_{an-1} \ne 0.

The algorithms behind radix_sqrtrem, with its interface: conversion to binary, and the root from a reciprocal square root by one Karp–Markstein step, with a final correction.

Set (\mathit{res}, n) to the fraction limbs of an approximation of 1/\sqrt a for 2 \le a < B, with error less than 2 B^{-n}.

Given (a, \mathit{an}) read as a fixed-point number a \in [B^{-2}, 1) with \mathit{an} fraction limbs (one of the top two limbs nonzero), set (q, n + 2) to an approximation of 1/\sqrt a \in (1, B] with n fraction limbs and two integral limbs, with absolute error less than 4 B^{-n}/\sqrt a. The range of a is wider than for radix_inv_approx so that inputs with an odd number of limbs can be handled, the scale factor B^{1/2} not being a limb shift.

With a as for radix_rsqrt_approx, set (q, n + 2) to an approximation of \sqrt a \in [1/B, 1) with n fraction limbs and two integral limbs (normally zero, though the value can round up to 1), with absolute error less than 4 B^{-n}/\sqrt a. The error is proportional to 1/\sqrt a rather than to \sqrt a, because the corrections of the last iteration are multiplied by an approximation of 1/\sqrt a \le B.

If (x, \mathit{xn}) is invertible modulo B^n (that is, x_0 is coprime to b), sets (\mathit{res}, n) to its inverse and returns 1; otherwise returns 0. Newton iteration on windowed products.

Given (b, \mathit{bn}) invertible modulo B, sets (q, n) to the B-adic (Hensel) quotient with qb \equiv a \pmod{B^n} and returns 1; returns 0 if b is not invertible. If rem is not NULL, the remainder (a - qb)/B^n \bmod B^{bn} is written to (\mathit{rem}, \mathit{bn}), so that a \equiv qb + B^n \mathit{rem} \pmod{B^{n + bn}}.

The algorithms behind radix_divmod_bn, with its interface: a one-limb divisor (with rem of one limb), the blockwise algorithm (one low and one high product per block of \mathit{bn} limbs), and the Karp–Markstein form with a half-precision inverse.

For a prime digit radix p (with e \ge 3 if p = 2): if (x, \mathit{xn}) is a unit and a square modulo p (modulo 8 if p = 2), set (\mathit{res}, n) to a square root of x modulo B^n, respectively to the inverse of that square root, and return 1; otherwise return 0. The square root returned is the one whose residue modulo p is the smallest nonnegative one.

Convert (a, \mathit{an}) from radix B to the binary radix, writing the limbs to res and returning the normalised length of the result. Leading zero limbs are not counted and may or may not be written. res must have room for the largest integer with \mathit{an} limbs in radix B, plus one extra limb, and must not overlap a. The basecase is Horner's rule (used below 50 limbs).

Convert the binary (a, \mathit{an}) to radix B, writing the limbs to res and returning the normalised length of the result. res must have room for radix_set_mpn_need_alloc(an, radix) limbs and must not overlap a. The basecase is repeated division by B (used below 40 limbs). The divide and conquer version caches the most recently used table of powers (per thread), so that repeated conversions of inputs of the same size, as when printing many numbers, do not recompute it.

Returns a number of output limbs that is sufficient for radix_set_mpn with an input of n limbs.

Frees a table of powers.

Returns the decimal representation of (x, n) as a null-terminated string, with a leading minus sign if negative is 1; for n = 0 the string is 0. The top limb must be nonzero (if n > 0). If res is NULL, a new string is allocated with flint_malloc; otherwise res must have room for the result and the terminating null. If b = 10, this takes time O(n) (each limb is printed with table lookups for pairs of digits); otherwise x is converted to binary and then to decimal, in time O(M(n) \log n).

Returns (x, n) written as a sum over its limbs. For example, 12345678 in radix 10^3 is printed as 12 * 10^6 + 345 * 10^3 + 678 with negative and ascending both 0, and as -(678 + 345 * 10^3 + 12 * 10^6) with both 1. Memory as for radix_get_str_decimal.

Initialises ctx as the generic ring of integers in radix b^e (with e = 0 selecting the largest admissible exponent), whose elements are radix_integer_t. The context owns its radix_t.

Initialise to zero and free an integer, and make sure that it has room for nlimbs limbs, returning the pointer to its limbs (which the caller may then fill, setting size accordingly).

Returns whether all limbs of x are reduced modulo B and the top limb is nonzero. Operations on a non-normalised integer have undefined behaviour.

Assignment and conversion; radix_integer_set_fmpz and radix_integer_get_fmpz use the radix conversion functions.

Set res to a random signed integer of at most max_limbs limbs, with a uniform respectively a test distribution of the limbs.

Predicates and comparisons; cmp and cmpabs return the sign of x - y and |x| - |y|.

The number of limbs, respectively digits, of |x|; the ssize versions return it negated if x is negative.

Return the limb, respectively the digit, of |x| at the given nonnegative position (zero beyond the top).

Set res to x with the limb (digit) at the nonnegative position index of |x| replaced by c < B (c < b). The sign of x is kept; if x is zero and c nonzero, the result is positive.

Negation, absolute value, sum, difference and product.

Set res to x B^n (x b^n), respectively to x / B^n (x / b^n) truncated towards zero, for n \ge 0.

Set res to x + y b^v, respectively x - y b^v, for any digit shift v \ge 0, without forming the shifted operand: the shift is split into whole limbs (a pointer offset) and v \bmod e digits (the fused radix_addlsh family). Any aliasing is allowed.

Return the largest n such that B^n (respectively b^n) divides x, that is, the number of trailing zero limbs (digits); 0 for x = 0.

Reductions modulo B^n (respectively b^n): trunc reduces |x| and keeps the sign of x; mod returns the representative in [0, B^n - 1] ([0, b^n - 1]); smod returns the symmetric representative, in [-\lfloor B^n/2 \rfloor, \lfloor B^n/2 \rfloor] for odd B and in [-B^n/2, B^n/2 - 1] for even B. These take time O(n) plus the time to copy.

Sets res to xy with the absolute value reduced modulo B^n, using a low product.

Inverse, quotient a/b, inverse square root and square root modulo B^n, by radix_invmod_bn, radix_divmod_bn, radix_rsqrtmod_bn and radix_sqrtmod_bn applied to the absolute values; the sign of the inputs is carried over to the result. Return 1 on success and 0 if the operand is not invertible (respectively not a square).

If b divides a, radix_integer_div sets q to a/b and returns 1; otherwise it sets q to zero and returns 0. radix_integer_divexact sets q to a/b assuming that the division is exact (the result is undefined otherwise).

Division with remainder a = qb + r, where the quotient is rounded towards zero (tdiv), towards -\infty (fdiv) or towards +\infty (cdiv), as in GMP's mpz functions. b must be nonzero.

x is invertible in \mathbb{Z} if and only if x = \pm 1. radix_integer_inv sets res to 1/x = x and returns 1 in that case and returns 0 otherwise.

radix_integer_sqrt sets s to \lfloor \sqrt x \rfloor and returns 1 if x is a perfect square and 0 otherwise; for negative x it returns 0 without modifying s. radix_integer_sqrtrem sets s and the remainder r = x - s^2, and throws an exception for negative x. radix_integer_is_square returns whether x is a perfect square (0 for negative x).

Computes the data for division by d \ge 1.

Return \lfloor x/d \rfloor for any word x: n_div_precomp selects the method from pre, and the other functions are the individual methods for callers that know the case: x \gg \ell for a power of two (m = 0), \lfloor xm / 2^{\mathtt{FLINT\_BITS}} \rfloor \gg \ell for c = 0, and \lfloor (x + 1) m / 2^{\mathtt{FLINT\_BITS}} \rfloor \gg \ell for c = 1, where x + 1 saturates at \beta - 1. The bounded versions require x \ne \beta - 1 and skip the saturation.

The remainder x \bmod d, and the quotient with the remainder stored in *r, with the same variants as the quotient functions.

Returns x + 1, or x if x = \beta - 1 (saturating increment).