ulong_extras.h

Returns the high word $\lfloor ab / B \rfloor$ of the double-word product $ab$.

Sets *a to $b + c$, $b - c$ or $bc$ modulo $B$ and returns 1 if the exact result does not fit in a ulong, otherwise 0.

Returns $n^{\mathit{exp}}$, with $0^0 = 1$. Overflow is not detected.

Returns $\lfloor \log_b n \rfloor$ and $\lceil \log_b n \rceil$ respectively. Requires $n \ge 1$ and $b \ge 2$.

Returns $\lceil \log_b 2^n \rceil$. Requires $b \ge 2$.

Returns $n$ with its lowest $b$ bits in reverse order; for example the 6-bit number $010110_2$ becomes $011010_2$.

Returns the number of digits of $n$ in the given base ($2 \le \mathit{base} \le 36$); n_sizeinbase returns 1 for $n = 0$. n_nonzero_sizeinbase10 is a branch-free version for base 10 and requires $n \ne 0$.

Returns $\lfloor \sqrt a \rfloor$. n_sqrtrem also sets *r to $a - \lfloor \sqrt a \rfloor^2$.

Returns 1 if $x$ is a perfect square, otherwise 0.

Returns $r = \lfloor n^{1/k} \rfloor$ for $k = \mathit{root}$, or 0 if $n$ or $k$ is 0. n_rootrem also sets the remainder $n - r^k$.

Returns $\lfloor n^{1/3} \rfloor$; n_cbrtrem also sets the remainder.

Alternative cube root algorithms used by n_cbrt: a bit-manipulation estimate in the spirit of the fast inverse square root, a binary search, and a Chebyshev approximation (requires $n \ne 0$).

n_is_perfect_power235 returns whether $n$ is a square, cube or fifth power. n_factor_power235 returns the root $r$ and sets *exp to $k \in \{2, 3, 5\}$ if $n = r^k$, and otherwise returns 0.

If $n = r^k$ with $k \ge 2$, returns $k$ and sets *root to $r$; otherwise returns 0. Both 0 and 1 count as squares. The $k$ found is not necessarily the largest possible.

Returns the precomputed inverse of $n \ne 0$ used by all _preinv functions [GraMol2010]. n_preinvert_limb_prenorm is a faster version for normalised $n$ (top bit set).

Returns $a \bmod n$. No restrictions on $a$; $n \ne 0$.

Returns $a \bmod n$ and sets *q to $\lfloor a/n \rfloor$. The last two variants take the number of leading zero bits of $n$ as norm; n_divrem_preinv_unnorm requires $\mathit{norm} > 0$.

Returns the quotient $\lfloor a/n \rfloor$. No restrictions on $a$; $n \ne 0$.

Returns $a \bmod n$ for the two-limb integer $a = a_{hi} B + a_{lo}$. No restrictions on $a$; $n \ne 0$.

Returns $a \bmod n$ for the three-limb integer $a = a_{hi} B^2 + a_{mi} B + a_{lo}$. Requires $a_{hi} < n$.

Returns $a \bmod n$ and sets *q to the quotient (0 or 1), assuming $n$ is normalised. No inverse is needed: the quotient is a single comparison.

Returns $1/n$ as a double, the precomputed inverse used by all _precomp functions. Requires $n > 0$.

Returns $a \bmod n$. Requires $n < 2^{53}$, $a < n^2$ and $a < 2^{\mathtt{FLINT\_BITS} - 1}$.

Returns $a \bmod n$ with no restrictions on $a$ or $n$.

Returns $a \bmod n$ and sets *q to the quotient. No restrictions on $a$ or $n$.

Returns the precomputed inverse for n_mod_barrett. Requires $n \ge 2$.

Returns $x \bmod n$ using Barrett reduction [Bar1987]. The lazy version returns a residue in $[0, 2n)$ instead of $[0, n)$. No restrictions on $x$. Both support $n = 1$ if npre is set to UWORD_MAX.

Returns the precomputed inverse for n_mod_lemire. Requires $n \ge 1$.

Returns $x \bmod n$ using the method of Lemire, Kaser and Kurz [LKK2019]. Only correct for $n, x < 2^{\mathtt{FLINT\_BITS} / 2}$.

Returns 1 and sets *q to $n/p$ if $p$ divides $n$; otherwise returns 0 and sets *q to 0.

