singular-data 4.0.3+ds-1 / usr / share / singular / LIB / arr.lib

//////////////////////////////////////////////////////////////////////////////
version="version arr.lib 4.0.2.0 Feb 2015 "; // $Id: da6cd859e2bb646dfdaa5235b8fc68962c7378b8 $
category="Miscellaneous";

/*

** TOPICS ** (Ctrl+f to search)
#1 NEWSTRUCTS & OVERLOADS
#2 CONSTRUCTORS
#3 ACCESS & ASSIGNEMENT
#4 DELETION
#5 COMPERATORS
#6 TYPE CASTING
#7 INHERITED FUNCTIONS
#8 PRINTING
#9 MANIPULATING VARIABLES
#10 CENTER COMPUTATIONS
#11 GEOMETRIC CONSTRUCTIONS
#11 EXAMPLES OF ARRANGEMENTS
#13 ORLIK SOLOMON AND POINCARE POLYNOMIAL
#14 FREENESS
#15 MULTI-ARRANGEMENTS
#16 COMBINATORICS

*/

info="
LIBRARY:    arr.lib  a library of algorithms for arrangements of hyperplanes

AUTHORS:    Randolf Scholz (rscholz@rhrk.uni-kl.de),
            Patrick Serwene (serwene@mathematik.uni-kl.de),
            Lukas Kuehne (lf.kuehne@gmail.com)

OVERLOADS:

