Functions in this section cover 4 classical arithmetic topics
* factorization over the integers, including primality proving;
* arithmetic functions;
* finite fields;
* factorization in polynomial rings, in particular over finite fields and over the rationals;
Arithmetic functions are naturally defined over ℤ or ℤ > 0 and related to integer factorization. An integer argument N can be given in the following alternate formats:
* t_MAT: its factorization fa = factor(N),
* t_VEC: a pair [N, fa] giving both the integer and
its factorization.
This allows to compute different arithmetic functions at a given N while factoring the latter only once.
? N = 10!; faN = factor(N);
? eulerphi(N)
%2 = 829440
? eulerphi(faN)
%3 = 829440
? eulerphi(S = [N, faN])
%4 = 829440
? sigma(S)
%5 = 15334088
All arithmetic functions in the narrow sense of the word — Euler's
totient function, the Moebius function,
the sums over divisors or powers of divisors, Ramanujan τ
function, etc. — call the factoring machinery described
under factorint. It includes trial division by small primes,
Shanks SQUFOF, Pollard Rho,
ECM and MPQS stages, with early exit options for functions
such as moebius and issquarefree. This machinery relies on a
fairly strong probabilistic primality test but, by default, integer
factorizations are not fully proven, see ispseudoprime. You
may also set
default(factor_proven, 1)
to ensure that all tentative factorizations are proven. This should not slow down PARI too much, unless prime numbers with hundreds of decimal digits occur frequently in your application.
Adds the integers contained in the
vector x (or the single integer x) to a special table of
"user-defined primes", and returns that table. Whenever factor is
subsequently called, it will trial divide by the elements in this table.
If x is empty or omitted, just returns the current list of extra
primes.
? addprimes(37975227936943673922808872755445627854565536638199)
? factor(15226050279225333605356183781326374297180681149613806\
88657908494580122963258952897654000350692006139)
%2 =
[37975227936943673922808872755445627854565536638199 1]
[40094690950920881030683735292761468389214899724061 1]
? ##
*** last result computed in 0 ms.
The entries in x must be primes: there is no internal check, even if
the factor_proven default is set. To remove primes from the list use
removeprimes.
The library syntax is GEN addprimes(GEN x = NULL).
Number of prime divisors of the integer |x| counted with multiplicity:
? factor(392) %1 = [2 3] [7 2] ? bigomega(392) %2 = 5; \\ = 3+2 ? omega(392) %3 = 2; \\ without multiplicity
The library syntax is long bigomega(GEN x).
If n is an integer written as
n = df2 with d squarefree, returns d. If flag is nonzero,
returns the two-element row vector [d,f]. By convention, we write 0 = 0
x 12, so core(0, 1) returns [0,1].
The library syntax is GEN core0(GEN n, long flag).
Also available are GEN core(GEN n) (flag = 0) and
GEN core2(GEN n) (flag = 1)
A fundamental discriminant is an integer of the form t = 1 mod 4 or 4t = 8,12 mod 16, with t squarefree (i.e. 1 or the discriminant of a quadratic number field). Given a nonzero integer n, this routine returns the (unique) fundamental discriminant d such that n = df2, f a positive rational number. If flag is nonzero, returns the two-element row vector [d,f]. If n is congruent to 0 or 1 modulo 4, f is an integer, and a half-integer otherwise.
By convention, coredisc(0, 1)) returns [0,1].
Note that quaddisc(n) returns the same value as coredisc(n),
and also works with rational inputs n ∈ ℚ*.
The library syntax is GEN coredisc0(GEN n, long flag).
Also available are GEN coredisc(GEN n) (flag = 0) and
GEN coredisc2(GEN n) (flag = 1)
x and y being vectors of perhaps different lengths but with y[1] != 0 considered as Dirichlet series, computes the quotient of x by y, again as a vector.
The library syntax is GEN dirdiv(GEN x, GEN y).
Computes the Dirichlet series attached to the Euler product of expression expr as p ranges through the primes from a to b. expr must be a polynomial or rational function in another variable than p (say X) and expr(X) is understood as the local factor expr(p-s).
The series is output as a vector of coefficients. If c is omitted, output the first b coefficients of the series; otherwise, output the first c coefficients. The following command computes the sigma function, attached to ζ(s)ζ(s-1):
? direuler(p=2, 10, 1/((1-X)*(1-p*X))) %1 = [1, 3, 4, 7, 6, 12, 8, 15, 13, 18] ? direuler(p=2, 10, 1/((1-X)*(1-p*X)), 5) \\ fewer terms %2 = [1, 3, 4, 7, 6]
Setting c < b is useless (the same effect would be achieved by setting b = c). If c > b, the computed coefficients are "missing" Euler factors:
? direuler(p=2, 10, 1/((1-X)*(1-p*X)), 15) \\ more terms, no longer = sigma ! %3 = [1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 0, 28, 0, 24, 24]
The library syntax is direuler(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b)
x and y being vectors of perhaps different lengths representing the Dirichlet series ∑n xn n-s and ∑n yn n-s, computes the product of x by y, again as a vector.
? dirmul(vector(10,n,1), vector(10,n,moebius(n))) %1 = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
The product
length is the minimum of #x*v(y) and #y*v(x),
where v(x) is the index of the first nonzero coefficient.
? dirmul([0,1], [0,1]); %2 = [0, 0, 0, 1]
The library syntax is GEN dirmul(GEN x, GEN y).
For positive integer N and complex number x, return the sum
f(1)1× + f(2)2× +...+ f(N)N×, where f is a completely
multiplicative function. If f is omitted, return
1× +...+ N×. When N ≤ 0, the function returns 0.
If both is set, return the pair for arguments (x,f) and
(-1-x,f). If both = 2, assume in addition that f is
real-valued (which is true when f is omitted, i.e. represents the constant
function f(n) = 1).
Caveat. when both is set, the present implementation assumes that |f(n)| is either 0 or 1, which is the case for Dirichlet characters.
A vector-valued multiplicative function f is allowed, in which case the above conditions must be met componentwise and the vector length must be constant.
Unlike variants using dirpowers(N,x), this function uses O(sqrt{N})
memory instead of O(N). And it is faster for large N. The return value
is usually a floating point number, but it will be exact if the result
is an integer. On the other hand, rational numbers are converted to
floating point approximations, since they are likely to blow up for large N.
? dirpowers(5, 2) %1 = [1, 4, 9, 16, 25] ? vecsum(%) %2 = 55 ? dirpowerssum(5, 2) %3 = 55 ? dirpowerssum(5, -2) %4 = 1.4636111111111111111111111111111111111 ? \p200 ? s = 1/2 + I * sqrt(3); N = 10^7; ? dirpowerssum(N, s); time = 11,425 ms. ? vecsum(dirpowers(N, s)) time = 19,365 ms. ? dirpowerssum(N, s, n->kronecker(-23,n)) time = 10,981 ms.
The dirpowerssum commands work with default stack size,
the dirpowers one requires a stacksize of at least 5GB.
The library syntax is dirpowerssumfun(ulong N, GEN x, void *E, GEN (*f)(void*, ulong, long), long prec). When f = NULL, one may use
GEN dirpowerssum(ulong N, GEN x, long prec).
Creates a row vector whose components are the
divisors of x. The factorization of x (as output by factor) can
be used instead. If flag = 1, return pairs [d, factor(d)].
By definition, these divisors are the products of the irreducible
factors of n, as produced by factor(n), raised to appropriate
powers (no negative exponent may occur in the factorization). If n is
an integer, they are the positive divisors, in increasing order.
? divisors(12)
%1 = [1, 2, 3, 4, 6, 12]
? divisors(12, 1) \\ include their factorization
%2 = [[1, matrix(0,2)], [2, Mat([2, 1])], [3, Mat([3, 1])],
[4, Mat([2, 2])], [6, [2, 1; 3, 1]], [12, [2, 2; 3, 1]]]
? divisors(x^4 + 2*x^3 + x^2) \\ also works for polynomials
%3 = [1, x, x^2, x + 1, x^2 + x, x^3 + x^2, x^2 + 2*x + 1,
x^3 + 2*x^2 + x, x^4 + 2*x^3 + x^2]
This function requires a lot of memory if x has many divisors. The following idiom runs through all divisors using very little memory, in no particular order this time:
F = factor(x); P = F[,1]; E = F[,2]; forvec(e = vectorv(#E,i,[0,E[i]]), d = factorback(P,e); ...)
If the factorization of d is also desired, then [P,e] almost provides it but not quite: e may contain 0 exponents, which are not allowed in factorizations. These must be sieved out as in:
? tofact(P,E) = matreduce(Mat([P,E])); ? tofact([2,3,5,7]~, [4,0,2,0]~) %4 = [2 4] [5 2]
We can then run the above loop with tofact(P,e) instead of,
or together with, factorback.
The library syntax is GEN divisors0(GEN x, long flag).
The functions GEN divisors(GEN N) (flag = 0) and
GEN divisors_factored(GEN N) (flag = 1) are also available.
Given three integers N > s > r ≥ 0 such that (r,s) = 1 and s3 > N, find all divisors d of N such that d = r (mod s). There are at most 11 such divisors and they are found in polynomial time Õ(log N)^2. (Lenstra).
? N = 245784; r = 19; s = 65 ; ? divisorslenstra(N, r, s) %2 = [19, 84, 539, 1254, 3724, 245784] ? [ d | d <- divisors(N), d % s == r] %3 = [19, 84, 539, 1254, 3724, 245784]
When the preconditions are not met, the result is undefined:
? N = 4484075232; r = 7; s = 1303; s^3 > N %4 = 0 ? divisorslenstra(N, r, s) ? [ d | d <- divisors(N), d % s == r ] %6 = [7, 2613, 9128, 19552, 264516, 3407352, 344928864]
(Divisors were missing but s3 < N.)
The library syntax is GEN divisorslenstra(GEN N, GEN r, GEN s).
Euler's φ (totient) function of the integer |x|, in other words |(ℤ/xℤ)*|.
? eulerphi(40) %1 = 16
According to this definition we let φ(0) := 2, since ℤ* = {-1,1};
this is consistent with znstar(0): we have
znstar(n).no = eulerphi(n) for all n ∈ ℤ.
The library syntax is GEN eulerphi(GEN x).
Factor x over domain D; if D is omitted, it is determined from x. For instance, if x is an integer, it is factored in ℤ, if it is a polynomial with rational coefficients, it is factored in ℚ[x], etc., see below for details. The result is a two-column matrix: the first contains the irreducibles dividing x (rational or Gaussian primes, irreducible polynomials), and the second the exponents. By convention, 0 is factored as 01.
x ∈ ℚ.
See factorint for the algorithms used. The factorization includes the
unit -1 when x < 0 and all other factors are positive; a denominator is
factored with negative exponents. The factors are sorted in increasing order.
? factor(-7/106) %1 = [-1 1] [ 2 -1] [ 7 1] [53 -1]
By convention, 1 is factored as matrix(0,2)
(the empty factorization, printed as [;]).
Large rational "primes" > 264 in the factorization are in fact
pseudoprimes (see ispseudoprime), a priori not rigorously proven
primes. Use isprime to prove primality of these factors, as in
? fa = factor(2^2^7 + 1) %2 = [59649589127497217 1] [5704689200685129054721 1] ? isprime( fa[,1] ) %3 = [1, 1]~ \\ both entries are proven primes
Another possibility is to globally set the default factor_proven, which
will perform a rigorous primality proof for each pseudoprime factor but will
slow down PARI.
A t_INT argument D can be added, meaning that we only trial divide
by all primes p < D and the addprimes entries, then skip all
expensive factorization methods. The limit D must be nonnegative.
In this case, one entry in the factorization may be a composite number: all
factors less than D2 and primes from the addprimes table
are actual primes. But (at most) one entry may not verify this criterion,
and it may be prime or composite: it is only known to be coprime to all
other entries and not a pure power.
? factor(2^2^7 +1, 10^5) %4 = [340282366920938463463374607431768211457 1]
Deprecated feature. Setting D = 0 is the same
as setting it to factorlimit + 1.
This routine uses trial division and perfect power tests, and should not be
used for huge values of D (at most 109, say):
factorint(, 1 + 8) will in general be faster. The latter does not
guarantee that all small prime factors are found, but it also finds larger
factors and in a more efficient way.
? F = (2^2^7 + 1) * 1009 * (10^5+3); factor(F, 10^5) \\ fast, incomplete time = 0 ms. %5 = [1009 1] [34029257539194609161727850866999116450334371 1] ? factor(F, 10^9) \\ slow time = 3,260 ms. %6 = [1009 1] [100003 1] [340282366920938463463374607431768211457 1] ? factorint(F, 1+8) \\ much faster and all small primes were found time = 8 ms. %7 = [1009 1] [100003 1] [340282366920938463463374607431768211457 1] ? factor(F) \\ complete factorization time = 60 ms. %8 = [1009 1] [100003 1] [59649589127497217 1] [5704689200685129054721 1]
x ∈ ℚ(i). The factorization is performed with Gaussian primes in ℤ[i] and includes Gaussian units in {±1, ± i}. Except for a possible leading unit, the Gaussian factors are normalized (rational factors are positive and irrational factors have positive imaginary part) and sorted by increasing imaginary part.
Unless factor_proven is set, large factors are actually pseudoprimes,
not proven primes; a rational factor is prime if less than 264 and an
irrational one if its norm is less than 264.
? factor(5*I) %9 = [ 2 + I 1] [1 + 2*I 1]
One can force the factorization of a rational number by setting the domain D = I:
? factor(-5, I) %10 = [ I 1] [ 2 + I 1] [1 + 2*I 1] ? factorback(%) %11 = -5
Univariate polynomials and rational functions.
PARI can factor univariate polynomials in K[t]. The following base fields
K are currently supported: ℚ, ℝ, ℂ, ℚp, finite fields and
number fields. See factormod and factorff for the algorithms used
over finite fields and nffactor for the algorithms over number fields.
The irreducible factors are sorted by increasing degree and normalized: they
are monic except when K = ℚ where they are primitive in ℤ[t].
The content is not included in the factorization, in particular
factorback will in general recover the original x only up to
multiplication by an element of K*: when K != ℚ, this scalar is
pollead(x) (since irreducible factors are monic); and when K = ℚ
you can either ask for the ℚ-content explicitly of use factorback:
? P = t^2 + 5*t/2 + 1; F = factor(P) %12 = [t + 2 1] [2*t + 1 1] ? content(P, 1) \\ Q-content %13 = 1/2 ? pollead(factorback(F)) / pollead(P) %14 = 2
You can specify K using the optional "domain" argument D as follows
* K = ℚ : D a rational number (t_INT or t_FRAC),
* K = ℤ/pℤ with p prime : D a t_INTMOD modulo p;
factoring modulo a composite number is not supported.
* K = 𝔽q : D a t_FFELT encoding the finite field; you can also
use a t_POLMOD of t_INTMOD modulo a prime p but this is usualy
less convenient;
* K = ℚ[X]/(T) a number field : D a t_POLMOD modulo T,
* K = ℚ(i) (alternate syntax for special case): D = I,
* K = ℚ(w) a quadratic number field (alternate syntax for special
case): D a t_QUAD,
* K = ℝ : D a real number (t_REAL); truncate the factorization
at accuracy precision(D). If x is inexact and precision(x)
is less than precision(D), then the precision of x is used instead.
* K = ℂ : D a complex number with a t_REAL component, e.g.
I * 1.; truncate the factorization as for K = ℝ,
* K = ℚp : D a t_PADIC; truncate the factorization at
p-adic accuracy padicprec(D), possibly less if x is inexact
with insufficient p-adic accuracy;
? T = x^2+1; ? factor(T, 1); \\ over Q ? factor(T, Mod(1,3)) \\ over F3 ? factor(T, ffgen(ffinit(3,2,'t))^0) \\ over F3^2 ? factor(T, Mod(Mod(1,3), t^2+t+2)) \\ over F3^2, again ? factor(T, O(3^6)) \\ over Q3, precision 6 ? factor(T, 1.) \\ over R, current precision ? factor(T, I*1.) \\ over C ? factor(T, Mod(1, y^3-2)) \\ over Q(21/3)
In most cases, it is possible and simpler to call a specialized variant rather than use the above scheme:
? factormod(T, 3) \\ over F3 ? factormod(T, [t^2+t+2, 3]) \\ over F3^2 ? factormod(T, ffgen([3,2], 't)) \\ over F3^2 ? factorpadic(T, 3,6) \\ over Q3, precision 6 ? nffactor(y^3-2, T) \\ over Q(21/3) ? polroots(T) \\ over C ? polrootsreal(T) \\ over R (real polynomial)
It is also possible to let the routine use the smallest field containing all
coefficients, taking into account quotient structures induced by
t_INTMODs and t_POLMODs (e.g. if a coefficient in ℤ/nℤ is known,
all rational numbers encountered are first mapped to ℤ/nℤ; different
moduli will produce an error):
? T = x^2+1; ? factor(T); \\ over Q ? factor(T*Mod(1,3)) \\ over F3 ? factor(T*ffgen(ffinit(3,2,'t))^0) \\ over F3^2 ? factor(T*Mod(Mod(1,3), t^2+t+2)) \\ over F3^2, again ? factor(T*(1 + O(3^6)) \\ over Q3, precision 6 ? factor(T*1.) \\ over R, current precision ? factor(T*(1.+0.*I)) \\ over C ? factor(T*Mod(1, y^3-2)) \\ over Q(21/3)
Multiplying by a suitable field element equal to 1 ∈ K in this way is error-prone and is not recommanded. Factoring existing polynomials with obvious fields of coefficients is fine, the domain argument D should be used instead ad hoc conversions.
Note on inexact polynomials. Polynomials with inexact coefficients (e.g. floating point or p-adic numbers) are first rounded to an exact representation, then factored to (potentially) infinite accuracy and we return a truncated approximation of that virtual factorization. To avoid pitfalls, we advise to only factor exact polynomials:
? factor(x^2-1+O(2^2)) \\ rounded to x^2 + 3, irreducible in Q2 %1 = [(1 + O(2^2))*x^2 + O(2^2)*x + (1 + 2 + O(2^2)) 1] ? factor(x^2-1+O(2^3)) \\ rounded to x^2 + 7, reducible ! %2 = [ (1 + O(2^3))*x + (1 + 2 + O(2^3)) 1] [(1 + O(2^3))*x + (1 + 2^2 + O(2^3)) 1] ? factor(x^2-1, O(2^2)) \\ no ambiguity now %3 = [ (1 + O(2^2))*x + (1 + O(2^2)) 1] [(1 + O(2^2))*x + (1 + 2 + O(2^2)) 1]
Note about inseparable polynomials. Polynomials with inexact
coefficients are considered to be squarefree: indeed, there exist a
squarefree polynomial arbitrarily close to the input, and they cannot be
distinguished at the input accuracy. This means that irreducible factors are
repeated according to their apparent multiplicity. On the contrary, using a
specialized function such as factorpadic with an exact rational
input yields the correct multiplicity when the (now exact) input is not
separable. Compare:
? factor(z^2 + O(5^2))) %1 = [(1 + O(5^2))*z + O(5^2) 1] [(1 + O(5^2))*z + O(5^2) 1] ? factor(z^2, O(5^2)) %2 = [1 + O(5^2))*z + O(5^2) 2]
Multivariate polynomials and rational functions. PARI recursively factors multivariate polynomials in K[t1,..., td] for the same fields K as above and the argument D is used in the same way to specify K. The irreducible factors are sorted by their main variable (least priority first) then by increasing degree.
? factor(x^2 + y^2, Mod(1,5)) %1 = [ x + Mod(2, 5)*y 1] [Mod(1, 5)*x + Mod(3, 5)*y 1] ? factor(x^2 + y^2, O(5^2)) %2 = [ (1 + O(5^2))*x + (O(5^2)*y^2 + (2 + 5 + O(5^2))*y + O(5^2)) 1] [(1 + O(5^2))*x + (O(5^2)*y^2 + (3 + 3*5 + O(5^2))*y + O(5^2)) 1] ? lift(%) %3 = [ x + 7*y 1] [x + 18*y 1]
Note that the implementation does not really support inexact real fields (ℝ or ℂ) and usually misses factors even if the input is exact:
? factor(x^2 + y^2, I) \\ over Q(i) %4 = [x - I*y 1] [x + I*y 1] ? factor(x^2 + y^2, I*1.) \\ over C %5 = [x^2 + y^2 1]
The library syntax is GEN factor0(GEN x, GEN D = NULL).
GEN factor(GEN x)
GEN boundfact(GEN x, ulong lim).
Gives back the factored object corresponding to a factorization. The integer 1 corresponds to the empty factorization.
If e is present, e and f must be vectors of the same length (e being integral), and the corresponding factorization is the product of the f[i]e[i].
If not, and f is vector, it is understood as in the preceding case with e
a vector of 1s: we return the product of the f[i]. Finally, f can be a
regular factorization, as produced with any factor command. A few
examples:
? factor(12) %1 = [2 2] [3 1] ? factorback(%) %2 = 12 ? factorback([2,3], [2,1]) \\ 2^2 * 3^1 %3 = 12 ? factorback([5,2,3]) %4 = 30
The library syntax is GEN factorback2(GEN f, GEN e = NULL).
Also available is GEN factorback(GEN f) (case e = NULL).
This function is obsolete, use factormod.
The library syntax is GEN factmod(GEN x, GEN p).
Obsolete, kept for backward compatibility: use factormod.
The library syntax is GEN factorff(GEN x, GEN p = NULL, GEN a = NULL).
Factorial of x. The expression x! gives a result which is an integer,
while factorial(x) gives a real number.
The library syntax is GEN mpfactr(long x, long prec).
GEN mpfact(long x) returns x! as a t_INT.
Factors the integer n into a product of
pseudoprimes (see ispseudoprime), using a combination of the
Shanks SQUFOF and Pollard Rho method (with modifications due to
Brent), Lenstra's ECM (with modifications by Montgomery), and
MPQS (the latter adapted from the LiDIA code with the kind
permission of the LiDIA maintainers), as well as a search for pure powers.
The output is a two-column matrix as for factor: the first column
contains the "prime" divisors of n, the second one contains the
(positive) exponents.
By convention 0 is factored as 01, and 1 as the empty factorization;
also the divisors are by default not proven primes if they are larger than
264, they only failed the BPSW compositeness test (see
ispseudoprime). Use isprime on the result if you want to
guarantee primality or set the factor_proven default to 1.
Entries of the private prime tables (see addprimes) are also included
as is.
This gives direct access to the integer factoring engine called by most arithmetical functions. flag is optional; its binary digits mean 1: avoid MPQS, 2: skip first stage ECM (we may still fall back to it later), 4: avoid Rho and SQUFOF, 8: don't run final ECM (as a result, a huge composite may be declared to be prime). Note that a (strong) probabilistic primality test is used; thus composites might not be detected, although no example is known.
You are invited to play with the flag settings and watch the internals at
work by using gp's debug default parameter (level 3 shows
just the outline, 4 turns on time keeping, 5 and above show an increasing
amount of internal details).
The library syntax is GEN factorint(GEN x, long flag).
Factors the polynomial f over the finite field defined by the domain D as follows:
* D = p a prime: factor over 𝔽p;
* D = [T,p] for a prime p and T(y) an irreducible polynomial over 𝔽p: factor over 𝔽p[y]/(T) (as usual the main variable of T must have lower priority than the main variable of f);
* D a t_FFELT: factor over the attached field;
* D omitted: factor over the field of definition of f, which must be a finite field.
The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix, the first column being the irreducible polynomials dividing f, and the second the exponents. By convention, the 0 polynomial factors as 01; a nonzero constant polynomial has empty factorization, a 0 x 2 matrix. The irreducible factors are ordered by increasing degree and the result is canonical: it will not change across multiple calls or sessions.
? factormod(x^2 + 1, 3) \\ over F3 %1 = [Mod(1, 3)*x^2 + Mod(1, 3) 1] ? liftall( factormod(x^2 + 1, [t^2+1, 3]) ) \\ over F9 %2 = [ x + t 1] [x + 2*t 1] \\ same, now letting GP choose a model ? T = ffinit(3,2,'t) %3 = Mod(1, 3)*t^2 + Mod(1, 3)*t + Mod(2, 3) ? liftall( factormod(x^2 + 1, [T, 3]) ) %4 = \\ t is a root of T ! [ x + (t + 2) 1] [x + (2*t + 1) 1] ? t = ffgen(t^2+Mod(1,3)); factormod(x^2 + t^0) \\ same using t_FFELT %5 = [ x + t 1] [x + 2*t 1] ? factormod(x^2+Mod(1,3)) %6 = [Mod(1, 3)*x^2 + Mod(1, 3) 1] ? liftall( factormod(x^2 + Mod(Mod(1,3), y^2+1)) ) %7 = [ x + y 1] [x + 2*y 1]
If flag is nonzero, outputs only the degrees of the irreducible polynomials (for example to compute an L-function). By convention, a constant polynomial (including the 0 polynomial) has empty factorization. The degrees appear in increasing order but need not correspond to the ordering with flag = 0 when multiplicities are present.
? f = x^3 + 2*x^2 + x + 2; ? factormod(f, 5) \\ (x+2)^2 * (x+3) %1 = [Mod(1, 5)*x + Mod(2, 5) 2] [Mod(1, 5)*x + Mod(3, 5) 1] ? factormod(f, 5, 1) \\ (deg 1) * (deg 1)^2 %2 = [1 1] [1 2]
The library syntax is GEN factormod0(GEN f, GEN D = NULL, long flag).
Distinct-degree factorization of the squarefree polynomial f over the finite field defined by the domain D as follows:
* D = p a prime: factor over 𝔽p;
* D = [T,p] for a prime p and T an irreducible polynomial over 𝔽p: factor over 𝔽p[x]/(T);
* D a t_FFELT: factor over the attached field;
* D omitted: factor over the field of definition of f, which must be a finite field.
If f is not squarefree, the result is undefined. The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix:
* the first column contains monic (squarefree, pairwise coprime) polynomials dividing f, all of whose irreducible factors have the same degree d;
* the second column contains the degrees of the irreducible factors.
The factorization is ordered by increasing degree d of irreducible factors, and the result is obviously canonical. This function is somewhat faster than full factorization.
? f = (x^2 + 1) * (x^2-1); ? factormodSQF(f,3) \\ squarefree over F3 %2 = [Mod(1, 3)*x^4 + Mod(2, 3) 1] ? factormodDDF(f, 3) %3 = [Mod(1, 3)*x^2 + Mod(2, 3) 1] \\ two degree 1 factors [Mod(1, 3)*x^2 + Mod(1, 3) 2] \\ irred of degree 2 ? for(i=1,10^5,factormodDDF(f,3)) time = 424 ms. ? for(i=1,10^5,factormod(f,3)) \\ full factorization is a little slower time = 464 ms. ? liftall( factormodDDF(x^2 + 1, [3, t^2+1]) ) \\ over F9 %6 = [x^2 + 1 1] \\ product of two degree 1 factors ? t = ffgen(t^2+Mod(1,3)); factormodDDF(x^2 + t^0) \\ same using t_FFELT %7 = [x^2 + 1 1] ? factormodDDF(x^2-Mod(1,3)) %8 = [Mod(1, 3)*x^2 + Mod(2, 3) 1]
The library syntax is GEN factormodDDF(GEN f, GEN D = NULL).
Squarefree factorization of the polynomial f over the finite field defined by the domain D as follows:
* D = p a prime: factor over 𝔽p;
* D = [T,p] for a prime p and T an irreducible polynomial over 𝔽p: factor over 𝔽p[x]/(T);
* D a t_FFELT: factor over the attached field;
* D omitted: factor over the field of definition of f, which must be a finite field.
The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix:
* the first column contains monic squarefree pairwise coprime polynomials dividing f;
* the second column contains the power to which the polynomial in column 1 divides f;
This is somewhat faster than full factorization. The factors are ordered by increasing exponent and the result is obviously canonical.
? f = (x^2 + 1)^3 * (x^2-1)^2; ? factormodSQF(f, 3) \\ over F3 %1 = [Mod(1, 3)*x^2 + Mod(2, 3) 2] [Mod(1, 3)*x^2 + Mod(1, 3) 3] ? for(i=1,10^5,factormodSQF(f,3)) time = 192 ms. ? for(i=1,10^5,factormod(f,3)) \\ full factorization is slower time = 409 ms. ? liftall( factormodSQF((x^2 + 1)^3, [3, t^2+1]) ) \\ over F9 %4 = [x^2 + 1 3] ? t = ffgen(t^2+Mod(1,3)); factormodSQF((x^2 + t^0)^3) \\ same using t_FFELT %5 = [x^2 + 1 3] ? factormodSQF(x^8 + x^7 + x^6 + x^2 + x + Mod(1,2)) %6 = [ Mod(1, 2)*x + Mod(1, 2) 2] [Mod(1, 2)*x^2 + Mod(1, 2)*x + Mod(1, 2) 3]
The library syntax is GEN factormodSQF(GEN f, GEN D = NULL).
'x})
Factors n-th cyclotomic polynomial Φn(x) mod p,
where p is a prime number not dividing n.
Much faster than factormod(polcyclo(n), p); the irreducible
factors should be identical and given in the same order.
If single is set, return a single irreducible factor; else (default)
return all the irreducible factors. Note that repeated calls of this
function with the single flag set may return different results because
the algorithm is probabilistic. Algorithms used are as follows.
Let F = ℚ(ζn). Let K be the splitting field of p in F and e the conductor of K. Then Φn(x) and Φe(x) have the same number of irreducible factors mod p and there is a simple algorithm constructing irreducible factors of Φn(x) from irreducible factors of Φe(x). So we may assume n is equal to the conductor of K. Let d be the order of p in (ℤ/nℤ)× and ϕ(n) = df. Then Φn(x) has f irreducible factors gi(x) (1 ≤ i ≤ f) of degree d over 𝔽p or ℤp.
* If d is small, then we factor gi(x) into d linear factors gij(x), 1 ≤ j ≤ d in 𝔽q[x] (q = pd) and construct Gi(x) = ∏j = 1d gij(x) ∈ 𝔽q[x]. Then Gi(x) ∈ 𝔽p[x] and gi(x) = Gi(x).
* If f is small, then we work in K, which is a Galois extension of degree f over ℚ. The Gaussian period θk = TrF/K(ζnk) is a sum of k-th power of roots of gi(x) and K = ℚ(θ1).
Now, for each k, there is a polynomial Tk(x) ∈ ℚ[x] satisfying θk = Tk(θ1) because all θk are in K. Let T(x) ∈ ℤ[x] be the minimal polynomial of θ1 over ℚ. We get θ1 mod p from T(x) and construct θ1,...,θd mod p using Tk(x). Finally we recover gi(x) from θ1,...,θd by Newton's formula.
? lift(factormodcyclo(15, 11)) %1 = [x^2 + 9*x + 4, x^2 + 4*x + 5, x^2 + 3*x + 9, x^2 + 5*x + 3] ? factormodcyclo(15, 11, 1) \\ single %2 = Mod(1, 11)*x^2 + Mod(5, 11)*x + Mod(3, 11) ? z1 = lift(factormod(polcyclo(12345),11311)[,1]); time = 32,498 ms. ? z2 = factormodcyclo(12345,11311); time = 47 ms. ? z1 == z2 %4 = 1
The library syntax is GEN factormodcyclo(long n, GEN p, long single, long v = -1) where v is a variable number.
Let k, l, m be three finite fields and f a (partial) map from l to m and g a (partial) map from k to l, return the (partial) map f o g from k to m.
a = ffgen([3,5],'a); b = ffgen([3,10],'b); c = ffgen([3,20],'c); m = ffembed(a, b); n = ffembed(b, c); rm = ffinvmap(m); rn = ffinvmap(n); nm = ffcompomap(n,m); ffmap(n,ffmap(m,a)) == ffmap(nm, a) %5 = 1 ffcompomap(rm, rn) == ffinvmap(nm) %6 = 1
The library syntax is GEN ffcompomap(GEN f, GEN g).
Given two finite fields elements a and b, return a map embedding the definition field of a to the definition field of b. Assume that the latter contains the former.
? a = ffgen([3,5],'a); ? b = ffgen([3,10],'b); ? m = ffembed(a, b); ? A = ffmap(m, a); ? minpoly(A) == minpoly(a) %5 = 1
The library syntax is GEN ffembed(GEN a, GEN b).
Extend the field K of definition of a by a root of the polynomial
P ∈ K[X] assumed to be irreducible over K. Return [r, m] where r
is a root of P in the extension field L and m is a map from K to L,
see ffmap.
If v is given, the variable name is used to display the generator of L,
else the name of the variable of P is used.
A generator of L can be recovered using b = ffgen(r).
The image of P in L[X] can be recovered using PL = ffmap(m,P).
? a = ffgen([3,5],'a); ? P = x^2-a; polisirreducible(P) %2 = 1 ? [r,m] = ffextend(a, P, 'b); ? r %3 = b^9+2*b^8+b^7+2*b^6+b^4+1 ? subst(ffmap(m, P), x, r) %4 = 0 ? ffgen(r) %5 = b
The library syntax is GEN ffextend(GEN a, GEN P, long v = -1) where v is a variable number.
Return the n-th power of the Frobenius map over the field of definition of m.
? a = ffgen([3,5],'a); ? f = fffrobenius(a); ? ffmap(f,a) == a^3 %3 = 1 ? g = fffrobenius(a, 5); ? ffmap(g,a) == a %5 = 1 ? h = fffrobenius(a, 2); ? h == ffcompomap(f,f) %7 = 1
The library syntax is GEN fffrobenius(GEN m, long n).
'x})
Return a generator for the finite field k as a t_FFELT.
The field k can be given by
* its order q
* the pair [p,f] where q = pf
* a monic irreducible polynomial with t_INTMOD coefficients modulo a
prime.
* a t_FFELT belonging to k.
If v is given, the variable name is used to display g, else the
variable of the polynomial or the t_FFELT is used, else x is used.
For efficiency, the characteristic is not checked to be prime; similarly
if a polynomial is given, we do not check whether it is irreducible.
When only the order is specified, the function uses the polynomial generated
by ffinit and is deterministic: two calls to the function with the
same parameters will always give the same generator.
To obtain a multiplicative generator, call ffprimroot on the result
(which is randomized). Its minimal polynomial then gives a primitive
polynomial, which can be used to redefine the finite field so that all
subsequent computations use the new primitive polynomial:
? g = ffgen(16, 't); ? g.mod \\ recover the underlying polynomial. %2 = t^4 + t^3 + t^2 + t + 1 ? g.pol \\ lift g as a t_POL %3 = t ? g.p \\ recover the characteristic %4 = 2 ? fforder(g) \\ g is not a multiplicative generator %5 = 5 ? a = ffprimroot(g) \\ recover a multiplicative generator %6 = t^3 + t^2 + t ? fforder(a) %7 = 15 ? T = minpoly(a) \\ primitive polynomial %8 = Mod(1, 2)*x^4 + Mod(1, 2)*x^3 + Mod(1, 2) ? G = ffgen(T); \\ is now a multiplicative generator ? fforder(G) %10 = 15
The library syntax is GEN ffgen(GEN k, long v = -1) where v is a variable number.
To create a generator for a prime finite field, the function
GEN p_to_GEN(GEN p, long v) returns ffgen(p,v)^0.
'x})Computes a monic polynomial of degree n which is irreducible over 𝔽p, where p is assumed to be prime. This function uses a fast variant of Adleman and Lenstra's algorithm.
It is useful in conjunction with ffgen; for instance if
P = ffinit(3,2), you can represent elements in 𝔽32 in term of
g = ffgen(P,'t). This can be abbreviated as
g = ffgen(3^2, 't), where the defining polynomial P can be later
recovered as g.mod.
The library syntax is GEN ffinit(GEN p, long n, long v = -1) where v is a variable number.
m being a map from K to L two finite fields, return the partial map p from L to K such that for all k ∈ K, p(m(k)) = k.
? a = ffgen([3,5],'a); ? b = ffgen([3,10],'b); ? m = ffembed(a, b); ? p = ffinvmap(m); ? u = random(a); ? v = ffmap(m, u); ? ffmap(p, v^2+v+2) == u^2+u+2 %7 = 1 ? ffmap(p, b) %8 = []
The library syntax is GEN ffinvmap(GEN m).
Discrete logarithm of the finite field element x in base g,
i.e. an e in ℤ such that ge = o. If
present, o represents the multiplicative order of g, see
Section se:DLfun; the preferred format for
this parameter is [ord, factor(ord)], where ord is the
order of g. It may be set as a side effect of calling ffprimroot.
The result is undefined if e does not exist. This function uses
* a combination of generic discrete log algorithms (see znlog)
* a cubic sieve index calculus algorithm for large fields of degree at least 5.
* Coppersmith's algorithm for fields of characteristic at most 5.
? t = ffgen(ffinit(7,5)); ? o = fforder(t) %2 = 5602 \\ not a primitive root. ? fflog(t^10,t) %3 = 10 ? fflog(t^10,t, o) %4 = 10 ? g = ffprimroot(t, &o); ? o \\ order is 16806, bundled with its factorization matrix %6 = [16806, [2, 1; 3, 1; 2801, 1]] ? fforder(g, o) %7 = 16806 ? fflog(g^10000, g, o) %8 = 10000
The library syntax is GEN fflog(GEN x, GEN g, GEN o = NULL).
Given a (partial) map m between two finite fields, return the image of x by m. The function is applied recursively to the component of vectors, matrices and polynomials. If m is a partial map that is not defined at x, return [].
? a = ffgen([3,5],'a); ? b = ffgen([3,10],'b); ? m = ffembed(a, b); ? P = x^2+a*x+1; ? Q = ffmap(m,P); ? ffmap(m,poldisc(P)) == poldisc(Q) %6 = 1
The library syntax is GEN ffmap(GEN m, GEN x).
Given a (partial) map m between two finite fields, express x as an algebraic element over the codomain of m in a way which is compatible with m. The function is applied recursively to the component of vectors, matrices and polynomials.
? a = ffgen([3,5],'a); ? b = ffgen([3,10],'b); ? m = ffembed(a, b); ? mi= ffinvmap(m); ? R = ffmaprel(mi,b) %5 = Mod(b,b^2+(a+1)*b+(a^2+2*a+2))
In particular, this function can be used to compute the relative minimal polynomial, norm and trace:
? minpoly(R) %6 = x^2+(a+1)*x+(a^2+2*a+2) ? trace(R) %7 = 2*a+2 ? norm(R) %8 = a^2+2*a+2
The library syntax is GEN ffmaprel(GEN m, GEN x).
Computes the number of monic irreducible polynomials over 𝔽q of degree exactly n (flag = 0 or omitted) or at most n (flag = 1).
The library syntax is GEN ffnbirred0(GEN q, long n, long flag).
Also available are
GEN ffnbirred(GEN q, long n) (for flag = 0)
and GEN ffsumnbirred(GEN q, long n) (for flag = 1).
Multiplicative order of the finite field element x. If o is
present, it represents a multiple of the order of the element,
see Section se:DLfun; the preferred format for
this parameter is [N, factor(N)], where N is the cardinality
of the multiplicative group of the underlying finite field.
? t = ffgen(ffinit(nextprime(10^8), 5)); ? g = ffprimroot(t, &o); \\ o will be useful! ? fforder(g^1000000, o) time = 0 ms. %5 = 5000001750000245000017150000600250008403 ? fforder(g^1000000) time = 16 ms. \\ noticeably slower, same result of course %6 = 5000001750000245000017150000600250008403
The library syntax is GEN fforder(GEN x, GEN o = NULL).
Return a primitive root of the multiplicative
group of the definition field of the finite field element x (not necessarily
the same as the field generated by x). If present, o is set to
a vector [ord, fa], where ord is the order of the group
and fa its factorization factor(ord). This last parameter is
useful in fflog and fforder, see Section se:DLfun.
? t = ffgen(ffinit(nextprime(10^7), 5)); ? g = ffprimroot(t, &o); ? o[1] %3 = 100000950003610006859006516052476098 ? o[2] %4 = [2 1] [7 2] [31 1] [41 1] [67 1] [1523 1] [10498781 1] [15992881 1] [46858913131 1] ? fflog(g^1000000, g, o) time = 1,312 ms. %5 = 1000000
The library syntax is GEN ffprimroot(GEN x, GEN *o = NULL).
Hilbert symbol of x and y modulo the prime p, p = 0 meaning the place at infinity (the result is undefined if p != 0 is not prime).
It is possible to omit p, in which case we take p = 0 if both x
and y are rational, or one of them is a real number. And take p = q
if one of x, y is a t_INTMOD modulo q or a q-adic. (Incompatible
types will raise an error.)
The quaternion algebra (x,y)_ℚ can be created with
alginit(nfinit('t),[x,y]) (for more details see alginit).
The library syntax is long hilbert(GEN x, GEN y, GEN p = NULL).
True (1) if D is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise. D can be input in factored form as for arithmetic functions:
? isfundamental(factor(-8)) %1 = 1 \\ count fundamental discriminants up to 10^8 ? c = 0; forfactored(d = 1, 10^8, if (isfundamental(d), c++)); c time = 40,840 ms. %2 = 30396325 ? c = 0; for(d = 1, 10^8, if (isfundamental(d), c++)); c time = 1min, 33,593 ms. \\ slower ! %3 = 30396325
The library syntax is long isfundamental(GEN D).
True (1) if the integer x is an s-gonal number, false (0) if not.
The parameter s > 2 must be a t_INT. If N is given, set it to n
if x is the n-th s-gonal number.
? ispolygonal(36, 3, &N) %1 = 1 ? N
The library syntax is long ispolygonal(GEN x, GEN s, GEN *N = NULL).
If k is given, returns true (1) if x is a k-th power, false
(0) if not. What it means to be a k-th power depends on the type of
x; see issquare for details.
If k is omitted, only integers and fractions are allowed for x and the
function returns the maximal k ≥ 2 such that x = nk is a perfect
power, or 0 if no such k exist; in particular ispower(-1),
ispower(0), and ispower(1) all return 0.
If a third argument &n is given and x is indeed a k-th power, sets n to a k-th root of x.
For a t_FFELT x, instead of omitting k (which is
not allowed for this type), it may be natural to set
k = (x.p ^ x.f - 1) / fforder(x)
The library syntax is long ispower(GEN x, GEN k = NULL, GEN *n = NULL).
Also available is
long gisanypower(GEN x, GEN *pty) (k omitted).
True (1) if x is a powerful integer, false (0) if not; an integer is powerful if and only if its valuation at all primes dividing x is greater than 1.
? ispowerful(50) %1 = 0 ? ispowerful(100) %2 = 1 ? ispowerful(5^3*(10^1000+1)^2) %3 = 1
The library syntax is long ispowerful(GEN x).
True (1) if x is a prime number, false (0) otherwise. A prime number is a positive integer having exactly two distinct divisors among the natural numbers, namely 1 and itself.
This routine proves or disproves rigorously that a number is prime, which can
be very slow when x is indeed a large prime integer. For instance
a 1000 digits prime should require 15 to 30 minutes with default algorithms.
Use ispseudoprime to quickly check for compositeness. Use
primecert in order to obtain a primality proof instead of a yes/no
answer; see also factor.
The function accepts vector/matrices arguments, and is then applied componentwise.
If flag = 0, use a combination of
* Baillie-Pomerance-Selfridge-Wagstaff compositeness test
(see ispseudoprime),
* Selfridge "p-1" test if x-1 is smooth enough,
* Adleman-Pomerance-Rumely-Cohen-Lenstra (APRCL) for general medium-sized x (less than 1500 bits),
* Atkin-Morain's Elliptic Curve Primality Prover (ECPP) for general large x.
If flag = 1, use Selfridge-Pocklington-Lehmer "p-1" test; this requires partially factoring various auxilliary integers and is likely to be very slow.
If flag = 2, use APRCL only.
If flag = 3, use ECPP only.
The library syntax is GEN gisprime(GEN x, long flag).
If x = pk is a prime power (p prime, k > 0), return k, else return 0. If a second argument &n is given and x is indeed the k-th power of a prime p, sets n to p.
The library syntax is long isprimepower(GEN x, GEN *n = NULL).
True (1) if x is a strong pseudo
prime (see below), false (0) otherwise. If this function returns false, x
is not prime; if, on the other hand it returns true, it is only highly likely
that x is a prime number. Use isprime (which is of course much
slower) to prove that x is indeed prime.
The function accepts vector/matrices arguments, and is then applied
componentwise.
If flag = 0, checks whether x has no small prime divisors (up to 101 included) and is a Baillie-Pomerance-Selfridge-Wagstaff pseudo prime. Such a pseudo prime passes a Rabin-Miller test for base 2, followed by a Lucas test for the sequence (P,1), where P ≥ 3 is the smallest odd integer such that P2 - 4 is not a square mod x. (Technically, we are using an "almost extra strong Lucas test" that checks whether Vn is ± 2, without computing Un.)
There are no known composite numbers passing the above test, although it is
expected that infinitely many such numbers exist. In particular, all
composites ≤ 264 are correctly detected (checked using
https://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html).
If flag > 0, checks whether x is a strong Miller-Rabin pseudo prime for flag randomly chosen bases (with end-matching to catch square roots of -1).
The library syntax is GEN gispseudoprime(GEN x, long flag).
If x = pk is a pseudo-prime power (p pseudo-prime as per
ispseudoprime, k > 0), return k, else
return 0. If a second argument &n is given and x is indeed
the k-th power of a prime p, sets n to p.
More precisely, k is always the largest integer such that x = nk for
some integer n and, when n ≤ 264 the function returns k > 0 if and
only if n is indeed prime. When n > 264 is larger than the threshold,
the function may return 1 even though n is composite: it only passed
an ispseudoprime(n) test.
The library syntax is long ispseudoprimepower(GEN x, GEN *n = NULL).
True (1) if x is a square, false (0)
if not. What "being a square" means depends on the type of x: all
t_COMPLEX are squares, as well as all nonnegative t_REAL; for
exact types such as t_INT, t_FRAC and t_INTMOD, squares are
numbers of the form s2 with s in ℤ, ℚ and ℤ/Nℤ respectively.
? issquare(3) \\ as an integer %1 = 0 ? issquare(3.) \\ as a real number %2 = 1 ? issquare(Mod(7, 8)) \\ in Z/8Z %3 = 0 ? issquare( 5 + O(13^4) ) \\ in Q_13 %4 = 0
If n is given, a square root of x is put into n.
? issquare(4, &n) %1 = 1 ? n %2 = 2
For polynomials, either we detect that the characteristic is 2 (and check directly odd and even-power monomials) or we assume that 2 is invertible and check whether squaring the truncated power series for the square root yields the original input.
For t_POLMOD x, we only support t_POLMODs of t_INTMODs
encoding finite fields, assuming without checking that the intmod modulus
p is prime and that the polmod modulus is irreducible modulo p.
? issquare(Mod(Mod(2,3), x^2+1), &n) %1 = 1 ? n %2 = Mod(Mod(2, 3)*x, Mod(1, 3)*x^2 + Mod(1, 3))
The library syntax is long issquareall(GEN x, GEN *n = NULL).
Also available is long issquare(GEN x). Deprecated
GP-specific functions GEN gissquare(GEN x) and
GEN gissquareall(GEN x, GEN *pt) return gen0 and gen1
instead of a boolean value.
True (1) if x is squarefree, false (0) if not. Here x can be an integer or a polynomial with coefficients in an integral domain.
? issquarefree(12) %1 = 0 ? issquarefree(6) %2 = 1 ? issquarefree(x^3+x^2) %3 = 0 ? issquarefree(Mod(1,4)*(x^2+x+1)) \\ Z/4Z is not a domain ! *** at top-level: issquarefree(Mod(1,4)*(x^2+x+1)) *** ^ — — — — — — — — — — -- *** issquarefree: impossible inverse in Fp_inv: Mod(2, 4).
A polynomial is declared squarefree if gcd(x,x') is
1. In particular a nonzero polynomial with inexact coefficients is
considered to be squarefree. Note that this may be inconsistent with
factor, which first rounds the input to some exact approximation before
factoring in the apropriate domain; this is correct when the input is not
close to an inseparable polynomial (the resultant of x and x' is not
close to 0).
An integer can be input in factored form as in arithmetic functions.
? issquarefree(factor(6)) %1 = 1 \\ count squarefree integers up to 10^8 ? c = 0; for(d = 1, 10^8, if (issquarefree(d), c++)); c time = 3min, 2,590 ms. %2 = 60792694 ? c = 0; forfactored(d = 1, 10^8, if (issquarefree(d), c++)); c time = 45,348 ms. \\ faster ! %3 = 60792694
The library syntax is long issquarefree(GEN x).
True (1) if x = φ(n) for some integer n, false (0) if not.
? istotient(14) %1 = 0 ? istotient(100) %2 = 0
If N is given, set N = n as well.
? istotient(4, &n) %1 = 1 ? n %2 = 10
The library syntax is long istotient(GEN x, GEN *N = NULL).
Return the largest non-negative integer e so that be ≤ x, where b > 1 is an integer and x ≥ 1 is a real number. If the parameter z is present, set it to be.
? logint(1000, 2) %1 = 9 ? 2^9 %2 = 512 ? logint(1000, 2, &z) %3 = 9 ? z %4 = 512 ? logint(Pi^2, 2, &z) %5 = 3 ? z %6 = 8
The number of digits used to write x in base b is
1 + logint(x,b):
? #digits(1000!, 10) %5 = 2568 ? logint(1000!, 10) %6 = 2567
This function may conveniently replace
floor( log(x) / log(b) )
which may not give the correct answer since PARI does not guarantee exact rounding.
The library syntax is long logint0(GEN x, GEN b, GEN *z = NULL).
Moebius μ-function of |x|; x must be a nonzero integer.
The library syntax is long moebius(GEN x).
Finds the smallest pseudoprime (see
ispseudoprime) greater than or equal to x. x can be of any real
type. Note that if x is a pseudoprime, this function returns x and not
the smallest pseudoprime strictly larger than x. To rigorously prove that
the result is prime, use isprime.
? nextprime(2) %1 = 2 ? nextprime(Pi) %2 = 5 ? nextprime(-10) %3 = 2 \\ primes are positive
If the optional parameter q is an integer, return the next prime congruent to 1 mod q; if q is an intmod, return the next prime in the given congruence class. The function will raise an exception if the class is not invertible.
? nextprime(2, Mod(4,7)) %4 = 11 ? nextprime(90, 25) %5 = 101
Despite the name, please note that the function is not guaranteed to return
a prime number, although no counter-example is known at present. The return
value is a guaranteed prime if x ≤ 264. To rigorously prove
that the result is prime in all cases, use isprime.
The library syntax is GEN nextprime0(GEN x, GEN q = NULL).
Also available is GEN nextprime(GEN N).
Number of divisors of |x|. x must be of type integer.
The library syntax is GEN numdiv(GEN x).
Number of distinct prime divisors of |x|. x must be of type integer.
? factor(392) %1 = [2 3] [7 2] ? omega(392) %2 = 2; \\ without multiplicity ? bigomega(392) %3 = 5; \\ = 3+2, with multiplicity
The library syntax is long omega(GEN x).
Finds the largest pseudoprime (see ispseudoprime) less than or equal
to x; the input x can be of any real type.
Returns 0 if x ≤ 1. Note that if x is a prime, this function returns x
and not the largest prime strictly smaller than x.
? precprime(2) %1 = 2 ? precprime(Pi) %2 = 3 ? precprime(-10) %3 = 0 \\ primes are positive
The function name comes from preceding prime.
Despite the name, please note that the function is not guaranteed to return
a prime number (although no counter-example is known at present); the return
value is a guaranteed prime if x ≤ 264. To rigorously prove
that the result is prime in all cases, use isprime.
The library syntax is GEN precprime(GEN x).
The n-th prime number
? prime(10^9) %1 = 22801763489
Uses checkpointing and a naive O(n) algorithm. Will need
several minutes for n up to 1011; make sure to start gp with
primelimit at least sqrt{pn}, e.g. the value
sqrt{nlog (nlog n)} is guaranteed to be sufficient; failing this condition
can result in a slowdown of a factor at least 30.
The running time varies a lot depending on the distance to the nearest
checkpoint; the largest checkpoint is currently 1011.
The library syntax is GEN prime(long n).
If N is a prime, return a PARI Primality Certificate for the prime N,
as described below. Otherwise, return 0. A Primality Certificate
c can be checked using primecertisvalid(c).
If flag = 0 (default), return an ECPP certificate (Atkin-Morain)
If flag = 0 and partial > 0, return a (potentially) partial ECPP certificate.
A PARI ECPP Primality Certificate for the prime N is either a prime
integer N < 264 or a vector C of length ℓ whose ith
component C[i] is a vector [Ni, ti, si, ai, Pi]
of length 5
where N1 = N. It is said to be valid if for each
i = 1,..., ℓ, all of the following conditions are satisfied
* Ni is a positive integer
* ti is an integer such that ti2 < 4Ni
* si is a positive integer which divides mi where mi = Ni + 1 - ti
* If we set qi = (mi)/(si), then
* qi > (Ni1/4+1)2
* qi = Ni+1 if 1 ≤ i < l
* qℓ ≤ 264 is prime
* ai is an integer
* P[i] is a vector of length 2 representing the affine
point Pi = (xi, yi) on the elliptic curve
E: y2 = x3 + aix + bi modulo Ni where
bi = yi2 - xi3 - aixi satisfying the following:
* mi Pi = oo
* si Pi != oo
Using the following theorem, the data in the vector C allows to
recursively certify the primality of N (and all the qi) under the single
assumption that qℓ be prime.
Theorem. If N is an integer and there exist positive integers m, q and a point P on the elliptic curve E: y2 = x3 + ax + b defined modulo N such that q > (N1/4 + 1)2, q is a prime divisor of m, mP = oo and (m)/(q)P != oo , then N is prime.
A partial certificate is identical except that the condition qℓ ≤ 264 is replaced by qℓ ≤ 2partial. Such partial certificate C can be extended to a full certificate by calling C = primecert(C), or to a longer partial certificate by calling C = primecert(C,,b) with b < partial.
? primecert(10^35 + 69) %1 = [[100000000000000000000000000000000069, 5468679110354 52074, 2963504668391148, 0, [60737979324046450274283740674 208692, 24368673584839493121227731392450025]], [3374383076 4501150277, -11610830419, 734208843, 0, [26740412374402652 72 4, 6367191119818901665]], [45959444779, 299597, 2331, 0 , [18022351516, 9326882 51]]] ? primecert(nextprime(2^64)) %2 = [[18446744073709551629, -8423788454, 160388, 1, [1059 8342506117936052, 2225259013356795550]]] ? primecert(6) %3 = 0 ? primecert(41) %4 = 41 ? N = 2^2000+841; ? Cp1 = primecert(N,,1500); \\ partial certificate time = 16,018 ms. ? Cp2 = primecert(Cp1,,1000); \\ (longer) partial certificate time = 5,890 ms. ? C = primecert(Cp2); \\ full certificate for N time = 1,777 ms. ? primecertisvalid(C) %9 = 1 ? primecert(N); time = 23,625 ms.
As the last command shows, attempting a succession of partial certificates should be about as fast as a direct computation.
If flag = 1 (very slow), return an N-1 certificate (Pocklington Lehmer)
A PARI N-1 Primality Certificate for the prime N is either a prime integer N < 264 or a pair [N, C], where C is a vector with ℓ elements which are either a single integer pi < 264 or a triple [pi,ai,Ci] with pi > 264 satisfying the following properties:
* pi is a prime divisor of N - 1;
* ai is an integer such that aiN-1 = 1 (mod N) and ai(N-1)/pi - 1 is coprime with N;
* Ci is an N-1 Primality Certificate for pi
* The product F of the pivp_{i(N-1)} is strictly larger than N1/3. Provided that all pi are indeed primes, this implies that any divisor of N is congruent to 1 modulo F.
* The Brillhart-Lehmer-Selfridge criterion is satisfied: when we write N = 1 + c1 F + c2 F2 in base F the polynomial 1 + c1 X + c2 X2 is irreducible over ℤ, i.e. c12 - 4c2 is not a square. This implies that N is prime.
This algorithm requires factoring partially p-1 for various prime integers p with an unfactored parted ≤ p2/3 and this may be exceedingly slow compared to the default.
The algorithm fails if one of the pseudo-prime factors is not prime, which is
exceedingly unlikely and well worth a bug report. Note that if you monitor
the algorithm at a high enough debug level, you may see warnings about
untested integers being declared primes. This is normal: we ask for partial
factorizations (sufficient to prove primality if the unfactored part is not
too large), and factor warns us that the cofactor hasn't been tested.
It may or may not be tested later, and may or may not be prime. This does
not affect the validity of the whole Primality Certificate.
The library syntax is GEN primecert0(GEN N, long flag, long partial).
Also available is
GEN ecpp0(GEN N, long partial) (flag = 0).
Returns a string suitable for print/write to display a primality certificate
from primecert, the format of which depends on the value of format:
* 0 (default): Human-readable format. See ??primecert for the
meaning of the successive N, t, s, a, m, q, E, P. The integer D is the
negative fundamental discriminant coredisc(t2 - 4N).
* 1: Primo format 4.
* 2: MAGMA format.
Currently, only ECPP Primality Certificates are supported.
? cert = primecert(10^35+69);
? s = primecertexport(cert); \\ Human-readable
? print(s)
[1]
N = 100000000000000000000000000000000069
t = 546867911035452074
s = 2963504668391148
a = 0
D = -3
m = 99999999999999999453132088964547996
q = 33743830764501150277
E = [0, 1]
P = [21567861682493263464353543707814204,
49167839501923147849639425291163552]
[2]
N = 33743830764501150277
t = -11610830419
s = 734208843
a = 0
D = -3
m = 33743830776111980697
q = 45959444779
E = [0, 25895956964997806805]
P = [29257172487394218479, 3678591960085668324]
\\ Primo format
? s = primecertexport(cert,1); write("cert.out", s);
\\ Magma format, write to file
? s = primecertexport(cert,2); write("cert.m", s);
? cert = primecert(10^35+69, 1); \\ N-1 certificate
? primecertexport(cert)
*** at top-level: primecertexport(cert)
*** ^ — — — — — — —
*** primecertexport: sorry, N-1 certificate is not yet implemented.
The library syntax is GEN primecertexport(GEN cert, long format).
Verifies if cert is a valid PARI ECPP Primality certificate, as described
in ??primecert.
? cert = primecert(10^35 + 69) %1 = [[100000000000000000000000000000000069, 5468679110354 52074, 2963504668391148, 0, [60737979324046450274283740674 208692, 24368673584839493121227731392450025]], [3374383076 4501150277, -11610830419, 734208843, 0, [26740412374402652 72 4, 6367191119818901665]], [45959444779, 299597, 2331, 0 , [18022351516, 9326882 51]]] ? primecertisvalid(cert) %2 = 1 ? cert[1][1]++; \\ random perturbation ? primecertisvalid(cert) %4 = 0 \\ no longer valid ? primecertisvalid(primecert(6)) %5 = 0
The library syntax is long primecertisvalid(GEN cert).
True (1) if the prime p is regular, false (0) otherwise. Raise an exception if p is not prime.
A regular prime (in the sense of Kummer) is such that p does not divide the class number of the p-th cyclotomic field ℚ(ζp), in particular 2 is considered to be regular. This function runs in time Õ(p).
The first irregular primes:
? [ p | p <- primes([2, 500]), !primeisregular(p) ] %1 = [37, 59, 67, 101, 103, 131, 149, 157, 233, 257, 263, 271, 283, 293, 307, 311, 347, 353, 379, 389, 401, 409, 421, 433, 461, 463, 467, 491]
The library syntax is long primeisregular(ulong p).
The prime counting function. Returns the number of primes p, p ≤ x.
? primepi(10) %1 = 4; ? primes(5) %2 = [2, 3, 5, 7, 11] ? primepi(10^11) %3 = 4118054813
Uses checkpointing and a naive O(x) algorithm;
make sure to start gp with primelimit at least sqrt{x}; failing this
condition can result in a slowdown of a factor at least 30.
The running time varies a lot depending on the distance to the nearest
checkpoint; the largest checkpoint is currently 2760727302517 ~ 241.
The library syntax is GEN primepi(GEN x).
Creates a row vector whose components are the first n prime numbers.
(Returns the empty vector for n ≤ 0.) A t_VEC n = [a,b] is also
allowed, in which case the primes in [a,b] are returned
? primes(10) \\ the first 10 primes %1 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] ? primes([0,29]) \\ the primes up to 29 %2 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] ? primes([15,30]) %3 = [17, 19, 23, 29]
The library syntax is GEN primes0(GEN n).
Compute the value of Ramanujan's tau function at an individual n,
assuming the truth of the GRH (to compute quickly class numbers of imaginary
quadratic fields using quadclassunit). If ell is 16, 18, 20, 22,
or 26, same for the newform of level 1 and corresponding weight. Otherwise,
compute the coefficient of the trace form at n.
The complexity is in Õ(n1/2) using O(log n) space.
If all values up to N are required, then
∑ τ(n)qn = q ∏n ≥ 1 (1-qn)24
and more generally, setting u = ℓ - 13 and C = 2/ζ(-u) for ℓ
> 12,
∑τℓ(n)qn = q ∏n ≥ 1
(1-qn)24 ( 1 + C∑n ≥ 1nu qn / (1-qn))
produces them in time Õ(N), against Õ(N3/2) for
individual calls to ramanujantau; of course the space complexity then
becomes Õ(N). For other values of ℓ,
mfcoefs(mftraceform([1,ell]),N) is much faster.
? tauvec(N) = Vec(q*eta(q + O(q^N))^24); ? N = 10^4; v = tauvec(N); time = 26 ms. ? ramanujantau(N) %3 = -482606811957501440000 ? w = vector(N, n, ramanujantau(n)); \\ much slower ! time = 13,190 ms. ? v == w %4 = 1
The library syntax is GEN ramanujantau(GEN n, long ell).
Returns a strong pseudo prime (see ispseudoprime) in [2,N-1].
A t_VEC N = [a,b] is also allowed, with a ≤ b in which case a
pseudo prime a ≤ p ≤ b is returned; if no prime exists in the
interval, the function will run into an infinite loop. If the upper bound
is less than 264 the pseudo prime returned is a proven prime.
? randomprime(100) %1 = 71 ? randomprime([3,100]) %2 = 61 ? randomprime([1,1]) *** at top-level: randomprime([1,1]) *** ^ — — — — — — *** randomprime: domain error in randomprime: *** floor(b) - max(ceil(a),2) < 0 ? randomprime([24,28]) \\ infinite loop
If the optional parameter q is an integer, return a prime congruent to 1 mod q; if q is an intmod, return a prime in the given congruence class. If the class contains no prime in the given interval, the function will raise an exception if the class is not invertible, else run into an infinite loop
? randomprime(100, 4) \\ 1 mod 4
%1 = 71
? randomprime(100, 4)
%2 = 13
? randomprime([10,100], Mod(2,5))
%3 = 47
? randomprime(100, Mod(0,2)) \\ silly but works
%4 = 2
? randomprime([3,100], Mod(0,2)) \\ not invertible
*** at top-level: randomprime([3,100],Mod(0,2))
*** ^ — — — — — — — — — --
*** randomprime: elements not coprime in randomprime:
0
2
? randomprime(100, 97) \\ infinite loop
The library syntax is GEN randomprime0(GEN N = NULL, GEN q = NULL).
Also available is GEN randomprime(GEN N = NULL).
Removes the primes listed in x from
the prime number table. In particular removeprimes(addprimes()) empties
the extra prime table. x can also be a single integer. List the current
extra primes if x is omitted.
The library syntax is GEN removeprimes(GEN x = NULL).
Sum of the k-th powers of the positive divisors of |x|. x and k must be of type integer.
The library syntax is GEN sumdivk(GEN x, long k).
Also available is GEN sumdiv(GEN n), for k = 1.
Returns the integer square root of x, i.e. the largest integer y such that y2 ≤ x, where x a nonnegative real number. If r is present, set it to the remainder r = x - y2, which satisfies 0 ≤ r < 2y + 1. Further, when x is an integer, r is an integer satisfying 0 ≤ r ≤ 2y.
? x = 120938191237; sqrtint(x) %1 = 347761 ? sqrt(x) %2 = 347761.68741970412747602130964414095216 ? y = sqrtint(x, &r); r %3 = 478116 ? x - y^2 %4 = 478116 ? sqrtint(9/4, &r) \\ not 3/2 ! %5 = 1 ? r %6 = 5/4
The library syntax is GEN sqrtint0(GEN x, GEN *r = NULL).
Also available is GEN sqrtint(GEN a).
Returns the integer n-th root of x, i.e. the largest integer y such that yn ≤ x, where x is a nonnegative real number.
? N = 120938191237; sqrtnint(N, 5) %1 = 164 ? N^(1/5) %2 = 164.63140849829660842958614676939677391 ? sqrtnint(Pi^2, 3) %3 = 2
The special case n = 2 is sqrtint
The library syntax is GEN sqrtnint(GEN x, long n).
Returns the Dedekind sum attached to the integers h and k, corresponding to a fast implementation of
s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
The library syntax is GEN sumdedekind(GEN h, GEN k).
Sum of digits in the integer n, when written in base B.
? sumdigits(123456789) %1 = 45 ? sumdigits(123456789, 2) %2 = 16 ? sumdigits(123456789, -2) %3 = 15
Note that the sum of bits in n is also returned by
hammingweight. This function is much faster than
vecsum(digits(n,B)) when B is 10 or a power of 2, and only
slightly faster in other cases.
The library syntax is GEN sumdigits0(GEN n, GEN B = NULL).
Also available is GEN sumdigits(GEN n), for B = 10.