padic_radix.h

Initialises ctx for $\mathbf{Q}_p$ with default relative precision prec_rel and absolute precision prec_abs, either of which may be PADIC_RADIX_PREC_INF, and the given combination of the flags below. Negative relative precisions are raised to 0, and finite precisions are clamped to $\pm$PADIC_RADIX_ERR_MAX. That $p$ is prime is only checked in builds with assertions. Returns GR_SUCCESS. The context is a field in the sense of gr_ctx_is_field, and is exact if and only if both precisions are infinite.

Initialises ctx with a random prime and a random configuration for testing: absolute precision only, relative precision only, both (with precisions in $[1, \mathtt{maxprec}]$) or neither, and random flags PADIC_RADIX_SIGNED, PADIC_RADIX_TEST_LIMITS and PADIC_RADIX_FAST_LINALG.

Implementations of gr_ctx_clear and gr_ctx_write; the description written has the form Radix 7-adic numbers (rel prec 10, abs prec inf).

Initialise to the exact zero, clear, swap in constant time, and copy the structure without copying the limbs.

Return the precision $N$ of x (PADIC_RADIX_EXACT for an exact element), respectively whether x is exact.

Sets res to a random element with small valuation, exact or with a random error term (never with an error term in an exact context).

Writes x in the form (u) * p^v + O(p^N), omitting the power of $p$ when $v = 0$ and the error term for exact elements, for example 188316833 + O(7^10). The unit is written by radix_get_str_sum, or as a decimal integer if the context has the flag PADIC_RADIX_DECIMAL. Strings of this form, including O(p^n), are accepted by gr_set_str.

Set res to the exact zero, respectively to one, which becomes $0 + O(p^{t})$ if $t = \min(a, r)$ is not positive for the context precisions $a$ and $r$.

Sets res to x, truncated to the context precision.

Set res to the integer $x$, truncated to the context precision. A negative $x$ in a context without PADIC_RADIX_SIGNED is stored as an inexact residue, or gives GR_UNABLE if the context precision is infinite.

Set res to the exact integer $x$, ignoring the context precision; the unit may be negative regardless of the flags.

Sets res to the value of x if x is an exact integer. Returns GR_DOMAIN if $v < 0$ and GR_UNABLE if x is inexact or its valuation exceeds the size limit of the context.

Set res to $-x$, $x + y$, $x - y$ and $xy$. The error term of a sum is $O(p^{\min(N_x, N_y)})$ and that of a product $O(p^{\min(N_x + v_y, N_y + v_x)})$; the result is then truncated to the context precision. An exact zero factor gives an exact zero. An inexact zero $O(p^{N})$ enters the error term of a product with valuation $N$, so that for example $O(p^{N_x}) \cdot O(p^{N_y}) = O(p^{N_x + N_y})$.

Set res to $xc$, with the result of converting $c$ with padic_radix_set_si (etc.) followed by padic_radix_mul.

Set res to $\mathrm{res} \pm xy$ (respectively $\pm xc$), with exactly the result of a multiplication into a temporary followed by an addition or subtraction. Product digits that cannot reach the result (above its precision) are not computed, and sums of at most two limbs are formed inline without temporaries.

Set res to $x^n$, computing the power of the unit with truncated multiplications and rounding once. For $x = u p^v + O(p^N)$ and $n \ge 1$ the error term is $O(p^{(n - 1)v + N + v_p(n)})$, which is valid since every term $\binom{n}{k} a^{n - k} h^k$, $k \ge 1$, of $(a + h)^n - a^n$ has valuation at least $(n - 1)v + N + v_p(n)$; for an inexact zero $O(p^N)$ it is $O(p^{nN})$. Negative exponents invert first. An exact input keeps an exact result when it fits the context precision.

Set res to $1/x$, respectively $x/y$. The relative precision of the result is the smaller relative precision of the operands, capped by the context. Exact operands give an exact result when the quotient terminates (for example $1/49$ or $6/3$). Return GR_DOMAIN when dividing by an exact zero, and GR_UNABLE when the divisor is not known to be nonzero or when the context is exact and the quotient has an infinite expansion.

Returns whether x is nonzero: T_TRUE if it is exact and nonzero or its leading digit is known, T_FALSE for an exact zero and T_UNKNOWN otherwise.