// OPERATORS
=           arrAdd         assignment
+           arrAdd         union of two arrs
[           arrGet         access to a single/multiple hyperplane(s)
-           arrMinus       deletes given hyperplanes from the arr
<=          arrLEQ         comparison
>=          arrGEQ         comparison
==          arrEQ          comparison
!=          arrNEQ         comparison
<           arrLNEQ        comparison
>           arrGNEQ        comparison

// TYPECASTING
matrix      arr2mat         coeff matrix
poly        arr2poly        defining polynomial

// OTHER
variables   arrVariables    ideal generated by the variables the arr depends on
nvars       arrNvars        number of variables the arr depends on
delete      arrDelete       deletes hyperplanes by indices
print       arrPrint        prints the arr on the screen

// IDEAL INHERITED FUNCTIONS
homog       arrHomog        checks if arrangement is homogeneous
simplify    arrSimplify     simplifies arrangement
size        arrSize         number of planes
subst       arrSubst        substitute variables

// MULTI-ARRANGEMENTS
=       multarrAdd          assignement of multarr
+       multarrAdd          union of multarr
poly    multarr2poly        defining polynomial
size    multarrSize         number of hyperplanes with mult.
print   multarrPrint        displays multiarr
delete  multarrDelete       deletes hyperplane

PROCEDURES:
arrSet(arr A, int k, poly p)    replaces the k-th Hyperplane with poly p

type2arr(#)                     converts general input to 'arr' using arrAdd.
mat2arr( matrix M)              affine  arrangement from coeff matrix
mat2carr(matrix M)              central arrangement from coeff matrix
arrPrintMatrix(arr A)           readable output as a coeff matrix
varMat(intvec v)                matrix of the corresponding ring_variables
varNum(def u)                   number of given variable (enh. version of varNum in dmod.lib)
arrSwapVar(arr A, i, j)         swaps two variables in the arrangement
arrLastVar(arr A)               ring_variable of largest index used in arrangement
arrCenter(arr A)                computes center of an arrangement
arrCentral(arr A)               checks if arrangement is central
arrCentered(arr A)              checks if arrangement is centered
arrCentralize(arr A)            makes centered arrangement central
arrCoordChange(A, T, #)         performs coordinate change
arrCoordNormalize(A, v)         performs projection onto coordinate hyperplane
arrCone(arr A, var)             coned arrangement
arrDecone(arr A, int k)         deconed arrangement
arrLocalize(arr A, intvec v)    localization of an arrangement onto a flat
arrRestrict(arr A, intvec v)    restricted arrangement onto a flat
arrIsEssential(arr A)           checks if arrangement is essential
arrEssentialize(arr A)          essentialized arragnement
arrBoolean(int v)               boolean arrangement
arrBraid(int v)                 braid arrangement
arrTypeB(int v)                 type B arrangement
arrTypeD(int v)                 type D arrangement
arrRandom(d,m,n)                random (affine) arrangement
arrRandomCentral(d,m,n)         random central arrangement
arrEdelmanReiner()              Edelman-Reiner arrangement
arrOrlikSolomon(arr A)          Orlik-Solomon algebra of the arrangement
arrDer(A)                       module of derivation
arrIsFree(A)                    checks if arrangement is free
arrExponents(A)                 exponents of a (free) arrangement
arr2multarr(arr A, intvec v)    converts normal arrangement to multiarrangement
multarr2arr(multarr A)          converts multiarrangement to normal arrangement
multarrRestrict(arr A, v)       restriction of A (as arr) to a flat with multiplicities
multarrMultRestrict(A, int k)   restriction of A (as multarr) to a hyperplane with multiplicities
arrFlats(arr A)                     intersection lattice
arrLattice(arr A)               computes the intersection lattice / poset
moebius(arrposet P)             computes moebius values
arrCharPoly(arr A)                 characteristic polynomial
arrPoincare(arr A)              poincare polynomial of the arrangement
arrChambers(arr A)              number of chambers of the arrangement
arrBoundedChambers(arr A)       number of bounded chambers of the arrangement
printMoebius(arr A)             displays the moebius values of all the flats in the poset
";


//============================================================================//
//-------------------------- #1 NEWSTRUCTS & OVERLOADS -----------------------//
//============================================================================//

// initialization of the library
static proc mod_init()
{
    // NEWSTRUCTS
    newstruct("arr","ideal l");
    newstruct("flat", "intvec REL, int moebius, intvec parents, int flag");
    newstruct("arrposet","arr A, list r");
    newstruct("arrflats","arr A, list r");
    newstruct("multarr","ideal l, intvec m");    // intvec: multiplicities of hyperplanes

    // OPERATORS
    system("install","arr","="          ,arrAdd         ,1);    // assignment
    system("install","arr","+"          ,arrAdd         ,2);    // union of arrs
    system("install","arr","["          ,arrGet         ,2);    // access
    system("install","arr","-"          ,arrMinus       ,2);    // delete plane
    system("install","arr","<="         ,arrLEQ         ,2);    // comparison
    system("install","arr",">="         ,arrGEQ         ,2);    // comparison
    system("install","arr","=="         ,arrEQ          ,2);    // comparison
    system("install","arr","!="         ,arrNEQ         ,2);    // comparison
    system("install","arr","<"          ,arrLNEQ        ,2);    // comparison
    system("install","arr",">"          ,arrGNEQ        ,2);    // comparison

    // TYPECASTING
    system("install","arr","matrix"     ,arr2mat        ,1);    // coeff matrix
    system("install","arr","poly"       ,arr2poly       ,1);    // defining polynomial
    system("install","arr","list"       ,arr2list       ,1);    // list of defining polynomials
    system("install","arr","ideal"      ,arr2ideal      ,1);    // list of defining polynomials

    // OTHER
    system("install","arr","variables"  ,arrVariables   ,1);
    system("install","arr","nvars"      ,arrNvars       ,1);
    system("install","arr","delete"     ,arrDelete      ,2);
    system("install","arr","print"      ,arrPrint       ,1);

    // IDEAL INHERITED FUNCTIONS
    system("install","arr","homog"      ,arrHomog       ,1);    // checks if homogeneous
    system("install","arr","homog"      ,arrHomog       ,2);    // checks if homogeneous
    system("install","arr","simplify"   ,arrSimplify    ,2);    // simplifies arrangement
    system("install","arr","size"       ,arrSize        ,1);    // number of planes
    system("install","arr","subst"      ,arrSubst       ,4);    // substitute variables

    // MULTI-ARRANGEMENTS
    system("install","multarr","="      ,multarrAdd     ,1);    // assignement of multarr
    system("install","multarr","+"      ,multarrAdd     ,2);    // union of multarr
    system("install","multarr","poly"   ,multarr2poly   ,1);    // defining polynomial
    system("install","multarr","size"   ,multarrSize    ,1);    // number of hyperplanes with mult.
    system("install","multarr","print"  ,multarrPrint   ,1);    // displays multiarr
    system("install","multarr","delete" ,multarrDelete  ,2);    // deletes hyperplane

    // COMBINATORICS
    system("install","arr","rank"           ,arrRank            ,1);
    system("install","arrflats","print"     ,arrPrintFlats      ,1);
    system("install","arrposet","print"     ,arrPrintPoset      ,1);

    // NEEDED LIBRARIES
    LIB "general.lib";
    LIB "monomialideal.lib";

}


//============================================================================//
//--------------------------- #2 CONSTRUCTORS --------------------------------//
//============================================================================//


// general method for creating arrangements
static proc arrAdd
"USAGE:     A = #; A +#;  # list containing arr/ideal/list/matrix/poly
RETURN:     [arr] Arrangement constructed by input parameters.
REMARKS:    The algorithm splits up the list # and uses the appropiate procedure
            to handle the input.
NOTE:       arrAdd simplifies certain inputs, for example A = (x,y,2x); gives the same arrangement
            as A = (x,y);
SEE ALSO:   arrAdd, type2arr
KEYWORDS:   arrangement; equal; constructor; operator
EXAMPLE:    example arrAdd; shows an example"

{
    arr A;
    for(int k=1; k<=size(#); k++){
        while(1){   //simulates switch, which singular doesn't offer
            if(typeof(#[k]) == "arr"   ){A = arrAddArr   (A, #[k]);break;}
            if(typeof(#[k]) == "poly"  ){A = arrAddPoly  (A, #[k]);break;}
            if(typeof(#[k]) == "ideal" ){A = arrAddIdeal (A, #[k]);break;}
            if(typeof(#[k]) == "matrix"){A = arrAddMatrix(A, #[k]);break;}
            if(typeof(#[k]) == "intmat"){A = arrAddMatrix(A, #[k]);break;}
            if(typeof(#[k]) == "list"  ){A = arrAdd      (   #[k]);break;}

            ERROR("bad input type");
        }
    }

    return (A);
}

example
{
"EXAMPLE: Creating a few arrangements"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A= ideal(x,y,z);        A;
arr B = A + ideal(x+1, x-1);     B;
arr C = list(A, x+1, x-1);  C;
arr D  = x2 - y2;           D;
}

// union of two arrangements
static proc arrAddArr(arr A, arr   B){
    return (arrAddIdeal(A, B.l));
}

// adds a single poly to the arrangement
// if the poly is linear, it is just added, otherwise Singular factorizes
static proc arrAddPoly(arr A, poly p){
    if(deg(p) == 0){
        ERROR("Given poly is not linear or Singluar is not able to factorize it");}
    else{
        if(deg(p) == 1){
            A.l = A.l + p;
            return (A);
        }
        else{
        ideal I = factorize(p,1);
            if(size(I) == 1){ERROR("Given poly is not a hyperplane");}
            else{ return (arrAdd(A,I)); }
        }
    }
    return(A);
}

// adds defining polys to the arrangement
static proc arrAddIdeal(arr A, ideal I){
    for(int k=1; k<=size(I); k++){
        A = arrAddPoly(A,I[k]);
    }

    return (A);
}

// adds defining polys to the arrangement
static proc arrAddMatrix(arr A, matrix M){
        return ( arrAddArr(A,mat2arr(M)) );
}


//============================================================================//
//--------------------------- #3 ACCESS & ASSIGNEMENT ------------------------//
//============================================================================//


// access to hyperplanes, overloads [] operator
static proc arrGet(arr A, intvec v)
"USAGE:     A[v]; v int/intvec
RETURN:     [poly] if v is [int] The defining poly of the the v-th hyperplane+
            [arr] if v is [intvec] The corresponing subarrangement
SEE ALSO:   arrGet, arrSet
KEYWORDS:   arrangement; get; operator
EXAMPLE:    example arrGet; shows an example"

{
    if(size(v) == 1){ return (A.l[v[1]]); }     //returns poly if v is integer
    ideal I = A.l;
    A = ideal(I[v]);
    return (A);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
A[2];
intvec v = 1,3;
A[v];
}

// replaces the k-th plane with poly p
proc arrSet(arr A, int k, poly p)
"USAGE:     arrSet(A, k, p);    arr A, int k, poly p;
RETURN:     [arr] Arrangement where the k-th hyperplane is replaced by p.
NOTE:       p must be linear
KEYWORDS:   arrangement; hyperplane; assign; set
EXAMPLE:    example arrSet; shows an example"

{
    if(deg(p) != 1){ERROR("Given poly is not a hyperplane");}
    else{                   // looks akward but needs to be done this way
        ideal I = A.l;
        I[k] = p;
        A = I;
    }
    return (A);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
arrSet(A,1,x+1);
}


//============================================================================//
//------------------------------- #4 DELETION --------------------------------//
//============================================================================//


// deletes all hyperplanes of the given indices
static proc arrDelete(arr A, intvec v)
"USAGE:     delete(A, v); v integer/intvec
RETURN:     [arr] Arrangement A without the hyperplanes given by v;
NOTE:       for deleting hyperplanes via polynomials, use arrMinus instead
SEE ALSO:   arrDelete, arrMinus
KEYWORDS:   arrangement; delete
EXAMPLE:    example arrDelete; shows an example"

{
    int i = 0; int k;
    int n = size(A);
    intvec u = 1..n;

// puts 0 in u for every element that needs to be deleted
// is done this way to deal with the case that the user gives the same index multiple times.
    for(k=1; k<=size(v); k++){
        if(u[v[k]] != 0){ u[v[k]] = 0; i++; }   // i = #elts that need to be deleted
    }

    if( i == n){return (emptyArr); }

// create intvec of the remaining indices
    v = 1..(n-i); i=1;
    for(k=1; k<=n; k++){
        if(u[k] != 0)   { v[i] = u[k]; i++; }
    }

    arr A' = arrGet(A,v);
    return (A');
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
delete(A,2);
intvec v = (1,3);
delete(A,v);
}


// deletes hyperplanes, overloads - operator
static proc arrMinus
"USAGE:     A - #; # list containing arr/ideal/list/matrix/poly
RETURN:     [arr] arrangement A without the hyperplanes of the arrangement defined by #.
REMARKS:    algorithm creates an arrangement B from # using arrAdd and
            then deletes hyperplanes which occur in both A and B.
NOTE:       The alorithm does not simplify by scalars, i.e. some hyperplanes might
            not be deleted. See example.
SEE ALSO:   arrDelete, arrMinus
KEYWORDS:   arrangement; delete; operator
EXAMPLE:    example arrMinus; shows an example"

{
    if(typeof(#[1]) != "arr"){ERROR("First input must be arr!");}
    arr A = #[1];
    arr B = #[2..size(#)];          // collects hyperplanes to be deleted
    list L;
    int k, l;
    for(k=1; k<=size(A); k++){      // create list of hyperplanes to be deleted
        for(l=1; l<=size(B); l++){
            if(A[k] == B[l]){L = insert(L,k);}
        }
    }

    l = size(L);                    // transforms list to intvec
    if(l != 0){
        intvec v = 1..l;
        for(k=1; k<=l; k++){v[k] = L[k];}
        A = delete(A,v);
    }

    return (A);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
A - ideal(x,y);
A - poly(2x);
}


//============================================================================//
//------------------------------- #5 COMPERATORS -----------------------------//
//============================================================================//


// returns logical value of 'A<=B'
static proc arrLEQ(arr A, arr B)
"USAGE:     A<=B; A,B arr
RETURN:     [0,1] true if A is a subarrangement of B
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; subarrangement; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    arr C = A - B;
    if(C[1] == 0){ return (1); }
    return (0);
}

example{example arrEQ;}


// returns logical value of 'A>=B'
static proc arrGEQ(arr A, arr B)
"USAGE:     A>=B; A,B arr
RETURN:     [0,1] true if B is a subarrangement of A
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; subarrangement; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    arr C = B - A;
    if(C[1] == 0){ return (1); }
    return (0);
}

example{example arrEQ;}


// returns logical value of 'A==B'
static proc arrEQ(arr A, arr B)
"USAGE:     A==B; A,B arr
RETURN:     [0,1] true if A equals B
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; equal; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    return ((A<=B) & (A>=B));
}

example
{
"EXAMPLE: Relationships between a few arrangements."; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
arr B = ideal(x,y);
A<=B;
A>=B;
A==B;
A!=B;
A<B;
A>B;
}


// returns logical value of 'A!=B'
static proc arrNEQ(arr A, arr B)
"USAGE:     A!=B; A,B arr
RETURN:     [0,1] true if A is not equal to B
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; equal; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    return (!(A==B));
}

example{example arrEQ;}


// returns logical value of 'A<B'
static proc arrLNEQ(arr A, arr B)
"USAGE:     A<B; A,B arr
RETURN:     [0,1] true if A is a proper subarrangement of B
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; subarrangement; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    return ((A<=B) & (A!=B));
}

example{example arrEQ;}


// returns logical value of 'A>B'
static proc arrGNEQ(arr A, arr B)
"USAGE:     A>B; A,B arr
RETURN:     [0,1] true if B is a proper subarrangement of A
NOTE:       algorithm is based on arrMinus and does not simplify by scalars, hence some
            technically equal hyperplanes might not be detected. See example.
SEE ALSO:   arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ
KEYWORDS:   comparison; subarrangement; operator
EXAMPLE:    example arrLEQ; shows an example"

{
    return ((A>=B) & (A!=B));
}

example{example arrEQ;}


//============================================================================//
//------------------------------ #6 TYPE CASTING -----------------------------//
//============================================================================//


// TYPE => ARRANGEMENT
proc type2arr(list #)
"USAGE:     type2arr(#); # def
RETURN:     [arr] Arrangement defined by the input
NOTE:       The procedure tries to cast the input to [arr] using arrAdd
KEYWORDS:   typecasting; arrangement
EXAMPLE:    example tye2arr; shows an example"

{
    return (arrAdd(#));
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
ideal I = x,y,z;
typeof(type2arr(I));
}


// ARRANGEMENT => POLY
static proc arr2poly(arr A)
"USAGE:     poly(A);  arr A
RETURN:     [poly] The defining polynomial which is the product of polynomials occuring in arr
NOTE:       The procedure will automatically simplify the polynomial by scalar multiplication.
SEE ALSO:   arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr
KEYWORDS:   typecasting; defining polynomial
EXAMPLE:    example arr2poly; shows an example"

{
    poly f = simplify(product(A.l),1);
    if(f == 0){ return (1); }
    return (f);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x+1, 2x-2, y);
arr2poly(A);
}


// ARRANGEMENT => LIST
static proc arr2list(arr A)
"USAGE:     arr2list(A); A arr
RETURN:     [list] containing all generators of A
SEE ALSO:   arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr
KEYWORDS:   typecasting; list
EXAMPLE:    example arr2list; shows an example"

{
    int n = size(A);
    list L;
    for(int k=1; k<=n; k++){L[k] = A[k];}

    return (L);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A= ideal(x,y,z);
arr2list(A);
}


// ARRANGEMENT => IDEAL
static proc arr2ideal(arr A)
"USAGE:     arr2ideal(A); A arr
RETURN:     [ideal], the internal ideal of A
NOTE:       arr2ideal(A); is the same as A.l - which may become private
SEE ALSO:   arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr
KEYWORDS:   typecasting; ideal
EXAMPLE:    example arr2ideal; shows an example"

{
    return (A.l);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
arr2ideal(A);
}


// ARRANGEMENT => MATRIX
static proc arr2mat(arr A, list #)
"USAGE:     matrix(A); A arr
            matrix(A, 'c') for central arrangements (shorter matrix)
RETURN:     [matrix] M = [T|b] representing the arrangement
NOTE:       If the arrangement is central or one is not interested in the const
            terms one can use "matrix(A, 'c')" instead to get the same matrix without
            the last column.
SEE ALSO:   arr2mat, mat2arr, mat2carr
KEYWORDS:   typecasting; matrix; coefficient
EXAMPLE:    example arr2mat;"

{
    matrix M = jacob(A.l);
    if(size(#) == 0){return ( concat(M, transpose(jet(A.l,0))) );}
    if(#[1] == 'c'){return (M);}
    ERROR("Bad optional input parameter!")
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
arr D = ideal(x+1, y-2, z, x+y+4);
print(arr2mat(A));
print(arr2mat(D));
}


// MATRIX => ARRANGEMENT
proc mat2arr(matrix M)
"USAGE:     mat2arr(M); matrix (M|b)
RETURN:     [arr] interprets the rows of the matrix as the defining polynomial equations
            of the arrangement where the last column will be considered as the constant
            terms, i.e. if M is an m*(n+1) matrix we have
            H_i = Ker( M_i1*x_1 +...+ M_in*x_n + M_i(n+1) ) for i=1...m and
            A = {H_1,...,H_m} the resulting arrangement.
SEE ALSO:   mat2carr
KEYWORDS:   typecasting; matrix; coefficient
EXAMPLE:    example mat2arr;"

{
    if(ncols(M) > nvars(basering)+1)
    {
        ERROR("Matrix too big! Please add variables to basering.");
    }

    int n = ncols(M)-1;
    matrix X[n+1][1];
    X[1..n,1] = varMat(1..n);
    X[n+1 ,1] = 1;
    arr A = ideal(M*X);

    return(A);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
matrix M[4][4] = 1,0,1,1,1,1,0,2,0,1,1,3,2,1,1,4;
print(M);
mat2arr(M);
}


// MATRIX => ARRANGEMENT (central)
proc mat2carr(matrix M)
"USAGE:     mat2carr(M); matrix M
RETURN:     [arr] interprets the rows of the matrix as the defining polynomial equations
            of the arrangement. I.e. if M is an m*n matrix we have
            H_i = Ker( M_i1*x_1 +...+ M_in*x_n) for i=1...m and
            A = {H_1,...,H_m} the resulting arrangement.
SEE ALSO:   mat2arr
KEYWORDS:   typecasting; matrix; coefficient; central
EXAMPLE:    example mat2carr;"

{
    if( ncols(M) > nvars(basering) ){
        ERROR("Error! not enough variables in the basering.");
    }

    arr A = ideal(M*varMat(1..ncols(M)));

    return(A);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
matrix M[4][3] = 1,0,1,1,1,0,0,1,1,2,1,1;
print(M);
mat2carr(M);
}


//============================================================================//
//---------------------------- #7 IDEAL FUNCTIONS ----------------------------//
//============================================================================//


// checks if A is central, homogenizes
static proc arrHomog(arr A, list #)
"USAGE:     homog(A); arr A
            homog(A, p); arr A, ring_variable p
RETURN:     [0,1] homog(A) is the same as arrCentral(A)
            [arr] homog(A,p) homogenizes A with respect to p
NOTE:       homog(A,p) is not the same as arrCone(A,p) as it does not add the additional hyperplane
SEE ALSO:   arrHomog, arrCentral, arrCone
KEYWORDS:   central; homogenize
EXAMPLE:    example arrHomog; shows an example"

{
    if(size(#) == 0){
        return (homog(A.l));
    }

    if(size(#) == 1){
        A = homog(A.l, #[1]);
        return (A);
    }

    ERROR("Too many innput arguments!");

}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y);
arr B = ideal(x,y,x+y+1);
homog(A);
homog(B);
homog(B,z);
homog(_);
}



// simplified arrangement
static proc arrSimplify(arr A, list #)
"USAGE:     arrSimplify(A);
            simplify(A, n); arr A, int n
RETURN:     [arr] simplified arrangement.
NOTE:       arrSimplify(A) is the same as simplify(A, 1), simplify with higher ints
            is not needed
SEE ALSO:   arrSimplify
KEYWORDS:   simplify
EXAMPLE:    example arrSimplify; shows an example"

{
    if(size(#) == 0){
        A = simplify(A.l, 1);
        return (A);
    }

    if(size(#) == 1){
        A = simplify(A.l, #[1]);
        return (A);
    }

    ERROR("Too many input arguments!");
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(2x,2y-1,2z-2);
arrSimplify(A);
}


// number of hyperplanes
static proc arrSize(arr A)
"USAGE:     size(A); arr A;
RETURN:     [int] number of hyperplanes in the arangement
NOTE:       size(A) also useable for multi-arrangements
SEE ALSO:   arrSize, multarrSize
KEYWORDS:   hyperplanes; size; number
EXAMPLE:    example arrSize; shows an example"

{
    return (size(A.l));
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
size(A);
}


// substitute variables
static proc arrSubst(arr A, list #)
"USAGE:     arrSubst(A, #); arr A, ring_variables/integers i,j;
RETURN:     [arr] with the corresponding substitutions made
NOTE:       applies 'subst' on the arrangement
SEE ALSO:   arrSubst
KEYWORDS:   variables; ring_variable; substitute
EXAMPLE:    example arrSubst; shows an example"

{
    if(size(#)%2 != 0){
        ERROR("Odd number of parameter inputs!");
    }

    for(int i=1; i<size(#); i=i+2){
        A = subst(A.l, #[i], #[i+1]);
    }

    return (A);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z,x+z);
arrSubst(A,x,y);
arrSubst(A,x,y,y,z);
}


//============================================================================//
//------------------------------- #8 PRINTING --------------------------------//
//============================================================================//


// prints arrangement in the console
static proc arrPrint(arr A)
"USAGE:     A; A arr
RETURN:     [] better readable output in the console as the newstruct print
SEE ALSO:   arrPrint, arrPrintMatrix
KEYWORDS:   print
EXAMPLE:    example arrPrint;"

{
    A.l;
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x,y,z);
A;
}


// readable as a coeff matrix
proc arrPrintMatrix(arr A)
"USAGE:     arPrintMatrix(arr A)
RETURN:     [] prints arr in matrix form
NOTE:       differs print(matrix(arr A)) since variables included
KEYWORDS:   print; matrix; coefficient
EXAMPLE:    example arrPrintMatrix;"

{
    matrix M = matrix(A);
    ideal I = variables(A);
    for(int k=1; k<=size(I); k++){
        M[1..nrows(M),k] = (M[1..nrows(M),k])*(I[k]);
    }

    print(M);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = mat2arr(random(20,5,4));
A;
arrPrintMatrix(A);
}


//============================================================================//
//------------------------- #9 MANIPULATING VARIABLES ------------------------//
//============================================================================//


// matrix of the corresponding ring_variables
proc varMat(intvec v)
"USAGE:     varMat(v); v intvec
RETURN:     [matrix] M containing the corresponding ring_variables
SEE ALSO:   varMat, varNum, arrSwapVar, arrLastVar
KEYWORDS:   variables; ring_variable
EXAMPLE:    example varMat; shows an example"

{
    matrix M[size(v)][1];
    for(int k=1; k<=size(v); k++){
        M[k,1] = var(v[k]);
    }
    return (M);
}

example
{
"EXAMPLE: 'even' ringvariables"; echo = 2;
ring R = 0,(x(1..6)),dp;
intvec v = 2,4,6;
varMat(v);
}


// number of given variable (enh. version of varNum in dmod.lib)
proc varNum(def u)
"USAGE:     varnum(string s);
            varnum(ring_variable);
RETURN:     [int] number of given ring variable, or 0 if it does not appear
NOTE:       This procedure has the same functionality as varNum from the dmod.lib
            package, but also accepts polys as input.
SEE ALSO:   varMat, varNum, arrSwapVar, arrLastVar
KEYWORDS:   variables; ring_variable; number
EXAMPLE:    example varNum; shows an example"

{
    if(typeof(u) == "string"){
        for(int i=1; i<=nvars(basering); i++){
            if(string(var(i)) == u){return (i);}
        }
        return (0);
    }

    if(typeof(u) == "poly"){
        for(int i=1; i<=nvars(basering); i++){
            if(var(i) == u){return (i);}
        }
        return (0);
    }

    ERROR("Wrong input type, expected string or ring_variable (poly)!");
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),dp;
varNum(y);
ring S = 0,(x(1..5),y(1..5)),dp;
varNum("y(3)");
}


// ideal of all variables the arrangement depends on
static proc arrVariables(arr A)
"USAGE:     variables(A); A arr
RETURN:     [ideal] whereas generators are variables A uses
NOTE:       inherited from the ideal class
SEE ALSO:   varMat, varNum, arrVariables ,arrNvars, arrSwapVar, arrLastVar
KEYWORDS:   variables; ring_variable
EXAMPLE:    example arrVariables; shows an example"

{
    return (variables(A.l));
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),lp;
arr A = ideal(x,y,z);
variables(A);
variables(A-y);
}


// number of variables the arrangement uses
static proc arrNvars(arr A)
"USAGE:     arrNvars(A); A arr
RETURN:     [int] number of variables A uses
NOTE:       inherited from the ideal class
SEE ALSO:   varMat, varNum, arrVariables ,arrNvars, arrSwapVar, arrLastVar
KEYWORDS:   variables; variables; ring_variable; number
EXAMPLE:    example arrNvars; shows an example"

{
    return (size(variables(A)));
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),lp;
arr A= ideal(x,y,z);
arrNvars(A);
arrNvars(A-y);
}


// swaps two variables in the arrangement
proc arrSwapVar(arr A, def i, def j)
"USAGE:     arrSwapVar(A, i, j); arr A, ring_variables/integers i,j
RETURN:     [arr] A where variables i and j are swapped.
NOTE:       if i and/or j are integers the algorithm considers the variables
            variables(A)[i] and/or variables(A)[j]
SEE ALSO:   varMat, varNum, arrSwapVar, arrLastVar
KEYWORDS:   swap; variables; ring_variable
EXAMPLE:    example arrSwapVar; shows an example"

{
    ideal I = variables(A);
    poly u,v;
    if(typeof(i) == "int"){ u = I[i]; } else{ u = i; }
    if(typeof(j) == "int"){ v = I[j]; } else{ v = j; }

    if(u == v){ return (A); }       // special case which messes up the rest

    // using the old trick on how to swap 2 cells without needing a third:
    // (a|b) =(a->a+b)=> (a+b|b) =(b->b-a)=> (b|b-a) =(a->b-a)=> (b|a)
    A = subst(A.l, u, u+v);
    A = subst(A.l, v, v-u);
    A = subst(A.l, u, v-u);

    return (A);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x,y,z),lp;
arr A  = ideal(x+1,x+y,z);
arrSwapVar(A,x,z);
}


//ring_variable of largest index used in arrangement
proc arrLastVar(arr A)
"USAGE:     arrLastVar(A); arr A
RETURN:     [int] number of the last variable A uses
NOTE:       useful if you want a list containing all variables x_1 ... x_k used in A,
            but you do not want to skip any like variables(A) does.
SEE ALSO:   varMat, varNum, arrSwapVar, arrLastVar
KEYWORDS:   variables; ring_variable
EXAMPLE:    example arrLastVar; shows an example"

{
    return ( rvar(variables(A)[arrNvars(A)]) );
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,x(1..10),dp;
arr A = ideal(x(1), x(2), x(3), x(6));
int n = arrLastVar(A);
varMat(1..n);
variables(A);
}


//============================================================================//
//--------------------------- #10 CENTER COMPUTATIONS ------------------------//
//============================================================================//


// checks if arr is central
proc arrCentral(arr A)
"USAGE:     arrCentral(A); arr A
RETURN:     [0,1] true if arr is central(i.e. all planes intersect in 0)
NOTE:       This is the same as homog(A)
SEE ALSO:   arrCentered, arrCentral, arrCenter, arrCentralize
KEYWORDS:   center; central
EXAMPLE:    example arrCentral;"

{
    return (homog(A));
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;

// centered and central
arr A = ideal(x,y,z);
arrCentered(A);
arrCentral(A);

// centered but not central (center: (-1,-1/2, 1))
arr B = ideal(x+1,2y+1,-z+1);
arrCentered(B);
arrCentral(B);
}


// checks wether arrangement has a center
proc arrCentered(arr A)
"USAGE:     arrCentered(A); arr A
RETURN:     [0,1] true if A is centered(i.e. intersection of all planes not empty)
NOTE:       The algorithm uses the rank of matrix: Ax=b has a solution iff
            rank(A) = rank(A|b)
SEE ALSO:   arrCentered, arrCentral, arrCenter, arrCentralize
KEYWORDS:   center
EXAMPLE:    example arrCentered;"

{
    matrix M = matrix(A);
    // classic test: Ax=b has a solution if and only if rank(A|b) == rank(A)
    if( rank(M) == rank(submat(M,1..nrows(M), 1..ncols(M)-1))){ return (1); }
    return (0);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A= ideal(x,y,x-y+1);        // centerless
arrCentral(A);
arr B= ideal(x,y,z);            // central with center being the origin
arrCentral(B);
arr C= ideal(x+1,2y+1,-z+1);    // central with center (-1,-1/2, 1)
arrCentral(C);
}


// computes center of an arrangement
proc arrCenter(arr A)
"USAGE:     arrCenter(A); arr A
RETURN:     [list] L entry 0 if A not centered or entries  1, x, H, where x is
        any particular point of the center and H is a matrix consisting of
        vectors which spanning linear intersection space.
            If there is exactly one solution, then H = 0.
NOTE:       The intersection of all hyperplanes can be expressed in forms of a
            linear system Ax=b, where (A|b) is the coeff. matrix of the arrange-
            ment, which is then solved using L-U decomposition
SEE ALSO:   arrCentered, arrCentral, arrCenter, arrCentralize
KEYWORDS:   center
EXAMPLE:    example arrCenter;"

{
    matrix M = matrix(A);       //return matrix (T|b)
    matrix T = submat(M, 1..nrows(M), 1..ncols(M)-1);
    matrix b = submat(M, 1..nrows(M), ncols(M))*(-1);
    list L = ludecomp(T);
    list Q = lusolve(L[1], L[2], L[3], b);
    return (Q);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;

arr A= ideal(x,y,x-y+1);    // centerless
arrCenter(A);

arr B= ideal(x,y,z);        // center is a single point
arrCenter(B);

arr C= ideal(x,z,x+z);      // center is a line
// here we get a wrong result because the matrix is simplified since A doesn't
// contain any "y" the matrix (A|b) will be 3x3 only.
arrCenter(C);
}


// makes centered arrangement central
proc arrCentralize(arr A)
"USAGE:     arrCenteralize(A); arr A
RETURN:     [arr] A after centralization via coordinate change
NOTE:       The coordinate change only does translation, vector of translation is the second
            output of arrCenter
SEE ALSO:   arrCentered, arrCentral, arrCenter, arrCentralize
KEYWORDS:   central; center; coordinate change
EXAMPLE:    example arrCenter;"

{
    if(arrCentral(A)){
        print("The arrangement is already central!");
        return ();
    }

    list L = arrCenter(A);

    if(L[1] == 0){
        print("The arrangement has no center and therefor cannot be centralized!");
        return ();
    }

    int n = nvars(basering);

    matrix T = diag(1,n);
    matrix c = L[2];

    A = arrCoordChange(A, T, c);

    return (A);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = ideal(x-1,y,x-z-1,x-z-1);
arrCentralize(A);
}


//============================================================================//
//------------------------ #11 GEOMETRIC CONSTRUCTIONS -----------------------//
//============================================================================//

// performs coordinate change
proc arrCoordChange(arr A, matrix T, list #)
"USAGE:     arrCoordChange(A, T);         arr A, (m*n mat) n*n or n*n+1 matrix T
            arrCoordChange(A, T , c);     arr A, n*n matrix T, n*1 matrix/vector
RETURN:     [arr]: Arrangement A [A|b] after a coordinate change f: x -> Tx + c with T invertible
            i.e. [A|b] => [AT^-1|b+AT^-1c] since we have
            f(H) = f(ker(a1*x1 + ... + an*xn - b)) = {f(x) : a'x -b = 0}
                        = {y :  a'f^-1(y) -b = 0)}
                        = {y :  a'(T^-1(y-c)) - b = 0}
                        = {y :  a'T^-1y  -(b + a'T^-1c) = 0}
NOTE:       There are 3 options how you can give the input (in each case n <= nvars(basering))
            1. Just a nxn matrix with
                => Will automatically complete T by a unit matrix and perform x -> Tx
            2. A nxn matrix T and a nx1 vector/matrix c with
                => Will automatically complete T and c and perform x -> Tx +c
            3. A nxn+1 matrix T with
                => will use last column as translation vector c
SEE ALSO:   arrCoordChange, arrCoordNormalize
KEYWORDS:   coordinate change
EXAMPLE:    example arrCoordChange; shows an example"

{
    int n = nvars(basering);
    int k = nrows(T);
    int l = ncols(T);
    matrix c;

    if( k>n || l-1>n ){
        ERROR("Matrix too big! (It has more cols/rows than there are variables.)");
        }

    // const vector integrated in matrix => split [T|c] into T and c
    if( l == k+1 ){
        if(size(#) > 0){
            ERROR("Bad input. If given a constant vector the matrix must be square!")
        }
        c = submat(T, 1..k, l);
        T = submat(T, 1..k, 1..k); l=k;
    }

    if( l != k || rank(T) != k){
        ERROR("Given matrix is not a base change matrix.")
    }

    if(size(#) > 0){
        if(size(#) >= 2){ERROR("Too many input arguments!");}
        c = matrix(#[1]);
        if(nrows(c) > n || ncols(c) > 1){
            ERROR("Constant vector maldefined. Dimension too big.")
        }

    }

    matrix M[n][n] = diag(1,n);
    M[1..k,1..k] = T;
    T = M;
    T = luinverse(T)[2];
    M = jacob(A.l);

    matrix b[n][1];
    b[1..nrows(c),1] = c;
    c = b;                          // gives c the right length.
    b = transpose(jet(A.l, 0));     // constants

    b = b - M*T*c;
    M = M*T;
    A = mat2arr(concat(M,b));

    return(A);
}

example
{
"EXAMPLE: ";    echo =2;
ring r = 0,(x,y,z),lp;
arr A = x,y,z;
arrCoordChange(A,1,[0,0,1]); //lifts z-hyperplane by 1 unit
matrix T[2][2] = [0,1,1,0];  // swaps x and y
arrCoordChange(A,T);
matrix c[2][1] = [1,0];
T = concat(T,c);            // now swap x and y and add 1 to x afterwards
arrCoordChange(A,T);
// Note how T doesn't even need to be a full 3x3 base change matrix.
}

// performs projection onto coordinate hyperplanes
proc arrCoordNormalize(arr A, intvec v)
"USAGE:     arrCoordChange(A, v);
RETURN:     [arr]: Arrangement after a coordinate change that transforms the arrangement such that
            after a tranformation x -> Tx + c we have the arrangement has the matrix representation
            [AT^-1|b+AT^-1c] such that [AT^-1]_v = I and [b+AT^-1c]_v = 0;
NOTE:       algorith performs a base change if H_k is homogenous (i.e. has no)
            constant term and an affine transformation otherwise
            Ax+b = 0, Transformation x = Ty+c: AT^-1y + AT^-1c + b = 0
            Now we want to have (AT^-1)_v = I and (AT^-1c +b)_v = AT^-1_v*c + b_v = 0
SEE ALSO:   arrCoordChange, arrCoordNormalize
KEYWORDS:   coordinate change
EXAMPLE:    example arrCoordChange; shows an example"

{
    int n = nvars(basering);
    if(size(v) > n){
        ERROR("Too many rows chosen, at max you can choose "+string(n));
    }

    matrix M = matrix(A);
    matrix Av[n][n];
    matrix bv[n][1];
    Av[1..size(v), 1..n] = submat(M,v,1..n);
    bv[1..size(v),1]     = submat(M,v, n+1);

    if(rank(Av) != n){
        if(rank(Av) != size(v)){
            ERROR("Normal vectors of the given hyperplanes are not linearily " +
            "independent. Cannot perform coordinate change!");
        }
        // Adds linearly independent lines to Av in order to make it invertible
        module F = freemodule(n);
        module Add = reduce(F,std(transpose(Av)));
        Add = transpose(simplify(transpose(Add),2));
        Av[(size(v)+1)..n, 1..n] = matrix(Add);
    }
    if(rank(Av) != n){ERROR("Av not invertible!");}

    A = arrCoordChange(A,Av,bv);

    return (A);
}
example
{
"EXAMPLE: "; echo=2;
ring r = 0,(x,y,z),lp;
arr A = ideal(x,y,z,x+z+4);
intvec v = 1,2,4;
arrCoordNormalize(A,v);
}

// coned arrangement
proc arrCone(arr A, list #)
"USAGE:     arrCone(A);
            arrCone(A, ring_variable); arr A arrangement in variables x_1...x_n;
RETURN:     arr, the coned hyperplane Arrangement cA with respect to the given
            ring_variable, or the last ring_variable if none was given.
NOTE:       The hyperplanes are homogenized w.r.t. v and a new hyperplane
            H = ker(x_n+1) is added.
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   cone; decone
EXAMPLE:    example arrCone; shows an example"

{
    poly p; int i;

    // case 1: no ring_variable given
    if(size(#) == 0){
        p = var(nvars(basering));
    }

    // case 2: a ring_variable is given
    else{ p =  #[1]; }

    // computation
    A = homog(A, p);
    A = A + p;

    return(A);
}

example
{
"EXAMPLE: Coning the arrangement A = (x+1, y) alongside x, y and z:"; echo=2;
ring R = 0,(x,y,z),dp;
arr  A = ideal(x+1, x,x-2,x-1);
arrCone(A, y);
arr B= ideal(x,y,x+y-1);
arrCone(B);
}


// deconed arrangement
proc arrDecone(arr A, int k)
"USAGE:     arrDecone(A, k); arrangement A, integer k;
RETURN:     arr: the deconed hyperplane Arrangement dA
NOTE:       A has to be non-empty and central. arrDecone is an inverse operation
            to arrCone since A == arrDecone(arrCone(A),size(A)+1) for any A.
            One can also decone a central arrangement with respect to any hyper-
            plane k, but than a coordinate change is necessary to make
            H_k = ker(x_k). Since such a coordinate change is not unique,
            use arrCoordchange to do so.
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   cone; decone
EXAMPLE:    example arrDecone; shows an example"

{
    if( !arrCentral(A) ){
        ERROR("Non-central arrangement can not be deconed!");}
    if( size(A) == 0){
        ERROR("Empty arrangement can not be deconed!");}
    if( size(A) < k){
        ERROR("There is no k-th hyperplane");}

    poly p = A[k];
    int n = rvar(p);

    if( n == 0 ){ ERROR("H_" + string(k) + " = " + string(p) +
        " is not of the form ker(x_i). Please do a coordinate change first." +
        "You can use arrCoordinateChange to transform the arrangement accordingly.");
        }

    A = A - p;
    A = arrSubst(A, p, 1);

    return(A);
}

example
{
"EXAMPLE: We decone arr consisting of (x,y,x+z) with respect to y";
echo = 2;
ring R = 0,(x,y,z),dp;
arr A= ideal(x,y,z,x+y-z);
arrDecone(A,3);
}

proc arrLocalize(arr A, intvec v)
"USAGE:     arrLocalize(A, v); arrangement A, intvec v;
RETURN:     arr: the localized arrangement A_X, i.e. A_X only contains the hyperplanes
            which contain the flat X, which is defined by the equations A[v]
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   localization
EXAMPLE:    example arrLocalize; shows an example"
{
    ideal I=std(arr2ideal(A[v]));
    arr L;
    int k;
    for(k=1; k<size(A); k++){
        if(NF(A[k],I)==0){
            L=L+A[k];
        }
    }
    return(L);
}

example
{
"EXAMPLE: We localize the Braid 3 arrangement at x and y";
echo = 2;
ring R = 0,(x,y,z),dp;
arr A=arrTypeB(3);
intvec v=5,8;
arr B=arrLocalize(A,v);
B;
}

// restricted arrangement onto a hyperplane
proc arrRestrict(arr A, intvec v, list #)
"USAGE:     arrRestrict(A, v); arrangement A, int/intvec v, optional argument "CC";
RETURN:     arr: the restricted hyperplane Arrangement (A^X)
NOTE:       A has to be non-empty.
REMARKS:    We restrict A to the flat X, defined by the equations in A[v].
            The restriction will only be performed, if the ideal defining
            the flat X is monomial (i.e. X is an intersection of coordinate planes).
            If the optional argument CC is given, the arrangement is transformed
            in such a way that X has the above form.
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   restriction
EXAMPLE:    example arrRestrict; shows an example"

{
    ideal I=A.l;
    I=I[v];
    option(redSB);
    I=std(I); // defining equations for flat X
    //option(none);
    ideal Al=arr2ideal(A);
    int i;

    if(isMonomial(I)==1){
        for(i=1;i<= size(I);i++){
            Al=subst(Al,I[i],0);
        }
        arr AR=Al;
        return(AR);
    }

    if(size(#) == 0){
        ERROR("The flat X is not defined by a monomial ideal. " +
        "Please do a coordinate change first. You can use arrCoordNormalize to " +
        "transform the arrangement accordingly by adding the argument CC.");
    }
    if(#[1]!="CC"){
        ERROR("The flat X is not defined by a monomial ideal. " +
        "Please do a coordinate change first. You can use arrCoordNormalize to " +
        "transform the arrangement accordingly by adding the argument CC.");
        }
    if(size(v)==0){return(A);}
    intvec w=v[1];
    intvec tmp;

    for(i=2;i<=size(v);i++){
        tmp=w,v[i];
        if(rank(matrix(A[w]))==size(tmp)){
            w=tmp;
        }
    }

    arr B=arrCoordNormalize(A,w);
    ideal Bl=arr2ideal(B);

    for(i=1;i<= size(w);i++){
        Bl=subst(Bl,Bl[w[i]],0);
    }
    arr BR=Bl;
    return(BR);
}

example
{
"EXAMPLE: We consider the possible restrictions of the type B3 arrangement ";
echo = 2;
ring S = 0,(x,y,z),dp;
arr A  = arrTypeB(3);
A;
arrRestrict(A,9);
arrRestrict(A,4,"CC");
intvec v=5,8;
arrRestrict(A,v);
}


// checks if arrangement is essential
proc arrIsEssential(arr A)
"USAGE:     arrIsEssential(A); arrangement A;
RETURN:     boolean: 1 if arr is essential, i.e. rank of maximal element of
            poset is dimension
NOTE:       A has to be non-empty.
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   essential
EXAMPLE:    example arrIsEssential; shows an example"

{
    ideal I = variables(A);

    if( var(size(I)) != I[size(I)] ){
        return(0);
    } if( rank(jacob(A.l)) == size(I) ){
        return(1);
    }

    return(0);
}

example
{
"EXAMPLE: We check wether theseare essential and a non-essential arrangement ";
echo = 2;
ring S = 0,(x,y,z),lp;
arr A = ideal(x,y,z);
arr B = ideal(x+y+z,x,y+z);
arrIsEssential(A);
arrIsEssential(B);
}


// essentialized arragnement
proc arrEssentialize(arr A)
"USAGE:     arrEssentialize(A); arrangement A;
RETURN:     essential arrangement by transformation
NOTE:       A has to be non-empty.
SEE ALSO:   arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize
KEYWORDS:   essential
EXAMPLE:    example arrEssentialize; shows an example"

{
    ideal I = variables(A);

    if( var(size(I)) != I[size(I)] ){
        for(int k=1; k<=size(I); k++){
            if( NF(var(k),std(I)) != 0 ){
                A = subst(A.l, I[size(I)], var(k));
                return (arrEssentialize(A));
            }
        }
    }

    while(arrIsEssential(A) == 0){
        A = subst(A.l, I[size(I)], 0);
    }

    return (A);
}

example
{
"EXAMPLE: We essentialize a non essential arrangement "; echo = 2;
ring S = 0,(x,y,z),dp;
arr A=arrBraid(3);
arrEssentialize(A);
arr B = ideal(x+y+z,x,y+z);
arrEssentialize(B);
}

//============================================================================//
//------------------------ #12 EXAMPLES OF ARRANGEMENTS ----------------------//
//============================================================================//


// boolean arrangement
proc arrBoolean(int v)
"USAGE:     arrBoolean(v); int v
RETURN:     arr, which uses the first v variables of ring for boolean arrangement
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; boolean
EXAMPLE:    example arrBoolean;"

{
    if(nvars(basering)<v){ return("not enough variables"); }

    arr A = mat2carr(unitmat(v));
    return(A);
}

example
{
"EXAMPLE:Boolean Arrangement, uses the seven first coordinate hyperplanes";
echo = 2;
ring R = 0,x(1..10),dp;
arrBoolean(7);
}


// braid arrangement
proc arrBraid(int v)
"USAGE:     arrBraid(v); int v
RETURN:     Type A (braid) arrangement of dimension v
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; braid
EXAMPLE:    example arrBraid;"

{
    if(nvars(basering)<v){
        return("not enough variables");
    }

    arr A;
    int i,j;
    for(i=1; i<=v; i++){
        for(j=i+1; j<=v; j++){
            A = A + (var(i) - var(j));
        }
    }
    return(A);
}

example
{
"EXAMPLE: The braid arrangement consists of the hyperplanes x(i)=x(j)";
echo = 2;
ring R = 0,x(1..10),dp;
arrBraid(7);
}


// type B arrangement
proc arrTypeB(int v)
"USAGE:     arrTypeB(v); int v
RETURN:     arrangement, which uses first v variables of ring for reflection
            arrangement of type B
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; type B
EXAMPLE:    example arrTypeB;"

{
    if(nvars(basering)<v){
        return("not enough variables");
    }

    arr A;
    int i,j;
    for(i=1; i<=v; i++){
        for(j=i+1; j<=v; j++){
            A = A + (var(i) - var(j));
            A = A + (var(i) + var(j));

        }
        A = A + var(i);
    }
    return(A);
}

example
{
"EXAMPLE: The type B reflection arrangement consists of the hyperplanes " +
"x(i)=x(j), x(i)=-x(j) and the coordinate hyperplanes";
echo = 2;
ring R = 0,x(1..10),dp;
arrTypeB(5);
}


// type D arrangement
proc arrTypeD(int v)
"USAGE:     arrTypeD(v); int v
RETURN:     arrangement, which uses first v variables of ring for reflection
            arrangement of type D
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; type D
EXAMPLE:    example arrTypeD;"

{
    if(nvars(basering)<v){
        return("not enough variables");
    }

    arr A;
    int i,j;
    for(i=1; i<=v; i++){
        for(j=i+1; j<=v; j++){
            A = A + (var(i) - var(j));
            A = A + (var(i) + var(j));
        }
    }
    return(A);
}

example
{
"EXAMPLE: The type D reflection arrangement consists of the hyperplanes " +
"x(i)=x(j) and x(i)=-x(j)";
echo = 2;
ring R = 0,x(1..10),dp;
arrTypeD(5);
}


// random (affine) arrangement
proc arrRandom(int d, int m, int n)
"USAGE:     arrRandom(n,v,N); int n,v,N
RETURN:     Random arrangement, where m is the number of hyperplanes, n the
            dimension, d the upper bound for absolute value of coefficients.
NOTE:       You can also write arr = random(d,m,n) to create random arrangements
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; random
EXAMPLE:    example arrRandom;"

{
    if(n > nvars(basering)){
        return("Error, too few variables or too high dimension");
    }

    intmat M = random(d,m,n+1);
    arr A = mat2arr(M);
    return(A);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,x(1..20),dp;
arrRandom(7,3,15);
}


// random central arrangement
proc arrRandomCentral(int d, int m, int n)
"USAGE:     arrRandomCentral(d,m,n); int d,m,n
RETURN:     Random central arrangement, where m is the number of hyperplanes, n
            the dimension, d the upper bound for absolute value of coefficients.
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; random; central
EXAMPLE:    example arrRandomCentral;"

{
    if(n > nvars(basering)){
        return("Error, too few variables or too high dimension");
    }

    intmat M = random(d,m,n);
    arr A = mat2carr(M);
    return(A);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,x(1..20),dp;
arrRandomCentral(7,3,15);
}


// Edelman-Reiner arrangement
proc arrEdelmanReiner()
"USAGE:     arrEdelmanReiner();
RETURN:     the Edelman-Reiner arrangement, which is a free arrangement but the
            restriction to the 6-th hyperplane is nonfree.
            (i.e. counterexample for Orlik-Conjecture)
NOTE:       the active ring must have at least five variables
SEE ALSO:   arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner
KEYWORDS:   example; Edelman-Reiner
EXAMPLE:    example arrEdelmanReiner;"

{
    if(nvars(basering) < 5){ ERROR("not enough variables"); }

    arr arrER = arrBoolean(5);
    int a,b,c,d;
    for(a=1; a<=2; a++){
    for(b=1; b<=2; b++){
    for(c=1; c<=2; c++){
    for(d=1; d<=2; d++){
        arrER = arrER + (var(1)+(-1)^a*var(2)+(-1)^b*var(3)+(-1)^c*var(4)+(-1)^d*var(5));
    }}}}

    return(arrER);
}

example
{
"EXAMPLE:"; echo = 2;
ring r=0,x(1..5),dp;
arrEdelmanReiner();
}


//============================================================================//
//--------------------- #13 Orlik Solomon and Poincare Poly ------------------//
//============================================================================//


// Orlik-Solomon algebra of the arrangement
proc arrOrlikSolomon(arr A)
"USAGE:     arrOrlikSolomon(A); arr A
RETURN:     [ring] exterior Algebra E as ring with Orlik-Solomon ideal as attribute I.
            The Orlik-Solomon ideal is generated by the differentials of dependent
            tuples of hyperplanes. For a complex arrangement the quotient E/I is
            isomorphic to the cohomology ring of the complement of the arrangement.
NOTE:       In order to access this ideal I activate this exterior algebra with setring.
SEE ALSO:   arrOrlikSolomon
KEYWORDS:   Orlik-Solomon
EXAMPLE:    example arrOrlikSolomon;"

{
    int central = arrCentral(A);  // 0: non-central, 1: central
    if(central == 0){
        A = arrCone(A);
    }

    module M = syz(A.l);
    M = jet(M,0); //Only use linear syzygies
    M = simplify(M,2); // drop zeros

    int n = ncols(A.l);
    def startRing = basering;
    ring R = 0,e(1..n),dp;
    def ER = Exterior();
    setring ER; // defines the Exterioralgebra
    ideal I; //final orlik solomon ideal

    matrix X = transpose(varMat(1..n));
    module M = fetch(startRing,M); //brings the module M to the Exterior Algebra
    ideal OSI = ideal(X*M);

    if(size(OSI) == 0){ // no relations among hyperplanes, I==0
        export(I);
        return(basering);
    }

    // monomial subideal procedure due to Sorin Popescu
    ideal K = OSI;
    ideal J;
    while(isMonomial(K) == 0){
        J = lead(std(K));
        K = intersect(J,OSI);
        K = interred(K);
    }

    OSI = K;
    poly diffOp; //now applying the differential operator
    ideal vars;
    int j;

    for(int i=1; i<=ncols(OSI); i++)
    {
        vars = variables(OSI[i]);
        diffOp = 0;
        for(j=1; j<=ncols(vars); j++)
        {
            diffOp = diffOp+(-1)^(j+1)*vars[j];
        }

        I = I + ideal( diff(diffOp,ideal(OSI[i])) );
    }

    if(central ==0 ){//project back in the non-central case
        n = nvars(basering)-1;
        ring R = 0,e(1..n),dp;
        def ERDC = Exterior();
        setring ERDC;
        ideal I = fetch(ER,I);
    }

    export(I);
    return (basering);
}

example
{
"EXAMPLE: Computing the Orlik-Solomon-Ideal for the D3-Arrangement"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = arrTypeB(3);
def E = arrOrlikSolomon(A);
setring E;
//The generators of the Orlik-Solomon-Ideal are:
I;
}


// The following 3 procedures were replaced by their faster combinatorical counterparts

/*

// poincare polynomial of the arrangement
proc arrPoincare(arr A)
"USAGE:     arrPoincare(A); arr A
RETURN:     [intvec] The Poincare polynomial as integer vector of the arrangement, which
            is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   Poincare polynomial
EXAMPLE:    example arrPoincare;"

{
    def startRing = basering;
    def Ext = arrOrlikSolomon(A);
    setring Ext;
    ideal OSI = lead(I);
    OSI = OSI + ideal(Ext);
    int n = nvars(Ext);
    ring suppRing = 0,x(1..n),dp;
    ideal I = fetch(Ext,OSI);
    intvec HP = hilb(std(I),2);

    return(HP);
}

example
{
"EXAMPLE: Computing the Poincare polynomial as intvec for the D3-Arrangement"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = arrTypeB(3);
//The coefficients of the Poincare polynomial are:
arrPoincare(A);
}


// number of chambers of the arrangement
proc arrChambers(arr A)
"USAGE:     arrChambers(A); arr A
RETURN:     [int] The number of chambers of an arrangement, which is equal to the
            evaluation of the Poincare polynomial at 1.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   chambers
EXAMPLE:    example arrChambers;"

{
    intvec HP = arrPoincare(A);
    intvec ones = 1:(size(HP));
    return ( transpose(ones)*HP );
}

example
{
"EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = arrTypeD(3);
//The number of chambers of the D3-Arrangement is:
arrChambers(A);
}


// number of bounded chambers of the arrangement
proc arrBoundedChambers(arr A)
"USAGE:     arrBoundedChambers(A); arr A
RETURN:     [int] The number of bounded chambers of an arrangement, which is equal to
            the evaluation of the Poincare polynomial at -1.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   chambers
EXAMPLE:    example arrBoundedChambers;"

{
    intvec HP = arrPoincare(A);
    intvec altOnes;
    for(int i=1; i<=size(HP); i++){
        altOnes[i] = (-1)^(i-1);
    }

    return ( transpose(altOnes)*HP );
}

example
{
"EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = arrTypeD(3);
// The number of bounded chambers of the D3-Arrangement is:
arrBoundedChambers(A);
}

*/

//============================================================================//
//------------------------------- #14 Freeness -------------------------------//
//============================================================================//


// module of derivation
proc arrDer
"USAGE:     arrDer(A); arr A , multarr A
RETURN:     [module] The module Der(A) of derivations of the (multi-)arrangement A, i.e.
            the derivations tangent to each hyperplane of A (resp. with multiplicities)
NOTE:       This is only defined for central (multi-)arrangements
SEE ALSO:   arrDer, arrIsFree, arrExponents
KEYWORDS:   derivation; multiarrangement
EXAMPLE:    example arrDer;"

{
    def A = #[1];
    if( (typeof(A) != "arr") && ( typeof(A) != "multarr") ){
        ERROR("bad input type!");
    } if(homog(A.l) == 0){
        ERROR("Arrangement not central!");
    }

    module fA = jacob(A.l);
    ideal J;

    if( typeof(A) == "arr" ){
        J = arr2ideal(A);
    } if(typeof(A) == "multarr"){
        J = multIdeal(A);
    }

    module tA = diag(matrix(J));
    module K  = modulo(fA,tA);
    module derivations = lessGenerators(K);

    return (derivations);
}

example
{
"EXAMPLE: Computing the derivation module of the boolean and braid arrangement";
echo = 2;
ring R = 0,(x,y,z),dp;
arr A3 = arrBoolean(3);
arr B3 = arrTypeB(3);
arr G = ideal(x,y,z,x+y+z);
//The derivation module of the Boolean 3-arrangement:
arrDer(A3);
//The derivation module of the Braid 3-arrangement:
arrDer(B3);
//The derivation module of the generic arrangement:
arrDer(G);
}


// checks if arrangement is free
proc arrIsFree(list #)
"USAGE:     arrIsFree(A); arr A, multarr A
RETURN:     [0,1] 1 if the (multi-)arrangement is free, i.e. Der(A) is a free module
NOTE:       only defined for central arrangements
SEE ALSO:   arrDer, arrIsFree, arrExponents
KEYWORDS:   free; multiarrangement
EXAMPLE:    example arrIsFree;"

{
    def A = #[1];
    module derivations = arrDer(A);
    return ( nvars(basering) == ncols(derivations) );
}

example
{
"EXAMPLE: checking freeness of the Edelman-Reiner arrangement and its restriction: ";
echo = 2;
ring R = 0,(x,y,z),dp;
arr A3 = arrBoolean(3);
arr B3 = arrTypeB(3);
arr G = ideal(x,y,z,x+y+z);

arrIsFree(A3);

arrIsFree(B3);

arrIsFree(G);
}


// exponents of a (free) arrangement
proc arrExponents
"USAGE:     arrExponents(A); arr A, multarr A
RETURN:     [intvec] The exponents of a free (multi-) arrangement, i.e. the degrees of a
            basis of D(A) the derivation module.
NOTE:       only defined for central arrangements
SEE ALSO:   arrDer, arrIsFree, arrExponents
KEYWORDS:   free; exponents; multiarrangement
EXAMPLE:    example arrExponents;"

{
    def A = #[1];

    if(arrIsFree(A) != 1){ ERROR("Arrangement is not free!"); }

    module der = arrDer(A);
    intvec exp;

    for(int i =1; i<=size(der); i++){
        exp[i] = deg(der[i]);
    }

    return(exp);
}

example
{
"EXAMPLE: computing the exponents of the Edelman-Reiner arrangement and its restriction: ";
echo = 2;
ring R = 0,(x,y,z),dp;
arr A3 = arrBoolean(3);
arr B3 = arrTypeB(3);
arr G = ideal(x,y,z,x+y+z);

arrExponents(A3);

arrExponents(B3);

}


// Tries to reduce the number of generators of a generating set for a module
static proc lessGenerators(module X){
    module Z = syz(X);
    matrix K = getColumnIndependentUnitPositions(Z);

    if( K == unitmat(nrows(Z)) ){
        return(X);
    }

    module Xnew = X*K;
    Xnew = simplify(Xnew,2);

    return(lessGenerators(Xnew));
}


// Looks for a generator which is redundant
static proc getColumnIndependentUnitPositions(module Z){
    int n = nrows(Z); // number of generators of D
    matrix K = unitmat(n);
    int i;
    for(int j=1; j<=ncols(Z); j++){
        for(i=1; i<=nrows(Z); i++){
            if(deg(Z[i,j]) == 0){
                K[i,i] = 0;
                return(K);
            }
        }
    }

    return(K);
}


// Outputs the ideal of powers of hyperplanes of a multiarrangement.
// Needed for the computation of arrDer.
static proc multIdeal(multarr A){
    ideal I;

    for(int i=1; i<=size(A.l); i++){
        I[i] = A.l[i]^A.m[i];
    }

    return(I);
}


//============================================================================//
//-------------------------- #15 MULTI-ARRANGEMENTS --------------------------//
//============================================================================//

//============================================================================//
//------------------------------- CONSTRUCTORS -------------------------------//
//============================================================================//


// general method for creating multiarrangements
static proc multarrAdd
"USAGE:     A = #; A +#;  # list containing arr/ideal/list/matrix/poly
RETURN:     [multarr] multiarrangement constructed by the input parameters.
NOTE:       algorithm splits up the list # and uses appropiate procedure
            to handle the input
KEYWORDS:   multiarrangement; equal; constructor; operator
EXAMPLE:    example multarrAdd; shows an example"

{
    multarr A;
    for(int k=1; k<=size(#); k++){
        while(1){   //simulates switch, which singular doesn't offer
            if(typeof(#[k]) == "poly" )  {A = multarrAddPoly (A, #[k]);break;}
            if(typeof(#[k]) == "ideal")  {A = multarrAddIdeal(A, #[k]);break;}
            if(typeof(#[k]) == "multarr"){A = multarrAddArr  (A, #[k]);break;}
            if(typeof(#[k]) == "list" )  {A = multarrAdd     (   #[k]);break;}
            ERROR("bad input type");
        }
    }
    return (A);
}

example
{
"EXAMPLE: Creating a few multiarrangements"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A = ideal(x,y,z,x);            A;
multarr B = A + ideal(x+1, x-1,y);     B;
multarr C = list(B, x+1, x-1);  C;
multarr D  = (x2 - y2)^2;           D;
}


// adds a single  poly  to the arrangement
// if the poly is linear, it is just added. If not Singular tries factorization
static proc multarrAddPoly(multarr A, poly p)
{
    p=simplify(p,1);
    if(deg(p) == 0){
        ERROR("Given poly is not linear or Singluar is not able to factorize it");}
    else{
        if(deg(p) == 1){
            int k;
            int b = 0;
            for(k=1; k<=size(A.l); k++){
                if(A.l[k] == p){
                    A.m[k] = A.m[k]+1;
                    b = 1;
                }
            }
            if(b == 0){
                A.l[size(A.l)+1] = p;
                A.m[size(A.l)] = 1;
            }
            return(A);
        }
        else{
            list I = factorize(p,2);
            if((size(I[1]) == 1) && (deg(I[1][1] > 1))){ERROR("Given poly is not a hyperplane");}
            else{
                int j,i;
                ideal J;
                for(i=1; i<=size(I[1]); i++)
                    {
                        for(j=1; j<=I[2][i]; j++){
                            J[size(J)+1] = I[1][i];
                        }
                    }
                return (multarrAdd(A,J));
            }
        }
    }
    return(A);
}


// adds defining polys to the arrangement
static proc multarrAddIdeal(multarr A, ideal I){
    for(int k=1; k<=size(I); k++){
        A = multarrAddPoly(A,I[k]);
    } return (A);
}


// union of two multarrangements
static proc multarrAddArr(multarr A, multarr   B){
    return (multarrAddIdeal(A, multIdeal(B)));
}


//============================================================================//
//------------------------------- TYPE CASTING -------------------------------//
//============================================================================//


// computes the defining polynomial
static proc multarr2poly(multarr A)
"USAGE:     multArrQPoly(multarr A);
RETURN:     [poly] q: defining polynomial of a multiarrangement with multiplicities
SEE ALSO:   arr2poly, multarr2poly
KEYWORDS:   multiarrangement; defining polynomial
EXAMPLE:    example multArrQPolys;"

{
    poly q=1;
    for(int i=1;i<=size(A.l);i=i+1)
    {
        q=q*A.l[i]^A.m[i];
    }
    return(q);
}

example
{
"EXAMPLE: Computing the Q-Poly for a multiarrangement"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A = ideal(x,y,z,x,y,x-y,x-z,x-z,y-z);            A;
poly q=multarr2poly(A);     q;
}


// converts simple arrangement to multiarrangement
proc arr2multarr(arr A, intvec v)
"USAGE:     multArrFromIntvec(arr A, intvec v);
RETURN:     [multarr] multiarrangement MA, which is the arrangement A with multiplicities v
NOTE:       the size of v must match the number of hyperplanes of the arrangement A
SEE ALSO:   arr2multarr, multarr2arr
KEYWORDS:   multiarrangement; arrangement; constructor
EXAMPLE:    example arr2multarr"

{
    if(ncols(A.l)!=size(v)){
        ERROR("Vector's size does not match the hyperplanes.");
    }
    multarr B;
    B.l=A.l;
    B.m=v;
    return(B);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
arr A = arrTypeB(3);            A;
intvec v=2:9;               v;
multarr MA=arr2multarr(A,v);
MA;
}


// converts multiarrangement to simple arrangement
proc multarr2arr(multarr A)
"USAGE:     multarr2arr(multarr A, intvec v);
RETURN:     [arr] arrangement A, with all multiplicities removed
SEE ALSO:   arr2multarr, multarr2arr
KEYWORDS:   multiarrangement; arrangement; constructor
EXAMPLE:    example multarr2arr"

{
    arr B;
    B.l=A.l;
    return(B);
}

example
{
"EXAMPLE:"; echo = 2;
ring r = 0,(x,y,z),dp;
multarr A=x2y3z5;   A;
arr AS = multarr2arr(A); AS;
}


//============================================================================//
//------------------------------- PRINTING -----------------------------------//
//============================================================================//


// prints arrangement in the console
static proc multarrPrint(multarr A)
"USAGE:     A; A arr
RETURN:     [] better readable output in the console as newstruct print.
SEE ALSO:   arrPrint, multarrPrint
KEYWORDS:   print
EXAMPLE:    example multarrPrint;"

{
    for(int j=1;j<=ncols(A.l);j++){
        print("_["+string(j)+"]=("+string(A.l[j])+")^"+ string(A.m[j]));
    }
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A = ideal(x2,y3,z);
A;
}


// number of hyperplanes with multiplicities
static proc multarrSize(multarr A)
"USAGE:     size(A); A multarr
RETURN:     [int] Number of hyperplanes with multiplicities
SEE ALSO:   arrSize
KEYWORDS:   multiarrangement; size; number; hyperplanes
EXAMPLE:    example multarrSize;"

{
    return (sum(A.m));
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A = ideal(x2,y3,z);
A;
}


//============================================================================//
//------------------------ RESTRICTION & DELETION ----------------------------//
//============================================================================//


// decrements the multiplicity of a hyperplane by one
static proc multarrDelete(multarr A, int k)
"USAGE:     multarrDelete(A, k); arrangement A, integer k;
RETURN:     [multarr] the hyperplane Multiarrangement A', i.e. the multiarrangement
            with multiplicity of H_k decremented by one. If m(H_k)=1, then the hyperplane
            H_k is deleted
SEE ALSO:   arrDelete, multarrDelete,
KEYWORDS:   multiarrangement; delete
EXAMPLE:    example multarrDelete;"

{
        if(k>ncols(A.l)){
            ERROR("There is no k-th hyperplane");
        }

        poly q = multarr2poly(A);
        q = q/(A.l[k]);
        multarr MA = q;
        return (MA);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A =ideal(x2,y3,z);
multarr AD = multarrDelete(A,2);  AD;
}


// restriction of A (as arr) to a hyperplane with multiplicities
proc multarrRestrict(arr A,intvec v, list #)
"USAGE:     multarrRestrict(A, v); arrangement A, int/intvec v, optional argument "CC";
RETURN:     [multarr] the restricted hyperplane Multi-Arrangement (A^X) with multiplicities
            i.e. counting how often one element of the restricted arrangement occurs
            as intersetion of hyperplane of the first arrangement. This definition
            is due to Guenter M. Ziegler.
NOTE:       A has to be non-empty.
REMARKS:    We restrict A to the flat X, defined by the equations in A[v].
            The restriction will only be performed, if the ideal defining
            the flat X is monomial (i.e. X is an intersection of coordinate planes).
            If the optional argument CC is given, the arrangement is transformed
            in such a way that X has the above form.
SEE ALSO:   multarrRestrict, multarrMultRestrict, arrRestrict
KEYWORDS:   multiarrangement; restriction
EXAMPLE:    example multarrRestrict; "

{
    ideal I=A.l;
    I=I[v];
    option(redSB);
    I=std(I); // defining equations for flat X
    //option(none);
    ideal Al=arr2ideal(A);
    int i;

    if(isMonomial(I)==1){
        for(i=1;i<= size(I);i++){
            Al=subst(Al,I[i],0);
        }
        multarr AR=simplify(Al,3);
        return(AR);
    }

    if(size(#) == 0){
        ERROR("The flat X is not defined by a monomial ideal. " +
        "Please do a coordinate change first. You can use arrCoordNormalize to " +
        "transform the arrangement accordingly by adding the argument CC.");
    }
    if(#[1]!="CC"){
        ERROR("The flat X is not defined by a monomial ideal. " +
        "Please do a coordinate change first. You can use arrCoordNormalize to " +
        "transform the arrangement accordingly by adding the argument CC.");
        }
    if(size(v)==0){return(A);}
    intvec w=v[1];
    intvec tmp;

    for(i=2;i<=size(v);i++){
        tmp=w,v[i];
        if(rank(matrix(A[w]))==size(tmp)){
            w=tmp;
        }
    }

    arr B=arrCoordNormalize(A,w);
    ideal Bl=arr2ideal(B);

    for(i=1;i<= size(w);i++){
        Bl=subst(Bl,Bl[w[i]],0);
    }
    multarr BR=simplify(Bl,3);
    return(BR);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,x(1..5),dp;
arr A = arrEdelmanReiner();   A;
multarr AR = multarrRestrict(A,6,"CC");  AR;
}


// restriction of A (as multarr) to a hyperplane with multiplicities
proc multarrMultRestrict(multarr A,int k)
"USAGE:     multarrMultRestrict(A, k); multiarrangement A, integer k;
RETURN:     [multarr] the restricted hyperplane Multi-Arrangement (A^H_k) with multiplicities, i.e.
            counting with multiplicities how often one element of the restricted arrangement occurs
            as intersetion of hyperplane of the first multiarrangement.
            This definition is due to Guenter M. Ziegler.
NOTE:       A has to be non-empty.
REMARKS:    The restriction will only be performed, if H_k = ker(x_i) for some i.
            One can also restrict an arrangement with respect to any hyper-
            plane k, but than a coordinate change is necessary first to make
            H_k = ker(x_k). Since such a coordinate change is not unique, please
            use arrCoordchange to do so.
SEE ALSO:   multarrRestrict, multarrMultRestrict, arrRestrict
KEYWORDS:   multiarrangement; restriction
EXAMPLE:    example multarrMultRestrict; "

{
    ideal I = variables(A.l);
    ideal J;
    int j,i;
    for(i=1;i<=size(A.l);i=i+1)
    {
        for(j=1;j<=A.m[i];j++){
            J[ncols(J)+1]=A.l[i];
        }
    }
    poly p = simplify(A.l[k],1);
    int n = rvar(p);
    if( n == 0 ){
        ERROR("H_" + string(k) + " = " + string(p) +
        " is not of the form ker(x_i). Please do a coordinate change first." +
        "You can use arrCoordinateChange to transform the arrangement accordingly.");
    }
    J = subst(J, p, 0);
    multarr MA = simplify(J,3);
    return(MA);
}

example
{
"EXAMPLE:"; echo = 2;
ring R = 0,(x,y,z),dp;
multarr A =ideal(x2,y2,z2,(x-y)^3,(x-z)^2,(y-z));    A;
//The restriction of the multiarrangement is:
multarr AR = multarrMultRestrict(A,1);  AR;
}


//============================================================================//
//----------------------------- #16  COMBINATORICS ---------------------------//
//============================================================================//

    /*
    RELATIONS can have 3 values:
    -1 : flat&hyperplane are parallel
    0 : hyperplane intersects the flat, but only further on in the lattice.
    +1 : hyperplane is part of the flat.

    */
    // intvec: multiplicities of hyperplanes


proc arrLattice(arr ARR)
"USAGE:     arrLattice(arr ARR)
RETURN:     [arrposet] intersection poset of the arrangement
NOTE:       The algorithm works by a bottom up approach, i.e. it calculates the
SEE ALSO:   arrLattice, arrFlats
KEYWORDS:   intersection lattice; poset; lattice
EXAMPLE:    example arrFlats;"

{

//Performance test
system("--ticks-per-sec",1);
timer = 0;
int start = timer;
// start
dbprint(3-voice,newline);
dbprint(3-voice,"=== Computing poset ===");
dbprint(3-voice,newline);

// Initialization
list L,L2;
intvec u,v, NFLAT;
int i,j,l,m,n,rk, rk2,rk3, tic, toc, numberOfFlats,numberOfHplanes, counter, ntests;
matrix NewRel, M;

// Initialization of matrix
ntests = 0;
M =  matrix(ARR);                  // m x n+1 matrix
m = size(ARR);                            // m = # hplanes
n = nvars(basering);                      // n = # variables
numberOfHplanes = m;

// Initialization of poset and flat
flat F;
arrposet P;         P.A = ARR;
F.REL = 0:m;        F.moebius = -1;
for(i=1; i<=n; i++){ P.r = insert(P.r,list());}
for(i=1; i<=m; i++){ F.REL[i] = 1;  L[i] = F; F.REL[i] = 0;}
F.moebius = 0;
P.r[1] = L;

// calculating r[i] step by step
for(rk=2; rk<=rank(ARR); rk++){                    // maximal flat possible has rank A.v ofc
    counter = 0;
    tic = timer;

    // Initialization of current round
    L = P.r[rk-1];                      // flats from the last iteration give rise to current ones.
    numberOfFlats = size(L);
    L2 = list();                        // new list for elts of rk=
    NewRel = getOldRel(L,m);            // contains all the information of the old arrangement

    // find the next leftmost child by looking for intersetion of F_j with hyperplanes
    for(j=1; j<=numberOfFlats; j++){           // goes over the elts of L.r[i-1]

        NFLAT = L[j].REL;

        for(i=1; i<=numberOfHplanes; i++){
        if (NewRel[i,j] == 0){
            //test this hyperplane
            NFLAT[i] = 1;
            v = collectHplanes(NFLAT);
            rk2 = rank(submat(M,v,1..n));
            ntests = ntests+2;

            // CHILD FOUND BY INTERSECTING WITH H_i
            if(rk2 == rank(submat(M,v,1..n+1))){
                NewRel[i,j] = 1; counter++;

                // potentially there exist more hyperplanes which intersect in the same space
                for(l=i+1; l<=m; l++){
                if(NewRel[l,j] == 0){
                    NFLAT[l] = 1;
                    u = collectHplanes(NFLAT);
                    rk3 = rank(submat(M, u, 1..n));
                    ntests = ntests +2;

                    // FLAT EXISTS
                    if(rk3 == rank(submat(M,u,1..n+1)) ){
                        if(rk2 == rk3){NewRel[l,j]   = 1;}
                        else{NFLAT[l] = 0;} //itersecting with both yields higher rank flat
                    }

                    // FLAT DOES NOT EXIST
                    else{NFLAT[l] = -1;}
                }
                }

                // Add new flat to list, reset NFLAT.
                F.REL = NFLAT;
                NFLAT = L[j].REL;                   // refresh.

                // Find Parents: (most expensive part it seems)
                F.parents = j:1;                      //parent in any case
                for(l=j+1; l<=numberOfFlats; l++){
                    if(isParent(L[l].REL, F.REL)){
                        F.parents = F.parents,l;
                        NewRel[1..m,l] = updateRel(F.REL, submat(NewRel,1..m,l));
                    }
                }

                L2 = insert(L2,F,size(L2));

            // means that the flat is parallel
            } else{
                L[j].REL[i] = -1;
                NewRel[i,j] = -1;
                NFLAT[i] =  -1;
            }

        }
        }
    }

    toc = timer;
    dbprint(3-voice,"rank "+string(rk)+": found "+ string(counter) +
            " flats in "+string(toc - tic)+"s");
    counter  = 0;
    P.r[rk] = L2;
}

dbprint(3-voice,newline);
//dbprint(3-voice,"Computation time: "+string(timer - start)+"s");
dbprint(3-voice,"Matrix tests: "+string(ntests));
return (P);
}

example
{
"EXAMPLE: Intersection lattice of the braid arrangement in 3 dimensions "; echo = 2;
ring r;
arrLattice(arrTypeB(3));
}


static proc arrRank(arr A){
    int n = rank(matrix(A));
    int m = nvars(A);
    if(n < m) {return (n);}
    return (m);
}

static proc getOldRel(list L, int m){
    int n = size(L);
    matrix NewRel[m][n];

    for(int i=1; i<=n; i++){
        NewRel[1..m,i] = L[i].REL;
    }

    return (NewRel);
}

static proc updateRel(intvec REL, matrix NewRel){
    for(int i=1; i<=size(REL); i++){
        if(NewRel[i,1] == 0 && REL[i] == 1){NewRel[i,1] = 1;}
    }

    return (NewRel);
}

// Flat F is a parent of Flat G if F[i] == 1 => G[i] == 1 for all i
static proc isParent(intvec parent, intvec child){
    int counter = 0;
    for(int i=1; i<= size(parent); i++){
        if(parent[i]==1){
            if(child[i] == 1){counter++;}
            else{ return (0); }
        }
    }
    if(counter == 0){ return (0); }
    return (1);
}

// transforms the "rel" intvec into an intvec which contains the indices of the hyperplanes.
static proc collectHplanes(intvec RELATIONS){
    intvec result;

    for(int i=1; i<=size(RELATIONS); i++){
        if(RELATIONS[i] > 0){result = result,i;}
    }

    result = result[2..size(result)];

    return (result);

}

static proc getFlat(arrposet P, int i, int j){
    list L = P.r[i];
    return (L[j]);
}

static proc setFlat(arrposet P, int i, int j, flat F){
    list L = P.r[i];
    L[j] = F;
    P.r[i] = L;
    return (P);
}

static proc getFlag(arrposet P, int i, int j){
    flat F = getFlat(P,i,j);
    return (F.flag);
}

static proc setFlag(arrposet P, int i, int j, int flag){
    flat F = getFlat(P,i,j);
    F.flag = flag;
    return (setFlat(P,i,j,F));
}

//============================================================================//
//------------------------- #16 INTERSECTION-LATTICE -------------------------//
//============================================================================//
//============================================================================//
//-------------------------- MOEBIUS    -------------------------------------//
//============================================================================//


// calculates the moebius values of the poset
proc moebius(arrposet P)
"USAGE:     moebius(arrposet P)
RETURN:     [arrposet] fills in the moebius values of the flats in the poset
SEE ALSO:   moebius
KEYWORDS:   moebius function
EXAMPLE:    example arrCharPoly; shows an example"

{
    list L;
    int i,j;

    for(i=2; i<=rank(P.A); i++){
        L = P.r[i];

        for(j=1; j<=size(L); j++){
            arrposet Q = P;
            export Q;
            L[j].moebius = -(1 + moebiusRecursion(i, j));
            kill Q;
        }
        P.r[i] = L;

    }

    return (P);
}

example
{
"EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2;
ring R = 0,(x,y,z,t),dp;
arr A = arrBraid(4);
arrposet P = arrLattice(A);
P;
//As you can see the values are not calculated yet:

printMoebius(P);
P = moebius(P);

//Now all entries are initialized:

printMoebius(P);
}


// moebiusrecursion(i,j) is the sum moebius(X) over {X | V>X>Flat_ij}
// Vss: all moebius values of rank < i are known and correctly saved in P
static proc moebiusRecursion(int i, int j){
    flat F = getFlat(Q, i, j);
    int m = F.moebius;

    if(getFlag(Q,i,j) == 1){return (0);}
    else {
        Q = setFlag(Q,i,j,1);
        if(i > 1){
            for(int k=1; k<=size(F.parents); k++){
                m = m + moebiusRecursion(i-1, F.parents[k]);
            }
        }
    }

    return (m);
}


// pi(A,t) = sum[X in L(A); moebius(X)*(-t)^r(X)]
proc arrPoincare(def input)
"USAGE:     arrPoincare(A); arr A
RETURN:     [intvec] The Poincare polynomial as integer vector of the arrangement, which
            is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   Poincare polynomial
EXAMPLE:    example arrPoincare;"

{
    while(1){
        if(typeof(input) == "arr"){
            arr A = input;
            arrposet P = arrLattice(A);
            break;
        }
        if(typeof(input) == "arrposet"){
            arr A = input.A;
            arrposet P = input;
            break;
        }
        ERROR("Bad input type");
    }

    P = moebius(P);
    int i, j, coeff, sign;
    list L;

    intvec v = 1;
    sign = 1;
    for(i=1; i<= rank(A); i++){
        L = P.r[i];
        sign = (-1)*sign;
        coeff = 0;
        for(j=1;j<=size(L);j++){
            coeff = coeff + L[j].moebius;
        }
        coeff = sign*coeff;
        v = v,coeff;
    }
    return (v);
}

example
{
"EXAMPLE: The poincare polynomial of the braid arrangement in k dimensions is given as: " +
"pi(A,t) = (1 + t)*...*(1 + (k-1)*t)"; echo = 2;
ring R = 0,(x,y,z,u,v),dp;
arr A = arrBraid(5);
intvec v = arrPoincare(A);
(1+x)*(1+2x)*(1+3x)*(1+4x);
v;
}


// X(A,t) = t^l * pi(A,-t^-1)
proc arrCharPoly(def input)
"USAGE:     arrCharPoly(arr A)
RETURN:     [intvec] coefficients of the characteristic polynomial of A in incresing order
REMARKS:    The algorithm only returns the coefficients of the characteristic polynomial since they
            are whole numbers but the basering could be something different.
SEE ALSO:   arrCharPoly, arrPoincare
KEYWORDS:   characteristic polynomial
EXAMPLE:    example arrCharPoly; shows an example"

{
    intvec v = arrPoincare(input);
    int l = size(v);
    v = v[l..1];            //reverses order
    for(int i=2; i<=l; i=i+2){
        v[i] = -v[i];
    }

    return (v);
}

example
{
"EXAMPLE: The characteristic polynomial of the Braid arrangement in k dimensions is given as: " +
"X(A,t) = t*(t-1)*...*(t-(k-1)))"; echo = 2;
ring R = 0,(x,y,z,u,v),dp;
arr A = arrBraid(5);
intvec v = arrCharPoly(A);
x*(x-1)*(x-2)*(x-3)*(x-4);
v;
}


// number of chambers of the arrangement
proc arrChambers(arr A)
"USAGE:     arrChambers(A); arr A
RETURN:     [int] The number of chambers of an arrangement, which is equal to the
            evaluation of the Poincare polynomial at 1.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   chambers
EXAMPLE:    example arrChambers;"

{
    intvec v = arrPoincare(A);
    return(sum(v));
}

example
{
"EXAMPLE: Computing the number of number of chambers in a simple arrangement"; echo = 2;
ring R = 0,(x,y),dp;
arr A = ideal(x,y,x+y-1);
arrChambers(A);
}


// number of bounded chambers of the arrangement
proc arrBoundedChambers(arr A)
"USAGE:     arrBoundedChambers(A); arr A
RETURN:     [int] The number of bounded chambers of an arrangement, which is equal to
            the evaluation of the Poincare polynomial at -1.
SEE ALSO:   arrPoincare, arrChambers, arrBoundedChambers
KEYWORDS:   chambers
EXAMPLE:    example arrBoundedChambers;"

{
    intvec v = arrPoincare(A);
    for(int i=2; i<=size(v); i=i+2){
        v[i] = -v[i];
    }
    return (sum(v));
}

example
{
"EXAMPLE: Computing the number of bounded chambers in a simple arrangement"; echo = 2;
ring R = 0,(x,y),dp;
arr A = ideal(x,y,x+y-1);
arrBoundedChambers(A);
}

//============================================================================//
//------------------------------- OUTPUT  ------------------------------------//
//============================================================================//

static proc arrPrintPoset(arrposet P){
    print("Given Arrangement:");
    P.A;
    print("Corresponding poset:");

    string s;
    int i,j;
    list L;
    for(i=1; i<= size(P.r); i++){
        print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======");
        s = "";
        L = P.r[i];
        for(j=1; j<=size(L); j++){
            s = s + " (" + string(collectHplanes(L[j].REL)) + "), ";
        }

        print(s);
        if(size(P.r[i]) == 0){break;}
    }
}


proc printMoebius(arrposet P)
"USAGE:     printMoebius(A); arr A
RETURN:     [] displays the moebius values of all the flats in the poset
REMARKS:    Mainly used for debugging.
EXAMPLE:    example printMoebius;"

{
    print("Moebius values: ");

    string s;
    int i,j;
    list L;
    for(i=1; i<= size(P.r); i++){
        print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======");
        s = "";
        L = P.r[i];
        for(j=1; j<=size(L); j++){
            s = s + " (" + string(L[j].moebius) + "), ";
        }

        print(s);
        if(size(P.r[i]) == 0){break;}
    }
}

example
{
"EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2;
ring R = 0,(x,y,z,t),dp;
arr A = arrBraid(4);
arrposet P = arrLattice(A);
P;
//As you can see the values are not calculated yet:

printMoebius(P);
P = moebius(P);

//Now all entries are initialized:

printMoebius(P);
}

//============================================================================//
//-------------------------------- arrFlats Stuff --------------------------------//
//============================================================================//
// test via hyperplanes.
proc arrFlats(arr ARR)
"USAGE:     size(A); A arr
RETURN:     [arrposet] Intersection lattice
SEE ALSO:   arrFlats
KEYWORDS:   intersection lattice
EXAMPLE:    example arrFlats;"

{
    print(newline);
    print("=== Computing poset ===");
    print(newline);

    //Performance test
    system("--ticks-per-sec",1000);
    timer = 0;
    int time = timer;

    // Initialization
    list L,L2;
    intvec u,v,w,src;
    int i,j,k,l,m,n,d, rk,rk2;
    int counter, ntests; ntests = 0;
    matrix S,T;
    intmat REL;

    // Initialization of matrix
    matrix M =  matrix(ARR);                  // m x n+1 matrix
    m = size(ARR);
    n = nvars(basering);

    // Initialization of P
    arrflats P;         P.A = ARR;
    for(i=1; i<=n; i++){ P.r = insert(P.r,list());}
    for(i=1; i<=m; i++){ L[i] = i;}
    P.r[1] = L;

    // calculating r[i] step by step
    for(i=2; i<=n; i++){                    // maximal flat possible has rank A.v ofc

        counter = 0;

        L = P.r[i-1];                   // flats from the last iteration give rise to current ones.
        if(size(L) <= 1){break;}            // finish early if no more intersections possible.
        L2 = list();                        // new list for elts of rk=i


        // REL is an intmat that contains information about the relations between
        // the i-1 flats and the Hyperplanes.
        // REL[i,j] =  1 means that H_i intersects F_j
        // REL[i,j] =  0 means that it has not been tested yet
        // REL[i,j] = -1 means that H_i does NOT intersect F_j, i.e. they are parallel
        // resetting the REL-matrix
        REL = intmat(intvec(0),m,size(L));
        // filles in the trivial entries, i.e. each flat intersects all hplanes it is composed of.
        for(j=1; j<=size(L); j++){
            u = L[j];
            for(k=1; k<= size(u); k++){
                REL[u[k],j] = 1;
            }
        }

        //find new flat of rank i
        for(j=1; j<=size(L); j++){           // goes over the elts of L.r[i-1]

            // Looking for intersetion of F_j with hyperplanes
            u = L[j];

            for(k=1; k<=m; k++){
            if (REL[k,j] == 0 ){
                v = insertVal(u,k);
                counter++;
                ntests = ntests +2;
                rk = rank(submat(M,v,1..n));
                if(rk == rank(submat(M,v,1..n+1))){
                    //INTERSECTION FOUND
                    // potentially there exist more hyperplanes which intersect in the same
                    // space. For example (x,y,x+y) all intersect in (0,0)

                    REL[k,j] = 1;

                    for(l=k+1; l<=m; l++){
                    if(REL[l,j] == 0){
                        u = insertVal(v,l);
                        ntests = ntests +2;
                        rk2 = rank(submat(M, u, 1..n));

                        if( rk == rk2 && rk == rank(submat(M,u,1..n+1)) ){
                            v = u;
                            REL[l,j] = 1;
                        }
                    }}

                // Add new flat to list, reset u.
                L2 = insert(L2,v,size(L2));
                u = L[j];                   //if k = m this is not necessary....
                }
            }}
        }

        print("rank: "+string(i)+", expected OPS: "+string(binomial(size(L),2))
                +", excecuted OPS: " + string(counter));
        counter  = 0;

        //cleanup => doing this during computation increases speed!!
        for(j=1; j<size(L2); j++){
            u = L2[j];
            for(k=j+1; k<= size(L2); k++){
                if(u == L2[k]){
                    L2  = delete(L2,k); k--; counter++;

                }
            }
        }

        P.r[i] = L2;
        print("Cleaned up: "+string(counter)+" hyperplanes");
        print(newline);
    }

    print(newline);
    //print("Computation time: "+string(timer - time));
    print("Matrix tests: "+string(ntests));

    return (P);
}

example
{
"EXAMPLE: "; echo = 2;
ring R = 0,(x(1..5)),dp;
arrFlats(arrBraid(5));
}

// inserst k in u in the correct lexicographical place
static proc insertVal(intvec u, int k){

    if(k<u[1]){
        u = k,u;
        return (u);
    }
    if(k>u[size(u)]){
        u = u,k;
        return (u);
    }

    int i = 2;
    while(u[i]<k){i++;}

    u = u[1..(i-1)],k,u[i..size(u)];


    return(u);

}

static proc arrPrintFlats(arrflats P){
    print("Given Arrangement:");
    P.A;
    print("Corresponding poset:");

    string s;
    int i,j;
    list L;
    for(i=1; i<= size(P.r); i++){
        print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======");
        s = "";
        L = P.r[i];
        for(j=1; j<=size(L); j++){
            s = s + " (" + string(L[j]) + "), ";
        }

        print(s);
    }
}


// Compares two posets
static proc cposet(arrflats P, arrflats Q){

list L = P.r;
list K = Q.r;
list L2, K2;
int i,j, isequal;

for(i =1; i<=size(L); i++){
    L2 = L[i];
    K2 = K[i];
    isequal = 1;
    for(j=1; j<= size(L2); j++){
        if(L2[j] != K2[j]){isequal = 0; break;}
    }
    print("@rank "+string(i)+" |P.r[i]| = "+string(size(L[i]))+", |Q.r[i]| = "+string(size(K[i])));
    print("Elements are the same: "+string(isequal));
}

}

//============================================================================//
//-------------------------------- END OF CODE -------------------------------//
//============================================================================//

//============================================================================//
//------------------------------- UNUSED PROCEDURES --------------------------//
//============================================================================//


static proc printParents(arrposet P)
"USAGE:     printParents(A); arr A
RETURN:     [] displays the parent-lists of all the flats in the poset
REMARKS:    Mainly used for debugging.
SEE ALSO:   printParents, printMoebius, printFlags
EXAMPLE:    example printParents;"

{
    print("Given Arrangement:");
    P.A;
    print("Corresponding poset:");

    string s;
    int i,j;
    list L;
    for(i=1; i<= size(P.r); i++){
        print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======");
        s = "";
        L = P.r[i];
        for(j=1; j<=size(L); j++){
            s = s + " (" + string(L[j].parents) + "), ";
        }

        print(s);
        if(size(P.r[i]) == 0){break;}
    }
}

static proc printFlags(arrposet P)
"USAGE:     printFlags(A); arr A
RETURN:     [] displays the flag bits of all the flats in the poset
REMARKS:    Mainly used for debugging.
SEE ALSO:   printParents, printMoebius, printFlags
EXAMPLE:    example printParents;"
{
    print("Flag values:");

    string s;
    int i,j;
    list L;
    for(i=1; i<= size(P.r); i++){
        print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======");
        s = "";
        L = P.r[i];
        for(j=1; j<=size(L); j++){
            s = s + " (" + string(L[j].flag) + "), ";
        }

        print(s);
        if(size(P.r[i]) == 0){break;}
    }
}

// expects 2 intvecs of stric ascending order
// returns merged intvec
static proc mergeIV(intvec u, intvec v){
    intvec res;
    int k;
    int m = 1;
    int n = 1;
    for(k=1; m<=size(u) && n<=size(v); k++){
        while(1){
            if(u[m] <  v[n]){res[k] = u[m]; m++;      break;}
            if(u[m] >  v[n]){res[k] = v[n];      n++; break;}
            if(u[m] == v[n]){res[k] = u[m]; m++; n++; break;}
            ERROR("Something went wrong in proc mergeIV");
        }
    }
    while(m <= size(u)){res[k] = u[m]; m++; k++;}
    while(n <= size(v)){res[k] = v[n]; n++; k++;}
    return (res);
}

// checks if v[i] < u[i] for all i
static proc isSmaller(intvec u, intvec v){
    for(int i= 1; i<=size(u); i++){
        if(u[i] > v[i]){
            return (0);
        }
    }

    return (1);
}

static proc new2old(arrposet P){
    arrflats Q;
    Q.A = P.A;
    int i,j;
    list L = P.r;
    for(i=1; i<=size(L); i++){
        list L2 = L[i];
        for(j=1; j<=size(L2); j++){
            L2[j] = collectHplanes(L2[j].REL);
        }
        L[i] = L2;
    }

    Q.r = L;

    return (Q);
}
/*********************************************************************
newstruct("arr","ideal l");

Defines a hyperplane arrangement by a list of linear polynomials such that the hyperplanes are the
varieties of those polynomials.

supported operators:
A = input               defines an arrangement by the input which may consist of the types arr/ideal/poly/matrix/list
A + input               adds arrangement defined by the right hand side to the left hand side
A - input               deletes hyperplanes from the arrangement
<, >, <=, >=, ==, !=    set theoretical comparison of two arrangements
A[int]                  access to a single hyperplane
A[intvec]               subarrangement defined by the indexed hyperplanes
A;                      prints the arrangement on the screen

supported type conversions:
matrix(A)               returns coefficient matrix of the defining polynomials
poly(A)                 returns the defining polynomial which is the product of all the single polynomials
list(A)                 returns the defining polynomials as a list
ideal(A)                returns the defining polynomials as an ideal
type2arr(input)         returns arrangement defined by the input which may consist of the types arr/ideal/poly/matrix/list

supported inherited functions:
delete(A, intvec)       deletes indexed hyperplanes from the arrangement
homog(A)                checks wether the defining polys are all homogenous <=> arr is central
homog(A, rvar)          homogenizes the def. polys with respect to the given ring variable
size(A)                 number of hyperplanes in the arrangement
subst(A,rvar,poly,...)  substitutes ringvariables with given polynomials, though they need to remain linear
variables(A)            ideal generated by the ring variables that A depends on
nvars(A)                number of ring variables that A depends on



newstruct("multarr","ideal l, intvec m");

Defines a hyperplane arrangement with multiplicities by a list of linear polynomials such that the hyperplanes are the
varieties of those polynomials and an intvec in which the multiplicities are saved.

supported operators:
M = input               defines an multarr by the input which may consist of the types multarr/ideal/poly/list
M + input               adds arrangement defined by the right hand side to the left hand side
M;                      prints the multarr on the screen

supported type conversions:
poly(M)                 returns defining polynomial with multiplicities
arr2multarr(A,intvec)   returns multarr with multiplicities defined by the intvec.
multarr2arr(M)          returns the internal arrangement (all multiplicities set to 1)
delete(M,int)           decrements the multiplicity of the hyperplane defined by the index by 1
size(M)                 returns number of hyperplanes counting multiplicities.
*-----------------------------------------------------------------------*/