radix.h – multiprecision arithmetic in general radix¶
This module implements multiprecision arithmetic in a general radix \(2 \le b < \beta\) where \(\beta = 2^{\mathtt{FLINT\_BITS}}\) is the machine word radix.
The user selects both a digit radix \(b\) and a limb radix \(B = b^e\) where \(B < \beta\), allowing several digits to be packed into each machine word. A limb is a chunk of \(e\) consecutive digits stored in one word. Efficiency is typically maximized by choosing the largest limb radix that fits in a machine word. For example, for decimal arithmetic one should typically select \(B = 10^{19}\) on a 64-bit machine and \(B = 10^9\) on a 32-bit machine. A smaller limb radix may however be more efficient in some circumstances.
Arithmetic in radix \(B\) is generally slower than arithmetic in the machine word radix \(\beta\), but has the advantage that base-related operations such as reduction modulo \(b^n\) and string input/output in base \(b\) can be done in time \(O(n)\) rather than \(O(M(n))\) or \(O(M(n) \log n)\). The two principal applications are decimal arithmetic (with \(b = 10\)) and \(p\)-adic arithmetic (with \(b = p\) a small odd prime number). This module can also be useful for experimenting with multi-limb algorithms where debugging and discovering corner cases is easier in a small base like 10 than with full-word limbs.
The functions in this module are not optimized for radix \(2^e\), where
the machine word radix is applicable (see GMP’s mpn and FLINT’s
flint_mpn routines).
Radix objects¶
-
type radix_struct¶
-
type radix_t¶
A context object defines the radix \(B = b^e\) and stores additional data such as various precomputed powers \(b^k\) and their inverses used for fast division and modular reduction.
-
void radix_init(radix_t radix, ulong b, unsigned int e)¶
Initializes context object representing radix \(B = b^e\) given \(b \ge 2\) and an exponent \(e \ge 1\) such that \(B < \beta\). Optionally, the user can input \(e = 0\), in which case \(e\) will automatically be set to the largest exponent such that \(b^e < \beta\).
-
void radix_init_randtest(radix_t radix, flint_rand_t state)¶
Initializes to a randomly chosen radix.
Low-level natural number arithmetic¶
A limb has type ulong. The type nn_ptr denotes a pointer
to an array of writable limbs and nn_srcptr denotes a pointer
to an array of readonly limbs.
We use the notation (x, n) for the natural number in the range \([0, B^n-1]\)
represented by an array of n limbs starting at address x.
Except where otherwise noted, the following rules apply:
Input limbs must be normalised to the range \([0, B - 1]\), and output limbs will also be normalised to this range.
Aliasing is permitted as long as aliased operands start at the same address.
-
void radix_rand_limbs(nn_ptr res, flint_rand_t state, slong n, const radix_t radix)¶
Sets (res, n) to a uniformly random integer in \([0, B^n-1]\).
-
void radix_rand_digits(nn_ptr res, flint_rand_t state, slong n, const radix_t radix)¶
Sets (res, rn) to a uniformly random integer in \([0, b^n-1]\) where \(rn = \lceil n / e \rceil\).
-
void radix_randtest_limbs(nn_ptr res, flint_rand_t state, slong n, const radix_t radix)¶
Sets (res, n) to a non-uniformly random integer in \([0, B^n-1]\). This will produce long strings of limbs with the value 0 or \(B-1\) with high probability.
-
void radix_randtest_digits(nn_ptr res, flint_rand_t state, slong n, const radix_t radix)¶
Sets (res, rn) to a non-uniformly random integer in \([0, b^n-1]\) where \(rn = \lceil n / e \rceil\).
-
ulong radix_neg(nn_ptr res, nn_srcptr x, slong n, const radix_t radix)¶
Sets (res, n) to the \(B\):s complement negation of (x, n), i.e., \((-x) \bmod B^n\), returning borrow.
-
ulong radix_add(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, const radix_t radix)¶
-
ulong radix_sub(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, const radix_t radix)¶
Sets (res, xn) to the sum or difference of (x, xn) and (y, yn), returning carry or borrow. Requires \(xn \ge yn \ge 1\).
-
void radix_mul(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, const radix_t radix)¶
Sets (res, xn + yn) to the product of (x, xn) and (y, yn). Requires \(xn \ge yn \ge 1\). Does not allow res to be aliased with either x or y.
This function is currently a wrapper of
radix_mulmid().
-
void radix_sqr(nn_ptr res, nn_srcptr x, slong xn, const radix_t radix)¶
Sets (res, 2 xn) to the square of (x, xn). Requires \(xn \ge 1\). Does not allow res to be aliased with x.
-
void radix_mulmid(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, slong lo, slong hi, const radix_t radix)¶
-
void radix_mulmid_classical(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, slong lo, slong hi, const radix_t radix)¶
-
void radix_mulmid_fft_small(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, slong lo, slong hi, const radix_t radix)¶
-
void radix_mulmid_KS(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, slong lo, slong hi, const radix_t radix)¶
-
void radix_mulmid_naive(nn_ptr res, nn_srcptr x, slong xn, nn_srcptr y, slong yn, slong lo, slong hi, const radix_t radix)¶
Short, truncated or middle product. Requires \(xn \ge yn \ge 1\) and \(0 \le lo < hi \le xn + yn\). Does not allow res to be aliased with either x or y.
Viewing as \(x\) and \(y\) as polynomials \(X, Y \in \mathbb{Z}[T]\) evaluated at \(T = B\), we implicitly compute the polynomial product \(Z = XY\) and extract the slice \(Z_{lo:hi} = \sum_{k=lo}^{hi-1} ([T^k] Z) T^{k-lo}\), and write \(Z_{lo:hi}(B) \bmod B^{hi-lo}\) to (res, hi - lo). We omit carries from lower order terms and we discard carries that would go into higher order terms.
In the special case \(lo = 0\), this sets (res, hi) to the low product \(xy \bmod B^{hi}\).
In the special case \(hi = xn + yn\), this sets (res, hi - lo) to an approximate high product \(r\) satisfying \(xy \ge r B^{lo} \ge xy - \min(xn,yn,lo) B^{lo+1}\).
In the special case \(lo = 0\), \(hi = xn + yn\), this computes the full product.
Various multiplication algorithms are implemented:
classical is an efficient implementation of the schoolbook algorithm, good when either yn or \(hi - lo\) is small.
fft_small multiplies using FFT. This is only available if FLINT is built with support for the fft_small module.
KS denotes Kronecker substitution, converting to a larger integer multiplication in the machine word radix.
naive is an unoptimized reference implementation which does a polynomial multiplication using
fmpz_poly.
Radix conversion¶
-
slong radix_get_mpn_basecase(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
-
slong radix_get_mpn_divconquer(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
-
slong radix_get_mpn(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
Convert (a, an) to the machine word radix, writing the output limbs to res and returning the exact size of the result in the machine word radix. Leading zero limbs are omitted from the output size and may or may not be written. Requires that res has enough space for the largest representable integer with an limbs in the given radix, plus one extra limb.
-
slong radix_set_mpn_basecase(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
-
slong radix_set_mpn_divconquer(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
-
slong radix_set_mpn(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)¶
Convert (a, an) from the machine word radix, writing the output limbs to res and returning the exact size of the result in the target radix. Leading zero limbs are omitted from the output size and may or may not be written. Requires that res has space for at least
radix_set_mpn_need_alloc(an, radix)limbs.
-
slong radix_set_mpn_need_alloc(slong n, const radix_t radix)¶
Return a number of output limbs for which
radix_set_mpn()is safe to call with input of length \(n\).