### Vectors, matrices, linear algebra and sets

#### algdep(x,k,{flag = 0})

x being real/complex, or p-adic, finds a polynomial of degree at most k with integer coefficients having x as approximate root. Note that the polynomial which is obtained is not necessarily the "correct" one. In fact it is not even guaranteed to be irreducible. One can check the closeness either by a polynomial evaluation (use subst), or by computing the roots of the polynomial given by algdep (use polroots).

Internally, lindep([1,x,...,x^k], flag) is used. If lindep is not able to find a relation and returns a lower bound for the sup norm of the smallest relation, algdep returns that bound instead. A non-zero value of flag may improve on the default behavior if the input number is known to a huge accuracy, and you suspect the last bits are incorrect (this truncates the number, throwing away the least significant bits), but default values are usually sufficient:

```\\\\\\\\\ LLL
? \p200
? algdep(2^(1/6)+3^(1/5), 30);      \\ wrong in 0.8s
? algdep(2^(1/6)+3^(1/5), 30, 100); \\ wrong in 0.4s
? algdep(2^(1/6)+3^(1/5), 30, 170); \\ right in 0.8s
? algdep(2^(1/6)+3^(1/5), 30, 200); \\ wrong in 1.0s
? \p250
? algdep(2^(1/6)+3^(1/5), 30);      \\ right in 1.0s
? algdep(2^(1/6)+3^(1/5), 30, 200); \\ right in 1.0s
? \p500
? algdep(2^(1/6)+3^(1/5), 30);      \\ right in 2.9s
? \p1000
? algdep(2^(1/6)+3^(1/5), 30);      \\ right in 10.6s
\\\\\\\\\ PSLQ
? \p200
? algdep(2^(1/6)+3^(1/5), 30, -3);  \\ failure in 15s
? \p250
? algdep(2^(1/6)+3^(1/5), 30, -3);  \\ right in 20s
? \p500
? algdep(2^(1/6)+3^(1/5), 30, -3);  \\ right in 52s
? \p1000
? algdep(2^(1/6)+3^(1/5), 30, -3);  \\ right in 164s
```

The changes in defaultprecision only affect the quality of the initial approximation to 2^{1/6} + 3^{1/5}, algdep itself uses exact operations (the size of its operands depend on the accuracy of the input of course: more accurate input means slower operations).

Proceeding by increments of 5 digits of accuracy, algdep with default flag produces its first correct result at 205 digits, and from then on a steady stream of correct results. Interestingly enough, our PSLQ also reliably succeeds from 205 digits on (and is 15 times slower at that accuracy).

The above example is the test case studied in a 2000 paper by Borwein and Lisonek: Applications of integer relation algorithms, Discrete Math., 217, p.65--82. The paper concludes in the superiority of the PSLQ algorithm, which either shows that PARI's implementation of PSLQ is lacking, or that its LLL is extremely good. The version of PARI tested there was 1.39, which succeeded reliably from precision 265 on, in about 200 as much time as the current version.

The library syntax is GEN algdep0(GEN x, long k, long flag). Also available is GEN algdep(GEN x, long k) (flag = 0).

#### charpoly(A,{v = x},{flag = 3})

characteristic polynomial of A with respect to the variable v, i.e.determinant of v*I-A if A is a square matrix. If A is not a square matrix, it returns the characteristic polynomial of the map "multiplication by A" if A is a scalar, in particular a polmod. E.g. charpoly(I) = x^2+1.

The value of flag is only significant for matrices. Let n be the dimension of A.

