Any nonzero value is interpreted as true and any zero as false (this includes empty vectors or matrices). The standard boolean operators || (inclusive or), && (and) and ! in prefix notation (not) are available. Their value is 1 (true) or 0 (false):

  ? a && b  \\ 1 iff a and b are nonzero
  ? a || b  \\ 1 iff a or b is nonzero
  ? !a      \\ 1 iff a is zero

The standard real comparison operators <= , < , >= , > , are available in GP. The result is 1 if the comparison is true, 0 if it is false. These operators allow to compare integers (t_INT), rational (t_FRAC) or real (t_REAL) numbers, real quadratic numbers (t_QUAD of positive discriminant) and infinity (oo, t_INFINITY).

By extension, two character strings (t_STR) are compared using the standard lexicographic order. Comparing a string to an object of a different type raises an exception. See also the cmp universal comparison function.

Two operators allow to test for equality: == (equality up to type coercion) and === (identity). The result is 1 if equality is decided, else 0.

The operator === is strict: objects of different type or length are never identical, polynomials in different variables are never identical, even if constant. On the contrary, == is very liberal: a  ==  b decides whether there is a natural map sending a to the domain of b or sending b to the domain of a, such that the comparison makes sense and equality holds. For instance

  ? 4 == Mod(1,3) \\ equal
  %1 = 1
  ? 4 === Mod(1,3) \\ but not identical
  %2 = 0
  
  ? 'x == 'y   \\ not equal (nonconstant and different variables)
  %3 = 0
  ? Pol(0,'x) == Pol(0,'y)  \\ equal (constant: ignore variable)
  %4 = 1
  ? Pol(0,'x) === Pol(0,'y)  \\ not identical
  %5 = 0
  
  ? 0 == Pol(0) \\ equal (not identical)
  %6 = 1
  ? [0] == 0     \\ equal (not identical)
  %7 = 1
  ? [0, 0] == 0  \\ equal (not identical)
  %8 = 1
  ? [0] == [0,0] \\ not equal
  %9 = 0

In particular == is not transitive in general. The operator === is transitive. The == operator allows two equivalent negated forms: != or <> ; there is no negated form for === .

Do not mistake = for == : the former is the assignment statement.

The expressions +x and -x refer to monadic operators: the first does nothing, the second negates x.

The library syntax is GEN gneg(GEN x) for -x.

The expression x + y is the sum of x and y. Addition between a scalar type x and a t_COL or t_MAT y returns respectively [y[1] + x, y[2],...] and y + x Id. Other additions between a scalar type and a vector or a matrix, or between vector/matrices of incompatible sizes are forbidden.

The library syntax is GEN gadd(GEN x, GEN y).

The expression x - y is the difference of x and y. Subtraction between a scalar type x and a t_COL or t_MAT y returns respectively [y[1] - x, y[2],...] and y - x Id. Other subtractions between a scalar type and a vector or a matrix, or between vector/matrices of incompatible sizes are forbidden.

The library syntax is GEN gsub(GEN x, GEN y) for x - y.

The expression x * y is the product of x and y. Among the prominent impossibilities are multiplication between vector/matrices of incompatible sizes, between a t_INTMOD or t_PADIC. Restricted to scalars, * is commutative; because of vector and matrix operations, it is not commutative in general.

Multiplication between two t_VECs or two t_COLs is not allowed; to take the scalar product of two vectors of the same length, transpose one of the vectors (using the operator ~ or the function mattranspose, see Section se:linear_algebra) and multiply a row vector by a column vector:

  ? a = [1,2,3];
  ? a * a
    ***   at top-level: a*a
    ***                  ^--
    *** _*_: forbidden multiplication t_VEC * t_VEC.
  ? a * a~
  %2 = 14

If x,y are binary quadratic forms, compose them; see also qfbnucomp and qfbnupow. If x,y are t_VECSMALL of the same length, understand them as permutations and compose them.

The library syntax is GEN gmul(GEN x, GEN y) for x * y. Also available is GEN gsqr(GEN x) for x * x.