Set res to a square root $\sqrt x$, respectively $1/\sqrt x$, with relative precision equal to that of $x$ (one digit less for $p = 2$), capped by the context. The root is a fixed canonical branch, so that a result at lower precision is the truncation of the result at higher precision, and rsqrt gives the reciprocal of sqrt. An exact perfect square gives its exact root. Return GR_DOMAIN if $x$ has odd valuation or its unit is not a square (and for rsqrt if $x$ is an exact zero), and GR_UNABLE if $x$ is not known to be nonzero, if $p = 2$ and the unit is known to fewer than 3 digits, or if the context is exact and the root is irrational. sqrt of an exact zero is zero.

Returns whether $x$ is a square: decided by the valuation and the unit modulo $p$ (modulo 8 for $p = 2$), and T_UNKNOWN for an inexact zero or, for $p = 2$, a unit known to fewer than 3 digits.

Elementwise $\mathrm{res}_i = \mathrm{vec}_i c$ and $\mathrm{res}_i = \mathrm{res}_i \pm \mathrm{vec}_i c$, with the results of the corresponding scalar functions; the scalar is converted once.

Ball predicates: T_TRUE only for exact elements with the property, T_FALSE if no value in the ball has it (for equality, if the balls are disjoint), and T_UNKNOWN otherwise.

Set res to $s \pm \sum_{i < len} a_i b_i$ with $a_i$ = vec1[i * stride1], $b_i$ = vec2[i * stride2], $s$ = initial (zero if initial is NULL) and the sign $-$ if subtract is nonzero. The naive version calls padic_radix_mul and padic_radix_add in a loop. The delayed version finds the minimum valuation and the precision of the result first, then accumulates the limb products of all terms in three-word accumulators, grouped by the residue of the digit offset modulo $e$ and by sign; at the end the residue classes are combined in one pass (each accumulator multiplied by $p^r$ for its residue $r$) and the carries are propagated once. The default uses the naive version for at most 4 terms at a precision of one limb, 3 terms at two limbs and 2 terms above (relative precision, or absolute if there is none), and, at precisions of up to two limbs, also for up to 7 terms whose valuations $v(a_i) + v(b_i)$ are not all equal.

When all terms are exact, the delayed version keeps the result exact if it fits in the context precision (as padic_radix_mul and padic_radix_add do), provided that the a priori extent of the terms, $\max_i (v(a_i) + e |a_i| + v(b_i) + e |b_i|)$ for units of $|a_i|$ and $|b_i|$ limbs, is at most one limb beyond the context precision; otherwise the sum is rounded to the context precision. In a context without finite precision, when the window of the delayed version (the extent of the exact terms, or the precision they determine) exceeds the size_limit of the context, it calls the naive version, which builds only the partial sums (and still gets exact cancellations of terms with huge valuations right).

The contiguous dot product and the dot product with vec2 reversed, as wrappers of padic_radix_dot_strided. These implement the generic dot product methods, which gr_mat_mul and the gr_poly products use.

Set $C = AB$. The packed version accumulates the limb products of every output directly. The fmpz version scales row $i$ of $A$ and column $j$ of $B$ to base frames $r_i = \min_t w(A_{it})$ and $c_j = \min_t w(B_{tj})$, truncates the residues to the largest output window, multiplies the resulting integer matrices with nmod_mat_mul (modulo $B$, when every window is a single limb) or fmpz_mat_mul (multimodular and blocked algorithms), and reduces every output modulo $B^s$ for its window of $s$ limbs before the conversion to radix $B$. The mpn_mod version is the same with the products done modulo $B^S$ for the largest window of $S$ limbs by mpn_mod_mat_mul (classical with delayed reduction, Waksman or multimodular), and returns GR_UNABLE if $B^S$ does not have between 2 and MPN_MOD_MAX_LIMBS binary limbs; padic_radix_mat_mul uses it in that range, the nmod version for one-limb windows and the fmpz version above. Outputs whose terms are all exact are summed exactly, in windows covering the a priori extents of their terms plus a carry and a sign limb, when padic_radix_dot_strided_delayed would keep them exact. Its precisions are computed with blockwise minima of the exponents, in $O(mn \cdot \text{blocks})$ operations for outputs whose row of $A$ or column of $B$ has exponents that are constant on blocks of the inner index (e.g. uniform valuations and precisions), and with the exact $O(k)$ minimum for the others; the results are identical to those of the packed version. The base frames $r_i + c_j$ can be far below the frame of an output whose terms all have higher valuations; the fmpz and mpn_mod versions return GR_UNABLE when the largest window in the base frames exceeds four times (plus 8 limbs) the largest window of an output in its own frame, since the integers would then be much larger than the precision requires (arbitrarily large for valuations spread beyond the precision).