Returns whether the odd number $d$ divides $n$, given $\mathit{inv1} = d^{-1} \bmod B$ and $\mathit{inv2} = \lfloor (B - 1)/d \rfloor$. Costs one multiplication. See Divisibility testing by multiplication.

Divides *n by the largest possible power of $p > 1$ and returns the exponent. Requires $n \ne 0$. n_remove2_precomp takes ppre = n_precompute_inverse(p).

Returns $\gcd(x, y)$, with $\gcd(x, 0) = \gcd(0, x) = x$.

Returns $g = \gcd(x, y)$ and sets *a, *b to unsigned values with $ax - by = g$. Requires $x \ge y$.

Returns $(a + b) \bmod n$, $(a - b) \bmod n$ and $-a \bmod n$ for reduced operands.

Returns $x^{-1} \bmod y$. Throws an exception if the inverse does not exist.

Returns $g = \gcd(x, y)$ and sets *a to the $0 \le a < y$ with $ax \equiv g \pmod y$. Requires $x < y$. When $y = 1$, returns 1 and sets *a to 0.

Returns the inverse of odd $n$ modulo $B$; the result is meaningless for even $n$.

Returns $ab \bmod n$. No restrictions on $a$, $b$ or $n$.

Returns $ab \bmod n$, computing the inverse of $n$ on the fly. Use n_mulmod2_preinv when $n$ is reused.

Returns $\lfloor a / 2^{\mathit{norm}} \rfloor \, b \bmod n$ for normalised $n$ and reduced $a, b$. If $a$, $b$ and $n$ are the operands and modulus shifted left by norm bits, this is their product modulo the unshifted modulus, shifted left by norm bits.

Returns $ab \bmod n$. Requires $n < 2^{53}$ and $a, b < n$.

Returns $\lfloor a B / n \rfloor$, the precomputed data for n_mulmod_shoup. Requires $a < n$.

Returns $ab \bmod n$ given a_precomp from n_mulmod_precomp_shoup. Requires $a < n$; $b$ is arbitrary.

Compute both the quotient $\lfloor aB/n \rfloor$ and the remainder $aB \bmod n$, as needed by n_mulmod_and_precomp_shoup. The second function recovers the remainder from a known quotient.

Sets *ab to $ab \bmod n$ and *ab_precomp to its precomputed data $\lfloor (ab \bmod n) B / n \rfloor$, without a division. Requires $a, b < n$. Useful for building tables of powers $a, a^2, a^3, \ldots$ together with their precomputed data.

Returns $a^{\mathit{exp}} \bmod n$ for $a < n$. Negative exponents are allowed; an exception is thrown if $a$ is then not invertible.

Returns $a^{\mathit{exp}} \bmod n$ for an unsigned exponent and $a < n$. n_powmod_ui_preinv requires normalised $n$, like n_mulmod_preinv.

Returns $a^{\mathit{exp}} \bmod n$ for a nonnegative multiprecision exponent and $a < n$.

Returns $a^{\mathit{exp}} \bmod n$ using a double inverse. Require $n < 2^{53}$ and $a < n$.

Returns the unique $0 \le x < m_1 m_2$ with $x \equiv r_1 \pmod{m_1}$ and $x \equiv r_2 \pmod{m_2}$. Requires coprime $m_1, m_2 > 1$ with $m_1 m_2$ fitting in a ulong, and reduced $r_1, r_2$. Nothing is checked.

Returns the Jacobi symbol $\left(\frac{x}{y}\right)$ for odd $y$ (not checked), with a binary algorithm [ShaSor1993]. n_jacobi_unsigned accepts a full-limb $x$.

Returns the smallest prime greater than $n$. Throws an exception if it does not fit in a ulong. The flag proved has no effect: the result is always a proved prime.

Return the number $\pi(n)$ of primes $p \le n$, or set $\mathit{lo} \le \pi(n) \le \mathit{hi}$.

Return the $n$-th prime $p_n$ (with $p_1 = 2$; $n = 0$ throws an exception), or set $\mathit{lo} \le p_n \le \mathit{hi}$.

Initialise and clear a prime iterator.

Returns the next prime and advances the iterator. The first call returns 2.

Makes the iterator continue with the first prime greater than $n$.

Maintain a per-thread cache of the first primes and their double inverses. The _readonly functions return pointers into the cache, which stay valid until n_cleanup_primes is called in the same thread.