If flag = 0, same method (Le Verrier's) as for computing the adjoint matrix, i.e.using the traces of the powers of A. Assumes that n! is invertible; uses O(n^4) scalar operations.

If flag = 1, uses Lagrange interpolation which is usually the slowest method. Assumes that n! is invertible; uses O(n^4) scalar operations.

If flag = 2, uses the Hessenberg form. Assumes that the base ring is a field. Uses O(n^3) scalar operations, but suffers from coefficient explosion unless the base field is finite or R.

If flag = 3, uses Berkowitz's division free algorithm, valid over any ring (commutative, with unit). Uses O(n^4) scalar operations.

If flag = 4, x must be integral. Uses a modular algorithm.

In practice one should use the default (Berkowitz) unless the base ring is Z (use flag = 4) or a field where coefficient explosion does not occur, e.g.a finite field or the reals (use flag = 2).

The library syntax is GEN charpoly0(GEN A, long v = -1, long flag), where v is a variable number. Also available are GEN caract(GEN A, long v) (flag = 1), GEN carhess(GEN A, long v) (flag = 2), \fun{GEN}{carberkowitz}{GEN A, long v} (flag = 3) and GEN caradj(GEN A, long v, GEN *pt). In this last case, if pt is not NULL, *pt receives the address of the adjoint matrix of A (see matadjoint), so both can be obtained at once.

#### concat(x,{y})

concatenation of x and y. If x or y is not a vector or matrix, it is considered as a one-dimensional vector. All types are allowed for x and y, but the sizes must be compatible. Note that matrices are concatenated horizontally, i.e.the number of rows stays the same. Using transpositions, it is easy to concatenate them vertically.

To concatenate vectors sideways (i.e.to obtain a two-row or two-column matrix), use Mat instead (see the example there). Concatenating a row vector to a matrix having the same number of columns will add the row to the matrix (top row if the vector is x, i.e.comes first, and bottom row otherwise).

The empty matrix [;] is considered to have a number of rows compatible with any operation, in particular concatenation. (Note that this is definitely not the case for empty vectors [] or []~.)

If y is omitted, x has to be a row vector or a list, in which case its elements are concatenated, from left to right, using the above rules.

```? concat([1,2], [3,4])
%1 = [1, 2, 3, 4]
? a = [[1,2]~, [3,4]~]; concat(a)
%2 =
[1 3]

[2 4]

? concat([1,2; 3,4], [5,6]~)
%3 =
[1 2 5]

[3 4 6]
? concat([%, [7,8]~, [1,2,3,4]])
%5 =
[1 2 5 7]

[3 4 6 8]

[1 2 3 4]
```

The library syntax is GEN concat(GEN x, GEN y = NULL). GEN concat1(GEN x) is a shortcut for concat(x,NULL).

#### lindep(x,{flag = 0})

x being a vector with p-adic or real/complex coefficients, finds a small integral linear combination among these coefficients.

If x is p-adic, flag is meaningless and the algorithm LLL-reduces a suitable (dual) lattice.

Otherwise, the value of flag determines the algorithm used; in the current version of PARI, we suggest to use non-negative values, since it is by far the fastest and most robust implementation. See the detailed example in Section [Label: se:algdep] ( algdep).

If flag >= 0, uses a floating point (variable precision) LLL algorithm. This is in general much faster than the other variants. If flag = 0 the accuracy is chosen internally using a crude heuristic. If flag > 0 the computation is done with an accuracy of flag decimal digits. To get meaningful results in the latter case, the parameter flag should be smaller than the number of correct decimal digits in the input.

If flag = -1, uses a variant of the LLL algorithm due to Hastad, Lagarias and Schnorr (STACS 1986). If the precision is too low, the routine may enter an infinite loop. Faster than the alternatives if it converges, especially when the accuracy is much larger than what is really necessary; usually diverges, though.

If flag = -2, x is allowed to be (and in any case interpreted as) a matrix. Returns a non trivial element of the kernel of x, or 0 if x has trivial kernel. The element is defined over the field of coefficients of x, and is in general not integral.

If flag = -3, uses the PSLQ algorithm. This may return a real number B, indicating that the input accuracy was exhausted and that no relation exist whose sup norm is less than B.

If flag = -4, uses an experimental 2-level PSLQ, which does not work at all. Don't use it!

The library syntax is GEN lindep0(GEN x, long flag). Also available are GEN lindep(GEN x) (flag = 0) GEN lindep2(GEN x, long bit) (flag >= 0, bypasses the check for p-adic inputs) and GEN deplin(GEN x) (flag = -2).

#### listcreate()

creates an empty list. This routine used to have a mandatory argument, which is now ignored (for backward compatibility). In fact, this function has become redundant and obsolete; it will disappear in future versions of PARI: just use List()

#### listinsert(L,x,n)

inserts the object x at position n in L (which must be of type t_LIST). This has complexity O(#L - n + 1): all the remaining elements of list (from position n+1 onwards) are shifted to the right.

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

#### listkill(L)

obsolete, retained for backward compatibility. Just use L = List() instead of listkill(L). In most cases, you won't even need that, e.g. local variables are automatically cleared when a user function returns.

The library syntax is void listkill(GEN L).

#### listpop(list,{n})

removes the n-th element of the list list (which must be of type t_LIST). If n is omitted, or greater than the list current length, removes the last element. This runs in time O(#L - n + 1).

The library syntax is void listpop(GEN list, long n).

#### listput(list,x,{n})

sets the n-th element of the list list (which must be of type t_LIST) equal to x. If n is omitted, or greater than the list length, appends x. You may put an element into an occupied cell (not changing the list length), but it is easier to use the standard list[n] = x construct. This runs in time O(#L) in the worst case (when the list must be reallocated), but in time O(1) on average: any number of successive listputs run in time O(#L), where #L denotes the list final length.

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

#### listsort(L,{flag = 0})

sorts the t_LIST list in place. If flag is non-zero, suppresses all repeated coefficients. This is faster than the vecsort command since no copy has to be made. No value returned.

The library syntax is void listsort(GEN L, long flag).

#### matadjoint(x,{flag = 0})

adjoint matrix of x, i.e.the matrix y of cofactors of x, satisfying x*y = det(x)*\Id. x must be a (non-necessarily invertible) square matrix of dimension n. If flag is 0 or omitted, use a fast algorithm which assumes that n! is invertible. If flag is 1, use a slower division-free algorithm.

```? a = [1,2,3;3,4,5;6,7,8] * Mod(1,2);
***   at top-level: matadjoint([1,2,3;3,
***                 ^--------------------
*** matadjoint: impossible inverse modulo: Mod(0, 2).
? matadjoint(a, 1)  \\ use safe algorithm
%2 =
[Mod(1, 2) Mod(1, 2) Mod(0, 2)]

[Mod(0, 2) Mod(0, 2) Mod(0, 2)]

[Mod(1, 2) Mod(1, 2) Mod(0, 2)]
```

Both algorithms use O(n^4) operations in the base ring.

The library syntax is GEN matadjoint0(GEN x, long flag). Also available are GEN adj(GEN x) (flag = 0) and GEN adjsafe(GEN x) (flag = 1).

#### matcompanion(x)

the left companion matrix to the polynomial x.

The library syntax is GEN matcompanion(GEN x).

#### matdet(x,{flag = 0})

determinant of x. x must be a square matrix.

If flag = 0, uses Gauss-Bareiss.

If flag = 1, uses classical Gaussian elimination, which is better when the entries of the matrix are reals or integers for example, but usually much worse for more complicated entries like multivariate polynomials.

The library syntax is GEN det0(GEN x, long flag). Also available are GEN det(GEN x) (flag = 0) and GEN det2(GEN x) (flag = 1).

#### matdetint(x)

x being an m x n matrix with integer coefficients, this function computes a non-zero multiple of the determinant of the lattice generated by the columns of x if it has maximal rank m, and returns zero otherwise, using the Gauss-Bareiss algorithm. When x is square, the exact determinant is obtained.

This function is useful in conjunction with mathnfmod, which needs to know such a multiple. If the rank is maximal and the matrix non-square, you can obtain the exact determinant using

```  matdet( mathnfmod(x, matdetint(x)) )
```

Note that as soon as one of the dimensions gets large (m or n is larger than 20, say), it will often be much faster to use mathnf(x, 1) or mathnf(x, 4) directly.

The library syntax is GEN detint(GEN x).

#### matdiagonal(x)

x being a vector, creates the diagonal matrix whose diagonal entries are those of x.

The library syntax is GEN diagonal(GEN x).

#### mateigen(x)

gives the eigenvectors of x as columns of a matrix.

The library syntax is GEN eigen(GEN x, long prec).

#### matfrobenius(M,{flag},{v = x})

returns the Frobenius form of the square matrix M. If flag = 1, returns only the elementary divisors as a vector of polynomials in the variable v. If flag = 2, returns a two-components vector [F,B] where F is the Frobenius form and B is the basis change so that M = B^{-1}FB.

The library syntax is GEN matfrobenius(GEN M, long flag, long v = -1), where v is a variable number.

#### mathess(x)

returns a matrix similar to the square matrix x, which is in upper Hessenberg form (zero entries below the first subdiagonal).

The library syntax is GEN hess(GEN x).

#### mathilbert(n)

x being a long, creates the Hilbert matrixof order x, i.e.the matrix whose coefficient (i,j) is 1/ (i+j-1).

The library syntax is GEN mathilbert(long n).

#### mathnf(x,{flag = 0})

if x is a (not necessarily square) matrix with integer entries, finds the upper triangular Hermite normal form of x. If the rank of x is equal to its number of rows, the result is a square matrix. In general, the columns of the result form a basis of the lattice spanned by the columns of x.

If flag = 0, uses the naive algorithm. This is in general fastest but may require too much memory as the dimension gets large (bigger than 100, say), in which case you may try mathnfmod(x, matdetint(x)) when x has maximal rank, and mathnf(x, 4) otherwise.

If flag = 1, outputs a two-component row vector [H,U], where H is the Hermite normal form of x defined as above, and U is the unimodular transformation matrix such that xU = [0|H]. When the kernel is large, U has in general huge coefficients. In the worst case, the running time is exponential with respect to the dimension, but the routine behaves well in small dimension (less than 50 or 100, say).

If flag = 3, uses Batut's algorithm and output [H,U,P], such that H and U are as before and P is a permutation of the rows such that P applied to xU gives H. This is in general slower than flag = 1 but the matrix U is smaller; it may still be large.

If flag = 4, as in case 1 above, but uses a variant of LLL reduction along the way. The matrix U is in general close to optimal (in terms of smallest L_2 norm), but the reduction is in general slow, although provably polynomial-time.

The library syntax is GEN mathnf0(GEN x, long flag). Also available are GEN hnf(GEN x) (flag = 0) and GEN hnfall(GEN x) (flag = 1). To reduce huge (say 400 x 400 and more) relation matrices (sparse with small entries), you can use the pair hnfspec / hnfadd. Since this is quite technical and the calling interface may change, they are not documented yet. Look at the code in basemath/alglin1.c.

#### mathnfmod(x,d)

if x is a (not necessarily square) matrix of maximal rank with integer entries, and d is a multiple of the (non-zero) determinant of the lattice spanned by the columns of x, finds the upper triangular Hermite normal form of x.

If the rank of x is equal to its number of rows, the result is a square matrix. In general, the columns of the result form a basis of the lattice spanned by the columns of x. Even when d is known, this is in general slower than mathnf but uses much less memory.

The library syntax is GEN hnfmod(GEN x, GEN d).

#### mathnfmodid(x,d)

outputs the (upper triangular) Hermite normal form of x concatenated with d times the identity matrix. Assumes that x has integer entries.

The library syntax is GEN hnfmodid(GEN x, GEN d).

#### matid(n)

creates the n x n identity matrix.

The library syntax is GEN matid(long n).

#### matimage(x,{flag = 0})

gives a basis for the image of the matrix x as columns of a matrix. A priori the matrix can have entries of any type. If flag = 0, use standard Gauss pivot. If flag = 1, use matsupplement (much slower: keep the default flag!).

The library syntax is GEN matimage0(GEN x, long flag). Also available is GEN image(GEN x) (flag = 0).

#### matimagecompl(x)

gives the vector of the column indices which are not extracted by the function matimage. Hence the number of components of matimagecompl(x) plus the number of columns of matimage(x) is equal to the number of columns of the matrix x.

The library syntax is GEN imagecompl(GEN x).

#### matindexrank(x)

x being a matrix of rank r, returns a vector with two t_VECSMALL components y and z of length r giving a list of rows and columns respectively (starting from 1) such that the extracted matrix obtained from these two vectors using vecextract(x,y,z) is invertible.

The library syntax is GEN indexrank(GEN x).

#### matintersect(x,y)

x and y being two matrices with the same number of rows each of whose columns are independent, finds a basis of the Q-vector space equal to the intersection of the spaces spanned by the columns of x and y respectively. The faster function idealintersect can be used to intersect fractional ideals (projective Z_K modules of rank 1); the slower but much more general function nfhnf can be used to intersect general Z_K-modules.

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

#### matinverseimage(x,y)

given a matrix x and a column vector or matrix y, returns a preimage z of y by x if one exists (i.e such that x z = y), an empty vector or matrix otherwise. The complete inverse image is z + {Ker} x, where a basis of the kernel of x may be obtained by matker.

```? M = [1,2;2,4];
? matinverseimage(M, [1,2]~)
%2 = [1, 0]~
? matinverseimage(M, [3,4]~)
%3 = []~    \\ no solution
? matinverseimage(M, [1,3,6;2,6,12])
%4 =
[1 3 6]

[0 0 0]
? matinverseimage(M, [1,2;3,4])
%5 = [;]    \\ no solution
? K = matker(M)
%6 =
[-2]

[1]
```

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

#### matisdiagonal(x)

returns true (1) if x is a diagonal matrix, false (0) if not.

The library syntax is GEN isdiagonal(GEN x).

#### matker(x,{flag = 0})

gives a basis for the kernel of the matrix x as columns of a matrix. The matrix can have entries of any type, provided they are compatible with the generic arithmetic operations (+, x and /).

If x is known to have integral entries, set flag = 1.

The library syntax is GEN matker0(GEN x, long flag). Also available are GEN ker(GEN x) (flag = 0), GEN keri(GEN x) (flag = 1).

#### matkerint(x,{flag = 0})

gives an LLL-reduced Z-basis for the lattice equal to the kernel of the matrix x as columns of the matrix x with integer entries (rational entries are not permitted).

If flag = 0, uses an integer LLL algorithm.

If flag = 1, uses matrixqz(x,-2). Many orders of magnitude slower than the default: never use this.

The library syntax is GEN matkerint0(GEN x, long flag). See also GEN kerint(GEN x) (flag = 0), which is a trivial wrapper around

```ZM_lll(ZM_lll(x, 0.99, LLL_KER), 0.99, LLL_INPLACE);
```
Remove the outermost ZM_lll if LLL-reduction is not desired (saves time).

#### matmuldiagonal(x,d)

product of the matrix x by the diagonal matrix whose diagonal entries are those of the vector d. Equivalent to, but much faster than x* matdiagonal(d).

The library syntax is GEN matmuldiagonal(GEN x, GEN d).

#### matmultodiagonal(x,y)

product of the matrices x and y assuming that the result is a diagonal matrix. Much faster than x*y in that case. The result is undefined if x*y is not diagonal.

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

#### matpascal(n,{q})

creates as a matrix the lower triangular Pascal triangle of order x+1 (i.e.with binomial coefficients up to x). If q is given, compute the q-Pascal triangle (i.e.using q-binomial coefficients).

The library syntax is GEN matqpascal(long n, GEN q = NULL). Also available is GEN matpascal(GEN x).

#### matrank(x)

rank of the matrix x.

The library syntax is long rank(GEN x).

#### matrix(m,n,{X},{Y},{expr = 0})

creation of the m x n matrix whose coefficients are given by the expression expr. There are two formal parameters in expr, the first one (X) corresponding to the rows, the second (Y) to the columns, and X goes from 1 to m, Y goes from 1 to n. If one of the last 3 parameters is omitted, fill the matrix with zeroes.

#### matrixqz(A,{p = 0})

A being an m x n matrix in M_{m,n}(Q), let {Im}_Q A (resp.{Im}_Z A) the Q-vector space (resp.the Z-module) spanned by the columns of A. This function has varying behavior depending on the sign of p:

If p >= 0, A is assumed to have maximal rank n <= m. The function returns a matrix B belongs to M_{m,n}(Z), with {Im}_Q B = {Im}_Q A, such that the GCD of all its n x n minors is coprime to p; in particular, if p = 0 (default), this GCD is 1.

```? minors(x) = vector(#x[,1], i, matdet( vecextract(x, Str("^",i), "..") ));
? A = [3,1/7; 5,3/7; 7,5/7]; minors(A)
%1 = [4/7, 8/7, 4/7]   \\ determinants of all 2x2 minors
? B = matrixqz(A)
%2 =
[3 1]

[5 2]

[7 3]
? minors(%)
%3 = [1, 2, 1]   \\ B integral with coprime minors
```

If p = -1, returns the HNF basis of the lattice Z^n cap {Im}_Z A.

If p = -2, returns the HNF basis of the lattice Z^n cap {Im}_Q A.

```? matrixqz(A,-1)
%4 =
[8 5]

[4 3]

[0 1]

? matrixqz(A,-2)
%5 =
[2 -1]

[1 0]

[0 1]
```

The library syntax is GEN matrixqz0(GEN A, GEN p = NULL).

#### matsize(x)

x being a vector or matrix, returns a row vector with two components, the first being the number of rows (1 for a row vector), the second the number of columns (1 for a column vector).

The library syntax is GEN matsize(GEN x).

#### matsnf(X,{flag = 0})

if X is a (singular or non-singular) matrix outputs the vector of elementary divisors of X, i.e.the diagonal of the Smith normal form of X, normalized so that d_n | d_{n-1} | ... | d_1.

The binary digits of flag mean:

1 (complete output): if set, outputs [U,V,D], where U and V are two unimodular matrices such that UXV is the diagonal matrix D. Otherwise output only the diagonal of D. If X is not a square matrix, then D will be a square diagonal matrix padded with zeros on the left or the top.

2 (generic input): if set, allows polynomial entries, in which case the input matrix must be square. Otherwise, assume that X has integer coefficients with arbitrary shape.

4 (cleanup): if set, cleans up the output. This means that elementary divisors equal to 1 will be deleted, i.e.outputs a shortened vector D' instead of D. If complete output was required, returns [U',V',D'] so that U'XV' = D' holds. If this flag is set, X is allowed to be of the form `vector of elementary divisors' or [U,V,D] as would normally be output with the cleanup flag unset.

The library syntax is GEN matsnf0(GEN X, long flag).

#### matsolve(M,B)

M being an invertible matrix and B a column vector, finds the solution X of MX = B, using Gaussian elimination. This has the same effect as, but is a bit faster, than M^{-1}*B.

The library syntax is GEN gauss(GEN M, GEN B).

#### matsolvemod(M,D,B,{flag = 0})

M being any integral matrix, D a column vector of non-negative integer moduli, and B an integral column vector, gives a small integer solution to the system of congruences sum_i m_{i,j}x_j = b_i (mod d_i) if one exists, otherwise returns zero. Shorthand notation: B (resp.D) can be given as a single integer, in which case all the b_i (resp.d_i) above are taken to be equal to B (resp.D).

```? M = [1,2;3,4];
? matsolvemod(M, [3,4]~, [1,2]~)
%2 = [-2, 0]~
? matsolvemod(M, 3, 1) \\ M X = [1,1]~ over F_3
%3 = [-1, 1]~
? matsolvemod(M, [3,0]~, [1,2]~) \\ x + 2y = 1 (mod 3), 3x + 4y = 2 (in Z)
%4 = [6, -4]~
```

If flag = 1, all solutions are returned in the form of a two-component row vector [x,u], where x is a small integer solution to the system of congruences and u is a matrix whose columns give a basis of the homogeneous system (so that all solutions can be obtained by adding x to any linear combination of columns of u). If no solution exists, returns zero.

The library syntax is GEN matsolvemod0(GEN M, GEN D, GEN B, long flag). Also available are GEN gaussmodulo(GEN M, GEN D, GEN B) (flag = 0) and GEN gaussmodulo2(GEN M, GEN D, GEN B) (flag = 1).

#### matsupplement(x)

assuming that the columns of the matrix x are linearly independent (if they are not, an error message is issued), finds a square invertible matrix whose first columns are the columns of x, i.e.supplement the columns of x to a basis of the whole space.

The library syntax is GEN suppl(GEN x).

#### mattranspose(x)

transpose of x (also x~). This has an effect only on vectors and matrices.

The library syntax is GEN gtrans(GEN x).

#### minpoly(A,{v = x})

minimal polynomial of A with respect to the variable v., i.e. the monic polynomial P of minimal degree (in the variable v) such that P(A) = 0.

The library syntax is GEN minpoly(GEN A, long v = -1), where v is a variable number.

#### qfgaussred(q)

decomposition into squares of the quadratic form represented by the symmetric matrix q. The result is a matrix whose diagonal entries are the coefficients of the squares, and the off-diagonal entries on each line represent the bilinear forms. More precisely, if (a_{ij}) denotes the output, one has q(x) = sum_i a_{ii} (x_i + sum_{j != i} a_{ij} x_j)^2

```? qfgaussred([0,1;1,0])
%1 =
[1/2 1]

[-1 -1/2]
```
This means that 2xy = (1/2)(x+y)^2 - (1/2)(x-y)^2.

The library syntax is GEN qfgaussred(GEN q). GEN qfgaussred_positive(GEN q) assumes that q is positive definite and is a little faster; returns NULL if a vector with negative norm occurs (non positive matrix or too many rounding errors).

#### qfjacobi(x)

x being a real symmetric matrix, this gives a vector having two components: the first one is the vector of (real) eigenvalues of x, sorted in increasing order, the second is the corresponding orthogonal matrix of eigenvectors of x. The method used is Jacobi's method for symmetric matrices.

The library syntax is GEN jacobi(GEN x, long prec).

#### qflll(x,{flag = 0})

LLL algorithm applied to the columns of the matrix x. The columns of x may be linearly dependent. The result is a unimodular transformation matrix T such that x .T is an LLL-reduced basis of the lattice generated by the column vectors of x. Note that if x is not of maximal rank T will not be square. The LLL parameters are (0.51,0.99), meaning that the Gram-Schmidt coefficients for the final basis satisfy mu_{i,j} <= |0.51|, and the Lov\'{a}sz's constant is 0.99.

If flag = 0 (default), assume that x has either exact (integral or rational) or real floating point entries. The matrix is rescaled, converted to integers and the behavior is then as in flag = 1.

If flag = 1, assume that x is integral. Computations involving Gram-Schmidt vectors are approximate, with precision varying as needed (Lehmer's trick, as generalized by Schnorr). Adapted from Nguyen and Stehlé's algorithm and Stehlé's code ( fplll-1.3).

If flag = 2, x should be an integer matrix whose columns are linearly independent. Returns a partially reduced basis for x, using an unpublished algorithm by Peter Montgomery: a basis is said to be partially reduced if |v_i ± v_j| >= |v_i| for any two distinct basis vectors v_i, v_j.

This is faster than flag = 1, esp. when one row is huge compared to the other rows (knapsack-style), and should quickly produce relatively short vectors. The resulting basis is not LLL-reduced in general. If LLL reduction is eventually desired, avoid this partial reduction: applying LLL to the partially reduced matrix is significantly slower than starting from a knapsack-type lattice.

If flag = 4, as flag = 1, returning a vector [K, T] of matrices: the columns of K represent a basis of the integer kernel of x (not LLL-reduced in general) and T is the transformation matrix such that x.T is an LLL-reduced Z-basis of the image of the matrix x.

If flag = 5, case as case 4, but x may have polynomial coefficients.

If flag = 8, same as case 0, but x may have polynomial coefficients.

The library syntax is GEN qflll0(GEN x, long flag). Also available are GEN lll(GEN x) (flag = 0), GEN lllint(GEN x) (flag = 1), and GEN lllkerim(GEN x) (flag = 4).

#### qflllgram(G,{flag = 0})

same as qflll, except that the matrix G = x~ * x is the Gram matrix of some lattice vectors x, and not the coordinates of the vectors themselves. In particular, G must now be a square symmetric real matrix, corresponding to a positive quadratic form (not necessarily definite: x needs not have maximal rank). The result is a unimodular transformation matrix T such that x.T is an LLL-reduced basis of the lattice generated by the column vectors of x. See qflll for further details about the LLL implementation.

If flag = 0 (default), assume that G has either exact (integral or rational) or real floating point entries. The matrix is rescaled, converted to integers and the behavior is then as in flag = 1.

If flag = 1, assume that G is integral. Computations involving Gram-Schmidt vectors are approximate, with precision varying as needed (Lehmer's trick, as generalized by Schnorr). Adapted from Nguyen and Stehlé's algorithm and Stehlé's code ( fplll-1.3).

flag = 4: G has integer entries, gives the kernel and reduced image of x.

flag = 5: same as 4, but G may have polynomial coefficients.

The library syntax is GEN qflllgram0(GEN G, long flag). Also available are GEN lllgram(GEN G) (flag = 0), GEN lllgramint(GEN G) (flag = 1), and GEN lllgramkerim(GEN G) (flag = 4).

#### qfminim(x,{b},{m},{flag = 0})

x being a square and symmetric matrix representing a positive definite quadratic form, this function deals with the vectors of x whose norm is less than or equal to b, enumerated using the Fincke-Pohst algorithm, storing at most m vectors (no limit if m is omitted). The function searches for the minimal non-zero vectors if b is omitted. The behavior is undefined if x is not positive definite (a "precision too low" error is most likely, although more precise error messages are possible). The precise behavior depends on flag.

If flag = 0 (default), seeks at most 2m vectors. The result is a three-component vector, the first component being the number of vectors found, the second being the maximum norm found, and the last vector is a matrix whose columns are the vectors found, only one being given for each pair ± v (at most m such pairs, unless m was omitted). The vectors are returned in no particular order.

If flag = 1, ignores m and returns the first vector whose norm is less than b. In this variant, an explicit b must be provided.

In these two cases, x must have integral entries. The implementation uses low precision floating point computations for maximal speed, which gives incorrect result when x has large entries. (The condition is checked in the code and the routine raises an error if large rounding errors occur.) A more robust, but much slower, implementation is chosen if the following flag is used:

If flag = 2, x can have non integral real entries. In this case, if b is omitted, the "minimal" vectors only have approximately the same norm. If b is omitted, m is an upper bound for the number of vectors that will be stored and returned, but all minimal vectors are nevertheless enumerated. If m is omitted, all vectors found are stored and returned; note that this may be a huge vector!

```? x = matid(2);
? qfminim(x)  \\ 4 minimal vectors of norm 1: ±[0,1], ±[1,0]
%2 = [4, 1, [0, 1; 1, 0]]
? { x =
[4, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 1, 0,-1, 0, 0, 0,-2;
2, 4,-2,-2, 0,-2, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1, 0, 1,-1,-1;
0,-2, 4, 0,-2, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 1, 0, 0, 1,-1,-1, 0, 0;
0,-2, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 1,-1, 0, 1,-1, 1, 0;
0, 0,-2, 0, 4, 0, 0, 0, 1,-1, 0, 0, 1, 0, 0, 0,-2, 0, 0,-1, 1, 1, 0, 0;
-2, -2,0, 0, 0, 4,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,-1, 1, 1;
0, 0, 0, 0, 0,-2, 4,-2, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 1,-1, 0;
0, 0, 0, 0, 0, 0,-2, 4, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1,-1,-1, 0, 1, 0;
0, 0, 0, 0, 1,-1, 0, 0, 4, 0,-2, 0, 1, 1, 0,-1, 0, 1, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0,-1, 0, 0, 0, 0, 4, 0, 0, 1, 1,-1, 1, 0, 0, 0, 1, 0, 0, 1, 0;
0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 4,-2, 0,-1, 0, 0, 0,-1, 0,-1, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 4,-1, 1, 0, 0,-1, 1, 0, 1, 1, 1,-1, 0;
1, 0,-1, 1, 1, 0, 0,-1, 1, 1, 0,-1, 4, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1,-1;
-1,-1, 1,-1, 0, 0, 1, 0, 1, 1,-1, 1, 0, 4, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1;
0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 1, 1, 0, 4, 0, 0, 0, 0, 1, 1, 0, 0;
0, 0, 1, 0,-2, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 4, 1, 1, 1, 0, 0, 1, 1;
1, 0, 0, 1, 0, 0,-1, 0, 1, 0,-1, 1, 1, 0, 0, 0, 1, 4, 0, 1, 1, 0, 1, 0;
0, 0, 0,-1, 0, 1, 0,-1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 4, 0, 1, 1, 0, 1;
-1, -1,1, 0,-1, 1, 0,-1, 0, 1,-1, 1, 0, 1, 0, 0, 1, 1, 0, 4, 0, 0, 1, 1;
0, 0,-1, 1, 1, 0, 0,-1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 4, 1, 0, 1;
0, 1,-1,-1, 1,-1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 4, 0, 1;
0,-1, 0, 1, 0, 1,-1, 1, 0, 1, 0,-1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 4, 1;
-2,-1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 4]; }
? qfminim(x,,0)  \\ the Leech lattice has 196560 minimal vectors of norm 4
time = 648 ms.
%4 = [196560, 4, [;]]
? qfminim(x,,0,2); \\ safe algorithm. Slower and unnecessary here.
time = 18,161 ms.
%5 = [196560, 4.000061035156250000, [;]]
```

In the last example, we store 0 vectors to limit memory use. All minimal vectors are nevertheless enumerated. Provided parisize is about 50MB, qfminim(x) succeeds in 2.5 seconds.

The library syntax is GEN qfminim0(GEN x, GEN b = NULL, GEN m = NULL, long flag, long prec). Also available are GEN minim(GEN x, GEN b = NULL, GEN m = NULL) (flag = 0), GEN minim2(GEN x, GEN b = NULL, GEN m = NULL) (flag = 1).

#### qfperfection(G)

G being a square and symmetric matrix with integer entries representing a positive definite quadratic form, outputs the perfection rank of the form. That is, gives the rank of the family of the s symmetric matrices v_iv_i^t, where s is half the number of minimal vectors and the v_i (1 <= i <= s) are the minimal vectors.

Since this requires computing the minimal vectors, the computations can become very lengthy as the dimension of x grows.

The library syntax is GEN perf(GEN G).

#### qfrep(q,B,{flag = 0})

q being a square and symmetric matrix with integer entries representing a positive definite quadratic form, outputs the vector whose i-th entry, 1 <= i <= B is half the number of vectors v such that q(v) = i. This routine uses a naive algorithm based on qfminim, and will fail if any entry becomes larger than 2^{31}.

The binary digits of flag mean:

* 1: count vectors of even norm from 1 to 2B.

* 2: return a t_VECSMALL instead of a t_VEC

The library syntax is GEN qfrep0(GEN q, GEN B, long flag).

#### qfsign(x)

returns [p,m] the signature of the quadratic form represented by the symmetric matrix x. Namely, p (resp.m) is the number of positive (resp.negative) eigenvalues of x.The result is computed using Gaussian reduction.

The library syntax is GEN qfsign(GEN x).

#### setintersect(x,y)

intersection of the two sets x and y (see setisset). The function also works if both x and y are vectors of strictly increasing entries, according to < ); in that case we return a vector of strictly increasing entries, not a set. Otherwise, the result is undefined.

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

#### setisset(x)

returns true (1) if x is a set, false (0) if not. In PARI, a set is a row vector whose entries are strictly increasing t_STRs. To convert any object into a set (this is most useful for vectors, of course), use the function Set.

```? a = [3, 1, 1, 2];
? setisset(a)
%2 = 0
? Set(a)
%3 = ["1", "2", "3"]
```

The library syntax is long setisset(GEN x).

#### setminus(x,y)

difference of the two sets x and y (see setisset), i.e.set of elements of x which do not belong to y. The function also works if both x and y are vectors of strictly increasing entries, according to < ); in that case we return a vector of strictly increasing entries, not a set. Otherwise, the result is undefined.

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

#### setsearch(S,x,{flag = 0})

searches if x belongs to the set S (see setisset). A set is a vector of t_STR, but this function works also if S is a arbitrary sorted vector or list (see listsort): if x is not a t_STR, we first replace it by Str(x) unless the first element of S is also not a t_STR.

If x belongs to the set and flag is zero or omitted, returns the index j such that S[j] = x, otherwise returns 0. If flag is non-zero returns the index j where x should be inserted, and 0 if it already belongs to S (this is meant to be used in conjunction with listinsert, see below).

```? T = [2,3,5,7]; S = Set(T);
? setsearch(S, 2)      \\ search in a true set, t_INT 2 converted to string
%2 = 1
? setsearch(S, Str(2)) \\ search in a true set, no need for conversion
%3 = 1
? setsearch(T, 2)      \\ search in a sorted vector, no need for conversion
%4 = 1
? setsearch(T, Str(2)) \\ search in a sorted vector, t_STR "2" not found
%5 = 0
%6 = 0
? setsearch(S, 4, 1)   \\ should have been inserted at index 3
%7 = 3
```

The library syntax is long setsearch(GEN S, GEN x, long flag).

#### setunion(x,y)

union of the two sets x and y (see setisset). The function also works if both x and y are vectors of strictly increasing entries, according to < ); in that case we return a vector of strictly increasing entries, not a set. Otherwise, the result is undefined.

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

#### trace(x)

this applies to quite general x. If x is not a matrix, it is equal to the sum of x and its conjugate, except for polmods where it is the trace as an algebraic number.

For x a square matrix, it is the ordinary trace. If x is a non-square matrix (but not a vector), an error occurs.

The library syntax is GEN gtrace(GEN x).

#### vecextract(x,y,{z})

extraction of components of the vector or matrix x according to y. In case x is a matrix, its components are as usual the columns of x. The parameter y is a component specifier, which is either an integer, a string describing a range, or a vector.

If y is an integer, it is considered as a mask: the binary bits of y are read from right to left, but correspond to taking the components from left to right. For example, if y = 13 = (1101)_2 then the components 1,3 and 4 are extracted.

If y is a vector, which must have integer entries, these entries correspond to the component numbers to be extracted, in the order specified.

If y is a string, it can be

* a single (non-zero) index giving a component number (a negative index means we start counting from the end).

* a range of the form "a..b", where a and b are indexes as above. Any of a and b can be omitted; in this case, we take as default values a = 1 and b = -1, i.e.the first and last components respectively. We then extract all components in the interval [a,b], in reverse order if b < a.

In addition, if the first character in the string is ^, the complement of the given set of indices is taken.

If z is not omitted, x must be a matrix. y is then the line specifier, and z the column specifier, where the component specifier is as explained above.

```? v = [a, b, c, d, e];
? vecextract(v, 5)          \\ mask
%1 = [a, c]
? vecextract(v, [4, 2, 1])  \\ component list
%2 = [d, b, a]
? vecextract(v, "2..4")     \\ interval
%3 = [b, c, d]
? vecextract(v, "-1..-3")   \\ interval + reverse order
%4 = [e, d, c]
? vecextract(v, "^2")       \\ complement
%5 = [a, c, d, e]
? vecextract(matid(3), "2..", "..")
%6 =
[0 1 0]

[0 0 1]
```

The library syntax is GEN extract0(GEN x, GEN y, GEN z = NULL).

#### vecsort(x,{cmp},{flag = 0})

sorts the vector x in ascending order, using a mergesort method. x must be a list, vector or matrix (seen as a vector of its columns). Note that mergesort is stable, hence the initial ordering of "equal" entries (with respect to the sorting criterion) is not changed.

If cmp is omitted, we use the standard comparison function < , thereby restricting the possible types for the elements of x (integers, fractions or reals). If cmp is present, it is understood as a comparison function and we sort according to it. The following possibilities exist:

* an integer k: sort according to the value of the k-th subcomponents of the components ofx.

* a vector: sort lexicographically according to the components listed in the vector. For example, if cmp = [2,1,3], sort with respect to the second component, and when these are equal, with respect to the first, and when these are equal, with respect to the third.

* a comparison function ( t_CLOSURE), with two arguments x and y, and returning an integer which is < 0, > 0 or = 0 if x < y, x > y or x = y respectively. The sign function is very useful in this context:

```? vecsort([3,0,2; 1,0,2],lex) \\ sort columns according to lex order
%1 =
[2 3]

[2 1]
? vecsort(v, (x,y)->sign(y-x))            \\ reverse sort
? vecsort(v, (x,y)->sign(abs(x)-abs(y)))  \\ sort by increasing absolute value
? cmp(x,y) = my(dx = poldisc(x), dy = poldisc(y)); sign(abs(dx) - abs(dy))
? vecsort([x^2+1, x^3-2, x^4+5*x+1], cmp)
```

The last example used the named cmp instead of an anonymous function, and sorts polynomials with respect to the absolute value of their discriminant. A more efficient approach would use precomputations to ensure a given discriminant is computed only once:

```? DISC = vector(#v, i, abs(poldisc(v[i])));
? perm = vecsort(vector(#v,i,i), (x,y)->sign(DISC[x]-DISC[y]))
? vecextract(v, perm)
```
Similar ideas apply whenever we sort according to the values of a function which is expensive to compute.

The binary digits of flag mean:

* 1: indirect sorting of the vector x, i.e.if x is an n-component vector, returns a permutation of [1,2,...,n] which applied to the components of x sorts x in increasing order. For example, vecextract(x, vecsort(x,,1)) is equivalent to vecsort(x).

* 2: sorts x by ascending lexicographic order (as per the lex comparison function).

* 4: use descending instead of ascending order.

* 8: remove "duplicate" entries with respect to the sorting function (keep the first occurring entry). For example:

```  ? vecsort([Pi,Mod(1,2),z], (x,y)->0, 8)   \\ make everything compare equal
%1 = [3.141592653589793238462643383]
? vecsort([[2,3],[0,1],[0,3]], 2, 8)
%2 = [[0, 1], [2, 3]]
```

The library syntax is GEN vecsort0(GEN x, GEN cmp = NULL, long flag).

#### vector(n,{X},{expr = 0})

creates a row vector (type t_VEC) with n components whose components are the expression expr evaluated at the integer points between 1 and n. If one of the last two arguments is omitted, fill the vector with zeroes.

Avoid modifying X within expr; if you do, the formal variable still runs from 1 to n. In particular, vector(n,i,expr) is not equivalent to

```v = vector(n)
for (i = 1, n, v[i] = expr)
```

as the following example shows:

```n = 3
v = vector(n); vector(n, i, i++)            ----> [2, 3, 4]
v = vector(n); for (i = 1, n, v[i] = i++)   ----> [2, 0, 4]
```

#### vectorsmall(n,{X},{expr = 0})

creates a row vector of small integers (type t_VECSMALL) with n components whose components are the expression expr evaluated at the integer points between 1 and n. If one of the last two arguments is omitted, fill the vector with zeroes.

#### vectorv(n,{X},{expr = 0})

as vector, but returns a column vector (type t_COL).