padic_radix_mat_mul implements the generic method gr_mat_mul: the classical algorithm for inner dimension at most 3, the modular versions above tuned crossovers (estimated from an upper bound for the output windows: from dimension 8 when all windows are one limb, about 32 at two limbs, 16 at 3 to 8 limbs, 28 up to 15 limbs and 48 beyond, raised when the limb frames of the entries are spread out and when the exact precision bookkeeping is needed), and the packed version otherwise.

Middle product by Kronecker substitution in radix $B$: the coefficients, scaled to the base frames of the two polynomials and truncated to the largest output window of $S$ limbs, are placed in slots of $2S + O(1)$ limbs of a single radix integer, and one radix_mulmid (which uses the small-prime FFT directly in radix $B$ for large operands, without conversions to binary) gives every output in its own slot. The mod version instead multiplies the residues modulo $B^S$, avoiding the padding of the slots: with _nmod_poly_mulmid when $S = 1$ and with _mpn_mod_poly_mulmid (classical, Karatsuba or the small-prime FFT, as tuned for mpn_mod) when $B^S$ has at most MPN_MOD_MAX_LIMBS binary limbs, converting to binary and back; otherwise it returns GR_UNABLE. The output precisions come from min-plus convolutions of the exponent sequences (below), and coefficients whose terms are all exact are summed exactly (in windows covering the extents of the terms) when padic_radix_dot_strided_delayed would keep them exact. The results are those of _padic_radix_poly_mulmid_packed, except with the context flag PADIC_RADIX_FAST_BOUNDS, which allows smaller (but never larger) precisions. Both return GR_UNABLE in the same cases as the packed version. Outputs whose window in the base frames exceeds four times (plus 8 limbs) the largest window of an output in its own frame (all of whose terms have valuations far above the base frames, e.g. next to a few coefficients of huge valuation) are computed separately by padic_radix_dot_strided_delayed when there are few of them (costing at most about the transform and a sixteenth of the classical algorithm); otherwise, that is when the valuations vary beyond the precision along the outputs (alternating, trending or random valuations spread over more digits than the precision), both return GR_UNABLE. High-valuation coefficients that are beyond the precision of every output they contribute to cost nothing.

Set res to the coefficients $n_{lo}, \ldots, n_{hi} - 1$ (respectively $0, \ldots, n - 1$) of the product of the polynomials. Coefficients that are a single product are computed with padic_radix_mul. The unpacked versions implement the generic methods. For both lengths at least 4, from an estimate $S$ of the output windows in limbs (the size of the units, plus one if the valuations are not multiples of $e$, plus the spread of the limb frames): the mod version when $S = 1$; for $S = 2$ from a length of about 64, for $S \le 10$ from about 24, and the radix KS version for $S > 10$ from about 32 (these crossovers raised when the frames are spread out); the packed version otherwise. Before choosing a Kronecker substitution route when the frames are spread out, the windows of 16 sampled outputs are estimated from the context precision, and the packed version is used if the base frames inflate at least three of them as above; if a route still gives up, the packed version is used without trying the other. Shorter products use the classical algorithm.

Implements the generic method gr_mat_find_nonzero_pivot, used by the generic LU decomposition, solving, inversion and reduced row echelon form algorithms: chooses a row whose entry in the given column has minimal valuation, so that the elimination only involves multipliers of nonnegative valuation. Returns GR_UNABLE if there is no nonzero entry but some inexact zero, and GR_DOMAIN if all entries are exact zeros.

