// Last change 12.02.2001 (Eric Westenberger) /////////////////////////////////////////////////////////////////////////////// version="$Id: zeroset.lib,v 1.19 2009-04-06 09:17:01 seelisch Exp $"; category="Symbolic-numerical solving"; info=" LIBRARY: zeroset.lib Procedures for roots and factorization AUTHOR: Thomas Bayer, email: tbayer@mathematik.uni-kl.de,@* http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/@* Current address: Hochschule Ravensburg-Weingarten OVERVIEW: Algorithms for finding the zero-set of a zero-dim. ideal in Q(a)[x_1,..,x_n], roots and factorization of univariate polynomials over Q(a)[t] where a is an algebraic number. Written in the scope of the diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli spaces of semiquasihomogeneous singularities and an implementation in Singular'. This library is meant as a preliminary extension of the functionality of Singular for univariate factorization of polynomials over simple algebraic extensions in characteristic 0. NOTE: Subprocedures with postfix 'Main' require that the ring contains a variable 'a' and no parameters, and the ideal 'mpoly', where 'minpoly' from the basering is stored. PROCEDURES: quotient(f, g) quotient q of f w.r.t. g (in f = q*g + remainder) remainder(f,g) remainder of the division of f by g roots(f) computes all roots of f in an extension field of Q sqfrNorm(f) norm of f (f must be squarefree) zeroSet(I) zero-set of the 0-dim. ideal I AUXILIARY PROCEDURES: egcdMain(f, g) gcd over an algebraic extension field of Q factorMain(f) factorization of f over an algebraic extension field invertNumberMain(c) inverts an element of an algebraic extension field quotientMain(f, g) quotient of f w.r.t. g remainderMain(f,g) remainder of the division of f by g rootsMain(f) computes all roots of f, might extend the ground field sqfrNormMain(f) norm of f (f must be squarefree) containedQ(data, f) f in data ? sameQ(a, b) a == b (list a,b) "; LIB "primitiv.lib"; LIB "primdec.lib"; // note : return a ring : ring need not be exported !!! // Artihmetic in Q(a)[x] without built-in procedures // assume basering = Q[x,a] and minpoly is represented by mpoly(a). // the algorithms are taken from "Polynomial Algorithms in Computer Algebra", // F. Winkler, Springer Verlag Wien, 1996. // To do : // squarefree factorization // multiplicities // Improvement : // a main problem is the growth of the coefficients. Try roots(x7 - 1) // return ideal mpoly ! // mpoly is not monic, comes from primitive_extra // IMPLEMENTATION // // In procedures with name 'proc-name'Main a polynomial ring over a simple // extension field is represented as Q[x...,a] together with the ideal // 'mpoly' (attribute "isSB"). The arithmetic in the extension field is // implemented in the procedures in the procedures 'MultPolys' (multiplication) // and 'InvertNumber' (inversion). After addition and substraction one should // apply 'SimplifyPoly' to the result to reduce the result w.r.t. 'mpoly'. // This is done by reducing each coefficient seperately, which is more // efficient for polynomials with many terms. /////////////////////////////////////////////////////////////////////////////// proc roots(poly f) "USAGE: roots(f); where f is a polynomial PURPOSE: compute all roots of f in a finite extension of the ground field without multiplicities. RETURN: ring, a polynomial ring over an extension field of the ground field, containing a list 'roots' and polynomials 'newA' and 'f': @format - 'roots' is the list of roots of the polynomial f (no multiplicities) - if the ground field is Q(a') and the extension field is Q(a), then 'newA' is the representation of a' in Q(a). If the basering contains a parameter 'a' and the minpoly remains unchanged then 'newA' = 'a'. If the basering does not contain a parameter then 'newA' = 'a' (default). - 'f' is the polynomial f in Q(a) (a' being substituted by 'newA') @end format ASSUME: ground field to be Q or a simple extension of Q given by a minpoly EXAMPLE: example roots; shows an example " { int dbPrt = printlevel-voice+3; // create a new ring where par(1) is replaced by the variable // with the same name or, if basering does not contain a parameter, // with a new variable 'a'. def ROB = basering; def ROR = TransferRing(basering); setring ROR; export(ROR); // get the polynomial f and find the roots poly f = imap(ROB, f); list result = rootsMain(f); // find roots of f // store the roots and the new representation of 'a' and transform // the coefficients of f. list roots = result[1]; poly newA = result[2]; map F = basering, maxideal(1); F[nvars(basering)] = newA; poly fn = SimplifyPoly(F(f)); // create a new ring with minploy = mpoly[1] (from ROR) def RON = NewBaseRing(); setring(RON); list roots = imap(ROR, roots); poly newA = imap(ROR, newA); poly f = imap(ROR, fn); kill ROR; export(roots); export(newA); export(f); dbprint(dbPrt," // 'roots' created a new ring which contains the list 'roots' and // the polynomials 'f' and 'newA' // To access the roots, newA and the new representation of f, type def R = roots(f); setring R; roots; newA; f; "); return(RON); } example {"EXAMPLE:"; echo = 2; ring R = (0,a), x, lp; minpoly = a2+1; poly f = x3 - a; def R1 = roots(f); setring R1; minpoly; newA; f; roots; map F; F[1] = roots[1]; F(f); } /////////////////////////////////////////////////////////////////////////////// proc rootsMain(poly f) "USAGE: rootsMain(f); where f is a polynomial PURPOSE: compute all roots of f in a finite extension of the ground field without multiplicities. RETURN: list, all entries are polynomials @format _[1] = roots of f, each entry is a polynomial _[2] = 'newA' - if the ground field is Q(b) and the extension field is Q(a), then 'newA' is the representation of b in Q(a) _[3] = minpoly of the algebraic extension of the ground field @end format ASSUME: basering = Q[x,a] ideal mpoly must be defined, it might be 0! NOTE: might change the ideal mpoly!! EXAMPLE: example roots; shows an example " { int i, linFactors, nlinFactors, dbPrt; intvec wt = 1,0; // deg(a) = 0 list factorList, nlFactors, nlMult, roots, result; poly fa, lc; dbPrt = printlevel-voice+3; // factor f in Q(a)[t] to obtain the roots lying in Q(a) // firstly, find roots of the linear factors, // nonlinear factors are processed later dbprint(dbPrt, "roots of " + string(f) + ", minimal polynomial = " + string(mpoly[1])); factorList = factorMain(f); // Factorize f dbprint(dbPrt, (" prime factors of f are : " + string(factorList[1]))); linFactors = 0; nlinFactors = 0; for(i = 2; i <= size(factorList[1]); i = i + 1) { // find linear and nonlinear factors fa = factorList[1][i]; if(deg(fa, wt) == 1) { linFactors++; // get the root from the linear factor lc = LeadTerm(fa, 1)[3]; fa = MultPolys(invertNumberMain(lc), fa); // make factor monic roots[linFactors] = var(1) - fa; // fa is monic !! } else { // ignore nonlinear factors nlinFactors++; nlFactors[nlinFactors] = factorList[1][i]; nlMult[nlinFactors] = factorList[2][i]; } } if(linFactors == size(factorList[1]) - 1) { // all roots of f are contained in the ground field result[1] = roots; result[2] = var(2); result[3] = mpoly[1]; return(result); } // process the nonlinear factors, i.e., extend the ground field // where a nonlinear factor (irreducible) is a minimal polynomial // compute the primitive element of this extension ideal primElem, minPolys, Fid; list partSol; map F, Xchange; poly f1, newA, mp, oldMinPoly; Fid = mpoly; F[1] = var(1); Xchange[1] = var(2); // the variables have to be exchanged Xchange[2] = var(1); // for the use of 'primitive' if(nlinFactors == 1) { // one nl factor // compute the roots of the nonlinear (irreducible, monic) factor f1 of f // by extending the basefield by a' with minimal polynomial f1 // Then call roots(f1) to find the roots of f1 over the new base field f1 = nlFactors[1]; if(mpoly[1] != 0) { mp = mpoly[1]; minPolys = Xchange(mp), Xchange(f1); primElem = primitive_extra(minPolys); // no random coord. change mpoly = std(primElem[1]); F = basering, maxideal(1); F[2] = primElem[2]; // transfer all to the new representation newA = primElem[2]; // new representation of a f1 = SimplifyPoly(F(f1)); //reduce(F(f1), mpoly); if(size(roots) > 0) {roots = SimplifyData(F(roots));} } else { mpoly = std(Xchange(f1)); newA = var(2); } result[3] = mpoly[1]; oldMinPoly = mpoly[1]; partSol = rootsMain(f1); // find roots of f1 over extended field if(oldMinPoly != partSol[3]) { // minpoly has changed ? // all previously computed roots must be transformed // because the minpoly has changed result[3] = partSol[3]; // new minpoly F[2] = partSol[2]; // new representation of algebraic number if(size(roots) > 0) {roots = SimplifyData(F(roots)); } newA = SimplifyPoly(F(newA)); // F(newA); } roots = roots + partSol[1]; // add roots result[2] = newA; result[1] = roots; } else { // more than one nonlinear (irreducible) factor (f_1,...,f_r) // solve each of them by rootsMain(f_i), append their roots // change the minpoly and transform all previously computed // roots if necessary. // Note that the for-loop is more or less book-keeping newA = var(2); result[2] = newA; for(i = 1; i <= size(nlFactors); i = i + 1) { oldMinPoly = mpoly[1]; partSol = rootsMain(nlFactors[i]); // main work nlFactors[i] = 0; // delete factor result[3] = partSol[3]; // store minpoly // book-keeping starts here as in the case 1 nonlinear factor if(oldMinPoly != partSol[3]) { // minpoly has changed F = basering, maxideal(1); F[2] = partSol[2]; // transfer all to the new representation newA = SimplifyPoly(F(newA)); // F(newA); new representation of a result[2] = newA; if(i < size(nlFactors)) { nlFactors = SimplifyData(F(nlFactors)); } // transform remaining factors if(size(roots) > 0) {roots = SimplifyData(F(roots));} } roots = roots + partSol[1]; // transform roots result[1] = roots; } // end more than one nl factor } return(result); } /////////////////////////////////////////////////////////////////////////////// proc zeroSet(ideal I, list #) "USAGE: zeroSet(I [,opt] ); I=ideal, opt=integer PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension of the ground field. RETURN: ring, a polynomial ring over an extension field of the ground field, containing a list 'zeroset', a polynomial 'newA', and an ideal 'id': @format - 'zeroset' is the list of the zeros of the ideal I, each zero is an ideal. - if the ground field is Q(b) and the extension field is Q(a), then 'newA' is the representation of b in Q(a). If the basering contains a parameter 'a' and the minpoly remains unchanged then 'newA' = 'a'. If the basering does not contain a parameter then 'newA' = 'a' (default). - 'id' is the ideal I in Q(a)[x_1,...] (a' substituted by 'newA') @end format ASSUME: dim(I) = 0, and ground field to be Q or a simple extension of Q given by a minpoly. OPTIONS: opt = 0: no primary decomposition (default) opt > 0: primary decomposition NOTE: If I contains an algebraic number (parameter) then I must be transformed w.r.t. 'newA' in the new ring. EXAMPLE: example zeroSet; shows an example " { int primaryDecQ, dbPrt; list rp; dbPrt = printlevel-voice+2; if(size(#) > 0) { primaryDecQ = #[1]; } else { primaryDecQ = 0; } // create a new ring 'ZSR' with one additional variable instead of the // parameter // if the basering does not contain a parameter then 'a' is used as the // additional variable. def RZSB = basering; def ZSR = TransferRing(RZSB); setring ZSR; // get ideal I and find the zero-set ideal id = std(imap(RZSB, I)); // print(dim(id)); if(dim(id) > 1) { // new variable adjoined to ZSR ERROR(" ideal not zerodimensional "); } list result = zeroSetMain(id, primaryDecQ); // store the zero-set, minimal polynomial and the new representative of 'a' list zeroset = result[1]; poly newA = result[2]; poly minPoly = result[3][1]; // transform the generators of the ideal I w.r.t. the new representation // of 'a' map F = basering, maxideal(1); F[nvars(basering)] = newA; id = SimplifyData(F(id)); // create a new ring with minpoly = minPoly def RZBN = NewBaseRing(); setring RZBN; list zeroset = imap(ZSR, zeroset); poly newA = imap(ZSR, newA); ideal id = imap(ZSR, id); kill ZSR; export(id); export(zeroset); export(newA); dbprint(dbPrt," // 'zeroSet' created a new ring which contains the list 'zeroset', the ideal // 'id' and the polynomial 'newA'. 'id' is the ideal of the input transformed // w.r.t. 'newA'. // To access the zero-set, 'newA' and the new representation of the ideal, type def R = zeroSet(I); setring R; zeroset; newA; id; "); setring RZSB; return(RZBN); } example {"EXAMPLE:"; echo = 2; ring R = (0,a), (x,y,z), lp; minpoly = a2 + 1; ideal I = x2 - 1/2, a*z - 1, y - 2; def T = zeroSet(I); setring T; minpoly; newA; id; zeroset; map F1 = basering, zeroset[1]; map F2 = basering, zeroset[2]; F1(id); F2(id); } /////////////////////////////////////////////////////////////////////////////// proc invertNumberMain(poly f) "USAGE: invertNumberMain(f); where f is a polynomial PURPOSE: compute 1/f if f is a number in Q(a), i.e., f is represented by a polynomial in Q[a]. RETURN: poly 1/f ASSUME: basering = Q[x_1,...,x_n,a], ideal mpoly must be defined and != 0 ! NOTE: outdated, use / instead " { if(diff(f, var(1)) != 0) { ERROR("number must not contain variable !");} int n = nvars(basering); def RINB = basering; string ringSTR = "ring RINR = 0, " + string(var(n)) + ", dp;"; execute(ringSTR); // new ring = Q[a] list gcdList; poly f, g, inv; f = imap(RINB, f); g = imap(RINB, mpoly)[1]; if(diff(f, var(1)) != 0) { inv = extgcd(f, g)[2]; } // f contains var(1) else { inv = 1/f;} // f element in Q setring(RINB); return(imap(RINR, inv)); } /////////////////////////////////////////////////////////////////////////////// proc MultPolys(poly f, poly g) "USAGE: MultPolys(f, g); poly f,g PURPOSE: multiply the polynomials f and g and reduce them w.r.t. mpoly RETURN: poly f*g ASSUME: basering = Q[x,a], ideal mpoly must be defined, it might be 0 ! " { return(SimplifyPoly(f * g)); } /////////////////////////////////////////////////////////////////////////////// proc LeadTerm(poly f, int i) "USAGE: LeadTerm(f); poly f, int i PURPOSE: compute the leading coef and term of f w.r.t var(i), where the last ring variable is treated as a parameter. RETURN: list of polynomials _[1] = leading term _[2] = leading monomial _[3] = leading coefficient ASSUME: basering = Q[x_1,...,x_n,a] " { list result; matrix co = coef(f, var(i)); result[1] = co[1, 1]*co[2, 1]; result[2] = co[1, 1]; result[3] = co[2, 1]; return(result); } /////////////////////////////////////////////////////////////////////////////// proc quotient(poly f, poly g) "USAGE: quotient(f, g); where f,g are polynomials; PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g) RETURN: list of polynomials @format _[1] = quotient q _[2] = remainder r @end format ASSUME: basering = Q[x] or Q(a)[x] NOTE: outdated, use div/mod instead EXAMPLE: example quotient; shows an example " { def QUOB = basering; def QUOR = TransferRing(basering); // new ring with parameter 'a' replaced by a variable setring QUOR; export(QUOR); poly f = imap(QUOB, f); poly g = imap(QUOB, g); list result = quotientMain(f, g); setring(QUOB); list result = imap(QUOR, result); kill QUOR; return(result); } example {"EXAMPLE:"; echo = 2; ring R = (0,a), x, lp; minpoly = a2+1; poly f = x4 - 2; poly g = x - a; list qr = quotient(f, g); qr; qr[1]*g + qr[2] - f; } proc quotientMain(poly f, poly g) "USAGE: quotientMain(f, g); where f,g are polynomials PURPOSE: compute the quotient q and remainder r s.th. f = g*q + r, deg(r) < deg(g) RETURN: list of polynomials @format _[1] = quotient q _[2] = remainder r @end format ASSUME: basering = Q[x,a] and ideal mpoly is defined (it might be 0), this represents the ring Q(a)[x] together with its minimal polynomial. NOTE: outdated, use div/mod instead EXAMPLE: example quotient; shows an example " { if(g == 0) { ERROR("Division by zero !");} def QMB = basering; def QMR = NewBaseRing(); setring QMR; poly f, g, h; h = imap(QMB, f) / imap(QMB, g); setring QMB; return(list(imap(QMR, h), 0)); } /////////////////////////////////////////////////////////////////////////////// proc remainder(poly f, poly g) "USAGE: remainder(f, g); where f,g are polynomials PURPOSE: compute the remainder of the division of f by g, i.e. a polynomial r s.t. f = g*q + r, deg(r) < deg(g). RETURN: poly ASSUME: basering = Q[x] or Q(a)[x] NOTE: outdated, use mod/reduce instead " { def REMB = basering; def REMR = TransferRing(basering); // new ring with parameter 'a' replaced by a variable setring(REMR); export(REMR); poly f = imap(REMB, f); poly g = imap(REMB, g); poly h = remainderMain(f, g); setring(REMB); poly r = imap(REMR, h); kill REMR; return(r); } example {"EXAMPLE:"; echo = 2; ring R = (0,a), x, lp; minpoly = a2+1; poly f = x4 - 1; poly g = x3 - 1; remainder(f, g); } proc remainderMain(poly f, poly g) "USAGE: remainderMain(f, g); where f,g are polynomials PURPOSE: compute the remainder r s.t. f = g*q + r, deg(r) < deg(g) RETURN: poly ASSUME: basering = Q[x,a] and ideal mpoly is defined (it might be 0), this represents the ring Q(a)[x] together with its minimal polynomial. NOTE: outdated, use mod/reduce instead " { int dg; intvec wt = 1,0;; poly lc, g1, r; if(deg(g, wt) == 0) { return(0); } lc = LeadTerm(g, 1)[3]; g1 = MultPolys(invertNumberMain(lc), g); // make g monic return(SimplifyPoly(reduce(f, std(g1)))); } /////////////////////////////////////////////////////////////////////////////// proc egcdMain(poly f, poly g) "USAGE: egcdMain(f, g); where f,g are polynomials PURPOSE: compute the polynomial gcd of f and g over Q(a)[x] RETURN: poly ASSUME: basering = Q[x,a] and ideal mpoly is defined (it might be 0), this represents the ring Q(a)[x] together with its minimal polynomial. NOTE: outdated, use gcd instead EXAMPLE: example EGCD; shows an example " { // might be extended to return s1, s2 s.t. f*s1 + g*s2 = gcd int i = 1; poly r1, r2, r; r1 = f; r2 = g; while(r2 != 0) { r = remainderMain(r1, r2); r1 = r2; r2 = r; } return(r1); } /////////////////////////////////////////////////////////////////////////////// proc MEGCD(poly f, poly g, int varIndex) "USAGE: MEGCD(f, g, i); poly f, g; int i PURPOSE: compute the polynomial gcd of f and g in the i'th variable RETURN: poly ASSUME: f, g are polynomials in var(i), last variable is the algebraic number EXAMPLE: example MEGCD; shows an example " // might be extended to return s1, s2 s.t. f*s1 + g*s2 = gc // not used ! { string @str, @sf, @sg, @mp, @parName; def @RGCDB = basering; @sf = string(f); @sg = string(g); @mp = string(minpoly); if(npars(basering) == 0) { @parName = "0";} else { @parName = "(0, " + parstr(basering) + ")"; } @str = "ring @RGCD = " + @parName + ", " + string(var(varIndex)) + ", dp;"; execute(@str); if(@mp != "0") { execute ("minpoly = " + @mp + ";"); } execute("poly @f = " + @sf + ";"); execute("poly @g = " + @sg + ";"); export(@RGCD); poly @h = gcd(@f, @g); setring(@RGCDB); poly h = imap(@RGCD, @h); kill @RGCD; return(h); } /////////////////////////////////////////////////////////////////////////////// proc sqfrNorm(poly f) "USAGE: sqfrNorm(f); where f is a polynomial PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x]. RETURN: list with 3 entries @format _[1] = squarefree norm of g (poly) _[2] = g (= f(x - s*a)) (poly) _[3] = s (int) @end format ASSUME: f must be squarefree, basering = Q(a)[x] and minpoly != 0. NOTE: the norm is an element of Q[x] EXAMPLE: example sqfrNorm; shows an example " { def SNB = basering; def SNR = TransferRing(SNB); // new ring with parameter 'a' // replaced by a variable setring SNR; poly f = imap(SNB, f); list result = sqfrNormMain(f); // squarefree norm of f setring SNB; list result = imap(SNR, result); kill SNR; return(result); } example {"EXAMPLE:"; echo = 2; ring R = (0,a), x, lp; minpoly = a2+1; poly f = x4 - 2*x + 1; sqfrNorm(f); } proc sqfrNormMain(poly f) "USAGE: sqfrNorm(f); where f is a polynomial PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x]. RETURN: list with 3 entries @format _[1] = squarefree norm of g (poly) _[2] = g (= f(x - s*a)) (poly) _[3] = s (int) @end format ASSUME: f must be squarefree, basering = Q[x,a] and ideal mpoly is equal to 'minpoly', this represents the ring Q(a)[x] together with 'minpoly'. NOTE: the norm is an element of Q[x] EXAMPLE: example SqfrNorm; shows an example " { def SNRMB = basering; int s = 0; intvec wt = 1,0; ideal mapId; // list result; poly g, N, N1, h; string ringSTR; mapId[1] = var(1) - var(2); // linear transformation mapId[2] = var(2); map Fs = SNRMB, mapId; N = resultant(f, mpoly[1], var(2)); // norm of f N1 = diff(N, var(1)); g = f; ringSTR = "ring SNRM1 = 0, " + string(var(1)) + ", dp;"; // univariate ring execute(ringSTR); poly N, N1, h; // N, N1 do not contain 'a', use built-in gcd h = gcd(imap(SNRMB, N), imap(SNRMB, N1)); setring(SNRMB); h = imap(SNRM1, h); while(deg(h, wt) != 0) { // while norm is not squarefree s = s + 1; g = reduce(Fs(g), mpoly); N = reduce(resultant(g, mpoly[1], var(2)), mpoly); // norm of g N1 = reduce(diff(N, var(1)), mpoly); setring(SNRM1); h = gcd(imap(SNRMB, N), imap(SNRMB, N1)); setring(SNRMB); h = imap(SNRM1, h); } return(list(N, g, s)); } /////////////////////////////////////////////////////////////////////////////// proc factorMain(poly f) "USAGE: factorMain(f); where f is a polynomial PURPOSE: compute the factorization of the squarefree poly f over Q(a)[t], minpoly = p(a). RETURN: list with 2 entries @format _[1] = factors, first is a constant _[2] = multiplicities (not yet implemented) @end format ASSUME: basering = Q[x,a], representing Q(a)[x]. An ideal mpoly must be defined, representing the minimal polynomial (it might be 0!). NOTE: outdated, use factorize instead EXAMPLE: example Factor; shows an example " { // extend this by a squarefree factorization !! // multiplicities are not valid !! int i, s; list normList, factorList, quo_rem; poly f1, h, h1, H, g, leadCoef, invCoeff; ideal fac1, fac2; map F; // if no minimal polynomial is defined then use 'factorize' // FactorOverQ is wrapped around 'factorize' if(mpoly[1] == 0) { // print(" factorize : deg = " + string(deg(f, intvec(1,0)))); factorList = factorize(f); // FactorOverQ(f); return(factorList); } // if mpoly != 0 and f does not contain the algebraic number, a root of // f might be contained in Q(a). Hence one must not use 'factorize'. fac1[1] = 1; fac2[1] = 1; normList = sqfrNormMain(f); // print(" factorize : deg = " + string(deg(normList[1], intvec(1,0)))); factorList = factorize(normList[1]); // factor squarefree norm of f over Q[x] g = normList[2]; s = normList[3]; F[1] = var(1) + s*var(2); // inverse transformation F[2] = var(2); fac1[1] = factorList[1][1]; fac2[1] = factorList[2][1]; for(i = 2; i <= size(factorList[1]); i = i + 1) { H = factorList[1][i]; h = egcdMain(H, g); quo_rem = quotientMain(g, h); g = quo_rem[1]; fac1[i] = SimplifyPoly(F(h)); fac2[i] = 1; // to be changed later } return(list(fac1, fac2)); } /////////////////////////////////////////////////////////////////////////////// proc zeroSetMain(ideal I, int primDecQ) "USAGE: zeroSetMain(ideal I, int opt); ideal I, int opt PURPOSE: compute the zero-set of the zero-dim. ideal I, in a simple extension of the ground field. RETURN: list - 'f' is the polynomial f in Q(a) (a' being substituted by newA) _[1] = zero-set (list), is the list of the zero-set of the ideal I, each entry is an ideal. _[2] = 'newA'; if the ground field is Q(a') and the extension field is Q(a), then 'newA' is the representation of a' in Q(a). If the basering contains a parameter 'a' and the minpoly remains unchanged then 'newA' = 'a'. If the basering does not contain a parameter then 'newA' = 'a' (default). _[3] = 'mpoly' (ideal), the minimal polynomial of the simple extension of the ground field. ASSUME: basering = K[x_1,x_2,...,x_n] where K = Q or a simple extension of Q given by a minpoly; dim(I) = 0. NOTE: opt = 0 no primary decomposition opt > 0 use a primary decomposition EXAMPLE: example zeroSet; shows an example " { // main work is done in zeroSetMainWork, here the zero-set of each ideal from the // primary decompostion is coputed by menas of zeroSetMainWork, and then the // minpoly and the parameter representing the algebraic extension are // transformed according to 'newA', i.e., only bookeeping is done. def altring=basering; int i, j, n, noMP, dbPrt; intvec w; list currentSol, result, idealList, primDecList, zeroSet; ideal J; map Fa; poly newA, oldMinPoly; dbPrt = printlevel-voice+2; dbprint(dbPrt, "zeroSet of " + string(I) + ", minpoly = " + string(minpoly)); n = nvars(basering) - 1; for(i = 1; i <= n; i++) { w[i] = 1;} w[n + 1] = 0; if(primDecQ == 0) { return(zeroSetMainWork(I, w, 0)); } newA = var(n + 1); if(mpoly[1] == 0) { noMP = 1;} else {noMP = 0;} primDecList = primdecGTZ(I); // primary decomposition dbprint(dbPrt, "primary decomposition consists of " + string(size(primDecList)) + " primary ideals "); // idealList = PDSort(idealList); // high degrees first for(i = 1; i <= size(primDecList); i = i + 1) { idealList[i] = primDecList[i][2]; // use prime component dbprint(dbPrt, string(i) + " " + string(idealList[i])); } // compute the zero-set of each primary ideal and join them. // If necessary, change the ground field and transform the zero-set dbprint(dbPrt, " find the zero-set of each primary ideal, form the union and keep track of the minimal polynomials "); for(i = 1; i <= size(idealList); i = i + 1) { J = idealList[i]; idealList[i] = 0; oldMinPoly = mpoly[1]; dbprint(dbPrt, " ideal#" + string(i) + " of " + string(size(idealList)) + " = " + string(J)); currentSol = zeroSetMainWork(J, w, 0); if(oldMinPoly != currentSol[3]) { // change minpoly and transform solutions dbprint(dbPrt, " change minpoly to " + string(currentSol[3][1])); dbprint(dbPrt, " new representation of algebraic number = " + string(currentSol[2])); if(!noMP) { // transform the algebraic number a Fa = basering, maxideal(1); Fa[n + 1] = currentSol[2]; newA = SimplifyPoly(Fa(newA)); // new representation of a if(size(zeroSet) > 0) {zeroSet = SimplifyZeroset(Fa(zeroSet)); } if(i < size(idealList)) { idealList = SimplifyZeroset(Fa(idealList)); } } else { noMP = 0;} } zeroSet = zeroSet + currentSol[1]; // add new elements } return(list(zeroSet, newA, mpoly)); } /////////////////////////////////////////////////////////////////////////////// proc zeroSetMainWork(ideal id, intvec wt, int sVars) "USAGE: zeroSetMainWork(I, wt, sVars); PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension of the ground field (without multiplicities). RETURN: list, all entries are polynomials _[1] = zeros, each entry is an ideal _[2] = newA; if the ground field is Q(a') this is the rep. of a' w.r.t. a _[3] = minpoly of the algebraic extension of the ground field (ideal) _[4] = name of algebraic number (default = 'a') ASSUME: basering = Q[x_1,x_2,...,x_n,a] ideal mpoly must be defined, it might be 0! NOTE: might change 'mpoly' !! EXAMPLE: example IdealSolve; shows an example " { def altring=basering; int i, j, k, nrSols, n, noMP; ideal I, generators, gens, solid, partsolid; list linSol, linearSolution, nLinSol, nonlinSolutions, partSol, sol, solutions, result; list linIndex, nlinIndex, index; map Fa, Fsubs; poly oldMinPoly, newA; if(mpoly[1] == 0) { noMP = 1;} else { noMP = 0;} n = nvars(basering) - 1; newA = var(n + 1); I = std(id); // find linear solutions of univariate generators linSol = LinearZeroSetMain(I, wt); generators = linSol[3]; // they are a standardbasis linIndex = linSol[2]; linearSolution = linSol[1]; if(size(linIndex) + sVars == n) { // all variables solved solid = SubsMapIdeal(linearSolution, linIndex, 0); result[1] = list(solid); result[2] = var(n + 1); result[3] = mpoly; return(result); } // find roots of the nonlinear univariate polynomials of generators // if necessary, transform linear solutions w.r.t. newA oldMinPoly = mpoly[1]; nLinSol = NonLinearZeroSetMain(generators, wt); // find solutions of univariate generators nonlinSolutions = nLinSol[1]; // store solutions nlinIndex = nLinSol[4]; // and index of solved variables generators = nLinSol[5]; // new generators // change minpoly if necessary and transform the ideal and the partial solutions if(oldMinPoly != nLinSol[3]) { newA = nLinSol[2]; if(!noMP && size(linearSolution) > 0) { // transform the algebraic number a Fa = basering, maxideal(1); Fa[n + 1] = newA; linearSolution = SimplifyData(Fa(linearSolution)); // ... } } // check if all variables are solved. if(size(linIndex) + size(nlinIndex) == n - sVars) { solutions = MergeSolutions(linearSolution, linIndex, nonlinSolutions, nlinIndex, list(), n); } else { // some variables are not solved. // substitute each partial solution in generators and find the // zero set of the resulting ideal by recursive application // of zeroSetMainWork ! index = linIndex + nlinIndex; nrSols = 0; for(i = 1; i <= size(nonlinSolutions); i = i + 1) { sol = linearSolution + nonlinSolutions[i]; solid = SubsMapIdeal(sol, index, 1); Fsubs = basering, solid; gens = std(SimplifyData(Fsubs(generators))); // substitute partial solution oldMinPoly = mpoly[1]; partSol = zeroSetMainWork(gens, wt, size(index) + sVars); if(oldMinPoly != partSol[3]) { // minpoly has changed Fa = basering, maxideal(1); Fa[n + 1] = partSol[2]; // a -> p(a), representation of a w.r.t. new minpoly newA = reduce(Fa(newA), mpoly); generators = std(SimplifyData(Fa(generators))); if(size(linearSolution) > 0) { linearSolution = SimplifyData(Fa(linearSolution));} if(size(nonlinSolutions) > 0) { nonlinSolutions = SimplifyZeroset(Fa(nonlinSolutions)); } sol = linearSolution + nonlinSolutions[i]; } for(j = 1; j <= size(partSol[1]); j++) { // for all partial solutions partsolid = partSol[1][j]; for(k = 1; k <= size(index); k++) { partsolid[index[k]] = sol[k]; } nrSols++; solutions[nrSols] = partsolid; } } } // end else return(list(solutions, newA, mpoly)); } /////////////////////////////////////////////////////////////////////////////// proc LinearZeroSetMain(ideal I, intvec wt) "USAGE: LinearZeroSetMain(I, wt) PURPOSE: solve the univariate linear polys in I ASSUME: basering = Q[x_1,...,x_n,a] RETURN: list _[1] = partial solution of I _[2] = index of solved vars _[3] = new generators (standardbasis) " { def altring=basering; int i, ok, n, found, nrSols; ideal generators, newGens; list result, index, totalIndex, vars, sol, temp; map F; poly f; result[1] = index; // sol[1] should be the empty list n = nvars(basering) - 1; generators = I; // might be wrong, use index ! ok = 1; nrSols = 0; while(ok) { found = 0; for(i = 1; i <= size(generators); i = i + 1) { f = generators[i]; vars = Variables(f, n); if(size(vars) == 1 && deg(f, wt) == 1) { // univariate,linear nrSols++; found++; index[nrSols] = vars[1]; sol[nrSols] = var(vars[1]) - MultPolys(invertNumberMain(LeadTerm(f, vars[1])[3]), f); } } if(found > 0) { F = basering, SubsMapIdeal(sol, index, 1); newGens = std(SimplifyData(F(generators))); // substitute, simplify alg. number if(size(newGens) == 0) {ok = 0;} generators = newGens; } else { ok = 0; } } if(nrSols > 0) { result[1] = sol;} result[2] = index; result[3] = generators; return(result); } /////////////////////////////////////////////////////////////////////////////// proc NonLinearZeroSetMain(ideal I, intvec wt) "USAGE: NonLinearZeroSetMain(I, wt); PURPOSE: solves the (nonlinear) univariate polynomials in I of the ground field (without multiplicities). RETURN: list, all entries are polynomials _[1] = list of solutions _[2] = newA _[3] = minpoly _[4] - index of solved variables _[5] = new representation of I ASSUME: basering = Q[x_1,x_2,...,x_n,a], ideal 'mpoly' must be defined, it might be 0 ! NOTE: might change 'mpoly' !! " { int i, nrSols, ok, n; ideal generators; list result, sols, index, vars, partSol; map F; poly f, newA; string ringSTR; def NLZR = basering; export(NLZR); n = nvars(basering) - 1; generators = I; newA = var(n + 1); result[2] = newA; // default nrSols = 0; ok = 1; i = 1; while(ok) { // test if the i-th generator of I is univariate f = generators[i]; vars = Variables(f, n); if(size(vars) == 1) { generators[i] = 0; generators = simplify(generators, 2); // remove 0 nrSols++; index[nrSols] = vars[1]; // store index of solved variable // create univariate ring ringSTR = "ring RIS1 = 0, (" + string(var(vars[1])) + ", " + string(var(n+1)) + "), lp;"; execute(ringSTR); ideal mpoly = std(imap(NLZR, mpoly)); list roots; poly f = imap(NLZR, f); export(RIS1); export(mpoly); roots = rootsMain(f); // get "old" basering with new minpoly setring(NLZR); partSol = imap(RIS1, roots); kill RIS1; if(mpoly[1] != partSol[3]) { // change minpoly mpoly = std(partSol[3]); F = NLZR, maxideal(1); F[n + 1] = partSol[2]; if(size(sols) > 0) {sols = SimplifyZeroset(F(sols)); } newA = reduce(F(newA), mpoly); // normal form result[2] = newA; generators = SimplifyData(F(generators)); // does not remove 0's } sols = ExtendSolutions(sols, partSol[1]); } // end univariate else { i = i + 1; } if(i > size(generators)) { ok = 0;} } result[1] = sols; result[3] = mpoly; result[4] = index; result[5] = std(generators); kill NLZR; return(result); } /////////////////////////////////////////////////////////////////////////////// static proc ExtendSolutions(list solutions, list newSolutions) "USAGE: ExtendSolutions(sols, newSols); list sols, newSols; PURPOSE: extend the entries of 'sols' by the entries of 'newSols', each entry of 'newSols' is a number. RETURN: list ASSUME: basering = Q[x_1,...,x_n,a], ideal 'mpoly' must be defined, it might be 0 ! NOTE: used by 'NonLinearZeroSetMain' " { int i, j, k, n, nrSols; list newSols, temp; nrSols = size(solutions); if(nrSols > 0) {n = size(solutions[1]);} else { n = 0; nrSols = 1; } k = 0; for(i = 1; i <= nrSols; i++) { for(j = 1; j <= size(newSolutions); j++) { k++; if(n == 0) { temp[1] = newSolutions[j];} else { temp = solutions[i]; temp[n + 1] = newSolutions[j]; } newSols[k] = temp; } } return(newSols); } /////////////////////////////////////////////////////////////////////////////// static proc MergeSolutions(list sol1, list index1, list sol2, list index2) "USAGE: MergeSolutions(sol1, index1, sol2, index2); all parameters are lists RETURN: list PURPOSE: create a list of solutions of size n, each entry of 'sol2' must have size n. 'sol1' is one partial solution (from 'LinearZeroSetMain') 'sol2' is a list of partial solutions (from 'NonLinearZeroSetMain') ASSUME: 'sol2' is not empty NOTE: used by 'zeroSetMainWork' { int i, j, k, m; ideal sol; list newSols; m = 0; for(i = 1; i <= size(sol2); i++) { m++; newSols[m] = SubsMapIdeal(sol1 + sol2[i], index1 + index2, 0); } return(newSols); } /////////////////////////////////////////////////////////////////////////////// static proc SubsMapIdeal(list sol, list index, int opt) "USAGE: SubsMapIdeal(sol,index,opt); list sol, index; int opt; PURPOSE: built an ideal I as follows. if i is contained in 'index' then set I[i] = sol[i] if i is not contained in 'index' then - opt = 0: set I[i] = 0 - opt = 1: set I[i] = var(i) if opt = 1 and n = nvars(basering) then set I[n] = var(n). RETURN: ideal ASSUME: size(sol) = size(index) <= nvars(basering) " { int k = 0; ideal I; for(int i = 1; i <= nvars(basering) - 1; i = i + 1) { // built subs. map if(containedQ(index, i)) { k++; I[index[k]] = sol[k]; } else { if(opt) { I[i] = var(i); } else { I[i] = 0; } } } if(opt) {I[nvars(basering)] = var(nvars(basering));} return(I); } /////////////////////////////////////////////////////////////////////////////// proc SimplifyZeroset(data) "USAGE: SimplifyZeroset(data); list data PURPOSE: reduce the entries of the elements of 'data' w.r.t. the ideal 'mpoly' 'data' is a list of ideals/lists. RETURN: list ASSUME: basering = Q[x_1,...,x_n,a], order = lp 'data' is a list of ideals ideal 'mpoly' must be defined, it might be 0 ! " { int i; list result; for(i = 1; i <= size(data); i++) { result[i] = SimplifyData(data[i]); } return(result); } /////////////////////////////////////////////////////////////////////////////// proc Variables(poly f, int n) "USAGE: Variables(f,n); poly f; int n; PURPOSE: list of variables among var(1),...,var(n) which occur in f. RETURN: list ASSUME: n <= nvars(basering) " { int i, nrV; list index; nrV = 0; for(i = 1; i <= n; i = i + 1) { if(diff(f, var(i)) != 0) { nrV++; index[nrV] = i; } } return(index); } /////////////////////////////////////////////////////////////////////////////// proc containedQ(data, f, list #) "USAGE: containedQ(data, f [, opt]); data=list; f=any type; opt=integer PURPOSE: test if f is an element of data. RETURN: int 0 if f not contained in data 1 if f contained in data OPTIONS: opt = 0 : use '==' for comparing f with elements from data@* opt = 1 : use @code{sameQ} for comparing f with elements from data " { int opt, i, found; if(size(#) > 0) { opt = #[1];} else { opt = 0; } i = 1; found = 0; while((!found) && (i <= size(data))) { if(opt == 0) { if(f == data[i]) { found = 1;} else {i = i + 1;} } else { if(sameQ(f, data[i])) { found = 1;} else {i = i + 1;} } } return(found); } ////////////////////////////////////////////////////////////////////////////// proc sameQ(a, b) "USAGE: sameQ(a, b); a,b=list/intvec PURPOSE: test a == b elementwise, i.e., a[i] = b[i]. RETURN: int 0 if a != b 1 if a == b " { if(typeof(a) == typeof(b)) { if(typeof(a) == "list" || typeof(a) == "intvec") { if(size(a) == size(b)) { int i = 1; int ok = 1; while(ok && (i <= size(a))) { if(a[i] == b[i]) { i = i + 1;} else {ok = 0;} } return(ok); } else { return(0); } } else { return(a == b);} } else { return(0);} } /////////////////////////////////////////////////////////////////////////////// static proc SimplifyPoly(poly f) "USAGE: SimplifyPoly(f); poly f PURPOSE: reduces the coefficients of f w.r.t. the ideal 'moly' if they contain the algebraic number 'a'. RETURN: poly ASSUME: basering = Q[x_1,...,x_n,a] ideal mpoly must be defined, it might be 0 ! NOTE: outdated, use reduce instead " { matrix coMx; poly f1, vp; vp = 1; for(int i = 1; i < nvars(basering); i++) { vp = vp * var(i);} coMx = coef(f, vp); f1 = 0; for(i = 1; i <= ncols(coMx); i++) { f1 = f1 + coMx[1, i] * reduce(coMx[2, i], mpoly); } return(f1); } /////////////////////////////////////////////////////////////////////////////// static proc SimplifyData(data) "USAGE: SimplifyData(data); ideal/list data; PURPOSE: reduces the entries of 'data' w.r.t. the ideal 'mpoly' if they contain the algebraic number 'a' RETURN: ideal/list ASSUME: basering = Q[x_1,...,x_n,a] ideal 'mpoly' must be defined, it might be 0 ! " { def altring=basering; int n; poly f; if(typeof(data) == "ideal") { n = ncols(data); } else { n = size(data);} for(int i = 1; i <= n; i++) { f = data[i]; data[i] = SimplifyPoly(f); } return(data); } /////////////////////////////////////////////////////////////////////////////// static proc TransferRing(R) "USAGE: TransferRing(R); PURPOSE: creates a new ring containing the same variables as R, but without parameters. If R contains a parameter then this parameter is added as the last variable and 'minpoly' is represented by the ideal 'mpoly' If the basering does not contain a parameter then 'a' is added and 'mpoly' = 0. RETURN: ring ASSUME: R = K[x_1,...,x_n] where K = Q or K = Q(a). NOTE: Creates the ring needed for all prodecures with name 'proc-name'Main " { def altring=basering; string ringSTR, parName, minPoly; setring(R); if(npars(basering) == 0) { parName = "a"; minPoly = "0"; } else { parName = parstr(basering); minPoly = string(minpoly); } ringSTR = "ring TR = 0, (" + varstr(basering) + "," + parName + "), lp;"; execute(ringSTR); execute("ideal mpoly = std(" + minPoly + ");"); export(mpoly); setring altring; return(TR); } /////////////////////////////////////////////////////////////////////////////// static proc NewBaseRing() "USAGE: NewBaseRing(); PURPOSE: creates a new ring, the last variable is added as a parameter. minpoly is set to mpoly[1]. RETURN: ring ASSUME: basering = Q[x_1,...,x_n, a], 'mpoly' must be defined " { int n = nvars(basering); int MP; string ringSTR, parName, varString; def BR = basering; if(mpoly[1] != 0) { parName = "(0, " + string(var(n)) + ")"; MP = 1; } else { parName = "0"; MP = 0; } for(int i = 1; i < n - 1; i++) { varString = varString + string(var(i)) + ","; } varString = varString + string(var(n-1)); ringSTR = "ring TR = " + parName + ", (" + varString + "), lp;"; execute(ringSTR); if(MP) { minpoly = number(imap(BR, mpoly)[1]); } setring BR; return(TR); } /////////////////////////////////////////////////////////////////////////////// /* Examples: // order = 20; ring S1 = 0, (s(1..3)), lp; ideal I = s(2)*s(3), s(1)^2*s(2)+s(1)^2*s(3)-1, s(1)^2*s(3)^2-s(3), s(2)^4-s(3)^4+s(1)^2, s(1)^4+s(2)^3-s(3)^3, s(3)^5-s(1)^2*s(3); ideal mpoly = std(0); // order = 10 ring S2 = 0, (s(1..5)), lp; ideal I = s(2)+s(3)-s(5), s(4)^2-s(5), s(1)*s(5)+s(3)*s(4)-s(4)*s(5), s(1)*s(4)+s(3)-s(5), s(3)^2-2*s(3)*s(5), s(1)*s(3)-s(1)*s(5)+s(4)*s(5), s(1)^2+s(4)^2-2*s(5), -s(1)+s(5)^3, s(3)*s(5)^2+s(4)-s(5)^3, s(1)*s(5)^2-1; ideal mpoly = std(0); //order = 126 ring S3 = 0, (s(1..5)), lp; ideal I = s(3)*s(4), s(2)*s(4), s(1)*s(3), s(1)*s(2), s(3)^3+s(4)^3-1, s(2)^3+s(4)^3-1, s(1)^3-s(4)^3, s(4)^4-s(4), s(1)*s(4)^3-s(1), s(5)^7-1; ideal mpoly = std(0); // order = 192 ring S4 = 0, (s(1..4)), lp; ideal I = s(2)*s(3)^2*s(4)+s(1)*s(3)*s(4)^2, s(2)^2*s(3)*s(4)+s(1)*s(2)*s(4)^2, s(1)*s(3)^3+s(2)*s(4)^3, s(1)*s(2)*s(3)^2+s(1)^2*s(3)*s(4), s(1)^2*s(3)^2-s(2)^2*s(4)^2, s(1)*s(2)^2*s(3)+s(1)^2*s(2)*s(4), s(1)^3*s(3)+s(2)^3*s(4), s(2)^4-s(3)^4, s(1)*s(2)^3+s(3)*s(4)^3, s(1)^2*s(2)^2-s(3)^2*s(4)^2, s(1)^3*s(2)+s(3)^3*s(4), s(1)^4-s(4)^4, s(3)^5*s(4)-s(3)*s(4)^5, s(3)^8+14*s(3)^4*s(4)^4+s(4)^8-1, 15*s(2)*s(3)*s(4)^7-s(1)*s(4)^8+s(1), 15*s(3)^4*s(4)^5+s(4)^9-s(4), 16*s(3)*s(4)^9-s(3)*s(4), 16*s(2)*s(4)^9-s(2)*s(4), 16*s(1)*s(3)*s(4)^8-s(1)*s(3), 16*s(1)*s(2)*s(4)^8-s(1)*s(2), 16*s(1)*s(4)^10-15*s(2)*s(3)*s(4)-16*s(1)*s(4)^2, 16*s(1)^2*s(4)^9-15*s(1)*s(2)*s(3)-16*s(1)^2*s(4), 16*s(4)^13+15*s(3)^4*s(4)-16*s(4)^5; ideal mpoly = std(0); ring R = (0,a), (x,y,z), lp; minpoly = a2 + 1; ideal I1 = x2 - 1/2, a*z - 1, y - 2; ideal I2 = x3 - 1/2, a*z2 - 3, y - 2*a; */