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.