Internal helpers of the prime iterator: extend its table of small primes until the second-largest entry is at least bound, or sieve the odd numbers in $[a, b]$ and point the iterator at the start of that block.

Returns 1 if $n$ is prime and 0 otherwise. Deterministic for all $n$.

Like n_is_prime for odd $n$, but skips trial division. Faster when $n$ is already known to have no small factors.

Primality test for the two-limb integer $n = n_{hi} B + n_{lo}$ with $n_{hi} \ne 0$. Returns 1 (prime), 0 (composite) or $-1$ (unknown). The answer is always definite below about $2^{81}$ on 64-bit machines.

Test whether odd $n$ is a strong probable prime to base $a$ [Mil1976; Rab1980], given the odd part $d$ of $n - 1$ and $0 < a < n$. The _precomp version requires $n < 2^{53}$. See Strong probable prime test.

Returns whether $i^{n-1} \equiv 1 \pmod n$. Requires $1 < i < n$ and $i \nmid n$.

Fibonacci probable prime test: checks $F_{n - (n/5)} \equiv 0 \pmod n$, where $(n/5)$ is a Jacobi symbol [CraPom2005]. Requires $\gcd(n, 10) = 1$.

Lucas probable prime test, similar to the one of Baillie and Wagstaff [BaiWag1980].

Baillie–PSW probable prime test [PSW1980; BaiWag1980]. No composite passing it is known, and it has been checked against tables of pseudoprimes up to $2^{64}$, so it is a primality test for word-size $n$.

Pocklington–Lehmer test [Poc1914; BLS1975]. Returns 1 if $n$ is proved prime, 0 if it is composite, and $-1$ if nothing was proved because the bases $2, \ldots, \mathit{iterations} - 1$ were not enough. If $\mathit{iterations}$ is 0 or 1, the number of bases is not limited and $-1$ is never returned for this reason.

Primality test based on pseudosquares, following Theorem 2.7 of [LukPatWil1996].

Obsolete; calls n_is_prime.

Initialises factors to the empty factorisation.

Sets factors to the complete factorisation of $n$ (empty for $n \le 1$); factors need not be initialised. The flag proved has no effect, as all primality tests used are deterministic for word-size integers.

Returns $p_1^{e_1} \cdots p_k^{e_k}$, or 0 if it does not fit in a ulong.

Multiplies the factorisation by $p^{\mathit{exp}}$, adding to the exponent of $p$ if it is already present.

Trial division by the first num_primes primes (starting at the prime with index start). Factors found are added to factors, which must be initialised; the unfactored cofactor is returned.

Partial factorisation: like n_factor_trial and n_factor, but stop as soon as the product of the prime factors found exceeds limit. Return the cofactor; n_factor_trial_partial also sets *prod to the product of the factors found.

Hart's one-line factoring algorithm [Har2012], a variant of Fermat's method. Returns a factor of $n$ or 0. Faster than SQUFOF below about $2^{40}$.

Lehman's algorithm [Leh1974]. Always returns a factor of $n$, but works only up to $10^{16}$ and is not used by n_factor.

Shanks' square form factorisation as described by Gower and Wagstaff [GowWag2008]. Returns a factor of $n$ that fits in one limb, or 0 if none was found within iters iterations. The ll version takes a two-limb $n = n_{hi} B + n_{lo}$.

Williams' $p + 1$ method [Wil1982] with prime bound $B_1$ and a random starting value $2 < c < n$. Returns a factor, or 0 or 1 on failure. The wrapper chooses parameters tuned for 31 to 64-bit $n$.

Pollard's rho method [Pol1975] with Brent's cycle detection [Bre1980], for composite $n$. Return 1 and set *factor to a nontrivial (not necessarily prime) factor on success, otherwise 0. The first function uses the polynomial $x^2 + a_i$ and starting value $x_i$; the second tries up to max_tries random choices.

Pollard–Brent rho for an odd two-limb $n$ with $n_{hi} \ne 0$. Returns 1 and writes a nontrivial two-limb factor on success, otherwise 0.

The elliptic curve method for an odd two-limb $n$ with $n_{hi} \ne 0$, trying up to curves curves with bounds $B_1$, $B_2$ (if $B_2 < B_1$, then $100 B_1$ is used). Returns 1 and writes a factor on success, otherwise 0.

