Pari/GP Reference Documentation | Contents - Global index - GP keyboard shortcuts |
O bezoutres deriv diffop eval factorpadic intformal padicappr padicfields polchebyshev polclass polcoeff polcyclo polcyclofactors poldegree poldisc poldiscreduced polgraeffe polhensellift polhermite polinterpolate poliscyclo poliscycloprod polisirreducible pollead pollegendre polmodular polrecip polresultant polresultantext polroots polrootsmod polrootspadic polrootsreal polsturm polsubcyclo polsylvestermatrix polsym poltchebi polzagier serconvol serlaplace serreverse subst substpol substvec sumformal taylor thue thueinit | |
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.
| |
O | |
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
| |
bezoutres | |
Deprecated alias for
The library syntax is
| |
deriv | |
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.
By definition, the main variable of a
The library syntax is
| |
diffop | |
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 Some examples: This function can be used to differentiate formal expressions: If E = exp(X^2) then we have E' = 2*X*E. We can derivate X*exp(X^2) as follow:
? diffop(E*X,[X,E],[1,2*X*E]) %1 = (2*X^2 + 1)*E
Let
? 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(var(i)=eval(Str("X",i))); my(x,v,dv); 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),polcoeff(x,n,'X)) }
The library syntax is
For n = 1, the function
| |
eval | |
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
If x is a character string,
? 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
? 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
? 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
| |
factorpadic | |
p-adic factorization
of the polynomial pol to precision r, the result being a
two-column matrix as in
? 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
The library syntax is
| |
intformal | |
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
? intformal(Mod(1,x^2+1), 'x) *** intformal: incorrect priority in intformal: variable x = x
The library syntax is
| |
padicappr | |
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
The library syntax is
| |
padicfields | |
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 non-isomorphic extensions; in particular, the number of polynomials is the number of classes of non-isomorphic 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
| |
polchebyshev | |
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 (
The library syntax is
| |
polclass | |
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
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 non-inert, 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 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
| |
polcoeff | |
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.
For greater flexibility, x can be a vector or matrix type and the
function then returns
The library syntax is
| |
polcyclo | |
n-th cyclotomic polynomial, evaluated at a (
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
The library syntax is
| |
polcyclofactors | |
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
| |
poldegree | |
Degree of the polynomial x in the main variable if v is omitted, in the variable v otherwise.
The degree of 0 is
The library syntax is
| |
poldisc | |
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
? z = 3*quadgen(8) + 4; ? poldisc(z) %2 = 8 ? q = Qfb(1,2,3); ? poldisc(q) %4 = -8
The library syntax is
| |
poldiscreduced | |
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
| |
polgraeffe | |
Returns the Graeffe transform g of f, such that g(x^2) = f(x) f(-x).
The library syntax is
| |
polhensellift | |
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
| |
polhermite | |
n-th Hermite polynomial H_n evaluated at a
(
The library syntax is
| |
polinterpolate | |
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]). If present, e will contain an error estimate on the returned value.
The library syntax is
| |
poliscyclo | |
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
| |
poliscycloprod | |
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
| |
polisirreducible | |
pol being a polynomial (univariate in the present version 2.9.2), returns 1 if pol is non-constant and irreducible, 0 otherwise. Irreducibility is checked over the smallest base field over which pol seems to be defined.
The library syntax is
| |
pollead | |
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
| |
pollegendre | |
n-th Legendre polynomial evaluated at a (
The library syntax is
| |
polmodular | |
Return the modular polynomial of prime level L in variables x and y
for the modular function specified by
? 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
| |
polrecip | |
Reciprocal polynomial of pol, i.e. the coefficients are in reverse order. pol must be a polynomial.
The library syntax is
| |
polresultant | |
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
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.
The library syntax is
| |
polresultantext | |
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
| |
polroots | |
Complex roots of the polynomial
x, given as a column vector where each root is repeated according to
its multiplicity. The precision is given as for transcendental functions: in
GP it is kept in the variable The algorithm used is a modification of A. Sch¨nhage's root-finding algorithm, due to and originally implemented by X. Gourdon. Barring bugs, it is guaranteed to converge and to give the roots to the required accuracy.
The library syntax is
| |
polrootsmod | |
Row vector of roots modulo p of the polynomial pol. Multiple roots are not repeated.
? polrootsmod(x^2-1,2) %1 = [Mod(1, 2)]~ If p is very small, you may set flag = 1, which uses a naive search.
The library syntax is
| |
polrootspadic | |
Vector of p-adic roots of the polynomial pol, given to p-adic precision r 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, 5) %1 = [2^3 + O(2^5), 2^3 + 2^4 + O(2^5), 1 + O(2^5)]~
If pol has inexact
The library syntax is
| |
polrootsreal | |
Real roots of the polynomial T with rational coefficients, multiple
roots being included according to their multiplicity. The roots are given
to a relative accuracy of
? \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"
( 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 non-negative 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
| |
polsturm | |
Number of real roots of the real squarefree polynomial T. If
the argument ab is present, it must be a vector [a,b] with
two real components (of type
If possible, you should stick to exact inputs, that is avoid
? 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 *** at top-level: polsturm(T^2,[0,4]) *** ^------------------- *** polsturm: domain error in polsturm: issquarefree(pol) = 0 ? polsturm((T*1.)^2, [0, 4]) \\ not squarefree AND inexact *** at top-level: polsturm((T*1.)^2,[0 *** ^-------------------- *** polsturm: precision too low in polsturm. In the last example, the input polynomial is not squarefree but there is no way to ascertain it from the given floating point approximation: we get a precision error in this case.
The library syntax is
| |
polsubcyclo | |
Gives polynomials (in variable v) defining the sub-Abelian extensions 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
The function
The library syntax is
| |
polsylvestermatrix | |
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
| |
polsym | |
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
| |
poltchebi | |
Deprecated alias for
The library syntax is
| |
polzagier | |
Creates Zagier's polynomial P_n^{(m)} used in
the functions 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
| |
serconvol | |
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
The library syntax is
| |
serlaplace | |
x must be a power series with non-negative exponents or a polynomial. If x = ∑ (a_k/k!)*X^k then the result is ∑ a_k*X^k.
The library syntax is
| |
serreverse | |
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
| |
subst | |
Replace the simple variable y by the argument z in the "polynomial" expression x. 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. Finally, if x is a vector, matrix or list, the substitution is applied to each individual entry.
Use the function
The library syntax is
| |
substpol | |
Replace the "variable" y by the argument z in the "polynomial"
expression x. Every type is allowed for x, but the same behavior
as
The difference with
? 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
| |
substvec | |
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]
The library syntax is
| |
sumformal | |
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
| |
taylor | |
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
? 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
| |
thue | |
Returns all solutions of the equation
P(x,y) = a in integers x and y, where tnf was created with
It is allowed to input directly the polynomial P instead of a tnf,
in which case, the function first performs
If tnf was computed without assuming GRH (flag 1 in * 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
? 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 non-monic 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.
The library syntax is
| |
thueinit | |
Initializes the tnf corresponding to P, a non-constant
univariate polynomial with integer coefficients.
The result is meant to be used in conjunction with
? 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 non-zero, certify results unconditionally. Otherwise, assume
GRH, this being much faster of course. In the latter case, the result
may still be unconditionally correct, see
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 use
? 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
? 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.
The library syntax is
| |