We group here all functions which are specific to polynomials or power series. Many other functions which can be applied on these objects are described in the other sections. Also, some of the functions described here can be applied to other types.

If p is an integer greater than 2, returns a p-adic 0 of precision e. In all other cases, returns a power series zero with precision given by e v, where v is the X-adic valuation of p with respect to its main variable.

The library syntax is `GEN `

.
**ggrando**()`GEN `

for a p-adic and
**zeropadic**(GEN p, long e)`GEN `

for a power series zero in variable v.**zeroser**(long v, long e)

Deprecated alias for `polresultantext`

The library syntax is `GEN `

where **polresultantext0**(GEN A, GEN B, long v = -1)`v`

is a variable number.

Derivative of x with respect to the main
variable if v is omitted, and with respect to v otherwise. The derivative
of a scalar type is zero, and the derivative of a vector or matrix is done
componentwise. One can use x' as a shortcut if the derivative is with
respect to the main variable of x; and also use x", etc., for multiple
derivatives altough `derivn`

is often preferrable.

By definition, the main variable of a `t_POLMOD`

is the main variable among
the coefficients from its two polynomial components (representative and
modulus); in other words, assuming a polmod represents an element of
R[X]/(T(X)), the variable X is a mute variable and the derivative is
taken with respect to the main variable used in the base ring R.

? f = (x/y)^5; ? deriv(f) %2 = 5/y^5*x^4 ? f' %3 = 5/y^5*x^4 ? deriv(f, 'x) \\ same since 'x is the main variable %4 = 5/y^5*x^4 ? deriv(f, 'y) %5 = -5/y^6*x^5

This function also operates on closures, in which case the variable
must be omitted. It returns a closure performing a numerical
differentiation as per `derivnum`

:

? f(x) = x^2; ? g = deriv(f) ? g(1) %3 = 2.0000000000000000000000000000000000000 ? f(x) = sin(exp(x)); ? deriv(f)(0) %5 = 0.54030230586813971740093660744297660373 ? cos(1) %6 = 0.54030230586813971740093660744297660373

The library syntax is `GEN `

where **deriv**(GEN x, long v = -1)`v`

is a variable number.

n-th derivative of x with respect to the main variable if v is omitted, and with respect to v otherwise; the integer n must be nonnegative. The derivative of a scalar type is zero, and the derivative of a vector or matrix is done componentwise. One can use x', x", etc., as a shortcut if the derivative is with respect to the main variable of x.

By definition, the main variable of a `t_POLMOD`

is the main variable among
the coefficients from its two polynomial components (representative and
modulus); in other words, assuming a polmod represents an element of
R[X]/(T(X)), the variable X is a mute variable and the derivative is
taken with respect to the main variable used in the base ring R.