The elliptic curve method [Len1987; ZD2006] for a single-limb $n$, trying up to curves random curves with stage bounds $B_1$ and $B_2$. Returns 1 or 2 if a factor was found in stage I or II, $-1$ if one was found while selecting a curve, and 0 on failure. The factor is stored in *f.

Compute the Möbius function $\mu(n)$, with $\mu(0) = 0$, or the values $\mu(0), \ldots, \mu(\mathit{len} - 1)$.

Returns whether $n$ is squarefree, that is, $\mu(n) \ne 0$. By convention 1 is squarefree and 0 is not.

Returns Euler's totient $\varphi(n)$, the number of $1 \le k \le n$ coprime to $n$.

Return $n! \bmod p$ for any modulus $p$ (no optimisations for composite $p$). The fast version uses fast multipoint evaluation in about $O(n^{1/2})$ operations up to logarithmic factors [Pol1974; BGS2007]; the other chooses between table lookup, a direct product and the fast algorithm depending on $n$.

Returns the square root of $a$ modulo the prime $p$ that lies in $[0, p/2]$, or 0 if $a$ is not a square. Requires $a < p$.

Compute all square roots of $a$ modulo $2^{\mathit{exp}}$, modulo $p^{\mathit{exp}}$, or modulo an integer $m$ given by its factorisation, and return their number. The roots are written to a newly allocated array stored in *sqrt, which the caller frees with flint_free. If $a$ is not a square, 0 is returned and *sqrt is set to NULL.

For odd non-square $n$, returns some $a$ with Jacobi symbol $\left(\frac{a}{n}\right) = -1$; for prime $n$ this is a quadratic nonresidue.

Return a primitive root modulo the prime $p$. The prefactor version takes the prime factors $q$ of $p - 1$ (exponents are ignored) and returns an $a$ with $a^{(p-1)/q} \ne 1$ for every $q$ given.

Returns $x$ with $a^x \equiv b \pmod n$ by baby-step giant-step [Sha1971], assuming $(\mathbb{Z}/n\mathbb{Z})^\times$ is cyclic, that is, $n$ is 2, 4, $p^k$ or $2p^k$ for an odd prime $p$.

Returns a uniformly random limb. _n_randlimb is an inline version for tight loops.

Returns a uniformly random integer with exactly bits bits (the top bit is set), or 0 if bits is 0.

Like n_randbits, but values with sparse binary representations are more likely.

Returns a uniformly random integer in $[0, \mathit{limit})$, or a random limb if limit is 0. The three functions are equivalent; _n_randint is inline.

Returns a random integer of random bit length between 0 and FLINT_BITS, with 0, 1, COEFF_MAX, WORD_MAX, UWORD_MAX and sparse values made more likely.

Like n_randtest, but never returns 0.

Returns a random prime with exactly bits bits, where $2 \le \mathit{bits} \le \mathtt{FLINT\_BITS}$. The result is always a proved prime; the flag proved has no effect.

Returns a random prime whose bit length is chosen at random between 2 and FLINT_BITS.

Computes the two-limb inverse $\lfloor 2^{3 \mathtt{FLINT\_BITS}} / m \rfloor$ of a two-limb $m > B$, as needed by n_ll_small_2_powmod and n_ll_small_powmod_triple.

Compute $b^e \bmod m$ for small single-limb bases $b$ and two-limb $e$ and $m > B$: the first function for $b = 2$, the second for three bases $b_1 \le b_2 \le b_3$ at once. Helpers for n_ll_is_prime.

Sets $(x : z) = 2 (x_0 : z_0)$, using $x = (x_0 + z_0)^2 (x_0 - z_0)^2$ and $z = 4 x_0 z_0 \left((x_0 - z_0)^2 + 4 a_{24} x_0 z_0\right)$.

Sets $(x : z) = (x_1 : z_1) + (x_2 : z_2)$ given their difference $(x_0 : z_0)$.

Sets $(x : z) = k (x_0 : z_0)$ with the Montgomery ladder.

Selects a curve and starting point from sigma with Suyama's parametrisation, storing $x_0$, $z_0$ and $(a + 2)/4$ in n_ecm_inf. Returns 1 if a factor of $n$ turned up along the way (stored in *f), otherwise 0.

Stage I with the num primes up to $B_1$ in prime_array, and stage II up to $B_2$ with primorial $P \approx \sqrt{B_2}$. Return 1 and set *f if a factor is found, otherwise 0.