Set res to the determinant of the square matrix A. The lu_full version uses Gaussian elimination with full pivoting on valuations (an entry of minimal valuation in the whole remaining submatrix). Then the multipliers and the pivot rows have valuations at least that of the pivot, the Schur complements keep the absolute precision of the input, and the pivots have the valuations $a_1 \le \cdots \le a_n$ of the elementary divisors. For a matrix with entries known to $O(p^N)$ (and $a_n < N$) the determinant is obtained to precision $N + a_1 + \cdots + a_{n-1}$, which is optimal: it is that of the differential $d \det(A) = \operatorname{tr}(\operatorname{adj}(A)\, dA)$. In comparison, division-free algorithms (Berkowitz, cofactor expansion) only give precision about $N$, losing up to $v(\det A) - a_n$ digits, and Gaussian elimination with the first nonzero pivot loses precision growing with the dimension for small $p$, even for generic matrices (about 19 digits at $n = 60$ for $p = 7$). If the remaining submatrix consists of zeros, the determinant of that block is bounded using the sums of the row and column minima of the precisions.

The levels version is a blocked form of full pivoting. At the minimal valuation $a$ of the current Schur complement $S$, a rank-revealing LU decomposition of $(S / p^a) \bmod p$ over $\mathbb{F}_p$ selects $r$ rows and columns such that $S_{11} = p^a U$ with $U$ invertible modulo $p$: these $r$ pivots all have valuation $a$, as in full pivoting, and are eliminated at once, with $\det S = \pm \det(S_{11}) \det(S_{22} - S_{21} S_{11}^{-1} S_{12})$, using the LU decomposition of $S_{11}$, triangular solves and one matrix product. The new Schur complement has entries of valuation at least $a + 1$. The precision is again optimal, and the cost is that of a blocked LU decomposition (a single level for a matrix invertible modulo $p$).

padic_radix_mat_det implements the generic method gr_mat_det: the cofactor formula for $n \le 2$, the lu_full version up to a dimension depending on the precision (40 at one limb, 24 at two limbs, 16 beyond), and the levels version above.

Set $X = A^{-1} B$ (respectively $X = A^{-1}$), returning GR_DOMAIN if $A$ is singular and GR_UNABLE if this cannot be decided. The algorithm is Gauss-Jordan elimination with full pivoting by valuation levels on $[A \mid B]$: at each level the pivot block (as in padic_radix_mat_det_levels) is applied to its rows, and its columns are eliminated from all the other rows, including the pivot rows of the previous levels, with one matrix product. This order of operations (that of Gauss-Jordan elimination with full pivoting) gives every entry of the solution with the optimal precision, that of the differential $dX = A^{-1}(dB - dA\, X)$; LU decompositions with column pivots and back substitution lose some digits when $A$ is not invertible modulo $p$. When $A / p^a$ is invertible modulo $p$ for the minimal valuation $a$ of its entries, all pivots have valuation $a$ and the generic LU solve (with the pivots of padic_radix_mat_find_nonzero_pivot) is used directly, as it is always with the context flag PADIC_RADIX_FAST_LINALG. These functions implement the generic methods gr_mat_nonsingular_solve and gr_mat_inv.

Set cp to the characteristic polynomial $\chi(x) = \det(xI - A)$ of the square matrix mat.

The berkowitz version is the division-free Berkowitz algorithm, with the matrix-vector products done on the matrix packed once (falling back to _gr_mat_charpoly_berkowitz for dimension at most 5, in contexts without a relative precision, and for matrices with exact nonzero entries). For integral entries known to $O(p^N)$ it gives every coefficient to precision about $N$, which is optimal when the matrix is cyclic modulo $p$ but can lose many digits otherwise (dozens for the low coefficients of moderately structured matrices of dimension 16).

The hessenberg version reduces the matrix to Hessenberg form by similarity transformations with pivots of minimal valuation (so that the transformations are in $\mathrm{GL}_n(\mathbb{Z}_p)$), in $O(n^3)$ operations; it loses a few digits.

The lift version gives the optimal precision. With the matrix scaled by a power of $p$ to make it integral, $M$ the exact lift of the entries and $\operatorname{adj}(xI - M) = \sum_k B_k x^k$, a perturbation $E$ with $v(E_{ij}) \ge N_{ij}$ changes $c_k$ by $-\operatorname{tr}(B_k E)$ plus terms of order at least two in $E$, whose coefficients (minors of $xI - M$) are integral. Hence $c_k(A) = c_k(M) + O(p^{N_k})$ with $N_k = \min(2 N_{\min}, \min_{ij} (N_{ij} + v((B_k)_{ji})))$, the precision of the differential of $\chi$. The coefficients $c_k(M)$ are computed modulo $p^{2 N_{\min}}$ by the hessenberg version at a raised working precision (retried with more digits if needed, with the berkowitz version as a last resort), and the $B_k$ by the recurrence $B_{n-1} = I$, $B_{k-1} = M B_k + c_k I$, modulo one word first and modulo $p^{N_{\min}}$ if valuations at that bound matter, with $n$ matrix products. Exact matrices give the exact characteristic polynomial (multimodular), rounded to the context precision.