? f = (x/y)^5; ? derivn(f, 2) %2 = 20/y^5*x^3 ? f'' %3 = 20/y^5*x^3 ? derivn(f, 2, 'x) \\ same since 'x is the main variable %4 = 20/y^5*x^3 ? derivn(f, 2, 'y) %5 = 30/y^7*x^5

This function also operates on closures, in which case the variable
must be omitted. It returns a closure performing a numerical
differentiation as per `derivnum`

:

? f(x) = x^10; ? g = derivn(f, 5) ? g(1) %3 = 30240.000000000000000000000000000000000 ? derivn(zeta, 2)(0) %4 = -2.0063564559085848512101000267299604382 ? zeta''(0) %5 = -2.0063564559085848512101000267299604382

The library syntax is `GEN `

where **derivn**(GEN x, long n, long v = -1)`v`

is a variable number.

Let v be a vector of variables, and d a vector of the same length,
return the image of x by the n-power (1 if n is not given) of the
differential operator D that assumes the value `d[i]`

on the variable
`v[i]`

. The value of D on a scalar type is zero, and D applies
componentwise to a vector or matrix. When applied to a `t_POLMOD`

, if no
value is provided for the variable of the modulus, such value is derived
using the implicit function theorem.

**Examples.**
This function can be used to differentiate formal expressions:
if E = exp(X^{2}) then we have E' = 2*X*E. We derivate X*exp(X^{2})
as follows:

? diffop(E*X,[X,E],[1,2*X*E]) %1 = (2*X^2 + 1)*E

Let `Sin`

and `Cos`

be two function such that
`Sin`

^{2}+`Cos`

^{2} = 1 and `Cos`

' = -`Sin`

.
We can differentiate `Sin`

/`Cos`

as follows,
PARI inferring the value of `Sin`

' from the equation:

? diffop(Mod('Sin/'Cos,'Sin^2+'Cos^2-1),['Cos],[-'Sin]) %1 = Mod(1/Cos^2, Sin^2 + (Cos^2 - 1))

Compute the Bell polynomials (both complete and partial) via the Faa di Bruno formula:

Bell(k,n=-1)= { my(x, v, dv, var = i->eval(Str("X",i))); v = vector(k, i, if (i==1, 'E, var(i-1))); dv = vector(k, i, if (i==1, 'X*var(1)*'E, var(i))); x = diffop('E,v,dv,k) / 'E; if (n < 0, subst(x,'X,1), polcoef(x,n,'X)); }

The library syntax is `GEN `

.**diffop0**(GEN x, GEN v, GEN d, long n)

For n = 1, the function `GEN `

is also
available.**diffop**(GEN x, GEN v, GEN d)

Replaces in x the formal variables by the values that
have been assigned to them after the creation of x. This is mainly useful
in GP, and not in library mode. Do not confuse this with substitution (see
`subst`

).

If x is a character string, `eval(x)`

executes x as a GP
command, as if directly input from the keyboard, and returns its
output.

? x1 = "one"; x2 = "two"; ? n = 1; eval(Str("x", n)) %2 = "one" ? f = "exp"; v = 1; ? eval(Str(f, "(", v, ")")) %4 = 2.7182818284590452353602874713526624978

Note that the first construct could be implemented in a
simpler way by using a vector `x = ["one","two"]; x[n]`

, and the second
by using a closure `f = exp; f(v)`

. The final example is more interesting:

? genmat(u,v) = matrix(u,v,i,j, eval( Str("x",i,j) )); ? genmat(2,3) \\ generic 2 x 3 matrix %2 = [x11 x12 x13] [x21 x22 x23]

A syntax error in the evaluation expression raises an `e_SYNTAX`

exception, which can be trapped as usual:

? 1a *** syntax error, unexpected variable name, expecting $end or ';': 1a *** ^- ? E(expr) = { iferr(eval(expr), e, print("syntax error"), errname(e) == "e_SYNTAX"); } ? E("1+1") %1 = 2 ? E("1a") syntax error

The library syntax is

.**geval**(GEN x)

p-adic factorization
of the polynomial *pol* to precision r, the result being a
two-column matrix as in `factor`

. Note that this is not the same
as a factorization over ℤ/p^{r}ℤ (polynomials over that ring do not form a
unique factorization domain, anyway), but approximations in ℚ/p^{r}ℤ of
the true factorization in ℚ_{p}[X].

? factorpadic(x^2 + 9, 3,5) %1 = [(1 + O(3^5))*x^2 + O(3^5)*x + (3^2 + O(3^5)) 1] ? factorpadic(x^2 + 1, 5,3) %2 = [ (1 + O(5^3))*x + (2 + 5 + 2*5^2 + O(5^3)) 1] [(1 + O(5^3))*x + (3 + 3*5 + 2*5^2 + O(5^3)) 1]

The factors are normalized so that their leading coefficient is a power of p. The method used is a modified version of the round 4 algorithm of Zassenhaus.

If *pol* has inexact `t_PADIC`

coefficients, this is not always
well-defined; in this case, the polynomial is first made integral by dividing
out the p-adic content, then lifted to ℤ using `truncate`

coefficientwise.
Hence we actually factor exactly a polynomial which is only p-adically
close to the input. To avoid pitfalls, we advise to only factor polynomials
with exact rational coefficients.

The library syntax is

. The function **factorpadic**(GEN f,GEN p, long r)`factorpadic0`

is
deprecated, provided for backward compatibility.

Let w = [1,z,...,z^{N-1}] from some primitive N-roots of unity z
where N is a power of 2, and P be a polynomial < N,
return the unnormalized discrete Fourier transform of P,
{ P(w[i]), 1 ≤ i ≤ N}. Also allow P to be a vector
[p_{0},...,p_{n}] representing the polynomial ∑_{i} p_{i} X^{i}.
Composing `fft`

and `fftinv`

returns N times the original input
coefficients.

? w = rootsof1(4); fft(w, x^3+x+1) %1 = [3, 1, -1, 1] ? fftinv(w, %) %2 = [4, 4, 0, 4] ? Polrev(%) / 4 %3 = x^3 + x + 1 ? w = powers(znprimroot(5),3); fft(w, x^3+x+1) %4 = [Mod(3,5),Mod(1,5),Mod(4,5),Mod(1,5)] ? fftinv(w, %) %5 = [Mod(4,5),Mod(4,5),Mod(0,5),Mod(4,5)]

The library syntax is `GEN `

.**FFT**(GEN w, GEN P)

Let w = [1,z,...,z^{N-1}] from some primitive N-roots of unity z
where N is a power of 2, and P be a polynomial < N,
return the unnormalized discrete Fourier transform of P,
{ P(1 / w[i]), 1 ≤ i ≤ N}. Also allow P to be a vector
[p_{0},...,p_{n}] representing the polynomial ∑_{i} p_{i} X^{i}.
Composing
`fft`

and `fftinv`

returns N times the original input coefficients.

? w = rootsof1(4); fft(w, x^3+x+1) %1 = [3, 1, -1, 1] ? fftinv(w, %) %2 = [4, 4, 0, 4] ? Polrev(%) / 4 %3 = x^3 + x + 1 ? N = 512; w = rootsof1(N); T = random(1000 * x^(N-1)); ? U = fft(w, T); time = 3 ms. ? V = vector(N, i, subst(T, 'x, w[i])); time = 65 ms. ? exponent(V - U) %7 = -97 ? round(Polrev(fftinv(w,U) / N)) == T %8 = 1

The library syntax is `GEN `

.**FFTinv**(GEN w, GEN P)

formal integration of x with respect to the variable v (wrt. the main variable if v is omitted). Since PARI cannot represent logarithmic or arctangent terms, any such term in the result will yield an error:

? intformal(x^2) %1 = 1/3*x^3 ? intformal(x^2, y) %2 = y*x^2 ? intformal(1/x) *** at top-level: intformal(1/x) *** ^ — — — — -- *** intformal: domain error in intformal: residue(series, pole) != 0

The argument x can be of any type. When x is a rational function, we assume that the base ring is an integral domain of characteristic zero.

By definition, the main variable of a `t_POLMOD`

is the main variable
among the coefficients from its two polynomial components
(representative and modulus); in other words, assuming a polmod represents an
element of R[X]/(T(X)), the variable X is a mute variable and the
integral is taken with respect to the main variable used in the base ring R.
In particular, it is meaningless to integrate with respect to the main
variable of `x.mod`

:

? intformal(Mod(1,x^2+1), 'x) *** intformal: incorrect priority in intformal: variable x = x

The library syntax is `GEN `

where **integ**(GEN x, long v = -1)`v`

is a variable number.

Vector of p-adic roots of the polynomial *pol* congruent to the
p-adic number a modulo p, and with the same p-adic precision as a.
The number a can be an ordinary p-adic number (type `t_PADIC`

, i.e. an
element of ℤ_{p}) or can be an integral element of a finite
*unramified* extension ℚ_{p}[X]/(T) of ℚ_{p}, given as a
`t_POLMOD`

`Mod`

(A,T) at least one of whose coefficients is a `t_PADIC`

and T
irreducible modulo p. In this case, the result is the vector of roots
belonging to the same extension of ℚ_{p} as a. The polynomial *pol*
should have exact coefficients; if not, its coefficients are first rounded
to ℚ or ℚ[X]/(T) and this is the polynomial whose roots we consider.

The library syntax is `GEN `

.
Also available is **padicappr**(GEN pol, GEN a)`GEN `

when a is a
**Zp_appr**(GEN f, GEN a)`t_PADIC`

.

Returns a vector of polynomials generating all the extensions of degree
N of the field ℚ_{p} of p-adic rational numbers; N is
allowed to be a 2-component vector [n,d], in which case we return the
extensions of degree n and discriminant p^{d}.

The list is minimal in the sense that two different polynomials generate
nonisomorphic extensions; in particular, the number of polynomials is the
number of classes of nonisomorphic extensions. If P is a polynomial in this
list, α is any root of P and K = ℚ_{p}(α), then α
is the sum of a uniformizer and a (lift of a) generator of the residue field
of K; in particular, the powers of α generate the ring of p-adic
integers of K.

If *flag* = 1, replace each polynomial P by a vector [P, e, f, d, c]
where e is the ramification index, f the residual degree, d the
valuation of the discriminant, and c the number of conjugate fields.
If *flag* = 2, only return the *number* of extensions in a fixed
algebraic closure (Krasner's formula), which is much faster.

The library syntax is `GEN `

.
Also available is
**padicfields0**(GEN p, GEN N, long flag)`GEN `

, which computes
extensions of ℚ**padicfields**(GEN p, long n, long d, long flag)_{p} of degree n and discriminant p^{d}.

Returns the n-th
Chebyshev polynomial of the first kind T_{n} (*flag* = 1) or the second
kind U_{n} (*flag* = 2), evaluated at a (`'x`

by default). Both series of
polynomials satisfy the 3-term relation
P_{n+1} = 2xP_{n} - P_{n-1},
and are determined by the initial conditions U_{0} = T_{0} = 1, T_{1} = x,
U_{1} = 2x. In fact T_{n}' = n U_{n-1} and, for all complex numbers z, we
have T_{n}(cos z) = cos (nz) and U_{n-1}(cos z) = sin(nz)/sin z.
If n ≥ 0, then these polynomials have degree n. For n < 0,
T_{n} is equal to T_{-n} and U_{n} is equal to -U_{-2-n}.
In particular, U_{-1} = 0.

The library syntax is `GEN `

.
Also available are
**polchebyshev_eval**(long n, long flag, GEN a = NULL)`GEN `

,
**polchebyshev**(long n, long flag, long v)`GEN `

and
**polchebyshev1**(long n, long v)`GEN `

for T**polchebyshev2**(long n, long v)_{n} and U_{n} respectively.

Return a polynomial in ℤ[x] generating the Hilbert class field for the
imaginary quadratic discriminant D. If inv is 0 (the default),
use the modular j-function and return the classical Hilbert polynomial,
otherwise use a class invariant. The following invariants correspond to
the different values of inv, where f denotes Weber's function
`weber`

, and w_{p,q} the double eta quotient given by
w_{p,q} = (η(x/p) η(x/q) )/(η(x) η(x/{pq}) )

The invariants w_{p,q} are not allowed unless they satisfy the following
technical conditions ensuring they do generate the Hilbert class
field and not a strict subfield:

***** if p ! = q, we need them both noninert, prime to the conductor of
ℤ[sqrt{D}]. Let P, Q be prime ideals above p and q; if both are
unramified, we further require that P^{± 1} Q^{± 1} be all distinct in
the class group of ℤ[sqrt{D}]; if both are ramified, we require that PQ
! = 1 in the class group.

***** if p = q, we want it split and prime to the conductor and
the prime ideal above it must have order ! = 1, 2, 4 in the class group.

Invariants are allowed under the additional conditions on D listed below.

***** 0 : j

***** 1 : f, D = 1 mod 8 and D = 1,2 mod 3;

***** 2 : f^{2}, D = 1 mod 8 and D = 1,2 mod 3;

***** 3 : f^{3}, D = 1 mod 8;

***** 4 : f^{4}, D = 1 mod 8 and D = 1,2 mod 3;

***** 5 : γ_{2} = j^{1/3}, D = 1,2 mod 3;

***** 6 : w_{2,3}, D = 1 mod 8 and D = 1,2 mod 3;

***** 8 : f^{8}, D = 1 mod 8 and D = 1,2 mod 3;

***** 9 : w_{3,3}, D = 1 mod 2 and D = 1,2 mod 3;

***** 10: w_{2,5}, D ! = 60 mod 80 and D = 1,2 mod 3;

***** 14: w_{2,7}, D = 1 mod 8;

***** 15: w_{3,5}, D = 1,2 mod 3;

***** 21: w_{3,7}, D = 1 mod 2 and 21 does not divide D

***** 23: w_{2,3}^{2}, D = 1,2 mod 3;

***** 24: w_{2,5}^{2}, D = 1,2 mod 3;

***** 26: w_{2,13}, D ! = 156 mod 208;

***** 27: w_{2,7}^{2}, D ! = 28 mod 112;

***** 28: w_{3,3}^{2}, D = 1,2 mod 3;

***** 35: w_{5,7}, D = 1,2 mod 3;

***** 39: w_{3,13}, D = 1 mod 2 and D = 1,2 mod 3;

The algorithm for computing the polynomial does not use the floating point approach, which would evaluate a precise modular function in a precise complex argument. Instead, it relies on a faster Chinese remainder based approach modulo small primes, in which the class invariant is only defined algebraically by the modular polynomial relating the modular function to j. So in fact, any of the several roots of the modular polynomial may actually be the class invariant, and more precise assertions cannot be made.

For instance, while `polclass(D)`

returns the minimal polynomial of
j(τ) with τ (any) quadratic integer for the discriminant D,
the polynomial returned by `polclass(D, 5)`

can be the minimal polynomial
of any of γ_{2} (τ), ζ_{3} γ_{2} (τ) or
ζ_{3}^{2} γ_{2} (τ), the three roots of the modular polynomial
j = γ_{2}^{3}, in which j has been specialised to j (τ).

The modular polynomial is given by
j = ((f^{24}-16)^{3} )/(f^{24}) for Weber's function f.

For the double eta quotients of level N = p q, all functions are covered
such that the modular curve X_{0}^{+} (N), the function field of which is
generated by the functions invariant under Γ^{0} (N) and the
Fricke-Atkin-Lehner involution, is of genus 0 with function field
generated by (a power of) the double eta quotient w.
This ensures that the full Hilbert class field (and not a proper subfield)
is generated by class invariants from these double eta quotients.
Then the modular polynomial is of degree 2 in j, and
of degree ψ (N) = (p+1)(q+1) in w.

? polclass(-163) %1 = x + 262537412640768000 ? polclass(-51, , 'z) %2 = z^2 + 5541101568*z + 6262062317568 ? polclass(-151,1) x^7 - x^6 + x^5 + 3*x^3 - x^2 + 3*x + 1

The library syntax is `GEN `

where **polclass**(GEN D, long inv, long x = -1)`x`

is a variable number.

Coefficient of degree n of the polynomial x, with respect to the main variable if v is omitted, with respect to v otherwise. If n is greater than the degree, the result is zero.

Naturally applies to scalars (polynomial of degree 0), as well as to rational functions whose denominator is a monomial. It also applies to power series: if n is less than the valuation, the result is zero. If it is greater than the largest significant degree, then an error message is issued.

The library syntax is `GEN `

where **polcoef**(GEN x, long n, long v = -1)`v`

is a variable number.

Deprecated alias for polcoef.

The library syntax is `GEN `

where **polcoef**(GEN x, long n, long v = -1)`v`

is a variable number.

n-th cyclotomic polynomial, evaluated at a (`'x`

by default). The
integer n must be positive.

Algorithm used: reduce to the case where n is squarefree; to compute the
cyclotomic polynomial, use Φ_{np}(x) = Φ_{n}(x^{p})/Φ(x); to compute
it evaluated, use Φ_{n}(x) = ∏_{d | n} (x^{d}-1)^{μ(n/d)}. In the
evaluated case, the algorithm assumes that a^{d} - 1 is either 0 or
invertible, for all d | n. If this is not the case (the base ring has
zero divisors), use `subst(polcyclo(n),x,a)`

.

The library syntax is `GEN `

.
The variant **polcyclo_eval**(long n, GEN a = NULL)`GEN `

returns the n-th
cyclotomic polynomial in variable v.**polcyclo**(long n, long v)

Returns a vector of polynomials, whose product is the product of distinct cyclotomic polynomials dividing f.

? f = x^10+5*x^8-x^7+8*x^6-4*x^5+8*x^4-3*x^3+7*x^2+3; ? v = polcyclofactors(f) %2 = [x^2 + 1, x^2 + x + 1, x^4 - x^3 + x^2 - x + 1] ? apply(poliscycloprod, v) %3 = [1, 1, 1] ? apply(poliscyclo, v) %4 = [4, 3, 10]

In general, the polynomials are products of cyclotomic polynomials and not themselves irreducible:

? g = x^8+2*x^7+6*x^6+9*x^5+12*x^4+11*x^3+10*x^2+6*x+3; ? polcyclofactors(g) %2 = [x^6 + 2*x^5 + 3*x^4 + 3*x^3 + 3*x^2 + 2*x + 1] ? factor(%[1]) %3 = [ x^2 + x + 1 1] [x^4 + x^3 + x^2 + x + 1 1]

The library syntax is `GEN `

.**polcyclofactors**(GEN f)

Degree of the polynomial x in the main variable if v is omitted, in the variable v otherwise.

The degree of 0 is `-oo`

. The degree of a nonzero scalar is 0.
Finally, when x is a nonzero polynomial or rational function, returns the
ordinary degree of x. Raise an error otherwise.

The library syntax is `GEN `

where **gppoldegree**(GEN x, long v = -1)`v`

is a variable number.
Also available is
`long `

, which returns **poldegree**(GEN x, long v)`-LONG_MAX`

if x = 0
and the degree as a `long`

integer.

Discriminant of the polynomial
*pol* in the main variable if v is omitted, in v otherwise. Uses a
modular algorithm over ℤ or ℚ, and the subresultant algorithm
otherwise.

? T = x^4 + 2*x+1; ? poldisc(T) %2 = -176 ? poldisc(T^2) %3 = 0

For convenience, the function also applies to types `t_QUAD`

and
`t_QFB`

:

? z = 3*quadgen(8) + 4; ? poldisc(z) %2 = 8 ? q = Qfb(1,2,3); ? poldisc(q) %4 = -8

The library syntax is `GEN `

where **poldisc0**(GEN pol, long v = -1)`v`

is a variable number.

Given a polynomial T with integer coefficients, return
[D, *faD*] where D is the discriminant of T and
*faD* is a cheap partial factorization of |D|: entries in its first
column are coprime and not perfect powers but need not be primes.
The factors are obtained by a combination of trial division, testing for
perfect powers, factorizations in coprimes, and computing Euclidean
remainder sequences for (T,T') modulo composite factors d of D
(which is likely to produce 0-divisors in ℤ/dℤ).
If *flag* is 1, finish the factorization using `factorint`

.

? T = x^3 - 6021021*x^2 + 12072210077769*x - 8092423140177664432; ? [D,faD] = poldiscfactors(T); print(faD); D [3, 3; 7, 2; 373, 2; 500009, 2; 24639061, 2] %2 = -27937108625866859018515540967767467 ? T = x^3 + 9*x^2 + 27*x - 125014250689643346789780229390526092263790263725; ? [D,faD] = poldiscfactors(T); print(faD) [2, 6; 3, 3; 125007125141751093502187, 4] ? [D,faD] = poldiscfactors(T, 1); print(faD) [2, 6; 3, 3; 500009, 12; 1000003, 4]

The library syntax is `GEN `

.**poldiscfactors**(GEN T, long flag)

Reduced discriminant vector of the (integral, monic) polynomial f. This is the vector of elementary divisors of ℤ[α]/f'(α)ℤ[α], where α is a root of the polynomial f. The components of the result are all positive, and their product is equal to the absolute value of the discriminant of f.

The library syntax is `GEN `

.**reduceddiscsmith**(GEN f)

Returns the monic polynomial in variable `v`

whose roots are the
components of the vector a with multiplicities, that is
∏_{i} (x - a_{i}).

? polfromroots([1,2,3]) %1 = x^3 - 6*x^2 + 11*x - 6 ? polfromroots([z, -z], 'y) %2 = y^2 - z^2

The library syntax is `GEN `

where **polfromroots**(GEN a, long v = -1)`v`

is a variable number.

Returns the Graeffe transform g of f, such that g(x^{2}) = f(x)
f(-x).

The library syntax is `GEN `

.**polgraeffe**(GEN f)

Given a prime p, an integral polynomial A whose leading coefficient
is a p-unit, a vector B of integral polynomials that are monic and
pairwise relatively prime modulo p, and whose product is congruent to
A/lc(A) modulo p, lift the elements of B to polynomials whose
product is congruent to A modulo p^{e}.

More generally, if T is an integral polynomial irreducible mod p, and
B is a factorization of A over the finite field 𝔽_{p}[t]/(T), you can
lift it to ℤ_{p}[t]/(T, p^{e}) by replacing the p argument with [p,T]:

? { T = t^3 - 2; p = 7; A = x^2 + t + 1; B = [x + (3*t^2 + t + 1), x + (4*t^2 + 6*t + 6)]; r = polhensellift(A, B, [p, T], 6) } %1 = [x + (20191*t^2 + 50604*t + 75783), x + (97458*t^2 + 67045*t + 41866)] ? liftall( r[1] * r[2] * Mod(Mod(1,p^6),T) ) %2 = x^2 + (t + 1)

The library syntax is `GEN `

.**polhensellift**(GEN A, GEN B, GEN p, long e)

n-th Hermite polynomial H_{n} evaluated at a
(`'x`

by default), i.e.
H_{n}(x) = (-1)^{n} e^{x^{2}} (d^{n})/(dx^{n})e^{-x^{2}}.
If *flag* is nonzero and n > 0, return [H_{n-1}(a), H_{n}(a)].

? polhermite(5) %1 = 32*x^5 - 160*x^3 + 120*x ? polhermite(5, -2) \\ H_{5}(-2) %2 = 16 ? polhermite(5,,1) %3 = [16*x^4 - 48*x^2 + 12, 32*x^5 - 160*x^3 + 120*x] ? polhermite(5,-2,1) %4 = [76, 16]

The library syntax is `GEN `

.
The variant **polhermite_eval0**(long n, GEN a = NULL, long flag)`GEN `

returns the n-th
Hermite polynomial in variable v. To obtain H**polhermite**(long n, long v)_{n}(a),
use `GEN `

.**polhermite_eval**(long n, GEN a)

Given the data vectors X and Y of the same length n (X containing the x-coordinates, and Y the corresponding y-coordinates), this function finds the interpolating polynomial P of minimal degree passing through these points and evaluates it at t. If Y is omitted, the polynomial P interpolates the (i,X[i]).

? v = [1, 2, 4, 8, 11, 13]; ? P = polinterpolate(v) \\ formal interpolation %1 = 7/120*x^5 - 25/24*x^4 + 163/24*x^3 - 467/24*x^2 + 513/20*x - 11 ? [ subst(P,'x,a) | a <- [1..6] ] %2 = [1, 2, 4, 8, 11, 13] ? polinterpolate(v,, 10) \\ evaluate at 10 %3 = 508 ? subst(P, x, 10) %4 = 508 ? P = polinterpolate([1,2,4], [9,8,7]) %5 = 1/6*x^2 - 3/2*x + 31/3 ? [subst(P, 'x, a) | a <- [1,2,4]] %6 = [9, 8, 7] ? P = polinterpolate([1,2,4], [9,8,7], 0) %7 = 31/3

If the goal is to extrapolate a function at a unique point, it is more efficient to use the t argument rather than interpolate formally then evaluate:

? x0 = 1.5; ? v = vector(20, i,random([-10,10])); ? for(i=1,10^3, subst(polinterpolate(v),'x, x0)) time = 352 ms. ? for(i=1,10^3, polinterpolate(v,,x0)) time = 111 ms. ? v = vector(40, i,random([-10,10])); ? for(i=1,10^3, subst(polinterpolate(v), 'x, x0)) time = 3,035 ms. ? for(i=1,10^3, polinterpolate(v,, x0)) time = 436 ms.

The threshold depends on the base field. Over small prime finite fields, interpolating formally first is more efficient

? bench(p, N, T = 10^3) = { my (v = vector(N, i, random(Mod(0,p)))); my (x0 = Mod(3, p), t1, t2); gettime(); for(i=1, T, subst(polinterpolate(v), 'x, x0)); t1 = gettime(); for(i=1, T, polinterpolate(v,, x0)); t2 = gettime(); [t1, t2]; } ? p = 101; ? bench(p, 4, 10^4) \\ both methods are equivalent %3 = [39, 40] ? bench(p, 40) \\ with 40 points formal is much faster %4 = [45, 355]

As the cardinality increases, formal interpolation requires more points to become interesting:

? p = nextprime(2^128); ? bench(p, 4) \\ formal is slower %3 = [16, 9] ? bench(p, 10) \\ formal has become faster %4 = [61, 70] ? bench(p, 100) \\ formal is much faster %5 = [1682, 9081] ? p = nextprime(10^500); ? bench(p, 4) \\ formal is slower %7 = [72, 354] ? bench(p, 20) \\ formal is still slower %8 = [1287, 962] ? bench(p, 40) \\ formal has become faster %9 = [3717, 4227] ? bench(p, 100) \\ faster but relatively less impressive %10 = [16237, 32335]

If t is a complex numeric value and e is present, e will contain an
error estimate on the returned value. More precisely, let P be the
interpolation polynomial on the given n points; there exist a subset
of n-1 points and Q the attached interpolation polynomial
such that e = `exponent`

(P(t) - Q(t)) (Neville's algorithm).

? f(x) = 1 / (1 + 25*x^2); ? x0 = 975/1000; ? test(X) = { my (P, e); P = polinterpolate(X, [f(x) | x <- X], x0, &e); [ exponent(P - f(x0)), e ]; } \\ equidistant nodes vs. Chebyshev nodes ? test( [-10..10] / 10 ) %4 = [6, 5] ? test( polrootsreal(polchebyshev(21)) ) %5 = [-15, -10] ? test( [-100..100] / 100 ) %7 = [93, 97] \\ P(x0) is way different from f(x0) ? test( polrootsreal(polchebyshev(201)) ) %8 = [-60, -55]

This is an example of Runge's phenomenon: increasing the number of equidistant nodes makes extrapolation much worse. Note that the error estimate is not a guaranteed upper bound (cf %4), but is reasonably tight in practice.

**Numerical stability.** The interpolation is performed in
a numerically stable way using ∏_{j ! = i} (X[i] - X[j]) instead of
Q'(X[i]) with Q = ∏_{i} (x - X[i]). Centering the interpolation
points X[i] around 0, thereby reconstructing P(x - m), for a suitable
m will further reduce the numerical error.

The library syntax is `GEN `

.**polint**(GEN X, GEN Y = NULL, GEN t = NULL, GEN *e = NULL)

P being a monic irreducible polynomial with integer coefficients,
return 0 if P is not a class polynomial for the j-invariant,
otherwise return the discriminant D < 0 such that `P = polclass(D)`

.

? polisclass(polclass(-47)) %1 = -47 ? polisclass(x^5+x+1) %2 = 0 ? apply(polisclass,factor(poldisc(polmodular(5)))[,1]) %3 = [-16,-4,-3,-11,-19,-64,-36,-24,-51,-91,-99,-96,-84]~

The library syntax is `long `

.**polisclass**(GEN P)

Returns 0 if f is not a cyclotomic polynomial, and n > 0 if f =
Φ_{n}, the n-th cyclotomic polynomial.

? poliscyclo(x^4-x^2+1) %1 = 12 ? polcyclo(12) %2 = x^4 - x^2 + 1 ? poliscyclo(x^4-x^2-1) %3 = 0

The library syntax is `long `

.**poliscyclo**(GEN f)

Returns 1 if f is a product of cyclotomic polynomial, and 0 otherwise.

? f = x^6+x^5-x^3+x+1; ? poliscycloprod(f) %2 = 1 ? factor(f) %3 = [ x^2 + x + 1 1] [x^4 - x^2 + 1 1] ? [ poliscyclo(T) | T <- %[,1] ] %4 = [3, 12] ? polcyclo(3) * polcyclo(12) %5 = x^6 + x^5 - x^3 + x + 1

The library syntax is `long `

.**poliscycloprod**(GEN f)

*pol* being a polynomial (univariate in the present version **2.17.0**),
returns 1 if *pol* is nonconstant and irreducible, 0 otherwise.
Irreducibility is checked over the smallest base field over which *pol*
seems to be defined.

The library syntax is `long `

.**polisirreducible**(GEN pol)

n-th Laguerre polynomial L^{(a)}_{n} of degree n and
parameter a evaluated at b (`'x`

by default), i.e.
L_{n}^{(a)}(x) =
(x^{-a}e^{×})/(n!) (d^{n})/(dx^{n})(e^{-x}x^{n+a}).
If *flag* is 1, return [L^{(a)}_{n-1}(b), L_{n}^{(a)}(b)].

The library syntax is `GEN `

.
To obtain the n-th Laguerre polynomial in variable v,
use **pollaguerre_eval0**(long n, GEN a = NULL, GEN b = NULL, long flag)`GEN `

. To obtain
L**pollaguerre**(long n, GEN a, GEN b, long v)^{(a)}_{n}(b), use `GEN `

.**pollaguerre_eval**(long n, GEN a, GEN b)

Leading coefficient of the polynomial or power series x. This is computed with respect to the main variable of x if v is omitted, with respect to the variable v otherwise.

The library syntax is `GEN `

where **pollead**(GEN x, long v = -1)`v`

is a variable number.

n-th Legendre polynomial P_{n} evaluated at a
(`'x`

by default), where
P_{n}(x) = (1)/(2^{n} n!) (d^{n})/(dx^{n})(x^{2}-1)^{n} .
If *flag* is 1, return [P_{n-1}(a), P_{n}(a)].

The library syntax is `GEN `

.
To obtain the n-th Legendre polynomial P**pollegendre_eval0**(long n, GEN a = NULL, long flag)_{n} in variable v,
use `GEN `

. To obtain P**pollegendre**(long n, long v)_{n}(a),
use `GEN `

.**pollegendre_eval**(long n, GEN a)

Return the modular polynomial of prime level L in variables x and y
for the modular function specified by `inv`

. If `inv`

is 0 (the
default), use the modular j function, if `inv`

is 1 use the
Weber-f function, and if `inv`

is 5 use γ_{2} =
sqrt[3]{j}.
See `polclass`

for the full list of invariants.
If x is given as `Mod(j, p)`

or an element j of
a finite field (as a `t_FFELT`

), then return the modular polynomial of
level L evaluated at j. If j is from a finite field and
`derivs`

is nonzero, then return a triple where the
last two elements are the first and second derivatives of the modular
polynomial evaluated at j.

? polmodular(3) %1 = x^4 + (-y^3 + 2232*y^2 - 1069956*y + 36864000)*x^3 + ... ? polmodular(7, 1, , 'J) %2 = x^8 - J^7*x^7 + 7*J^4*x^4 - 8*J*x + J^8 ? polmodular(7, 5, 7*ffgen(19)^0, 'j) %3 = j^8 + 4*j^7 + 4*j^6 + 8*j^5 + j^4 + 12*j^2 + 18*j + 18 ? polmodular(7, 5, Mod(7,19), 'j) %4 = Mod(1, 19)*j^8 + Mod(4, 19)*j^7 + Mod(4, 19)*j^6 + ... ? u = ffgen(5)^0; T = polmodular(3,0,,'j)*u; ? polmodular(3, 0, u,'j,1) %6 = [j^4 + 3*j^2 + 4*j + 1, 3*j^2 + 2*j + 4, 3*j^3 + 4*j^2 + 4*j + 2] ? subst(T,x,u) %7 = j^4 + 3*j^2 + 4*j + 1 ? subst(T',x,u) %8 = 3*j^2 + 2*j + 4 ? subst(T'',x,u) %9 = 3*j^3 + 4*j^2 + 4*j + 2

The library syntax is `GEN `

where **polmodular**(long L, long inv, GEN x = NULL, long y = -1, long derivs)`y`

is a variable number.

Reciprocal polynomial of *pol* with respect to its main variable,
i.e. the coefficients of the result are in reverse order; *pol* must be
a polynomial.

? polrecip(x^2 + 2*x + 3) %1 = 3*x^2 + 2*x + 1 ? polrecip(2*x + y) %2 = y*x + 2

The library syntax is `GEN `

.**polrecip**(GEN pol)

Resultant of the two
polynomials x and y with exact entries, with respect to the main
variables of x and y if v is omitted, with respect to the variable v
otherwise. The algorithm assumes the base ring is a domain. If you also need
the u and v such that x*u + y*v = Res(x,y), use the
`polresultantext`

function.

If *flag* = 0 (default), uses the algorithm best suited to the inputs,
either the subresultant algorithm (Lazard/Ducos variant, generic case),
a modular algorithm (inputs in ℚ[X]) or Sylvester's matrix (inexact
inputs).

If *flag* = 1, uses the determinant of Sylvester's matrix instead; this should
always be slower than the default.

If x or y are multivariate with a huge *polynomial* content, it
is advisable to remove it before calling this function. Compare:

? a = polcyclo(7) * ((t+1)/(t+2))^100; ? b = polcyclo(11)* ((t+2)/(t+3))^100); ? polresultant(a,b); time = 3,833 ms. ? ca = content(a); cb = content(b); \ polresultant(a/ca,b/cb)*ca^poldegree(b)*cb*poldegree(a); \\ instantaneous

The function only removes rational denominators and does
not compute automatically the content because it is generically small and
potentially *very* expensive (e.g. in multivariate contexts).
The choice is yours, depending on your application.

The library syntax is `GEN `

where **polresultant0**(GEN x, GEN y, long v = -1, long flag)`v`

is a variable number.

Finds polynomials U and V such that A*U + B*V = R, where R is the resultant of U and V with respect to the main variables of A and B if v is omitted, and with respect to v otherwise. Returns the row vector [U,V,R]. The algorithm used (subresultant) assumes that the base ring is a domain.

? A = x*y; B = (x+y)^2; ? [U,V,R] = polresultantext(A, B) %2 = [-y*x - 2*y^2, y^2, y^4] ? A*U + B*V %3 = y^4 ? [U,V,R] = polresultantext(A, B, y) %4 = [-2*x^2 - y*x, x^2, x^4] ? A*U+B*V %5 = x^4

The library syntax is `GEN `

where **polresultantext0**(GEN A, GEN B, long v = -1)`v`

is a variable number.
Also available is
`GEN `

.**polresultantext**(GEN x, GEN y)

Complex roots of the polynomial T, given as a column vector where each
root is repeated according to its multiplicity and given as floating point
complex numbers at the current `realprecision`

:

? polroots(x^2) %1 = [0.E-38 + 0.E-38*I, 0.E-38 + 0.E-38*I]~ ? polroots(x^3+1) %2 = [-1.00... + 0.E-38*I, 0.50... - 0.866...*I, 0.50... + 0.866...*I]~

The algorithm used is a modification of Schönhage's root-finding algorithm, due to and originally implemented by Gourdon. It runs in polynomial time in deg(T) and the precision. If furthermore T has rational coefficients, roots are guaranteed to the required relative accuracy. If the input polynomial T is exact, then the ordering of the roots does not depend on the precision: they are ordered by increasing |Im z|, then by increasing Re z; in case of tie (conjugates), the root with negative imaginary part comes first.

The library syntax is `GEN `

.**roots**(GEN T, long prec)

Return a sharp upper bound B for the modulus of
the largest complex root of the polynomial T with complex coefficients
with relative error τ. More precisely, we have |z| ≤ B for all roots
and there exist one root such that |z_{0}| ≥ B exp(-2τ). Much faster
than either polroots or polrootsreal.

? T=poltchebi(500); ? vecmax(abs(polroots(T))) time = 5,706 ms. %2 = 0.99999506520185816611184481744870013191 ? vecmax(abs(polrootsreal(T))) time = 1,972 ms. %3 = 0.99999506520185816611184481744870013191 ? polrootsbound(T) time = 217 ms. %4 = 1.0098792554165905155 ? polrootsbound(T, log(2)/2) \\ allow a factor 2, much faster time = 51 ms. %5 = 1.4065759938190154354 ? polrootsbound(T, 1e-4) time = 504 ms. %6 = 1.0000920717983847741 ? polrootsbound(T, 1e-6) time = 810 ms. %7 = 0.9999960628901692905 ? polrootsbound(T, 1e-10) time = 1,351 ms. %8 = 0.9999950652993869760

The library syntax is `GEN `

.**polrootsbound**(GEN T, GEN tau = NULL)

Obsolete, kept for backward compatibility: use factormod.

The library syntax is `GEN `

.**polrootsff**(GEN x, GEN p = NULL, GEN a = NULL)

Vector of roots 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(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.

Multiple roots are *not* repeated.

? polrootsmod(x^2-1,2) %1 = [Mod(1, 2)]~ ? polrootsmod(x^2+1,3) %2 = []~ ? polrootsmod(x^2+1, [y^2+1,3]) %3 = [Mod(Mod(1, 3)*y, Mod(1, 3)*y^2 + Mod(1, 3)), Mod(Mod(2, 3)*y, Mod(1, 3)*y^2 + Mod(1, 3))]~ ? polrootsmod(x^2 + Mod(1,3)) %4 = []~ ? liftall( polrootsmod(x^2 + Mod(Mod(1,3),y^2+1)) ) %5 = [y, 2*y]~ ? t = ffgen(y^2+Mod(1,3)); polrootsmod(x^2 + t^0) %6 = [y, 2*y]~

The library syntax is `GEN `

.**polrootsmod**(GEN f, GEN D = NULL)

Vector of p-adic roots of the polynomial *pol*, given to
p-adic precision r; the integer p is assumed to be a prime.
Multiple roots are
*not* repeated. Note that this is not the same as the roots in
ℤ/p^{r}ℤ, rather it gives approximations in ℤ/p^{r}ℤ of the true roots
living in ℚ_{p}:

? polrootspadic(x^3 - x^2 + 64, 2, 4) %1 = [2^3 + O(2^4), 2^3 + O(2^4), 1 + O(2^4)]~ ? polrootspadic(x^3 - x^2 + 64, 2, 5) %2 = [2^3 + O(2^5), 2^3 + 2^4 + O(2^5), 1 + O(2^5)]~

As the second commands show, the first two roots *are*
distinct in ℚ_{p}, even though they are equal modulo 2^{4}.

More generally, if T is an integral polynomial irreducible
mod p and f has coefficients in ℚ[t]/(T), the argument p
may be replaced by the vector [T,p]; we then return the roots of f in
the unramified extension ℚ_{p}[t]/(T).

? polrootspadic(x^3 - x^2 + 64*y, [y^2+y+1,2], 5) %3 = [Mod((2^3 + O(2^5))*y + (2^3 + O(2^5)), y^2 + y + 1), Mod((2^3 + 2^4 + O(2^5))*y + (2^3 + 2^4 + O(2^5)), y^2 + y + 1), Mod(1 + O(2^5), y^2 + y + 1)]~

If *pol* has inexact `t_PADIC`

coefficients, this need not
well-defined; in this case, the polynomial is first made integral by
dividing out the p-adic content, then lifted to ℤ using `truncate`

coefficientwise. Hence the roots given are approximations of the roots of an
exact polynomial which is p-adically close to the input. To avoid pitfalls,
we advise to only factor polynomials with exact rational coefficients.

The library syntax is `GEN `

.**polrootspadic**(GEN f, GEN p, long r)

Real roots of the polynomial T with real coefficients, multiple
roots being included according to their multiplicity. If the polynomial
does not have rational coefficients, it is first rescaled and rounded.
The roots are given to a relative accuracy of `realprecision`

.
If argument *ab* is
present, it must be a vector [a,b] with two components (of type
`t_INT`

, `t_FRAC`

or `t_INFINITY`

) and we restrict to roots belonging
to that closed interval.

? \p9 ? polrootsreal(x^2-2) %1 = [-1.41421356, 1.41421356]~ ? polrootsreal(x^2-2, [1,+oo]) %2 = [1.41421356]~ ? polrootsreal(x^2-2, [2,3]) %3 = []~ ? polrootsreal((x-1)*(x-2), [2,3]) %4 = [2.00000000]~

The algorithm used is a modification of Uspensky's method (relying on
Descartes's rule of sign), following Rouillier and Zimmerman's article
"Efficient isolation of a polynomial real roots"
(`https://hal.inria.fr/inria-00072518/`

). Barring bugs, it is guaranteed
to converge and to give the roots to the required accuracy.

**Remark.** If the polynomial T is of the
form Q(x^{h}) for some h ≥ 2 and *ab* is omitted, the routine will
apply the algorithm to Q (restricting to nonnegative roots when h is
even), then take h-th roots. On the other hand, if you want to specify
*ab*, you should apply the routine to Q yourself and a suitable
interval [a',b'] using approximate h-th roots adapted to your problem:
the function will not perform this change of variables if *ab* is present.

The library syntax is `GEN `

.**realroots**(GEN T, GEN ab = NULL, long prec)

Number of distinct real roots of the real polynomial *T*. If
the argument *ab* is present, it must be a vector [a,b] with
two real components (of type `t_INT`

, `t_REAL`

, `t_FRAC`

or `t_INFINITY`

) and we count roots belonging to that closed interval.

If possible, you should stick to exact inputs, that is avoid `t_REAL`

s in
T and the bounds a,b: the result is then guaranteed and we use a fast
algorithm (Uspensky's method, relying on Descartes's rule of sign, see
`polrootsreal`

). Otherwise, the polynomial is rescaled and rounded first
and the result may be wrong due to that initial error. If only a or b is
inexact, on the other hand, the interval is first thickened using rational
endpoints and the result remains guaranteed unless there exist a root
*very* close to a nonrational endpoint (which may be missed or unduly
included).

? T = (x-1)*(x-2)*(x-3); ? polsturm(T) %2 = 3 ? polsturm(T, [-oo,2]) %3 = 2 ? polsturm(T, [1/2,+oo]) %4 = 3 ? polsturm(T, [1, Pi]) \\ Pi inexact: not recommended ! %5 = 3 ? polsturm(T*1., [0, 4]) \\ T*1. inexact: not recommended ! %6 = 3 ? polsturm(T^2, [0, 4]) \\ not squarefree: roots are not repeated! %7 = 3

The library syntax is `long `

or
**RgX_sturmpart**(GEN T, GEN ab)`long `

(for the case **sturm**(GEN T)`ab = NULL`

). The function
`long `

is obsolete and deprecated.**sturmpart**(GEN T, GEN a, GEN b)

Gives polynomials (in variable v) defining the (Abelian) subextensions
of degree d of the cyclotomic field ℚ(ζ_{n}), where d | φ(n).

If there is exactly one such extension the output is a polynomial, else it is
a vector of polynomials, possibly empty. To get a vector in all cases,
use `concat([], polsubcyclo(n,d))`

.

Each such polynomial is the minimal polynomial for a Gaussian period
Tr_{ℚ(ζ_{f})/L} (ζ_{f}), where L is the degree d
subextension of ℚ(ζ_{n}) and f | n is its conductor. In
Galois-theoretic terms, L = ℚ(ζ_{n})^{H}, where H runs through all
index d subgroups of (ℤ/nℤ)^{*}.

The function `galoissubcyclo`

allows to specify exactly which
sub-Abelian extension should be computed by giving H.

**Complexity.** Ignoring logarithmic factors, `polsubcyclo`

runs
in time O(n). The function `polsubcyclofast`

returns different, less
canonical, polynomials but runs in time O(d^{4}), again ignoring logarithmic
factors; thus it can handle much larger values of n.

The library syntax is `GEN `

where **polsubcyclo**(long n, long d, long v = -1)`v`

is a variable number.

If 1 ≤ d ≤ 6 or a prime, finds an equation for the subfields of
ℚ(ζ_{n}) with Galois group C_{d}; the special value d = -4 provides
the subfields with group V_{4} = C_{2} x C_{2}. Contrary to
`polsubcyclo`

, the
output is always a (possibly empty) vector of polynomials. If s = 0 (default)
all signatures, otherwise s = 1 (resp., -1) for totally real
(resp., totally complex). Set `exact = 1`

for subfields of conductor n.

The argument n can be given as in arithmetic functions: as an integer, as a
factorization matrix, or (preferred) as a pair [N, `factor`

(N)].

**Comparison with polsubcyclo.** First

`polsubcyclofast`

does not usually return Gaussian periods, but ad hoc polynomials which do
generate the same field. Roughly speaking (ignoring
logarithmic factors), the complexity of `polsubcyclo`

is independent of
d and the complexity of `polsubcyclofast`

is independent of n.
Ignoring logarithmic factors, `polsubcylo`

runs in time O(n) and
`polsubcyclofast`

in time O(d`polsubcyclo`

if n is large,
but gets slower as d increases and becomes unusable for d ≥ 40 or so.

? polsubcyclo(10^7+19,7); time = 1,852 ms. ? polsubcyclofast(10^7+19,7); time = 15 ms. ? polsubcyclo(10^17+21,5); \\ won't finish *** polsubcyclo: user interrupt after 2h ? polsubcyclofast(10^17+21,5); time = 3 ms. ? polsubcyclofast(10^17+3,7); time = 26 ms. ? polsubcyclo(10^6+117,13); time = 193 ms. ? polsubcyclofast(10^6+117,13); time = 50 ms. ? polsubcyclofast(10^6+199,19); time = 202 ms. ? polsubcyclo(10^6+199,19); \\ about as fast time = 3191ms. ? polsubcyclo(10^7+271,19); time = 2,067 ms. ? polsubcyclofast(10^7+271,19); time = 201 ms.

The library syntax is `GEN `

.**polsubcyclofast**(GEN n, long d, long s, long exact)

Forms the Sylvester matrix corresponding to the two polynomials x and y, where the coefficients of the polynomials are put in the columns of the matrix (which is the natural direction for solving equations afterwards). The use of this matrix can be essential when dealing with polynomials with inexact entries, since polynomial Euclidean division doesn't make much sense in this case.

The library syntax is `GEN `

.**sylvestermatrix**(GEN x, GEN y)

Creates the column vector of the symmetric powers of the roots of the polynomial x up to power n, using Newton's formula.

The library syntax is `GEN `

.**polsym**(GEN x, long n)

Deprecated alias for `polchebyshev`

The library syntax is `GEN `

where **polchebyshev1**(long n, long v = -1)`v`

is a variable number.

Given T ∈ 𝔽_{p}[X] return the polynomial P ∈ ℤ_{p}[X] whose roots
(resp. leading coefficient) are the Teichmuller lifts of the roots
(resp. leading coefficient) of T, to p-adic precision r. If T is
monic, P is the reduction modulo p^{r} of the unique monic polynomial
congruent to T modulo p such that P(X^{p}) = 0 (mod P(X),p^{r}).

? T = ffinit(3, 3, 't) %1 = Mod(1,3)*t^3 + Mod(1,3)*t^2 + Mod(1,3)*t + Mod(2,3) ? P = polteichmuller(T,3,5) %2 = t^3 + 166*t^2 + 52*t + 242 ? subst(P, t, t^3) % (P*Mod(1,3^5)) %3 = Mod(0, 243) ? [algdep(a+O(3^5),2) | a <- Vec(P)] %4 = [x - 1, 5*x^2 + 1, x^2 + 4*x + 4, x + 1]

When T is monic and irreducible mod p, this provides
a model ℚ_{p}[X]/(P) of the unramified extension ℚ_{p}[X] / (T) where
the Frobenius has the simple form X mod P ` ⟼ `

X^{p} mod P.

The library syntax is `GEN `

.**polteichmuller**(GEN T, ulong p, long r)

Let T ∈ ℚ[x] be a nonzero polynomial; returns U monic in ℤ[x]
such that U(x) = C T(x/L) for some C,L ∈ ℚ. If the pointer argument
`&L`

is present, set `L`

to L.

? poltomonic(9*x^2 - 1/2) %1 = x^2 - 2 ? U = poltomonic(9*x^2 - 1/2, &L) %2 = x^2 - 2 ? L %3 = 6 ? U / subst(9*x^2 - 1/2, x, x/L) %4 = 4

This function does not compute discriminants or maximal orders and runs
with complexity almost linear in the input size. If T is already monic with
integer coefficient, `poltomonic`

may still transform it if ℤ[x]/(T)
is contained in a trivial subring of the maximal order, generated by L x:

? poltomonic(x^2 + 4, &L) %5 = x^2 + 1 ? L %6 = 1/2

If T is irreducible, the functions `polredabs`

(exponential time) and `polredbest`

(polynomial time) also find a monic
integral generating polynomial for the number field ℚ[x]/(T), with
explicit guarantees on its size, but are orders of magnitude slower.

The library syntax is `GEN `

.**poltomonic**(GEN T, GEN *L = NULL)

Creates Zagier's polynomial P_{n}^{(m)} used in
the functions `sumalt`

and `sumpos`

(with *flag* = 1), see
"Convergence acceleration of alternating series", Cohen et al.,
*Experiment. Math.*, vol. 9, 2000, pp. 3–12.

If m < 0 or m ≥ n, P_{n}^{(m)} = 0.
We have
P_{n} := P_{n}^{(0)} is T_{n}(2x-1), where T_{n} is the Legendre
polynomial of the second kind. For n > m > 0, P_{n}^{(m)} is the m-th
difference with step 2 of the sequence n^{m+1}P_{n}; in this case, it
satisfies
2 P_{n}^{(m)}(sin^{2} t)
= (d^{m+1})/(dt^{m+1}) (sin(2t)^{m} sin(2(n-m)t)).

The library syntax is `GEN `

.**polzag**(long n, long m)

finds a linear relation between powers (1,s,
..., s^{p}) of the series s, with polynomial coefficients of degree
≤ r. In case no relation is found, return 0.

? s = 1 + 10*y - 46*y^2 + 460*y^3 - 5658*y^4 + 77740*y^5 + O(y^6); ? seralgdep(s, 2, 2) %2 = -x^2 + (8*y^2 + 20*y + 1) ? subst(%, x, s) %3 = O(y^6) ? seralgdep(s, 1, 3) %4 = (-77*y^2 - 20*y - 1)*x + (310*y^3 + 231*y^2 + 30*y + 1) ? seralgdep(s, 1, 2) %5 = 0

The series main variable must not be x, so as to be able to express the result as a polynomial in x.

The library syntax is `GEN `

.**seralgdep**(GEN s, long p, long r)

Convolution (or Hadamard product) of the
two power series x and y; in other words if x = ∑ a_{k}*X^{k}
and y = ∑ b_{k}*X^{k} then `serconvol`

(x,y) = ∑ a_{k}*b_{k}*X^{k}.

The library syntax is `GEN `

.**convol**(GEN x, GEN y)

Find a linear relation between the derivatives (s, s',..., s^{p}) of
the series s and 1, with polynomial coefficients of degree ≤ r. In
case no relation is found, return 0, otherwise return [E,P] such that
E(d)(S) = P where d is the standard derivation.

? S = sum(i=0, 50, binomial(3*i,i)*T^i) + O(T^51); ? serdiffdep(S, 3, 3) %2 = [(27*T^2 - 4*T)*x^2 + (54*T - 2)*x + 6, 0] ? (27*T^2 - 4*T)*S'' + (54*T - 2)*S' + 6*S %3 = O(T^50) ? S = exp(T^2) + T^2; ? serdiffdep(S, 3, 3) %5 = [x-2*T, -2*T^3+2*T] ? S'-2*T*S %6 = 2*T-2*T^3+O(T^17)

The series main variable must not be x, so as to be able to express the result as a polynomial in x.

The library syntax is `GEN `

.**serdiffdep**(GEN s, long p, long r)

x must be a power series with nonnegative
exponents or a polynomial. If x = ∑ (a_{k}/k!)*X^{k} then the result isi
∑ a_{k}*X^{k}.

The library syntax is `GEN `

.**laplace**(GEN x)

Reverse power series of s, i.e. the series t such that t(s) = x; s must be a power series whose valuation is exactly equal to one.

? \ps 8 ? t = serreverse(tan(x)) %2 = x - 1/3*x^3 + 1/5*x^5 - 1/7*x^7 + O(x^8) ? tan(t) %3 = x + O(x^8)

The library syntax is `GEN `

.**serreverse**(GEN s)

Replace the simple variable y by the argument z in the "polynomial"
expression x. If z is a vector, return the vector of the evaluated
expressions `subst(x, y, z[i])`

.

Every type is allowed for x, but if it is not a genuine polynomial (or power series, or rational function), the substitution will be done as if the scalar components were polynomials of degree zero. In particular, beware that:

? subst(1, x, [1,2; 3,4]) %1 = [1 0] [0 1] ? subst(1, x, Mat([0,1])) *** at top-level: subst(1,x,Mat([0,1]) *** ^ — — — — — — -- *** subst: forbidden substitution by a non square matrix.

If x is a power series, z must be either a polynomial, a power series, or a rational function. If x is a vector, matrix or list, the substitution is applied to each individual entry.

Use the function `substvec`

to replace several variables at once,
or the function `substpol`

to replace a polynomial expression.

The library syntax is `GEN `

where **gsubst**(GEN x, long y, GEN z)`y`

is a variable number.

Replace the "variable" y by the argument z in the "polynomial"
expression x. Every type is allowed for x, but the same behavior
as `subst`

above apply.

The difference with `subst`

is that y is allowed to be any polynomial
here. The substitution is done moding out all components of x
(recursively) by y - t, where t is a new free variable of lowest
priority. Then substituting t by z in the resulting expression. For
instance

? substpol(x^4 + x^2 + 1, x^2, y) %1 = y^2 + y + 1 ? substpol(x^4 + x^2 + 1, x^3, y) %2 = x^2 + y*x + 1 ? substpol(x^4 + x^2 + 1, (x+1)^2, y) %3 = (-4*y - 6)*x + (y^2 + 3*y - 3)

The library syntax is `GEN `

.
Further, **gsubstpol**(GEN x, GEN y, GEN z)`GEN `

attempts to
write T(x) in the form t(x**gdeflate**(GEN T, long v, long d)^{d}), where x = `pol`

(v), and returns
_{x}`NULL`

if the substitution fails (for instance in the example `%2`

above).

v being a vector of monomials of degree 1 (variables),
w a vector of expressions of the same length, replace in the expression
x all occurrences of v_{i} by w_{i}. The substitutions are done
simultaneously; more precisely, the v_{i} are first replaced by new
variables in x, then these are replaced by the w_{i}:

? substvec([x,y], [x,y], [y,x]) %1 = [y, x] ? substvec([x,y], [x,y], [y,x+y]) %2 = [y, x + y] \\ not [y, 2*y]

As in `subst`

, variables may be replaced
by a vector of values, in which case the cartesian product is returned:

? substvec([x,y], [x,y], [[1,2], 3]) %3 = [[1, 3], [2, 3]] ? substvec([x,y], [x,y], [[1,2], [3,4]]) %4 = [[1, 3], [2, 3], [1, 4], [2, 4]]

The library syntax is `GEN `

.**gsubstvec**(GEN x, GEN v, GEN w)

formal sum of the polynomial expression f with respect to the
main variable if v is omitted, with respect to the variable v otherwise;
it is assumed that the base ring has characteristic zero. In other words,
considering f as a polynomial function in the variable v,
returns F, a polynomial in v vanishing at 0, such that F(b) - F(a)
= sum_{v = a+1}^{b} f(v):

? sumformal(n) \\ 1 + ... + n %1 = 1/2*n^2 + 1/2*n ? f(n) = n^3+n^2+1; ? F = sumformal(f(n)) \\ f(1) + ... + f(n) %3 = 1/4*n^4 + 5/6*n^3 + 3/4*n^2 + 7/6*n ? sum(n = 1, 2000, f(n)) == subst(F, n, 2000) %4 = 1 ? sum(n = 1001, 2000, f(n)) == subst(F, n, 2000) - subst(F, n, 1000) %5 = 1 ? sumformal(x^2 + x*y + y^2, y) %6 = y*x^2 + (1/2*y^2 + 1/2*y)*x + (1/3*y^3 + 1/2*y^2 + 1/6*y) ? x^2 * y + x * sumformal(y) + sumformal(y^2) == % %7 = 1

The library syntax is `GEN `

where **sumformal**(GEN f, long v = -1)`v`

is a variable number.

Taylor expansion around 0 of x with respect to
the simple variable t. x can be of any reasonable type, for example a
rational function. Contrary to `Ser`

, which takes the valuation into
account, this function adds O(t^{d}) to all components of x.

? taylor(x/(1+y), y, 5) %1 = (y^4 - y^3 + y^2 - y + 1)*x + O(y^5) ? Ser(x/(1+y), y, 5) *** at top-level: Ser(x/(1+y),y,5) *** ^ — — — — — - *** Ser: main variable must have higher priority in gtoser.

The library syntax is `GEN `

where **tayl**(GEN x, long t, long precdl)`t`

is a variable number.

Returns all solutions of the equation
P(x,y) = a in integers x and y, where *tnf* was created with
`thueinit`

(P). If present, *sol* must contain the solutions of
Norm(x) = a modulo units of positive norm in the number field
defined by P (as computed by `bnfisintnorm`

). If there are infinitely
many solutions, an error is issued.

It is allowed to input directly the polynomial P instead of a *tnf*,
in which case, the function first performs `thueinit(P,0)`

. This is
very wasteful if more than one value of a is required.

If *tnf* was computed without assuming GRH (flag 1 in `thueinit`

),
then the result is unconditional. Otherwise, it depends in principle of the
truth of the GRH, but may still be unconditionally correct in some
favorable cases. The result is conditional on the GRH if
a ! = ± 1 and P has a single irreducible rational factor, whose
attached tentative class number h and regulator R (as computed
assuming the GRH) satisfy

***** h > 1,

***** R/0.2 > 1.5.

Here's how to solve the Thue equation x^{13} - 5y^{13} = - 4:

? tnf = thueinit(x^13 - 5); ? thue(tnf, -4) %1 = [[1, 1]]

In this case, one checks that `bnfinit(x^13 -5).no`

is 1. Hence, the only solution is (x,y) = (1,1) and the result is
unconditional. On the other hand:

? P = x^3-2*x^2+3*x-17; tnf = thueinit(P); ? thue(tnf, -15) %2 = [[1, 1]] \\ a priori conditional on the GRH. ? K = bnfinit(P); K.no %3 = 3 ? K.reg %4 = 2.8682185139262873674706034475498755834

This time the result is conditional. All results computed using this
particular *tnf* are likewise conditional, *except* for a right-hand
side of ± 1.
The above result is in fact correct, so we did not just disprove the GRH:

? tnf = thueinit(x^3-2*x^2+3*x-17, 1 /*unconditional*/); ? thue(tnf, -15) %4 = [[1, 1]]

Note that reducible or nonmonic polynomials are allowed:

? tnf = thueinit((2*x+1)^5 * (4*x^3-2*x^2+3*x-17), 1); ? thue(tnf, 128) %2 = [[-1, 0], [1, 0]]

Reducible polynomials are in fact much easier to handle.

**Note.** When P is irreducible without a real root, the default
strategy is to use brute force enumeration in time |a|^{1/deg P} and
avoid computing a tough *bnf* attached to P, see `thueinit`

.
Besides reusing a quantity you might need for other purposes, the
default argument *sol* can also be used to use a different strategy
and prove that there are no solutions; of course you need to compute a
*bnf* on you own to obtain *sol*. If there *are* solutions
this won't help unless P is quadratic, since the enumeration will be
performed in any case.

The library syntax is `GEN `

.**thue**(GEN tnf, GEN a, GEN sol = NULL)

Initializes the *tnf* corresponding to P, a nonconstant
univariate polynomial with integer coefficients.
The result is meant to be used in conjunction with `thue`

to solve Thue
equations P(X / Y)Y^{deg P} = a, where a is an integer. Accordingly,
P must either have at least two distinct irreducible factors over ℚ,
or have one irreducible factor T with degree > 2 or two conjugate
complex roots: under these (necessary and sufficient) conditions, the
equation has finitely many integer solutions.

? S = thueinit(t^2+1); ? thue(S, 5) %2 = [[-2, -1], [-2, 1], [-1, -2], [-1, 2], [1, -2], [1, 2], [2, -1], [2, 1]] ? S = thueinit(t+1); *** at top-level: thueinit(t+1) *** ^ — — — — - *** thueinit: domain error in thueinit: P = t + 1

The hardest case is when deg P > 2 and P is irreducible with at least one real root. The routine then uses Bilu-Hanrot's algorithm.

If *flag* is nonzero, certify results unconditionally. Otherwise, assume
GRH, this being much faster of course. In the latter case, the result
may still be unconditionally correct, see `thue`

. For instance in most
cases where P is reducible (not a pure power of an irreducible), *or*
conditional computed class groups are trivial *or* the right hand side
is ±1, then results are unconditional.

**Note.** The general philosophy is to disprove the existence of large
solutions then to enumerate bounded solutions naively. The implementation
will overflow when there exist huge solutions and the equation has degree
> 2 (the quadratic imaginary case is special, since we can stick to
`bnfisintnorm`

, there are no fundamental units):

? thue(t^3+2, 10^30) *** at top-level: L=thue(t^3+2,10^30) *** ^ — — — — — -- *** thue: overflow in thue (SmallSols): y <= 80665203789619036028928. ? thue(x^2+2, 10^30) \\ quadratic case much easier %1 = [[-1000000000000000, 0], [1000000000000000, 0]]

**Note.** It is sometimes possible to circumvent the above, and in any
case obtain an important speed-up, if you can write P = Q(x^{d}) for some
d > 1 and Q still satisfying the `thueinit`

hypotheses. You can then
solve
the equation attached to Q then eliminate all solutions (x,y) such that
either x or y is not a d-th power.

? thue(x^4+1, 10^40); \\ stopped after 10 hours ? filter(L,d) = my(x,y); [[x,y] | v<-L, ispower(v[1],d,&x)&&ispower(v[2],d,&y)]; ? L = thue(x^2+1, 10^40); ? filter(L, 2) %4 = [[0, 10000000000], [10000000000, 0]]

The last 2 commands use less than 20ms.

**Note.** When P is irreducible without a real root, the equation
can be solved unconditionnally in time |a|^{1/deg P}. When this
latter quantity is huge and the equation has no solutions, this fact
may still be ascertained via arithmetic conditions but this now implies
solving norm equations, computing a *bnf* and possibly assuming the GRH.
When there is no real root, the code does not compute a *bnf*
(with certification if *flag* = 1) if it expects this to be an "easy"
computation (because the result would only be used for huge values of a).
See `thue`

for a way to compute an expensive *bnf* on your own and
still get a result where this default cheap strategy fails.

The library syntax is `GEN `

.**thueinit**(GEN P, long flag, long prec)