The expression x / y is the quotient of x and y. In addition to the impossibilities for multiplication, note that if the divisor is a matrix, it must be an invertible square matrix, and in that case the result is x*y^-1. Furthermore note that the result is as exact as possible: in particular, division of two integers always gives a rational number (which may be an integer if the quotient is exact) and not the Euclidean quotient (see x \ y for that), and similarly the quotient of two polynomials is a rational function in general. To obtain the approximate real value of the quotient of two integers, add 0. to the result; to obtain the approximate p-adic value of the quotient of two integers, add O(p^k) to the result; finally, to obtain the Taylor series expansion of the quotient of two polynomials, add O(X^k) to the result or use the taylor function (see Section se:taylor).

The library syntax is GEN gdiv(GEN x, GEN y) for x / y.

The expression x \y is the Euclidean quotient of x and y. If y is a real scalar, this is defined as floor(x/y) if y > 0, and ceil(x/y) if y < 0 and the division is not exact. Hence the remainder x - (x\y)*y is in [0, |y|[.

Note that when y is an integer and x a polynomial, y is first promoted to a polynomial of degree 0. When x is a vector or matrix, the operator is applied componentwise.

The library syntax is GEN gdivent(GEN x, GEN y) for x \ y.

The expression x \/ y evaluates to the rounded Euclidean quotient of x and y. This is the same as x \y except for scalar division: the quotient is such that the corresponding remainder is smallest in absolute value and in case of a tie the quotient closest to + oo is chosen (hence the remainder would belong to [{-}|y|/2, |y|/2[).

When x is a vector or matrix, the operator is applied componentwise.

The library syntax is GEN gdivround(GEN x, GEN y) for x \/ y.

The expression x % y evaluates to the modular Euclidean remainder of x and y, which we now define. When x or y is a nonintegral real number, x%y is defined as x - (x\y)*y. Otherwise, if y is an integer, this is the smallest nonnegative integer congruent to x modulo y. (This actually coincides with the previous definition if and only if x is an integer.) If y is a polynomial, this is the polynomial of smallest degree congruent to x modulo y. For instance:

  ? (1/2) % 3
  %1 = 2
  ? 0.5 % 3
  %2 = 0.5000000000000000000000000000
  ? (1/2) % 3.0
  %3 = 1/2

Note that when y is an integer and x a polynomial, y is first promoted to a polynomial of degree 0. When x is a vector or matrix, the operator is applied componentwise.

The library syntax is GEN gmod(GEN x, GEN y) for x % y.

The expression n! is the factorial of the non-negative integer n.

The library syntax is GEN mpfact(long n)

The expression n# is the primorial of the non-negative integer n, that is the product of all prime numbers less than or equal to x.

The library syntax is GEN mpprimorial(long n)

When op is a binary arithmetic operator among +, -, *, %, /, \ or \/, the construct x op =  y is a shortcut for x  =  x op y.

  ? v[1] += 10  \\ increment v[1] by 10
  ? a /= 2 \\ divide a by 2

x++ is a shortcut for x = x + 1 and for x += 1.

x-- is a shortcut for x = x - 1 and for x -= 1.

The expression x^n is powering.

* If the exponent n is an integer, then exact operations are performed using binary (left-shift) powering techniques. By definition, x⁰ is (an empty product interpreted as) an exact 1 in the underlying prime ring:

  ? 0.0 ^ 0
  %1 = 1
  ? (1 + O(2^3)) ^ 0
  %2 = 1
  ? (1 + O(x)) ^ 0
  %3 = 1
  ? Mod(2,4)^0
  %4 = Mod(1,4)
  ? Mod(x,x^2)^0
  %5 = Mod(1, x^2)

If x is a p-adic number, its precision will increase if v_p(n) > 0 and n != 0. Powering a binary quadratic form (type t_QFB) returns a representative of the class, which is reduced if the input was. (In particular, x^1 returns x itself, whether it is reduced or not.)

PARI rewrites the multiplication x * x of two identical objects as x². Here, identical means the operands are reference the same chunk of memory; no equality test is performed. This is no longer true when more than two arguments are involved.

  ? a = 1 + O(2); b = a;
  ? a * a  \\ = a^2, precision increases
  %2 = 1 + O(2^3)
  ? a * b \\ not rewritten as a^2
  %3 = 1 + O(2)
  ? a*a*a \\ not rewritten as a^3
  %4 = 1 + O(2)

* If the exponent is a rational number p/q the behaviour depends on x. If x is a complex number, return exp(n log x) (principal branch), in an exact form if possible:

  ? 4^(1/2)  \\ 4 being a square, this is exact
  %1 = 2
  ? 2^(1/2)  \\ now inexact
  %2 = 1.4142135623730950488016887242096980786
  ? (-1/4)^(1/2) \\ exact again
  %3 = 1/2*I
  ? (-1)^(1/3)
  %4 = 0.500...+ 0.866...*I

Note that even though -1 is an exact cube root of -1, it is not exp(log(-1)/3); the latter is returned.

Otherwise return a solution y of y^q = x^p if it exists; beware that this is defined up to q-th roots of 1 in the base field. Intmods modulo composite numbers are not supported.

  ? Mod(7,19)^(1/2)
  %1 = Mod(11, 19) \\ is any square root
  ? sqrt(Mod(7,19))
  %2 = Mod(8, 19)  \\ is the smallest square root
  ? Mod(1,4)^(1/2)
   ***   at top-level: Mod(1,4)^(1/2)
   ***                         ^ —  — 
   *** _^_: not a prime number in gpow: 4.

* If the exponent is a negative integer or rational number, an inverse must be computed. For noninvertible t_INTMOD x, this will fail and (for n an integer) implicitly exhibit a factor of the modulus:

  ? Mod(4,6)^(-1)
    ***   at top-level: Mod(4,6)^(-1)
    ***                         ^ — --
    *** _^_: impossible inverse modulo: Mod(2, 6).

Here, a factor 2 is obtained directly. In general, take the gcd of the representative and the modulus. This is most useful when performing complicated operations modulo an integer N whose factorization is unknown. Either the computation succeeds and all is well, or a factor d is discovered and the computation may be restarted modulo d or N/d.

For noninvertible t_POLMOD x, the behavior is the same:

  ? Mod(x^2, x^3-x)^(-1)
    ***   at top-level: Mod(x^2,x^3-x)^(-1)
    ***                               ^ — --
    *** _^_: impossible inverse in RgXQ_inv: Mod(x^2, x^3 - x).

Note that the underlying algorihm (subresultant) assumes that the base ring is a domain:

  ? a = Mod(3*y^3+1, 4); b = y^6+y^5+y^4+y^3+y^2+y+1; c = Mod(a,b);
  ? c^(-1)
    ***   at top-level: Mod(a,b)^(-1)
    ***                         ^ — --
    *** _^_: impossible inverse modulo: Mod(2, 4).

In fact c is invertible, but ℤ/4ℤ is not a domain and the algorithm fails. It is possible for the algorithm to succeed in such situations and any returned result will be correct, but chances are that an error will occur first. In this specific case, one should work with 2-adics. In general, one can also try the following approach

  ? inversemod(a, b) =
  { my(m, v = variable(b));
    m = polsylvestermatrix(polrecip(a), polrecip(b));
    m = matinverseimage(m, matid(#m)[,1]);
    Polrev(m[1..poldegree(b)], v);
  }
  ? inversemod(a,b)
  %2 = Mod(2,4)*y^5 + Mod(3,4)*y^3 + Mod(1,4)*y^2 + Mod(3,4)*y + Mod(2,4)

This is not guaranteed to work either since matinverseimage must also invert pivots. See Section se:linear_algebra.

For a t_MAT x, the matrix is expected to be square and invertible, except in the special case x^(-1) which returns a left inverse if one exists (rectangular x with full column rank).

  ? x = Mat([1;2])
  %1 =
  [1]
  
  [2]
  
  ? x^(-1)
  %2 =
  [1 0]

* Finally, if the exponent n is not a rational number, powering is treated as the transcendental function exp(nlog x), although it will be more precise than the latter when n and x are exact:

  ? s = 1/2 + 10^14 * I
  ? localprec(200); z = 2^s  \\ for reference
  ? exponent(2^s - z)
  %3 = -127  \\ perfect
  ? exponent(exp(s * log(2)) - z)
  %4 = -84 \\ not so good

The second computation is less precise because log(2) is first computed to 38 decimal digits, then multiplied by s, which has a huge imaginary part amplifying the error.

In this case, x ⟼ xⁿ is treated as a transcendental function and and in particular acts componentwise on vector or matrices, even square matrices ! (See Section se:trans.) If x is 0 and n is an inexact 0, this will raise an exception:

  ? 4 ^ 1.0
  %1 = 4.0000000000000000000000000000000000000
  ? 0^ 0.0
   ***   at top-level: 0^0.0
   ***                  ^ — -
   *** _^_: domain error in gpow(0,n): n <= 0

The library syntax is GEN gpow(GEN x, GEN n, long prec) for x^n.

Using variants of the extended Euclidean algorithm, returns a rational approximation a/b to x. If B is present, it must be a positive real scalar and it imposes 0 < q ≤ B. If B is omitted, returns a good approximation affordable given the input accuracy: for a t_REAL x, we choose a B whose bit size is half the accuracy of x. If you are looking for true rational numbers, presumably approximated to sufficient accuracy, you should first try that option.

* If x is a t_REAL or a t_FRAC, this function uses continued fractions.

  ? 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 a/b is either the last convergent p_n/q_n of the continued fraction of x with q_n ≤ B or given by (p_n-1 + tp_n) / (q_n-1 + tq_n) for the largest integer t such that the denominator is ≤ B.

* If x is a t_INTMOD modulo N or a t_PADIC of precision N = p^k, this function performs rational modular reconstruction modulo N. The routine then returns the unique rational number a/b in coprime integers |a| < N/2B and b ≤ B which is congruent to x modulo N. Omitting B amounts to choosing it of the order of sqrt{N/2}. If rational reconstruction is not possible (no suitable a/b exists), returns [].

  ? 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, returns [].

The library syntax is GEN bestappr(GEN x, GEN B = NULL).

Using variants of the extended Euclidean algorithm (Padé approximants), returns a rational function approximation a/b to x, whose denominator is limited by B, if present. If B is omitted, returns 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 nonnegative real (impose 0 ≤ degree(b) ≤ B).

* If x is a t_POLMOD modulo N this function performs rational modular reconstruction modulo N. The routine then returns the unique rational function a/b in coprime polynomials, with degree(b) ≤ B and degree(a) minimal, which is congruent to x modulo N. Omitting B amounts to choosing it equal to the floor of degree(N) / 2. If rational reconstruction is not possible (no suitable a/b exists), returns [].

  ? T = Mod(x^3 + x^2 + x + 3, x^4 - 2);
  ? bestapprPade(T)
  %2 = (2*x - 1)/(x - 1)
  ? U = Mod(1 + x + x^2 + x^3 + x^5, x^9);
  ? bestapprPade(U)  \\ internally chooses B = 4
  %3 = []
  ? bestapprPade(U, 5) \\ with B = 5, a solution exists
  %4 = (2*x^4 + x^3 - x - 1)/(-x^5 + x^3 + x^2 - 1)

* If x is a t_SER, we implicitly convert the input to a t_POLMOD modulo N = t^k where k is the series absolute precision.

  ? T = 1 + t + t^2 + t^3 + t^4 + t^5 + t^6 + O(t^7); \\ mod t^7
  ? bestapprPade(T)
  %1 = -1/(t - 1)

* If x is a t_SER and both B and Q are nonnegative, returns a rational function approximation a/b to x, with a of degree at most B and b of degree at most Q, with x-a/b = O(t^B+Q+1+v) if t is the variable, where v is the valuation of x, the empty vector if not possible.

* If x is a t_RFRAC, we implicitly convert the input to a t_POLMOD modulo N = t^k where k = 2B + 1. If B was omitted, we return x:

  ? T = (4*t^2 + 2*t + 3)/(t+1)^10;
  ? bestapprPade(T,1)
  %2 = [] \\ impossible
  ? bestapprPade(T,2)
  %3 = 27/(337*t^2 + 84*t + 9)
  ? bestapprPade(T,3)
  %4 = (4253*t - 3345)/(-39007*t^3 - 28519*t^2 - 8989*t - 1115)

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 GEN bestapprPade0(GEN x, long B, long Q).

GEN bestapprPade(GEN x, long B) as bestapprPade0 when Q is ommited.

T being an integral polynomial and V being a scalar, vector, or matrix with complex coefficients, return a reasonable approximation of V with polmods modulo T. T can also be any number field structure, in which case the minimal polynomial attached to the structure (T.pol) is used. The rootT argument, if present, must be an element of polroots(T) (or T.pol), i.e. a complex root of T fixing an embedding of ℚ[x]/(T) into ℂ.

  ? bestapprnf(sqrt(5), polcyclo(5))
  %1 = Mod(-2*x^3 - 2*x^2 - 1, x^4 + x^3 + x^2 + x + 1)
  ? bestapprnf(sqrt(5), polcyclo(5), exp(4*I*Pi/5))
  %2 = Mod(2*x^3 + 2*x^2 + 1, x^4 + x^3 + x^2 + x + 1)

When the output has huge rational coefficients, try to increase the working realbitprecision: if the answer does not stabilize, consider that the reconstruction failed. Beware that if T is not Galois over ℚ, some embeddings may not allow to reconstruct V:

  ? T = x^3-2; vT = polroots(T); z = 3*2^(1/3)+1;
  ? bestapprnf(z, T, vT[1])
  %2 = Mod(3*x + 1, x^3 - 2)
  ? bestapprnf(z, T, vT[2])
  %3 = 4213714286230872/186454048314072  \\ close to 3*2^(1/3) + 1

The library syntax is GEN bestapprnf(GEN V, GEN T, GEN rootT = NULL, long prec).

Deprecated alias for gcdext

The library syntax is GEN gcdext0(GEN x, GEN y).

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 p-adics p^v (u + O(p^d)) of non-negative valuation v, converted to the obvious intmod p^v u modulo p^v+d and O(p^v) is converted to 0 modulo p^v.

  ? chinese(1 + O(2), 2 + O(3))
  %3 = Mod(5, 6)

Finally, we allow vector and matrix arguments of same dimensions, 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)]
  ? M = mathilbert(3);
  ? chinese(Mod(M,7), Mod(M,11))
  %4 =
  [ Mod(1, 77) Mod(39, 77) Mod(26, 77)]
  
  [Mod(39, 77) Mod(26, 77) Mod(58, 77)]
  
  [Mod(26, 77) Mod(58, 77) Mod(31, 77)]

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 understood as 0 modulo N if the low degree terms are defined mod N:

   ? chinese((x+1)*Mod(1,2), (x^2+2*x+1)*Mod(1,3))
   %3 = Mod(4, 6)*x^2 + Mod(5, 6)*x + Mod(3, 6)

If y is omitted, and x is a vector, chinese is applied recursively to the components of x, yielding a residue belonging to the same class as all components of x.

Finally chinese(x,x) = x regardless of the type of x; this allows vector arguments to contain other data, so long as they are identical in both vectors.

The library syntax is GEN chinese(GEN x, GEN y = NULL). GEN chinese1(GEN x) is also available.

Gives the result of a comparison between arbitrary objects x and y (as -1, 0 or 1). The underlying order relation is transitive, the function returns 0 if and only if x  ===  y. It has no mathematical meaning but satisfies the following properties when comparing entries of the same type:

* two t_INTs compare as usual (i.e. cmp(x,y) < 0 if and only if x < y);

* two t_VECSMALLs of the same length compare lexicographically;

* two t_STRs compare lexicographically.

In case all components are equal up to the smallest length of the operands, the more complex is considered to be larger. More precisely, the longest is the largest; when lengths are equal, we have matrix > vector > scalar. For example:

  ? cmp(1, 2)
  %1 = -1
  ? cmp(2, 1)
  %2 = 1
  ? cmp(1, 1.0)   \\ note that 1 == 1.0, but (1===1.0) is false.
  %3 = -1
  ? cmp(x + Pi, [])
  %4 = -1

This function is mostly useful to handle sorted lists or vectors of arbitrary objects. For instance, if v is a vector, the construction vecsort(v, cmp) is equivalent to Set(v).

The library syntax is int cmp_universal(GEN x, GEN y).

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 (t_INT or t_FRAC), and either 1 (inexact input) or x (exact input) otherwise; the result should be identical to gcd(x, 0).

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 t_VECSMALL is computed assuming the entries are signed integers.

The optional argument D allows to control over which ring we compute and get a more predictable behaviour:

* 1: we only consider the underlying ℚ-structure and the denominator is a (positive) rational number

* a simple variable, say 'x: all entries are considered as rational functions in K(x) where K is the field of fraction of some factorial ring (for instance ℚ or polynomial rings) and the content is an element of K.

  ? f = x + 1/y + 1/2;
  ? content(f) \\ as a t_POL in x
  %2 = 1/(2*y)
  ? content(f, 1) \\ Q-content
  %3 = 1/2
  ? content(f, y) \\ as a rational function in y
  %4 = 1/2
  ? g = x^2*y + y^2*x;
  ? content(g, x)
  %6 = y
  ? content(g, y)
  %7 = x

The library syntax is GEN content0(GEN x, GEN D = NULL).

Returns the row vector whose components are the partial quotients of the continued fraction expansion of x. In other words, a result [a₀,...,a_n] means that x ~ a₀+1/(a₁+...+1/a_n). The output is normalized so that a_n != 1 (unless we also have n = 0).

If x is a real t_QFB q, x is understood as the largest root of q(X,1).

The number of partial quotients n+1 is limited by nmax. If nmax is omitted, the expansion stops at the last significant partial quotient.

If x is a real t_QUAD or a real t_QFB and no limits are specified, the function returns a vector [P,Q] where P is the preperiod and Q is the period.

  ? \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]
  ? w=quadgen(4*114);
  ? contfrac(w,,20)
  %4 = [10,1,2,10,2,1,20,1,2,10,2,1,20,1,2,10,2,1,20,1]
  ? contfrac(w)
  %5 = [[10],[1,2,10,2,1,20]]

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₀)(a₀+b₁/(a₁+...+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 t_REAL).

A direct implementation of the numerical continued fraction contfrac(x,b) described above would be

  \\ "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/(a₁a₂) +...), 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 contfrac(x,, b).

The library syntax is GEN contfrac0(GEN x, GEN b = NULL, long nmax). Also available are GEN gboundcf(GEN x, long nmax), GEN gcf(GEN x) and GEN gcf2(GEN b, GEN x).

When x is a vector or a one-row matrix, x is considered as the list of partial quotients [a₀,a₁,...,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₀+1/(a₁+...+1/a_n). If x is a matrix with two rows [b₀,b₁,...,b_n] and [a₀,a₁,...,a_n], this is then considered as a generalized continued fraction and we have similarly p_n/q_n = (1/b₀)(a₀+b₁/(a₁+...+b_n/a_n)). Note that in this case one usually has b₀ = 1.

If n ≥ 0 is present, returns all convergents from p₀/q₀ up to p_n/q_n. (All convergents if x is too small to compute the n+1 requested convergents.)

  ? a = contfrac(Pi,10)
  %1 = [3, 7, 15, 1, 292, 1, 1, 1, 3]
  ? allpnqn(x) = contfracpnqn(x,#x) \\ all convergents
  ? allpnqn(a)
  %3 =
  [3 22 333 355 103993 104348 208341 312689 1146408]
  
  [1  7 106 113  33102  33215  66317  99532  364913]
  ? contfracpnqn(a) \\ last two convergents
  %4 =
  [1146408 312689]
  
  [ 364913  99532]
  
  ? contfracpnqn(a,3) \\ first three convergents
  %5 =
  [3 22 333 355]
  
  [1  7 106 113]

The library syntax is GEN contfracpnqn(GEN x, long n). also available is GEN pnqn(GEN x) for n = -1.

Creates a column vector with two components, the first being the Euclidean quotient (x \y), the second the Euclidean remainder (x - (x\y)*y), of the division of x by y. This avoids the need to do two divisions if one needs both the quotient and the remainder. If v is present, and x, y are multivariate polynomials, divide with respect to the variable v.

Beware that divrem(x,y)[2] is in general not the same as x % y; no GP operator corresponds to it:

  ? divrem(1/2, 3)[2]
  %1 = 1/2
  ? (1/2) % 3
  %2 = 2
  ? divrem(Mod(2,9), 3)[2]
   ***   at top-level: divrem(Mod(2,9),3)[2
   ***                 ^ —  —  —  —  —  — --
   ***   forbidden division t_INTMOD \ t_INT.
  ? Mod(2,9) % 6
  %3 = Mod(2,3)

The library syntax is GEN divrem(GEN x, GEN y, long v = -1) where v is a variable number. Also available is GEN gdiventres(GEN x, GEN y) when v is not needed.

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 gcdext function. x and y can have rather quite general types, for instance both rational numbers. If y is omitted and x is a vector, returns the gcd of all components of x, i.e. this is equivalent to content(x).

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 gcd(x, y[i]), resp. gcd(x, y[,i]). Else if x is a vector/matrix the result has the same type as x and an analogous definition. Note that for these types, gcd is not commutative.

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" nontrivial, is to use polresultant: a value close to 0 means that a small deformation of the inputs has nontrivial gcd. You may also use gcdext, which does try to compute an approximate gcd d and provides u, v to check whether u x + v y is close to d.

The library syntax is GEN ggcd0(GEN x, GEN y = NULL). Also available are GEN ggcd(GEN x, GEN y), if y is not NULL, and GEN content(GEN x), if y = NULL.

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 gcd(x,y) which is defined to be a scalar in this case:

  ? 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 GEN gcdext0(GEN x, GEN y).

Let inputs x and y be both integers, or both polynomials in the same variable. Return a vector [M, [a,b]~], where M is an invertible 2 x 2 matrix such that M*[x,y]~ = [a,b]~, where b is small. More precisely,

* polynomial case: det M has degree 0 and we have deg a ≥ ceil{max(deg x,deg y))/2} > deg b.

* integer case: det M = ± 1 and we have a ≥ ceil{sqrt{max(|x|,|y|)}} > b. Assuming x and y are nonnegative, then M^-1 has nonnegative coefficients, and det M is equal to the sign of both main diagonal terms M[1,1] and M[2,2].

The library syntax is GEN ghalfgcd(GEN x, GEN y).

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 nonnegative 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 lcm(x, y[i]), resp. lcm(x, y[,i]). Else if x is a vector/matrix the result has the same type as x and an analogous definition. Note that for these types, lcm is not commutative.

Note that lcm(v) is quite different from

  l = v[1]; for (i = 1, #v, l = lcm(l, v[i]))

Indeed, lcm(v) is a scalar, but l may not be (if one of the v[i] is a vector/matrix). The computation uses a divide-conquer tree and should be much more efficient, especially when using the GMP multiprecision kernel (and more subquadratic algorithms become available):

  ? 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 GEN glcm0(GEN x, GEN y = NULL).

Gives the result of a lexicographic comparison between x and y (as -1, 0 or 1). This is to be interpreted in quite a wide sense: it is admissible to compare objects of different types (scalars, vectors, matrices), provided the scalars can be compared, as well as vectors/matrices of different lengths; finally, when comparing two scalars, a complex number a + I*b is interpreted as a vector [a,b] and a real number a as [a,0]. The comparison is recursive.

In case all components are equal up to the smallest length of the operands, the more complex is considered to be larger. More precisely, the longest is the largest; when lengths are equal, we have matrix > vector > scalar. For example:

  ? lex([1,3], [1,2,5])
  %1 = 1
  ? lex([1,3], [1,3,-1])
  %2 = -1
  ? lex([1], [[1]])
  %3 = -1
  ? lex([1], [1]~)
  %4 = 0
  ? lex(2 - I, 1)
  %5 = 1
  ? lex(2 - I, 2)
  %6 = -1

The library syntax is int lexcmp(GEN x, GEN y).

Creates the maximum of x and y when they can be compared.

The library syntax is GEN gmax(GEN x, GEN y).

Creates the minimum of x and y when they can be compared.

The library syntax is GEN gmin(GEN x, GEN y).

Shifts x componentwise left by n bits if n ≥ 0 and right by |n| bits if n < 0. May be abbreviated as x << n or x >> (-n). A left shift by n corresponds to multiplication by 2ⁿ. A right shift of an integer x by |n| corresponds to a Euclidean division of x by 2^|n| with a remainder of the same sign as x, hence is not the same (in general) as x \ 2ⁿ.

The library syntax is GEN gshift(GEN x, long n).

Multiplies x by 2ⁿ. The difference with shift is that when n < 0, ordinary division takes place, hence for example if x is an integer the result may be a fraction, while for shifts Euclidean division takes place when n < 0 hence if x is an integer the result is still an integer.

The library syntax is GEN gmul2n(GEN x, long n).

sign (0, 1 or -1) of x, which must be of type integer, real or fraction; t_QUAD with positive discriminants and t_INFINITY are also supported.

The library syntax is int gsigne(GEN x).

If x is a list, vector or matrix, returns the largest entry of x, otherwise returns a copy of x. Error if x is empty. Here, largest refers to the ordinary real ordering ( <= ).

If v is given, set it to the index of a largest entry (indirect maximum), when x is a vector or list. If x is a matrix, set v to coordinates [i,j] such that x[i,j] is a largest entry. This argument v is ignored for other types. When the vector has equal largest entries, the first occurence is chosen; in a matrix, the smallest j is chosen first, then the smallest i. vector or matrix.

  ? vecmax([10, 20, -30, 40])
  %1 = 40
  ? vecmax([10, 20, -30, 40], &v); v
  %2 = 4
  ? vecmax([10, 20; -30, 40], &v); v
  %3 = [2, 2]

The library syntax is GEN vecmax0(GEN x, GEN *v = NULL). When v is not needed, the function GEN vecmax(GEN x) is also available.

If x is a list, vector or matrix, returns the smallest entry of x, otherwise returns a copy of x. Error if x is empty. Here, smallest refers to the ordinary real ordering ( <= ).

If v is given, set it to the index of a smallest entry (indirect minimum), when x is a vector or list. If x is a matrix, set v to coordinates [i,j] such that x[i,j] is a smallest entry. This argument v is ignored for other types. When a vector has equal smallest entries, the first occurence is chosen; in a matrix, the smallest j is chosen first, then the smallest i.

  ? vecmin([10, 20, -30, 40])
  %1 = -30
  ? vecmin([10, 20, -30, 40], &v); v
  %2 = 3
  ? vecmin([10, 20; -30, 40], &v); v
  %3 = [2, 1]
  ? vecmin([1,0;0,0], &v); v
  %3 = [2, 1]

The library syntax is GEN vecmin0(GEN x, GEN *v = NULL). When v is not needed, the function GEN vecmin(GEN x) is also available.