The lifted entries are reduced modulo the largest working precision, so that entries of high valuation cost nothing. The working precision $2(N_{\min} + s)$, where $p^s$ is the scaling, is large when an entry has a much smaller valuation (or error term) than the others; the lift version returns GR_UNABLE when $N_{\min} + s > 64(R + n) + 1024$, where $R$ is the largest relative precision of an inexact entry. For exact matrices it returns GR_UNABLE when the exact characteristic polynomial would exceed the size_limit of the context, when the entries of $M$ have more than $4 r + 256$ digits for a finite context precision $r$, or when they are dominated by powers of $p$ (valuations above four times the size of the units plus 256 digits), which Berkowitz's algorithm in $\mathbf{Q}_p$ handles without expanding them.

padic_radix_mat_charpoly implements the generic method. For matrices with integral entries of uniform precision $N$ that are cyclic modulo $p$ (tested with a random Krylov basis over $\mathbb{F}_p$), the optimal precision is $N$ for every coefficient: the hessenberg version is tried first and accepted if it attains it (as for matrices that are generic modulo $p$), and otherwise the berkowitz version, which always does, is used. Other matrices use the lift version, together with Danilevsky's algorithm (keeping the more precise result for each coefficient) where the bound $2 N_{\min}$ limits it, and the berkowitz version when the lift version gives up. With the context flag PADIC_RADIX_FAST_LINALG, the hessenberg version is used for all matrices.

Returns the number $n$ of terms after which the exponential series for an argument of valuation $v$ can be truncated modulo $p^N$, $n = \max(1, \lceil ((p-1)N - 1)/((p-1)v - 1) \rceil)$. Assumes that the series converges ($v \ge 1$, or $v \ge 2$ for $p = 2$).

Set rop to the nonnegative residue of $\exp(u p^v)$ modulo $p^N$, for $v \ge 1$ ($v \ge 2$ if $p = 2$), not checked. The rectangular version uses rectangular splitting, the balanced version the $p$-adic bit-burst method with binary splitting. _padic_radix_exp uses the rectangular version for $N$ below a cutoff between 80 and 200 that decreases with the size of $p$.

Set res to $\exp(x)$, with precision $\min(N_x, a, r)$ for the context precisions $a$ and $r$ (the result is a unit). Return GR_DOMAIN if the series diverges, and GR_UNABLE if convergence cannot be decided (an inexact zero of too low precision) or the context precision is infinite.

Returns the least $b \ge 2$ such that $iv - \operatorname{ord}_p(i) \ge N$ for all $i \ge b$, so that the series for $\log(1 - y)$ with $\operatorname{ord}_p(y) = v$ can be truncated after $b - 1$ terms.

Set rop to the residue of $\log(1 - y)$ modulo $p^N$ for an integer $y$ of valuation at least 1 (not checked). _padic_radix_log uses the rectangular version for $N$ below a cutoff between 80 and 350 that decreases with the size of $p$.

Set res to $\log(x)$, with absolute precision $\min(N_x, a, r)$. Return GR_DOMAIN unless $x \equiv 1 \pmod p$ (for every $p$, including $p = 2$), and GR_UNABLE if this cannot be decided or if $x$ is exact and the context precision infinite.

Canonicalises a freshly computed element (moving factors of $p$ from the unit to the valuation), enforces the representability bounds of PADIC_RADIX_ERR_MAX and truncates to the context precision. Returns GR_UNABLE if the valuation or precision cannot be represented, in which case res is set to $O(p^{-\mathrm{ERR\_MAX}})$ so that further arithmetic on it stays within the representable range.

Simple reference implementations of addition ($x + y$, or $x - y$ if sub is nonzero) and multiplication, used to test the optimised functions.