Pari/GP Reference Documentation | Contents - Global index - GP keyboard shortcuts |
Arithmetic functions and the factoring engine Dirichlet characters Orders in finite groups and Discrete Logarithm functions addprimes bestappr bestapprPade bezout bigomega binomial charconj chardiv chareval charker charmul charorder chinese content contfrac contfracpnqn core coredisc dirdiv direuler dirmul divisors eulerphi factor factorback factorcantor factorff factorial factorint factormod ffgen ffinit fflog ffnbirred fforder ffprimroot fibonacci gcd gcdext hilbert isfundamental ispolygonal ispower ispowerful isprime isprimepower ispseudoprime ispseudoprimepower issquare issquarefree istotient kronecker lcm logint moebius nextprime numbpart numdiv omega partitions polrootsff precprime prime primepi primes qfbclassno qfbcompraw qfbhclassno qfbnucomp qfbnupow qfbpowraw qfbprimeform qfbred qfbredsl2 qfbsolve quadclassunit quaddisc quadgen quadhilbert quadpoly quadray quadregulator quadunit ramanujantau randomprime removeprimes sigma sqrtint sqrtnint stirling sumdedekind sumdigits zncharinduce zncharisodd znchartokronecker znconreychar znconreyconductor znconreyexp znconreylog zncoppersmith znlog znorder znprimroot znstar | |
These functions are by definition functions whose natural domain of definition is either ℤ (or ℤ_{ > 0}). The way these functions are used is completely different from transcendental functions in that there are no automatic type conversions: in general only integers are accepted as arguments. An integer argument N can be given in the following alternate formats:
*
* 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
| |
Arithmetic functions and the factoring engine | |
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 etc.--- call, after trial
division by small primes, the same versatile factoring machinery described
under
default(factor_proven, 1) to ensure that all tentative factorizations are fully proven. This should not slow down PARI too much, unless prime numbers with hundreds of decimal digits occur frequently in your application.
| |
Orders in finite groups and Discrete Logarithm functions | |
The following functions compute the order of an element in a finite group:
All such functions allow an optional argument specifying an integer N, representing the order of the group. (The order functions also allows any non-zero multiple of the order, with a minor loss of efficiency.) That optional argument follows the same format as given above:
*
*
* When the group is fixed and many orders or discrete logarithms will be computed, it is much more efficient to initialize this data once and for all and pass it to the relevant functions, as in
? p = nextprime(10^40); ? v = [p-1, factor(p-1)]; \\ data for discrete log & order computations ? znorder(Mod(2,p), v) %3 = 500000000000000000000000000028 ? g = znprimroot(p); ? znlog(2, g, v) %5 = 543038070904014908801878611374
| |
Dirichlet characters | |
The finite abelian group G = (ℤ/Nℤ)^* can be written G = ⨁ _{i ≤ n} (ℤ/d_iℤ) g_i, with d_n | ... | d_2 | d_1 (SNF condition), all d_i > 0, and ∏_i d_i = φ(N). The SNF condition makes the d_i unique, but the generators g_i, of respective order d_i, are definitely not unique. The ⨁ notation means that all elements of G can be written uniquely as ∏_i g_i^{n_i} where n_i ∈ ℤ/d_iℤ. The g_i are the so-called SNF generators of G. * a character on the abelian group ⨁ (ℤ/d_jℤ) g_j is given by a row vector χ = [a_1,...,a_n] of integers 0 ≤ a_i < d_i such that χ(g_j) = e(a_j / d_j) for all j, with the standard notation e(x) := exp(2iπ x). In other words, χ(∏ g_j^{n_j}) = e(∑ a_j n_j / d_j). This will be generalized to more general abelian groups in later sections (Hecke characters), but in the present case of (ℤ/Nℤ)^*, there is a useful alternate convention : namely, it is not necessary to impose the SNF condition and we can use Chinese reminders instead. If N = ∏ p^{e_p} is the factorization of N into primes, the so-called Conrey generators of G are the generators of the (ℤ/p^{e_p}ℤ)^* lifted to (ℤ/Nℤ)^* by requesting that they be congruent to 1 modulo N/p^{e_p} (for p odd we take the smallest positive primitive root, and for p = 2 we take -1 if e_2 > 1 and additionally 5 if e_2 > 2). We can again write G = ⨁ _{i ≤ n} (ℤ/D_iℤ) G_i, where again ∏_i D_i = φ(N). These generators don't satisfy the SNF condition in general since their orders are now (p-1)p^{e_p-1} for p odd; for p = 2, the generator -1 has order 2 and 5 has order 2^{e_2-2} (e_2 > 2). Nevertheless, any m ∈ (ℤ/Nℤ)^* can be uniquely decomposed as ∏ G_i^{m_i} for some m_i modulo D_i and we can define a character by χ(G_j) = e(m_j / D_j) for all j.
* The column vector of the m_j, 0 ≤ m_j < D_j is called the
Conrey logarithm of m (discrete logarithm in terms of the Conrey
generators). Note that discrete logarithms in PARI/GP are always expressed as
* The attached character is called the Conrey character attached to m.
To sum up a Dirichlet character can be defined by a Concretely, this works as follows:
Also available are
| |
addprimes | |
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
The entries in x must be primes: there is no internal check, even if
the
The library syntax is
| |
bestappr | |
Using variants of the extended Euclidean algorithm, returns a rational approximation a/b to x, whose denominator is limited by B, if present. If B is omitted, return the best approximation affordable given the input accuracy; if you are looking for true rational numbers, presumably approximated to sufficient accuracy, you should first try that option. Otherwise, B must be a positive real scalar (impose 0 < b ≤ B).
* If x is a
? bestappr(Pi, 100) %1 = 22/7 ? bestappr(0.1428571428571428571428571429) %2 = 1/7 ? bestappr([Pi, sqrt(2) + 'x], 10^3) %3 = [355/113, x + 1393/985] By definition, a/b is the best rational approximation to x if |b x - a| < |v x - u| for all integers (u,v) with 0 < v ≤ B. (Which implies that n/d is a convergent of the continued fraction of x.)
* If x is a
? bestappr(Mod(18526731858, 11^10)) %1 = 1/7 ? bestappr(Mod(18526731858, 11^20)) %2 = [] ? bestappr(3 + 5 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + O(5^7)) %2 = -1/3 In most concrete uses, B is a prime power and we performed Hensel lifting to obtain x. The function applies recursively to components of complex objects (polynomials, vectors,...). If rational reconstruction fails for even a single entry, return [].
The library syntax is
| |
bestapprPade | |
Using variants of the extended Euclidean algorithm, returns a rational function approximation a/b to x, whose denominator is limited by B, if present. If B is omitted, return the best approximation affordable given the input accuracy; if you are looking for true rational functions, presumably approximated to sufficient accuracy, you should first try that option. Otherwise, B must be a non-negative real (impose 0 ≤ degree(b) ≤ B).
* If x is a
? bestapprPade((1-x^11)/(1-x)+O(x^11)) %1 = 1/(-x + 1) ? bestapprPade([1/(1+x+O(x^10)), (x^3-2)/(x^3+1)], 1) %2 = [1/(x + 1), -2]
* If x is a
? bestapprPade(Mod(1+x+x^2+x^3+x^4, x^4-2)) %1 = (2*x - 1)/(x - 1) ? % * Mod(1,x^4-2) %2 = Mod(x^3 + x^2 + x + 3, x^4 - 2) ? bestapprPade(Mod(1+x+x^2+x^3+x^5, x^9)) %2 = [] ? bestapprPade(Mod(1+x+x^2+x^3+x^5, x^10)) %3 = (2*x^4 + x^3 - x - 1)/(-x^5 + x^3 + x^2 - 1) The function applies recursively to components of complex objects (polynomials, vectors,...). If rational reconstruction fails for even a single entry, return [].
The library syntax is
| |
bezout | |
Deprecated alias for
The library syntax is
| |
bigomega | |
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
| |
binomial | |
binomial coefficient binom{x}{y}. Here y must be an integer, but x can be any PARI object.
The library syntax is
| |
charconj | |
Let cyc represent a finite abelian group by its elementary
divisors, i.e. (d_j) represents ∑_{j ≤ k} ℤ/d_jℤ with d_k
| ... | d_1; any object which has a This function returns the conjugate character.
? cyc = [15,5]; chi = [1,1]; ? charconj(cyc, chi) %2 = [14, 4] ? bnf = bnfinit(x^2+23); ? bnf.cyc %4 = [3] ? charconj(bnf, [1]) %5 = [2]
For Dirichlet characters (when
? G = idealstar(,8); \\ (Z/8Z)^* ? charorder(G, 3) \\ Conrey label %2 = 2 ? chi = znconreylog(G, 3); ? charorder(G, chi) \\ Conrey logarithm %4 = 2
The library syntax is
| |
chardiv | |
Let cyc represent a finite abelian group by its elementary
divisors, i.e. (d_j) represents ∑_{j ≤ k} ℤ/d_jℤ with d_k
| ... | d_1; any object which has a Given two characters a and b, return the character a / b = a b.
? cyc = [15,5]; a = [1,1]; b = [2,4]; ? chardiv(cyc, a,b) %2 = [14, 2] ? bnf = bnfinit(x^2+23); ? bnf.cyc %4 = [3] ? chardiv(bnf, [1], [2]) %5 = [2]
For Dirichlet characters on (ℤ/Nℤ)^*, additional
representations are available (Conrey labels, Conrey logarithm),
see Section se:dirichletchar or
? G = idealstar(,100); ? G.cyc %2 = [20, 2] ? a = [10, 1]; \\ usual representation for characters ? b = 7; \\ Conrey label; ? c = znconreylog(G, 11); \\ Conrey log ? chardiv(G, b,b) %6 = 1 \\ Conrey label ? chardiv(G, a,b) %7 = [0, 5]~ \\ Conrey log ? chardiv(G, a,c) %7 = [0, 14]~ \\ Conrey log
The library syntax is
| |
chareval | |
Let G be an abelian group structure affording a discrete logarithm
method, e.g G = Note on characters. Let K be some field. If G is an abelian group, let χ: G → K^* be a character of finite order and let o be a multiple of the character order such that χ(n) = ζ^{c(n)} for some fixed ζ ∈ K^* of multiplicative order o and a unique morphism c: G → (ℤ/oℤ,+). Our usual convention is to write G = (ℤ/o_1ℤ) g_1 ⨁ ...⨁ (ℤ/o_dℤ) g_d for some generators (g_i) of respective order d_i, where the group has exponent o := lcm_i o_i. Since ζ^o = 1, the vector (c_i) in ∏ (ℤ/o_iℤ) defines a character χ on G via χ(g_i) = ζ^{c_i (o/o_i)} for all i. Classical Dirichlet characters have values in K = ℂ and we can take ζ = exp(2iπ/o).
Note on Dirichlet characters.
In the special case where bid is attached to G = (ℤ/qℤ)^*
(as per The character value is encoded as follows, depending on the optional argument z: * If z is omitted: return the rational number c(x)/o for x coprime to q, where we normalize 0 ≤ c(x) < o. If x can not be mapped to the group (e.g. x is not coprime to the conductor of a Dirichlet or Hecke character) we return the sentinel value -1. * If z is an integer o, then we assume that o is a multiple of the character order and we return the integer c(x) when x belongs to the group, and the sentinel value -1 otherwise. * z can be of the form [zeta, o], where zeta is an o-th root of 1 and o is a multiple of the character order. We return ζ^{c(x)} if x belongs to the group, and the sentinel value 0 otherwise. (Note that this coincides with the usual extension of Dirichlet characters to ℤ, or of Hecke characters to general ideals.) * Finally, z can be of the form [vzeta, o], where vzeta is a vector of powers ζ^0,..., ζ^{o-1} of some o-th root of 1 and o is a multiple of the character order. As above, we return ζ^{c(x)} after a table lookup. Or the sentinel value 0.
The library syntax is
| |
charker | |
Let cyc represent a finite abelian group by its elementary
divisors, i.e. (d_j) represents ∑_{j ≤ k} ℤ/d_jℤ with d_k
| ... | d_1; any object which has a
This function returns the kernel of χ, as a matrix K in HNF which is a
left-divisor of
? cyc = [15,5]; chi = [1,1]; ? charker(cyc, chi) %2 = [15 12] [ 0 1] ? bnf = bnfinit(x^2+23); ? bnf.cyc %4 = [3] ? charker(bnf, [1]) %5 = [3]
Note that for Dirichlet characters (when
? G = idealstar(,8); \\ (Z/8Z)^* ? charker(G, 1) \\ Conrey label for trivial character %2 = [1 0] [0 1]
The library syntax is
| |
charmul | |
Let cyc represent a finite abelian group by its elementary
divisors, i.e. (d_j) represents ∑_{j ≤ k} ℤ/d_jℤ with d_k
| ... | d_1; any object which has a Given two characters a and b, return the product character ab.
? cyc = [15,5]; a = [1,1]; b = [2,4]; ? charmul(cyc, a,b) %2 = [3, 0] ? bnf = bnfinit(x^2+23); ? bnf.cyc %4 = [3] ? charmul(bnf, [1], [2]) %5 = [0]
For Dirichlet characters on (ℤ/Nℤ)^*, additional
representations are available (Conrey labels, Conrey logarithm), see
Section se:dirichletchar or
? G = idealstar(,100); ? G.cyc %2 = [20, 2] ? a = [10, 1]; \\ usual representation for characters ? b = 7; \\ Conrey label; ? c = znconreylog(G, 11); \\ Conrey log ? charmul(G, b,b) %6 = 49 \\ Conrey label ? charmul(G, a,b) %7 = [0, 15]~ \\ Conrey log ? charmul(G, a,c) %7 = [0, 6]~ \\ Conrey log
The library syntax is
| |
charorder | |
Let cyc represent a finite abelian group by its elementary
divisors, i.e. (d_j) represents ∑_{j ≤ k} ℤ/d_jℤ with d_k
| ... | d_1; any object which has a
This function returns the order of the character
? cyc = [15,5]; chi = [1,1]; ? charorder(cyc, chi) %2 = 15 ? bnf = bnfinit(x^2+23); ? bnf.cyc %4 = [3] ? charorder(bnf, [1]) %5 = 3
For Dirichlet characters (when
? G = idealstar(,100); \\ (Z/100Z)^* ? charorder(G, 7) \\ Conrey label %2 = 4
The library syntax is
| |
chinese | |
If x and y are both intmods or both polmods, creates (with the same type) a z in the same residue class as x and in the same residue class as y, if it is possible.
? chinese(Mod(1,2), Mod(2,3)) %1 = Mod(5, 6) ? chinese(Mod(x,x^2-1), Mod(x+1,x^2+1)) %2 = Mod(-1/2*x^2 + x + 1/2, x^4 - 1) This function also allows vector and matrix arguments, in which case the operation is recursively applied to each component of the vector or matrix.
? chinese([Mod(1,2),Mod(1,3)], [Mod(1,5),Mod(2,7)]) %3 = [Mod(1, 10), Mod(16, 21)]
For polynomial arguments in the same variable, the function is applied to each
coefficient; if the polynomials have different degrees, the high degree terms
are copied verbatim in the result, as if the missing high degree terms in the
polynomial of lowest degree had been
? P = x+1; Q = x^2+2*x+1; ? chinese(P*Mod(1,2), Q*Mod(1,3)) %4 = Mod(1, 3)*x^2 + Mod(5, 6)*x + Mod(3, 6) ? chinese(Vec(P,3)*Mod(1,2), Vec(Q,3)*Mod(1,3)) %5 = [Mod(1, 6), Mod(5, 6), Mod(4, 6)] ? Pol(%) %6 = Mod(1, 6)*x^2 + Mod(5, 6)*x + Mod(4, 6)
If y is omitted, and x is a vector,
Finally
The library syntax is
| |
content | |
Computes the gcd of all the coefficients of x, when this gcd makes sense. This is the natural definition if x is a polynomial (and by extension a power series) or a vector/matrix. This is in general a weaker notion than the ideal generated by the coefficients:
? content(2*x+y) %1 = 1 \\ = gcd(2,y) over Q[y]
If x is a scalar, this simply returns the absolute value of x if x is
rational ( The content of a rational function is the ratio of the contents of the numerator and the denominator. In recursive structures, if a matrix or vector coefficient x appears, the gcd is taken not with x, but with its content:
? content([ [2], 4*matid(3) ]) %1 = 2
The content of a
The library syntax is
| |
contfrac | |
Returns the row vector whose components are the partial quotients of the continued fraction expansion of x. In other words, a result [a_0,...,a_n] means that x ~ a_0+1/(a_1+...+1/a_n). The output is normalized so that a_n != 1 (unless we also have n = 0).
The number of partial quotients n+1 is limited by
? \p19 realprecision = 19 significant digits ? contfrac(Pi) %1 = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2] ? contfrac(Pi,, 3) \\ n = 2 %2 = [3, 7, 15] x can also be a rational function or a power series.
If a vector b is supplied, the numerators are equal to the coefficients
of b, instead of all equal to 1 as above; more precisely, x ~
(1/b_0)(a_0+b_1/(a_1+...+b_n/a_n)); for a numerical continued fraction
(x real), the a_i are integers, as large as possible; if x is a
rational function, they are polynomials with deg a_i = deg b_i + 1.
The length of the result is then equal to the length of b, unless the next
partial quotient cannot be reliably computed, in which case the expansion
stops. This happens when a partial remainder is equal to zero (or too small
compared to the available significant digits for x a
A direct implementation of the numerical continued fraction
\\ "greedy" generalized continued fraction cf(x, b) = { my( a= vector(#b), t ); x *= b[1]; for (i = 1, #b, a[i] = floor(x); t = x - a[i]; if (!t || i == #b, break); x = b[i+1] / t; ); a; } There is some degree of freedom when choosing the a_i; the program above can easily be modified to derive variants of the standard algorithm. In the same vein, although no builtin function implements the related Engel expansion (a special kind of Egyptian fraction decomposition: x = 1/a_1 + 1/(a_1a_2) +... ), it can be obtained as follows:
\\ n terms of the Engel expansion of x engel(x, n = 10) = { my( u = x, a = vector(n) ); for (k = 1, n, a[k] = ceil(1/u); u = u*a[k] - 1; if (!u, break); ); a }
Obsolete hack. (don't use this): if b is an integer, nmax
is ignored and the command is understood as
The library syntax is
| |
contfracpnqn | |
When x is a vector or a one-row matrix, x is considered as the list of partial quotients [a_0,a_1,...,a_n] of a rational number, and the result is the 2 by 2 matrix [p_n,p_{n-1};q_n,q_{n-1}] in the standard notation of continued fractions, so p_n/q_n = a_0+1/(a_1+...+1/a_n). If x is a matrix with two rows [b_0,b_1,...,b_n] and [a_0,a_1,...,a_n], this is then considered as a generalized continued fraction and we have similarly p_n/q_n = (1/b_0)(a_0+b_1/(a_1+...+b_n/a_n)). Note that in this case one usually has b_0 = 1. If n ≥ 0 is present, returns all convergents from p_0/q_0 up to p_n/q_n. (All convergents if x is too small to compute the n+1 requested convergents.)
? a=contfrac(Pi,20) %1 = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2] ? contfracpnqn(a,3) %2 = [3 22 333 355] [1 7 106 113] ? contfracpnqn(a,7) %3 = [3 22 333 355 103993 104348 208341 312689] [1 7 106 113 33102 33215 66317 99532]
The library syntax is
| |
core | |
If n is an integer written as
n = df^2 with d squarefree, returns d. If flag is non-zero,
returns the two-element row vector [d,f]. By convention, we write 0 = 0
x 1^2, so
The library syntax is
| |
coredisc | |
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 non-zero integer n, this routine returns the (unique) fundamental discriminant d such that n = df^2, f a positive rational number. If flag is non-zero, 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,
Note that
The library syntax is
| |
dirdiv | |
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
| |
direuler | |
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
| |
dirmul | |
x and y being vectors of perhaps different lengths representing the Dirichlet series ∑_n x_n n^{-s} and ∑_n y_n 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
? dirmul([0,1], [0,1]); %2 = [0, 0, 0, 1]
The library syntax is
| |
divisors | |
Creates a row vector whose components are the
divisors of x. The factorization of x (as output by
By definition, these divisors are the products of the irreducible
factors of n, as produced by
The library syntax is
| |
eulerphi | |
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
The library syntax is
| |
factor | |
General factorization function, where x is a rational (including integers), a complex number with rational real and imaginary parts, or a rational function (including polynomials). 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 0^1.
ℚ and ℚ(i).
See
? fa = factor(2^2^7 + 1) %1 = [59649589127497217 1] [5704689200685129054721 1] ? isprime( fa[,1] ) %2 = [1, 1]~ \\ both entries are proven primes
Another possibility is to set the global default
A
? factor(2^2^7 +1, 10^5) %3 = [340282366920938463463374607431768211457 1]
Deprecated feature. Setting lim = 0 is the same
as setting it to
This routine uses trial division and perfect power tests, and should not be
used for huge values of lim (at most 10^9, say):
? F = (2^2^7 + 1) * 1009 * 100003; factor(F, 10^5) \\ fast, incomplete time = 0 ms. %4 = [1009 1] [34029257539194609161727850866999116450334371 1] ? factor(F, 10^9) \\ very slow time = 6,892 ms. %6 = [1009 1] [100003 1] [340282366920938463463374607431768211457 1] ? factorint(F, 1+8) \\ much faster, all small primes were found time = 12 ms. %7 = [1009 1] [100003 1] [340282366920938463463374607431768211457 1] ? factor(F) \\ complete factorisation time = 112 ms. %8 = [1009 1] [100003 1] [59649589127497217 1] [5704689200685129054721 1] Over ℚ, the prime factors are sorted in increasing order.
Rational functions.
The polynomials or rational functions to be factored must have scalar
coefficients. In particular PARI does not know how to factor
multivariate polynomials. The following domains are currently
supported: ℚ, ℝ, ℂ, ℚ_p, finite fields and number fields.
See
The routine guesses a sensible ring over which to factor: the
smallest ring containing all coefficients, taking into account quotient
structures induced by
? T = x^2+1; ? factor(T); \\ over Q ? factor(T*Mod(1,3)) \\ over F_3 ? factor(T*ffgen(ffinit(3,2,'t))^0) \\ over F_{3^2} ? factor(T*Mod(Mod(1,3), t^2+t+2)) \\ over F_{3^2}, again ? factor(T*(1 + O(3^6)) \\ over Q_3, 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(2^{1/3})
In most cases, it is clearer and simpler to call an
explicit variant than to rely on the generic
? factormod(T, 3) \\ over F_3 ? factorff(T, 3, t^2+t+2)) \\ over F_{3^2} ? factorpadic(T, 3,6) \\ over Q_3, precision 6 ? nffactor(y^3-2, T) \\ over Q(2^{1/3}) ? polroots(T) \\ over C Note that factorization of polynomials is done up to multiplication by a constant. In particular, the factors of rational polynomials will have integer coefficients, and the content of a polynomial or rational function is discarded and not included in the factorization. If needed, you can always ask for the content explicitly:
? factor(t^2 + 5/2*t + 1) %1 = [2*t + 1 1] [t + 2 1] ? content(t^2 + 5/2*t + 1) %2 = 1/2
The irreducible factors are sorted by increasing degree.
See also
The library syntax is
| |
factorback | |
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(12) %1 = [2 2] [3 1] ? factorback(%) %2 = 12 ? factorback([2,3], [2,1]) \\ 2^3 * 3^1 %3 = 12 ? factorback([5,2,3]) %4 = 30
The library syntax is
| |
factorcantor | |
Factors the polynomial x modulo the
prime p, using distinct degree plus
Cantor-Zassenhaus. The coefficients of x must be
operation-compatible with ℤ/pℤ. The result is a two-column matrix, the
first column being the irreducible polynomials dividing x, and the second
the exponents. If you want only the degrees of the irreducible
polynomials (for example for computing an L-function), use
The library syntax is
| |
factorff | |
Factors the polynomial x in the field 𝔽_q defined by the irreducible polynomial a over 𝔽_p. The coefficients of x must be operation-compatible with ℤ/pℤ. The result is a two-column matrix: the first column contains the irreducible factors of x, and the second their exponents. If all the coefficients of x are in 𝔽_p, a much faster algorithm is applied, using the computation of isomorphisms between finite fields.
Either a or p can omitted (in which case both are ignored) if x has
? factorff(x^2 + 1, 5, y^2+3) \\ over F_5[y]/(y^2+3) ~ F_25 %1 = [Mod(Mod(1, 5), Mod(1, 5)*y^2 + Mod(3, 5))*x + Mod(Mod(2, 5), Mod(1, 5)*y^2 + Mod(3, 5)) 1] [Mod(Mod(1, 5), Mod(1, 5)*y^2 + Mod(3, 5))*x + Mod(Mod(3, 5), Mod(1, 5)*y^2 + Mod(3, 5)) 1] ? t = ffgen(y^2 + Mod(3,5), 't); \\ a generator for F_25 as a t_FFELT ? factorff(x^2 + 1) \\ not enough information to determine the base field *** at top-level: factorff(x^2+1) *** ^--------------- *** factorff: incorrect type in factorff. ? factorff(x^2 + t^0) \\ make sure a coeff. is a t_FFELT %3 = [x + 2 1] [x + 3 1] ? factorff(x^2 + t + 1) %11 = [x + (2*t + 1) 1] [x + (3*t + 4) 1] Notice that the second syntax is easier to use and much more readable.
The library syntax is
| |
factorial | |
Factorial of x. The expression x! gives a result which is an integer,
while
The library syntax is
| |
factorint | |
Factors the integer n into a product of
pseudoprimes (see
By convention 0 is factored as 0^1, and 1 as the empty factorization;
also the divisors are by default not proven primes is they are larger than
2^{64}, they only failed the BPSW compositeness test (see
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
The library syntax is
| |
factormod | |
Factors the polynomial x modulo the prime integer p, using
Berlekamp. The coefficients of x must be operation-compatible with
ℤ/pℤ. The result is a two-column matrix, the first column being the
irreducible polynomials dividing x, and the second the exponents. If flag
is non-zero, outputs only the degrees of the irreducible polynomials
(for example, for computing an L-function). A different algorithm for
computing the mod p factorization is
The library syntax is
| |
ffgen | |
Return a
? g = ffgen(8, 't); ? g.mod %2 = t^3 + t^2 + 1 ? g.p %3 = 2 ? g.f %4 = 3 ? ffgen(6) *** at top-level: ffgen(6) *** ^-------- *** ffgen: not a prime number in ffgen: 6. Alternative syntax: instead of a prime power q = p^f, one may input the pair [p,f]:
? g = ffgen([2,4], 't); ? g.p %2 = 2 ? g.mod %3 = t^4 + t^3 + t^2 + t + 1
Finally, one may input
directly the polynomial P (monic, irreducible, with
The library syntax is
To create a generator for a prime finite field, the function
| |
ffinit | |
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
The library syntax is
| |
fflog | |
Discrete logarithm of the finite field element x in base g, i.e.
an e in ℤ such that g^e = o. If
present, o represents the multiplicative order of g, see
Section se:DLfun; the preferred format for
this parameter is If no o is given, assume that g is a primitive root. The result is undefined if e does not exist. This function uses
* a combination of generic discrete log algorithms (see * 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
| |
ffnbirred | |
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
| |
fforder | |
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
? 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
| |
ffprimroot | |
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
? 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
| |
fibonacci | |
x-th Fibonacci number.
The library syntax is
| |
gcd | |
Creates the greatest common divisor of x and y.
If you also need the u and v such that x*u + y*v = gcd(x,y),
use the
When x and y are both given and one of them is a vector/matrix type,
the GCD is again taken recursively on each component, but in a different way.
If y is a vector, resp. matrix, then the result has the same type as y,
and components equal to The algorithm used is a naive Euclid except for the following inputs: * integers: use modified right-shift binary ("plus-minus" variant). * univariate polynomials with coefficients in the same number field (in particular rational): use modular gcd algorithm. * general polynomials: use the subresultant algorithm if coefficient explosion is likely (non modular coefficients). If u and v are polynomials in the same variable with inexact coefficients, their gcd is defined to be scalar, so that
? a = x + 0.0; gcd(a,a) %1 = 1 ? b = y*x + O(y); gcd(b,b) %2 = y ? c = 4*x + O(2^3); gcd(c,c) %3 = 4
A good quantitative check to decide whether such a
gcd "should be" non-trivial, is to use
The library syntax is
| |
gcdext | |
Returns [u,v,d] such that d is the gcd of x,y, x*u+y*v = gcd(x,y), and u and v minimal in a natural sense. The arguments must be integers or polynomials.
? [u, v, d] = gcdext(32,102) %1 = [16, -5, 2] ? d %2 = 2 ? gcdext(x^2-x, x^2+x-2) %3 = [-1/2, 1/2, x - 1]
If x,y are polynomials in the same variable and inexact
coefficients, then compute u,v,d such that x*u+y*v = d, where d
approximately divides both and x and y; in particular, we do not obtain
? a = x + 0.0; gcd(a,a) %1 = 1 ? gcdext(a,a) %2 = [0, 1, x + 0.E-28] ? gcdext(x-Pi, 6*x^2-zeta(2)) %3 = [-6*x - 18.8495559, 1, 57.5726923] For inexact inputs, the output is thus not well defined mathematically, but you obtain explicit polynomials to check whether the approximation is close enough for your needs.
The library syntax is
| |
hilbert | |
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
The library syntax is
| |
isfundamental | |
True (1) if x is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise.
The library syntax is
| |
ispolygonal | |
True (1) if the integer x is an s-gonal number, false (0) if not.
The parameter s > 2 must be a
? ispolygonal(36, 3, &N) %1 = 1 ? N
The library syntax is
| |
ispower | |
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
If k is omitted, only integers and fractions are allowed for x and the
function returns the maximal k ≥ 2 such that x = n^k is a perfect
power, or 0 if no such k exist; in particular 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
k = (x.p ^ x.f - 1) / fforder(x)
The library syntax is
| |
ispowerful | |
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
| |
isprime | |
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 prime and has more than 1000 digits, say.
Use
If flag = 0, use a combination of Baillie-PSW pseudo primality test (see
If flag = 1, use Selfridge-Pocklington-Lehmer "p-1" test and output a primality certificate as follows: return * 0 if x is composite, * 1 if x is small enough that passing Baillie-PSW test guarantees its primality (currently x < 2^{64}, as checked by Jan Feitsma), * 2 if x is a large prime whose primality could only sensibly be proven (given the algorithms implemented in PARI) using the APRCL test. * Otherwise (x is large and x-1 is smooth) output a three column matrix as a primality certificate. The first column contains prime divisors p of x-1 (such that ∏ p^{v_p(x-1)} > x^{1/3}), the second the corresponding elements a_p as in Proposition 8.3.1 in GTM 138 , and the third the output of isprime(p,1).
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
If flag = 2, use APRCL.
The library syntax is
| |
isprimepower | |
If x = p^k 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
| |
ispseudoprime | |
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 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), P smallest positive integer such that P^2 - 4 is not a square mod x).
There are no known composite numbers passing the above test, although it is
expected that infinitely many such numbers exist. In particular, all
composites ≤ 2^{64} are correctly detected (checked using
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
| |
ispseudoprimepower | |
If x = p^k is a pseudo-prime power (p pseudo-prime as per
More precisely, k is always the largest integer such that x = n^k for
some integer n and, when n ≤ 2^{64} the function returns k > 0 if and
only if n is indeed prime. When n > 2^{64} is larger than the threshold,
the function may return 1 even though n is composite: it only passed
an
The library syntax is
| |
issquare | |
True (1) if x is a square, false (0)
if not. What "being a square" means depends on the type of x: all
? 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
? 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
| |
issquarefree | |
True (1) if x is squarefree, false (0) if not. Here x can be an integer or a polynomial.
The library syntax is
| |
istotient | |
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
| |
kronecker | |
Kronecker symbol (x|y), where x and y must be of type integer. By definition, this is the extension of Legendre symbol to ℤ x ℤ by total multiplicativity in both arguments with the following special rules for y = 0, -1 or 2: * (x|0) = 1 if |x |= 1 and 0 otherwise. * (x|-1) = 1 if x ≥ 0 and -1 otherwise. * (x|2) = 0 if x is even and 1 if x = 1,-1 mod 8 and -1 if x = 3,-3 mod 8.
The library syntax is
| |
lcm | |
Least common multiple of x and y, i.e. such that lcm(x,y)*gcd(x,y) = x*y, up to units. If y is omitted and x is a vector, returns the lcm of all components of x. For integer arguments, return the non-negative lcm.
When x and y are both given and one of them is a vector/matrix type,
the LCM is again taken recursively on each component, but in a different way.
If y is a vector, resp. matrix, then the result has the same type as y,
and components equal to
Note that
l = v[1]; for (i = 1, #v, l = lcm(l, v[i]))
Indeed,
? v = vector(10^5, i, random); ? lcm(v); time = 546 ms. ? l = v[1]; for (i = 1, #v, l = lcm(l, v[i])) time = 4,561 ms.
The library syntax is
| |
logint | |
Return the largest integer e so that b^e ≤ x, where the parameters b > 1 and x > 0 are both integers. If the parameter z is present, set it to b^e.
? logint(1000, 2) %1 = 9 ? 2^9 %2 = 512 ? logint(1000, 2, &z) %3 = 9 ? z %4 = 512
The number of digits used to write b in base x is
? #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
| |
moebius | |
Moebius μ-function of |x|. x must be of type integer.
The library syntax is
| |
nextprime | |
Finds the smallest pseudoprime (see
The library syntax is
| |
numbpart | |
Gives the number of unrestricted partitions of
n, usually called p(n) in the literature; in other words the number of
nonnegative integer solutions to a+2b+3c+.. .= n. n must be of type
integer and n < 10^{15} (with trivial values p(n) = 0 for n < 0 and
p(0) = 1). The algorithm uses the Hardy-Ramanujan-Rademacher formula.
To explicitly enumerate them, see
The library syntax is
| |
numdiv | |
Number of divisors of |x|. x must be of type integer.
The library syntax is
| |
omega | |
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
| |
partitions | |
Returns the vector of partitions of the integer k as a sum of positive
integers (parts); for k < 0, it returns the empty set
? partitions(3) %1 = [Vecsmall([3]), Vecsmall([1, 2]), Vecsmall([1, 1, 1])]
correspond to 3, 1+2 and 1+1+1. The number
of (unrestricted) partitions of k is given
by
? #partitions(50) %1 = 204226 ? numbpart(50) %2 = 204226 Optional parameters n and a are as follows: * n = nmax (resp. n = [nmin,nmax]) restricts partitions to length less than nmax (resp. length between nmin and nmax), where the length is the number of nonzero entries. * a = amax (resp. a = [amin,amax]) restricts the parts to integers less than amax (resp. between amin and amax).
? partitions(4, 2) \\ parts bounded by 2 %1 = [Vecsmall([2, 2]), Vecsmall([1, 1, 2]), Vecsmall([1, 1, 1, 1])] ? partitions(4,, 2) \\ at most 2 parts %2 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])] ? partitions(4,[0,3], 2) \\ at most 2 parts %3 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])] By default, parts are positive and we remove zero entries unless amin ≤ 0, in which case nmin is ignored and X is of constant length nmax:
? partitions(4, [0,3]) \\ parts between 0 and 3 %1 = [Vecsmall([0, 0, 1, 3]), Vecsmall([0, 0, 2, 2]),\ Vecsmall([0, 1, 1, 2]), Vecsmall([1, 1, 1, 1])]
The library syntax is
| |
polrootsff | |
Returns the vector of distinct roots of the polynomial x in the field
𝔽_q defined by the irreducible polynomial a over 𝔽_p. The
coefficients of x must be operation-compatible with ℤ/pℤ.
Either a or p can omitted (in which case both are ignored) if x has
? polrootsff(x^2 + 1, 5, y^2+3) \\ over F_5[y]/(y^2+3) ~ F_25 %1 = [Mod(Mod(3, 5), Mod(1, 5)*y^2 + Mod(3, 5)), Mod(Mod(2, 5), Mod(1, 5)*y^2 + Mod(3, 5))] ? t = ffgen(y^2 + Mod(3,5), 't); \\ a generator for F_25 as a t_FFELT ? polrootsff(x^2 + 1) \\ not enough information to determine the base field *** at top-level: polrootsff(x^2+1) *** ^----------------- *** polrootsff: incorrect type in factorff. ? polrootsff(x^2 + t^0) \\ make sure one coeff. is a t_FFELT %3 = [3, 2] ? polrootsff(x^2 + t + 1) %4 = [2*t + 1, 3*t + 4] Notice that the second syntax is easier to use and much more readable.
The library syntax is
| |
precprime | |
Finds the largest pseudoprime (see
The library syntax is
| |
prime | |
The n-th prime number
? prime(10^9) %1 = 22801763489
Uses checkpointing and a naive O(n) algorithm. Will need
about 30 minutes for n up to 10^{11}; make sure to start gp with
The library syntax is
| |
primepi | |
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
The library syntax is
| |
primes | |
Creates a row vector whose components are the first n prime numbers.
(Returns the empty vector for n ≤ 0.) A
? 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
| |
qfbclassno | |
Ordinary class number of the quadratic order of discriminant D, for "small" values of D.
* if D > 0 or flag = 1, use a O(|D|^{1/2})
algorithm (compute L(1,χ_D) with the approximate functional equation).
This is slower than
* if D < 0 and flag = 0 (or omitted), use a O(|D|^{1/4})
algorithm (Shanks's baby-step/giant-step method). It should
be faster than
Important warning. In the latter case, this function only
implements part of Shanks's method (which allows to speed it up
considerably). It gives unconditionnally correct results for |D| < 2.
10^{10}, but may give incorrect results for larger values if the class
group has many cyclic factors. We thus recommend to double-check results
using the function Warning. Contrary to what its name implies, this routine does not compute the number of classes of binary primitive forms of discriminant D, which is equal to the narrow class number. The two notions are the same when D < 0 or the fundamental unit ϵ has negative norm; when D > 0 and Nϵ > 0, the number of classes of forms is twice the ordinary class number. This is a problem which we cannot fix for backward compatibility reasons. Use the following routine if you are only interested in the number of classes of forms:
QFBclassno(D) = qfbclassno(D) * if (D < 0 || norm(quadunit(D)) < 0, 1, 2) Here are a few examples:
? qfbclassno(400000028) time = 3,140 ms. %1 = 1 ? quadclassunit(400000028).no time = 20 ms. \\ { much faster} %2 = 1 ? qfbclassno(-400000028) time = 0 ms. %3 = 7253 \\ { correct, and fast enough} ? quadclassunit(-400000028).no time = 0 ms. %4 = 7253
See also
The library syntax is
Finally
| |
qfbcompraw | |
composition of the binary quadratic forms x and y, without reduction of the result. This is useful e.g. to compute a generating element of an ideal. The result is undefined if x and y do not have the same discriminant.
The library syntax is
| |
qfbhclassno | |
Hurwitz class number of x, where
x is non-negative and congruent to 0 or 3 modulo 4. For x > 5.
10^5, we assume the GRH, and use
The library syntax is
| |
qfbnucomp | |
composition of the primitive positive
definite binary quadratic forms x and y (type The current implementation is slower than the generic routine for small D, and becomes faster when D has about 45 bits.
The library syntax is
| |
qfbnupow | |
n-th power of the primitive positive definite
binary quadratic form x using Shanks's NUCOMP and NUDUPL algorithms;
if set, L should be equal to The current implementation is slower than the generic routine for small discriminant D, and becomes faster for D ~ 2^{45}.
The library syntax is
| |
qfbpowraw | |
n-th power of the binary quadratic form
x, computed without doing any reduction (i.e. using
The library syntax is
| |
qfbprimeform | |
Prime binary quadratic form of discriminant
x whose first coefficient is p, where |p| is a prime number.
By abuse of notation,
p = ± 1 is also valid and returns the unit form. Returns an
error if x is not a quadratic residue mod p, or if x < 0 and p < 0.
(Negative definite
The library syntax is
| |
qfbred | |
Reduces the binary quadratic form x (updating Shanks's distance function if x is indefinite). The binary digits of flag are toggles meaning 1: perform a single reduction step 2: don't update Shanks's distance The arguments d, isd, sd, if present, supply the values of the discriminant, floor{sqrt{d}}, and sqrt{d} respectively (no checking is done of these facts). If d < 0 these values are useless, and all references to Shanks's distance are irrelevant.
The library syntax is
and for indefinite forms:
| |
qfbredsl2 | |
Reduction of the (real or imaginary) binary quadratic form x, return
[y,g] where y is reduced and g in SL(2,ℤ) is such that
g.x = y; data, if
present, must be equal to [D,
The library syntax is
| |
qfbsolve | |
Solve the equation Q(x,y) = p over the integers, where Q is a binary quadratic form and p a prime number. Return [x,y] as a two-components vector, or zero if there is no solution. Note that this function returns only one solution and not all the solutions.
Let D = disc Q. The algorithm used runs in probabilistic polynomial time
in p (through the computation of a square root of D modulo p); it is
polynomial time in D if Q is imaginary, but exponential time if Q is
real (through the computation of a full cycle of reduced forms). In the
latter case, note that
The library syntax is
| |
quadclassunit | |
Buchmann-McCurley's sub-exponential algorithm for computing the class group of a quadratic order of discriminant D.
This function should be used instead of The result is a vector v whose components should be accessed using member functions:
*
*
*
*
The flag is obsolete and should be left alone. In older versions,
it supposedly computed the narrow class group when D > 0, but this did not
work at all; use the general function Optional parameter tech is a row vector of the form [c_1, c_2], where c_1 ≤ c_2 are non-negative real numbers which control the execution time and the stack size, see se:GRHbnf. The parameter is used as a threshold to balance the relation finding phase against the final linear algebra. Increasing the default c_1 means that relations are easier to find, but more relations are needed and the linear algebra will be harder. The default value for c_1 is 0 and means that it is taken equal to c_2. The parameter c_2 is mostly obsolete and should not be changed, but we still document it for completeness: we compute a tentative class group by generators and relations using a factorbase of prime ideals ≤ c_1 (log |D|)^2, then prove that ideals of norm ≤ c_2 (log |D|)^2 do not generate a larger group. By default an optimal c_2 is chosen, so that the result is provably correct under the GRH --- a famous result of Bach states that c_2 = 6 is fine, but it is possible to improve on this algorithmically. You may provide a smaller c_2, it will be ignored (we use the provably correct one); you may provide a larger c_2 than the default value, which results in longer computing times for equally correct outputs (under GRH).
The library syntax is
| |
quaddisc | |
Discriminant of the étale algebra ℚ(sqrt{x}), where x ∈ ℚ^*.
This is the same as
? quaddisc(7) %1 = 28 ? quaddisc(-7) %2 = -7
The library syntax is
| |
quadgen | |
Creates the quadratic number ω = (a+sqrt{D})/2 where a = 0 if D = 0 mod 4, a = 1 if D = 1 mod 4, so that (1,ω) is an integral basis for the quadratic order of discriminant D. D must be an integer congruent to 0 or 1 modulo 4, which is not a square.
The library syntax is
| |
quadhilbert | |
Relative equation defining the Hilbert class field of the quadratic field of discriminant D. If D < 0, uses complex multiplication (Schertz's variant).
If D > 0 Stark units are used and (in rare cases) a
vector of extensions may be returned whose compositum is the requested class
field. See
The library syntax is
| |
quadpoly | |
Creates the "canonical" quadratic
polynomial (in the variable v) corresponding to the discriminant D,
i.e. the minimal polynomial of
The library syntax is
| |
quadray | |
Relative equation for the ray
class field of conductor f for the quadratic field of discriminant D
using analytic methods. A For D < 0, uses the σ function and Schertz's method.
For D > 0, uses Stark's conjecture, and a vector of relative equations may be
returned. See
The library syntax is
| |
quadregulator | |
Regulator of the quadratic field of positive discriminant x. Returns
an error if x is not a discriminant (fundamental or not) or if x is a
square. See also
The library syntax is
| |
quadunit | |
Fundamental unit of the real quadratic field ℚ(sqrt D) where D is the positive discriminant of the field. If D is not a fundamental discriminant, this probably gives the fundamental unit of the corresponding order. D must be an integer congruent to 0 or 1 modulo 4, which is not a square; the result is a quadratic number (see Section se:quadgen).
The library syntax is
| |
ramanujantau | |
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
? 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
| |
randomprime | |
Returns a strong pseudo prime (see
The library syntax is
| |
removeprimes | |
Removes the primes listed in x from
the prime number table. In particular
The library syntax is
| |
sigma | |
Sum of the k-th powers of the positive divisors of |x|. x and k must be of type integer.
The library syntax is
| |
sqrtint | |
Returns the integer square root of x, i.e. the largest integer y such that y^2 ≤ x, where x a non-negative integer.
? N = 120938191237; sqrtint(N) %1 = 347761 ? sqrt(N) %2 = 347761.68741970412747602130964414095216
The library syntax is
| |
sqrtnint | |
Returns the integer n-th root of x, i.e. the largest integer y such that y^n ≤ x, where x is a non-negative integer.
? N = 120938191237; sqrtnint(N, 5) %1 = 164 ? N^(1/5) %2 = 164.63140849829660842958614676939677391
The special case n = 2 is
The library syntax is
| |
stirling | |
Stirling number of the first kind s(n,k) (flag = 1, default) or of the second kind S(n,k) (flag = 2), where n, k are non-negative integers. The former is (-1)^{n-k} times the number of permutations of n symbols with exactly k cycles; the latter is the number of ways of partitioning a set of n elements into k non-empty subsets. Note that if all s(n,k) are needed, it is much faster to compute ∑_k s(n,k) x^k = x(x-1)...(x-n+1). Similarly, if a large number of S(n,k) are needed for the same k, one should use ∑_n S(n,k) x^n = (x^k)/((1-x)...(1-kx)). (Should be implemented using a divide and conquer product.) Here are simple variants for n fixed:
/* list of s(n,k), k = 1..n */ vecstirling(n) = Vec( factorback(vector(n-1,i,1-i*'x)) ) /* list of S(n,k), k = 1..n */ vecstirling2(n) = { my(Q = x^(n-1), t); vector(n, i, t = divrem(Q, x-i); Q=t[1]; simplify(t[2])); }
The library syntax is
| |
sumdedekind | |
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
| |
sumdigits | |
Sum of digits in the integer n, when written in base B > 1.
? sumdigits(123456789) %1 = 45 ? sumdigits(123456789, 2) %1 = 16
Note that the sum of bits in n is also returned by
The library syntax is
| |
zncharinduce | |
Let G be attached to (ℤ/qℤ)^* (as per
* a
* a
Let N be a multiple of q, return the character modulo N induced by
? G = idealstar(,4); ? chi = znconreylog(G,1); \\ trivial character mod 4 ? zncharinduce(G, chi, 80) \\ now mod 80 %3 = [0, 0, 0]~ ? zncharinduce(G, 1, 80) \\ same using directly Conrey label %4 = [0, 0, 0]~ ? G2 = idealstar(,80); ? zncharinduce(G, 1, G2) \\ same %4 = [0, 0, 0]~ ? chi = zncharinduce(G, 3, G2) \\ induce the non-trivial character mod 4 %5 = [1, 0, 0]~ ? znconreyconductor(G2, chi, &chi0) %6 = [4, Mat([2, 2])] ? chi0 %7 = [1]~ Here is a larger example:
? G = idealstar(,126000); ? label = 1009; ? chi = znconreylog(G, label) %3 = [0, 0, 0, 14, 0]~ ? N0 = znconreyconductor(G, label, &chi0) %4 = [125, Mat([5, 3])] ? chi0 \\ primitive character mod 5^3 attached to chi %5 = [14]~ ? G0 = idealstar(,N0); ? zncharinduce(G0, chi0, G) \\ induce back %7 = [0, 0, 0, 14, 0]~ ? znconreyexp(G, %) %8 = 1009
The library syntax is
| |
zncharisodd | |
Let G be attached to (ℤ/Nℤ)^* (as per
* a
* a
Return 1 if and only if
? G = idealstar(,8); ? zncharisodd(G, 1) \\ trivial character %2 = 0 ? zncharisodd(G, 3) %3 = 1 ? chareval(G, 3, -1) %4 = 1/2
The library syntax is
| |
znchartokronecker | |
Let G be attached to (ℤ/Nℤ)^* (as per
* a
* a
If flag = 0, return the discriminant D if If flag = 1, return the fundamental discriminant attached to the corresponding primitive character.
? G = idealstar(,8); CHARS = [1,3,5,7]; \\ Conrey labels ? apply(t->znchartokronecker(G,t), CHARS) %2 = [4, -8, 8, -4] ? apply(t->znchartokronecker(G,t,1), CHARS) %3 = [1, -8, 8, -4]
The library syntax is
| |
znconreychar | |
Given a bid attached to (ℤ/qℤ)^* (as per
Let q = ∏_p p^{e_p} be the factorization of q into distinct primes. For all odd p with e_p > 0, let g_p be the element in (ℤ/qℤ)^* which is * congruent to 1 mod q/p^{e_p}, * congruent mod p^{e_p} to the smallest integer whose order is φ(p^{e_p}). For p = 2, we let g_4 (if 2^{e_2} ≥ 4) and g_8 (if furthermore (2^{e_2} ≥ 8) be the elements in (ℤ/qℤ)^* which are * congruent to 1 mod q/2^{e_2}, * g_4 = -1 mod 2^{e_2}, * g_8 = 5 mod 2^{e_2}.
Then the g_p (and the extra g_4 and g_8 if 2^{e_2} ≥ 2) are
independent
generators of (ℤ/qℤ)^*, i.e. every m in (ℤ/qℤ)^* can be written
uniquely as ∏_p g_p^{m_p}, where m_p is defined modulo the order
o_p of g_p
and p ∈ S_q, the set of prime divisors of q together with 4
if 4 | q and 8 if 8 | q. Note that the g_p are in general
not SNF
generators as produced by
The Conrey logarithm of m is the vector (m_p)_{p ∈ S_q}, obtained
via Note. It is useless to include the generators in the bid, except for debugging purposes: they are well defined from elementary matrix operations and Chinese remaindering, their explicit value as elements in (ℤ/qℤ)^* is never used.
? G = idealstar(,8,2); /*add generators for debugging:*/ ? G.cyc %2 = [2, 2] \\ Z/2 x Z/2 ? G.gen %3 = [7, 3] ? znconreychar(G,1) \\ 1 is always the trivial character %4 = [0, 0] ? znconreychar(G,2) \\ 2 is not coprime to 8 !!! *** at top-level: znconreychar(G,2) *** ^----------------- *** znconreychar: elements not coprime in Zideallog: 2 8 *** Break loop: type 'break' to go back to GP prompt break> ? znconreychar(G,3) %5 = [0, 1] ? znconreychar(G,5) %6 = [1, 1] ? znconreychar(G,7) %7 = [1, 0] We indeed get all 4 characters of (ℤ/8ℤ)^*. For convenience, we allow to input the Conrey logarithm of m instead of m:
? G = idealstar(,55); ? znconreychar(G,7) %2 = [7, 0] ? znconreychar(G, znconreylog(G,7)) %3 = [7, 0]
The library syntax is
| |
znconreyconductor | |
Let bid be attached to (ℤ/qℤ)^* (as per
* a
* a
Return the conductor of
If
? G = idealstar(,126000); ? znconreyconductor(G,11) \\ primitive %2 = 126000 ? znconreyconductor(G,1) \\ trivial character, not primitive! %3 = [1, matrix(0,2)] ? N0 = znconreyconductor(G,1009, &chi0) \\ character mod 5^3 %4 = [125, Mat([5, 3])] ? chi0 %5 = [14]~ ? G0 = idealstar(,N0); \\ format [N,factor(N)] accepted ? znconreyexp(G0, chi0) %7 = 9 ? znconreyconductor(G0, chi0) \\ now primitive, as expected %8 = 125
The group
The library syntax is
| |
znconreyexp | |
Given a bid attached to (ℤ/qℤ)^* (as per
The character chi is given either as a
*
*
? G = idealstar(,126000) ? znconreylog(G,1) %2 = [0, 0, 0, 0, 0]~ ? znconreyexp(G,%) %3 = 1 ? G.cyc \\ SNF generators %4 = [300, 12, 2, 2, 2] ? chi = [100, 1, 0, 1, 0]; \\ some random character on SNF generators ? znconreylog(G, chi) \\ in terms of Conrey generators %6 = [0, 3, 3, 0, 2]~ ? znconreyexp(G, %) \\ apply to a Conrey log %7 = 18251 ? znconreyexp(G, chi) \\ ... or a char on SNF generators %8 = 18251 ? znconreychar(G,%) %9 = [100, 1, 0, 1, 0]
The library syntax is
| |
znconreylog | |
Given a bid attached to (ℤ/qℤ)^* (as per
Let q = ∏_p p^{e_p} be the factorization of q into distinct primes, where we assume e_2 = 0 or e_2 ≥ 2. (If e_2 = 1, we can ignore 2 from the factorization, as if we replaced q by q/2, since (ℤ/qℤ)^* ~ (ℤ/(q/2)ℤ)^*.) For all odd p with e_p > 0, let g_p be the element in (ℤ/qℤ)^* which is * congruent to 1 mod q/p^{e_p}, * congruent mod p^{e_p} to the smallest integer whose order is φ(p^{e_p}) for p odd, For p = 2, we let g_4 (if 2^{e_2} ≥ 4) and g_8 (if furthermore (2^{e_2} ≥ 8) be the elements in (ℤ/qℤ)^* which are * congruent to 1 mod q/2^{e_2}, * g_4 = -1 mod 2^{e_2}, * g_8 = 5 mod 2^{e_2}.
Then the g_p (and the extra g_4 and g_8 if 2^{e_2} ≥ 2) are
independent
generators of ℤ/qℤ^*, i.e. every m in (ℤ/qℤ)^* can be written
uniquely as ∏_p g_p^{m_p}, where m_p is defined modulo the
order o_p of g_p
and p ∈ S_q, the set of prime divisors of q together with 4
if 4 | q and 8 if 8 | q.
Note that the g_p are in general not SNF
generators as produced by
The Conrey logarithm of m is the vector (m_p)_{p ∈ S_q}. The inverse
function
? G = idealstar(,126000); ? znconreylog(G,1) %2 = [0, 0, 0, 0, 0]~ ? znconreyexp(G, %) %3 = 1 ? znconreylog(G,2) \\ 2 is not coprime to modulus !!! *** at top-level: znconreylog(G,2) *** ^----------------- *** znconreylog: elements not coprime in Zideallog: 2 126000 *** Break loop: type 'break' to go back to GP prompt break> ? znconreylog(G,11) \\ wrt. Conrey generators %4 = [0, 3, 1, 76, 4]~ ? log11 = ideallog(,11,G) \\ wrt. SNF generators %5 = [178, 3, -75, 1, 0]~
For convenience, we allow to input the ordinary discrete log of m,
? znconreylog(G, log11) %7 = [0, 3, 1, 76, 4]~
We also allow a character (
? G.cyc %8 = [300, 12, 2, 2, 2] ? chi = [10,1,0,1,1]; ? znconreylog(G, chi) %10 = [1, 3, 3, 10, 2]~ ? n = znconreyexp(G, chi) %11 = 84149 ? znconreychar(G, n) %12 = [10, 1, 0, 1, 1]
The library syntax is
| |
zncoppersmith | |
N being an integer and P ∈ ℤ[X], finds all integers x with
|x| ≤ X such that
gcd(N, P(x)) ≥ B,
using Coppersmith's algorithm (a famous application of the LLL
algorithm). X must be smaller than exp(log^2 B / (deg(P) log N)):
for B = N, this means X < N^{1/deg(P)}. Some x larger than X may
be returned if you are very lucky. The smaller B (or the larger X), the
slower the routine will be. The strength of Coppersmith method is the
ability to find roots modulo a general composite N: if N is a prime
or a prime power, We shall now present two simple applications. The first one is finding non-trivial factors of N, given some partial information on the factors; in that case B must obviously be smaller than the largest non-trivial divisor of N.
setrand(1); \\ to make the example reproducible interval = [10^30, 10^31]; p = randomprime(interval); q = randomprime(interval); N = p*q; p0 = p % 10^20; \\ assume we know 1) p > 10^29, 2) the last 19 digits of p L = zncoppersmith(10^19*x + p0, N, 10^12, 10^29) \\ result in 10ms. %6 = [738281386540] ? gcd(L[1] * 10^19 + p0, N) == p %7 = 1 and we recovered p, faster than by trying all possibilities < 10^{12}. The second application is an attack on RSA with low exponent, when the message x is short and the padding P is known to the attacker. We use the same RSA modulus N as in the first example:
setrand(1); P = random(N); \\ known padding e = 3; \\ small public encryption exponent X = floor(N^0.3); \\ N^(1/e - epsilon) x0 = random(X); \\ unknown short message C = lift( (Mod(x0,N) + P)^e ); \\ known ciphertext, with padding P zncoppersmith((P + x)^3 - C, N, X) \\ result in 244ms. %14 = [2679982004001230401] ? %[1] == x0 %15 = 1 We guessed an integer of the order of 10^{18}, almost instantly.
The library syntax is
| |
znlog | |
This functions allows two distinct modes of operation depending on g:
* if g is the output of * else g is an explicit element in (ℤ/Nℤ)^*, we compute the discrete logarithm of x in (ℤ/Nℤ)^* in base g. The rest of this entry describes the latter possibility. The result is [] when x is not a power of g, though the function may also enter an infinite loop in this case.
If present, o represents the multiplicative order of g, see
Section se:DLfun; the preferred format for this parameter is
? p = nextprime(10^4); g = znprimroot(p); o = [p-1, factor(p-1)]; ? for(i=1,10^4, znlog(i, g, o)) time = 205 ms. ? for(i=1,10^4, znlog(i, g)) time = 244 ms. \\ a little slower The result is undefined if g is not invertible mod N or if the supplied order is incorrect. This function uses * a combination of generic discrete log algorithms (see below). * in (ℤ/Nℤ)^* when N is prime: a linear sieve index calculus method, suitable for N < 10^{50}, say, is used for large prime divisors of the order. The generic discrete log algorithms are: * Pohlig-Hellman algorithm, to reduce to groups of prime order q, where q | p-1 and p is an odd prime divisor of N, * Shanks baby-step/giant-step (q < 2^{32} is small), * Pollard rho method (q > 2^{32}). The latter two algorithms require O(sqrt{q}) operations in the group on average, hence will not be able to treat cases where q > 10^{30}, say. In addition, Pollard rho is not able to handle the case where there are no solutions: it will enter an infinite loop.
? g = znprimroot(101) %1 = Mod(2,101) ? znlog(5, g) %2 = 24 ? g^24 %3 = Mod(5, 101) ? G = znprimroot(2 * 101^10) %4 = Mod(110462212541120451003, 220924425082240902002) ? znlog(5, G) %5 = 76210072736547066624 ? G^% == 5 %6 = 1 ? N = 2^4*3^2*5^3*7^4*11; g = Mod(13, N); znlog(g^110, g) %7 = 110 ? znlog(6, Mod(2,3)) \\ no solution %8 = [] For convenience, g is also allowed to be a p-adic number:
? g = 3+O(5^10); znlog(2, g) %1 = 1015243 ? g^% %2 = 2 + O(5^10)
The library syntax is
| |
znorder | |
x must be an integer mod n, and the
result is the order of x in the multiplicative group (ℤ/nℤ)^*. Returns
an error if x is not invertible.
The parameter o, if present, represents a non-zero
multiple of the order of x, see Section se:DLfun; the preferred format for
this parameter is
The library syntax is
| |
znprimroot | |
Returns a primitive root (generator) of (ℤ/nℤ)^*, whenever this latter group is cyclic (n = 4 or n = 2p^k or n = p^k, where p is an odd prime and k ≥ 0). If the group is not cyclic, the result is undefined. If n is a prime power, then the smallest positive primitive root is returned. This may not be true for n = 2p^k, p odd. Note that this function requires factoring p-1 for p as above, in order to determine the exact order of elements in (ℤ/nℤ)^*: this is likely to be costly if p is large.
The library syntax is
| |
znstar | |
Gives the structure of the multiplicative group (ℤ/nℤ)^*. The output G depends on the value of flag:
* flag = 0 (default), an abelian group structure [h,d,g],
where h = φ(n) is the order (
* flag = 1 the result is a * flag = 2 same as flag = 1 with generators.
? G = znstar(40) %1 = [16, [4, 2, 2], [Mod(17, 40), Mod(21, 40), Mod(11, 40)]] ? G.no \\ eulerphi(40) %2 = 16 ? G.cyc \\ cycle structure %3 = [4, 2, 2] ? G.gen \\ generators for the cyclic components %4 = [Mod(17, 40), Mod(21, 40), Mod(11, 40)] ? apply(znorder, G.gen) %5 = [4, 2, 2]
According to the above definitions,
The library syntax is
| |