source: git/Singular/LIB/zeroset.lib @ 6fe3a0

fieker-DuValspielwiese
Last change on this file since 6fe3a0 was 6fe3a0, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes/lossen: new libs, new tests git-svn-id: file:///usr/local/Singular/svn/trunk@7891 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 46.8 KB
Line 
1// Last change 12.02.2001 (Eric Westenberger)
2///////////////////////////////////////////////////////////////////////////////
3version="$Id: zeroset.lib,v 1.11 2005-04-25 10:13:07 Singular Exp $";
4category="Symbolic-numerical solving";
5info="
6LIBRARY:  zeroset.lib      Procedures For Roots and Factorization
7AUTHOR:   Thomas Bayer,    email: tbayer@mathematik.uni-kl.de
8          http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
9          Current Adress: Institut fuer Informatik, TU Muenchen
10
11OVERVIEW:
12 Algorithms for finding the zero-set of a zero-dim. ideal in Q(a)[x_1,..,x_n],
13 Roots and Factorization of univariate polynomials over Q(a)[t]
14 where a is an algebraic number. Written in the frame of the
15 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli
16 spaces of semiquasihomogeneous singularities and an implementation in Singular'.
17 This library is meant as a preliminary extension of the functionality
18 of Singular for univariate factorization of polynomials over simple algebraic
19 extensions in characteristic 0.
20 Subprocedures with postfix 'Main' require that the ring contains a variable
21 'a' and no parameters, and the ideal 'mpoly', where 'minpoly' from the
22 basering is stored.
23
24PROCEDURES:
25 EGCD(f, g)    gcd over an algebraic extension field of Q
26 Factor(f)    factorization of f over an algebraic extension field
27 Quotient(f, g)    quotient q  of f w.r.t. g (in f = q*g + remainder)
28 Remainder(f,g)    remainder of the division of f by g
29 Roots(f)    computes all roots of f in an extension field of Q
30 SQFRNorm(f)    norm of f (f must be squarefree)
31 ZeroSet(I)    zero-set of the 0-dim. ideal I
32
33AUXILIARY PROCEDURES:
34 EGCDMain(f, g)    gcd over an algebraic extension field of Q
35 FactorMain(f)    factorization of f over an algebraic extension field
36 InvertNumberMain(c)  inverts an element of an algebraic extension field
37 QuotientMain(f, g)  quotient of f w.r.t. g
38 RemainderMain(f,g)  remainder of the division of f by g
39 RootsMain(f)    computes all roots of f, might extend the ground field
40 SQFRNormMain(f)  norm of f (f must be squarefree)
41 ContainedQ(data, f)  f in data ?
42 SameQ(a, b)    a == b (list a,b)
43";
44
45LIB "primitiv.lib";
46LIB "primdec.lib";
47
48// note : return a ring : ring need not be exported !!!
49
50// Artihmetic in Q(a)[x] without built-in procedures
51// assume basering = Q[x,a] and minpoly is represented by mpoly(a).
52// the algorithms are taken from "Polynomial Algorithms in Computer Algebra",
53// F. Winkler, Springer Verlag Wien, 1996.
54
55
56// To do :
57// squarefree factorization
58// multiplicities
59
60// Improvement :
61// a main problem is the growth of the coefficients. Try Roots(x7 - 1)
62// return ideal mpoly !
63// mpoly is not monic, comes from primitive_extra
64
65// IMPLEMENTATION
66//
67// In procedures with name 'proc-name'Main a polynomial ring over a simple
68// extension field is represented as Q[x...,a] together with the ideal
69// 'mpoly' (attribute "isSB"). The arithmetic in the extension field is
70// implemented in the procedures in the procedures 'MultPolys' (multiplication)
71// and 'InvertNumber' (inversion). After addition and substraction one should
72// apply 'SimplifyPoly' to the result to reduce the result w.r.t. 'mpoly'.
73// This is done by reducing each coefficient seperately, which is more
74// efficient for polynomials with many terms.
75
76
77///////////////////////////////////////////////////////////////////////////////
78
79proc Roots(poly f)
80"USAGE:   Roots(f); where f is a polynomial
81PURPOSE: compute all roots of f in a finite extension of the ground field
82         without multiplicities.
83RETURN:  ring, a polynomial ring over an extension field of the ground field,
84         containing a list 'roots' and polynomials 'newA' and 'f':
85  @format
86  - 'roots' is the list of roots of the polynomial f (no multiplicities)
87  - if the ground field is Q(a') and the extension field is Q(a), then
88    'newA' is the representation of a' in Q(a).
89    If the basering contains a parameter 'a' and the minpoly remains unchanged
90    then 'newA' = 'a'.
91    If the basering does not contain a parameter then 'newA' = 'a' (default).
92  - 'f' is the polynomial f in Q(a) (a' being substituted by 'newA')
93  @end format
94ASSUME:  ground field to be Q or a simple extension of Q given by a minpoly
95EXAMPLE: example  Roots; shows an example
96"
97{
98  int dbPrt = printlevel-voice+3;
99
100  // create a new ring where par(1) is replaced by the variable
101  // with the same name or, if basering does not contain a parameter,
102  // with a new variable 'a'.
103
104  def ROB = basering;
105  def ROR = TransferRing(basering);
106  setring ROR;
107  export(ROR);
108
109  // get the polynomial f and find the roots
110
111  poly f = imap(ROB, f);
112  list result = RootsMain(f);  // find roots of f
113
114  // store the roots and the the new representation of 'a' and transform
115  // the coefficients of f.
116
117  list roots = result[1];
118  poly newA = result[2];
119  map F = basering, maxideal(1);
120  F[nvars(basering)] = newA;
121  poly fn = SimplifyPoly(F(f));
122
123  // create a new ring with minploy = mpoly[1] (from ROR)
124
125  def RON = NewBaseRing();
126  setring(RON);
127  list roots = imap(ROR, roots);
128  poly newA = imap(ROR, newA);
129  poly f = imap(ROR, fn);
130  kill(ROR);
131  export(roots);
132  export(newA);
133  export(f); dbprint(dbPrt,"
134// 'Roots' created a new ring which contains the list 'roots' and
135// the polynomials 'f' and 'newA'
136// To access the roots, newA and the new representation of f, type
137   def R = Roots(f); setring R; roots; newA; f;
138");
139  return(RON);
140}
141example
142{"EXAMPLE:";  echo = 2;
143  ring R = (0,a), x, lp;
144  minpoly = a2+1;
145  poly f = x3 - a;
146  def R1 = Roots(f);
147  setring R1;
148  minpoly;
149  newA;
150  f;
151  roots;
152  map F;
153  F[1] = roots[1];
154  F(f);
155}
156
157///////////////////////////////////////////////////////////////////////////////
158
159proc RootsMain(poly f)
160"USAGE:   RootsMain(f); where f is a polynomial
161PURPOSE: compute all roots of f in a finite extension of the ground field
162         without multiplicities.
163RETURN:  list, all entries are polynomials
164  @format
165  _[1] = roots of f, each entry is a polynomial
166  _[2] = 'newA' - if the ground field is Q(a') and the extension field
167         is Q(a), then 'newA' is the representation of a' in Q(a)
168  _[3] = minpoly of the algebraic extension of the ground field
169  @end format
170ASSUME:  basering = Q[x,a] ideal mpoly must be defined, it might be 0!
171NOTE:    might change the ideal mpoly !!
172EXAMPLE: example  Roots; shows an example
173"
174{
175  int i, linFactors, nlinFactors, dbPrt;
176  intvec wt = 1,0;    // deg(a) = 0
177  list factorList, nlFactors, nlMult, roots, result;
178  poly fa, lc;
179
180  dbPrt = printlevel-voice+3;
181
182  // factor f in Q(a)[t] to obtain the roots lying in Q(a)
183  // firstly, find roots of the linear factors,
184  // nonlinear factors are processed later
185
186  dbprint(dbPrt, "Roots of " + string(f) +  ", minimal polynomial = " + string(mpoly[1]));
187  factorList = FactorMain(f);          // Factorize f
188  dbprint(dbPrt, (" prime factors of f are : " + string(factorList[1])));
189
190  linFactors = 0;
191  nlinFactors = 0;
192  for(i = 2; i <= size(factorList[1]); i = i + 1) {  // find linear and nonlinear factors
193    fa = factorList[1][i];
194    if(deg(fa, wt) == 1) {
195      linFactors++;        // get the root from the linear factor
196      lc = LeadTerm(fa, 1)[3];
197      fa = MultPolys(InvertNumberMain(lc), fa); // make factor monic
198      roots[linFactors] = var(1) - fa;  // fa is monic !!
199    }
200    else {            // ignore nonlinear factors
201      nlinFactors++;
202      nlFactors[nlinFactors] = factorList[1][i];
203      nlMult[nlinFactors] = factorList[2][i];
204    }
205  }
206  if(linFactors == size(factorList[1]) - 1) {    // all roots of f are contained in the ground field
207    result[1] = roots;
208    result[2] = var(2);
209    result[3] = mpoly[1];
210    return(result);
211  }
212
213  // process the nonlinear factors, i.e., extend the ground field
214  // where a nonlinear factor (irreducible) is a minimal polynomial
215  // compute the primitive element of this extension
216
217  ideal primElem, minPolys, Fid;
218  list partSol;
219  map F, Xchange;
220  poly f1, newA, mp, oldMinPoly;
221
222  Fid = mpoly;
223  F[1] = var(1);
224  Xchange[1] = var(2);      // the variables have to be exchanged
225  Xchange[2] = var(1);      // for the use of 'primitive'
226
227  if(nlinFactors == 1) {            // one nl factor
228
229    // compute the roots of the nonlinear (irreducible, monic) factor f1 of f
230    // by extending the basefield by a' with minimal polynomial f1
231    // Then call Roots(f1) to find the roots of f1 over the new base field
232
233    f1 = nlFactors[1];
234    if(mpoly[1] != 0) {
235      mp = mpoly[1];
236      minPolys = Xchange(mp), Xchange(f1);
237      primElem = primitive_extra(minPolys);  // no random coord. change
238      mpoly = std(primElem[1]);
239      F = basering, maxideal(1);
240      F[2] = primElem[2];      // transfer all to the new representation
241      newA = primElem[2];      // new representation of a
242      f1 = SimplifyPoly(F(f1));     //reduce(F(f1), mpoly);
243      if(size(roots) > 0) {roots = SimplifyData(F(roots));}
244    }
245    else {
246      mpoly = std(Xchange(f1));
247      newA = var(2);
248    }
249    result[3] = mpoly[1];
250    oldMinPoly = mpoly[1];
251    partSol = RootsMain(f1);    // find roots of f1 over extended field
252
253    if(oldMinPoly != partSol[3]) {    // minpoly has changed ?
254      // all previously computed roots must be transformed
255      // because the minpoly has changed
256      result[3] = partSol[3];    // new minpoly
257      F[2] = partSol[2];    // new representation of algebraic number
258      if(size(roots) > 0) {roots = SimplifyData(F(roots)); }
259      newA = SimplifyPoly(F(newA)); // F(newA);
260    }
261    roots = roots + partSol[1];  // add roots
262    result[2] = newA;
263    result[1] = roots;
264  }
265  else {  // more than one nonlinear (irreducible) factor (f_1,...,f_r)
266    // solve each of them by RootsMain(f_i), append their roots
267    // change the minpoly and transform all previously computed
268    // roots if necessary.
269    // Note that the for-loop is more or less book-keeping
270
271    newA = var(2);
272    result[2] = newA;
273    for(i = 1; i <= size(nlFactors); i = i + 1) {
274      oldMinPoly = mpoly[1];
275      partSol = RootsMain(nlFactors[i]);    // main work
276      nlFactors[i] = 0;        // delete factor
277      result[3] = partSol[3];        // store minpoly
278
279      // book-keeping starts here as in the case 1 nonlinear factor
280
281      if(oldMinPoly != partSol[3]) { // minpoly has changed
282        F = basering, maxideal(1);
283        F[2] = partSol[2];    // transfer all to the new representation
284        newA = SimplifyPoly(F(newA));    // F(newA); new representation of a
285        result[2] = newA;
286        if(i < size(nlFactors)) {
287          nlFactors = SimplifyData(F(nlFactors));
288        } // transform remaining factors
289        if(size(roots) > 0) {roots = SimplifyData(F(roots));}
290      }
291      roots = roots + partSol[1];    // transform roots
292      result[1] = roots;
293    }  // end more than one nl factor
294
295  }
296  return(result);
297}
298
299///////////////////////////////////////////////////////////////////////////////
300
301proc ZeroSet(ideal I, list #)
302"USAGE:   ZeroSet(I [,opt] ); I=ideal, opt=integer
303PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
304         of the ground field.
305RETURN:  ring, a polynomial ring over an extension field of the ground field,
306         containing a list 'zeroset', a polynomial 'newA', and an
307         ideal 'id':
308  @format
309  - 'zeroset' is the list of the zeros of the ideal I, each zero is an ideal.
310  - if the ground field is Q(a') and the extension field is Q(a), then
311    'newA' is the representation of a' in Q(a).
312    If the basering contains a parameter 'a' and the minpoly remains unchanged
313    then 'newA' = 'a'.
314    If the basering does not contain a parameter then 'newA' = 'a' (default).   
315  - 'id' is the ideal I in Q(a)[x_1,...] (a' substituted by 'newA')
316  @end format
317ASSUME:  dim(I) = 0, and ground field to be Q or a simple extension of Q given
318         by a minpoly.
319OPTIONS: opt = 0 no primary decomposition (default)
320         opt > 0 primary decomposition
321NOTE:    If I contains an algebraic number (parameter) then 'I' must be
322         transformed w.r.t. 'newA' in the new ring.
323EXAMPLE: example ZeroSet; shows an example
324"
325{
326  int primaryDecQ, dbPrt;
327  list rp;
328
329  dbPrt = printlevel-voice+2;
330
331  if(size(#) > 0) { primaryDecQ = #[1]; }
332  else { primaryDecQ = 0; }
333
334  // create a new ring 'ZSR' with one additional variable instead of the
335  // parameter
336  // if the basering does not contain a parameter then 'a' is used as the
337  // additional variable.
338
339  def RZSB = basering;
340  def ZSR = TransferRing(RZSB);
341  setring ZSR;
342
343  // get ideal I and find the zero-set
344
345  ideal id = std(imap(RZSB, I));
346//  print(dim(id));
347  if(dim(id) > 1) {       // new variable adjoined to ZSR
348    ERROR(" ideal not zerodimensional ");
349  }
350
351  list result = ZeroSetMain(id, primaryDecQ);
352
353  // store the zero-set, minimal polynomial and the new representative of 'a'
354
355  list zeroset = result[1];
356  poly newA = result[2];
357  poly minPoly = result[3][1];
358
359  // transform the generators of the ideal I w.r.t. the new representation
360  // of 'a'
361
362  map F = basering, maxideal(1);
363  F[nvars(basering)] = newA;
364  id = SimplifyData(F(id));
365
366  // create a new ring with minpoly = minPoly
367
368  def RZBN = NewBaseRing();
369  setring RZBN;
370
371  list zeroset = imap(ZSR, zeroset);
372  poly newA = imap(ZSR, newA);
373  ideal id = imap(ZSR, id);
374  kill(ZSR);
375
376  export(id);
377  export(zeroset);
378  export(newA);
379    dbprint(dbPrt,"
380// 'ZeroSet' created a new ring which contains the list 'zeroset', the ideal
381// 'id' and the polynomial 'newA'. 'id' is the ideal of the input transformed
382// w.r.t. 'newA'.
383// To access the zero-set, 'newA' and the new representation of the ideal, type
384   def R = ZeroSet(I); setring R; zeroset; newA; id;
385");
386  setring RZSB;
387  return(RZBN);
388}
389example
390{"EXAMPLE:";  echo = 2;
391  ring R = (0,a), (x,y,z), lp;
392  minpoly = a2 + 1;
393  ideal I = x2 - 1/2, a*z - 1, y - 2;
394  def T = ZeroSet(I);
395  setring T;
396  minpoly;
397  newA;
398  id;
399  zeroset;
400  map F1 = basering, zeroset[1];
401  map F2 = basering, zeroset[2];
402  F1(id);
403  F2(id);
404}
405
406///////////////////////////////////////////////////////////////////////////////
407
408proc InvertNumberMain(poly f)
409"USAGE:   InvertNumberMain(f); where f is a polynomial
410PURPOSE: compute 1/f if f is a number in Q(a) i.e., f is represented by a
411         polynomial in Q[a].
412RETURN:  poly 1/f
413ASSUME:  basering = Q[x_1,...,x_n,a], ideal mpoly must be defined and != 0 !
414"
415{
416  if(diff(f, var(1)) != 0) { ERROR("number must not contain variable !");}
417
418  int n = nvars(basering);
419  def RINB = basering;
420  string ringSTR = "ring RINR = 0, " + string(var(n)) + ", dp;";
421  execute(ringSTR);        // new ring = Q[a]
422
423  list gcdList;
424  poly f, g, inv;
425
426  f = imap(RINB, f);
427  g = imap(RINB, mpoly)[1];
428
429  if(diff(f, var(1)) != 0) { inv = extgcd(f, g)[2]; }  // f contains var(1)
430  else {  inv = 1/f;}          // f element in Q
431
432  setring(RINB);
433  return(imap(RINR, inv));
434}
435
436///////////////////////////////////////////////////////////////////////////////
437
438proc MultPolys(poly f, poly g)
439"USAGE:   MultPolys(f, g); poly f,g
440PURPOSE: multiply the polynomials f and g and reduce them w.r.t. mpoly
441RETURN:  poly f*g
442ASSUME:  basering = Q[x,a], ideal mpoly must be defined, it might be 0 !
443"
444{
445  return(SimplifyPoly(f * g));
446}
447
448///////////////////////////////////////////////////////////////////////////////
449
450proc LeadTerm(poly f, int i)
451"USAGE:   LeadTerm(f); poly f, int i
452PURPOSE: compute the leading coef and term of f w.r.t var(i), where the last
453         ring variable is treated as a parameter.
454RETURN:  list of polynomials
455         _[1] = leading term
456         _[2] = leading monomial
457         _[3] = leading coefficient
458ASSUME:  basering = Q[x_1,...,x_n,a]
459"
460{
461  list result;
462  matrix co = coef(f, var(i));
463  result[1] = co[1, 1]*co[2, 1];
464  result[2] = co[1, 1];
465  result[3] = co[2, 1];
466  return(result);
467}
468
469///////////////////////////////////////////////////////////////////////////////
470
471proc Quotient(poly f, poly g)
472"USAGE:   Quotient(f, g); where f,g are polynomials;
473PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g)
474RETURN:  list of polynomials
475  @format
476  _[1] = quotient  q
477  _[2] = remainder r
478  @end format
479ASSUME:  basering = Q[x] or Q(a)[x]
480EXAMPLE: example  Quotient; shows an example
481"
482{
483  def QUOB = basering;
484  def QUOR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
485  setring QUOR;
486  export(QUOR);
487  poly f = imap(QUOB, f);
488  poly g = imap(QUOB, g);
489  list result = QuotientMain(f, g);
490
491  setring(QUOB);
492  list result = imap(QUOR, result);
493  kill(QUOR);
494  return(result);
495}
496example
497{"EXAMPLE:";  echo = 2;
498 ring R = (0,a), x, lp;
499 minpoly = a2+1;
500 poly f =  x4 - 2;
501 poly g = x - a;
502 list qr = Quotient(f, g);
503 qr;
504 qr[1]*g + qr[2] - f;
505}
506
507proc QuotientMain(poly f, poly g)
508"USAGE:   QuotientMain(f, g); where f,g are polynomials
509PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g)
510RETURN:  list of polynomials
511  @format
512  _[1] = quotient  q
513  _[2] = remainder r
514  @end format
515ASSUME:  basering = Q[x,a] and ideal mpoly is defined (it might be 0),
516         this represents the ring Q(a)[x] together with its minimal polynomial.
517EXAMPLE: example  Quotient; shows an example
518"
519{
520  if(g == 0) { ERROR("Division by zero !");}
521
522  def QMB = basering;
523  def QMR = NewBaseRing();
524  setring QMR;
525  poly f, g, h;
526  h = imap(QMB, f) / imap(QMB, g);
527  setring QMB;
528  return(list(imap(QMR, h), 0));
529}
530
531///////////////////////////////////////////////////////////////////////////////
532
533proc Remainder(poly f, poly g)
534"USAGE:   Remainder(f, g); where f,g are polynomials
535PURPOSE: compute the remainder of the division of f by g, i.e. a polynomial r
536         s.t. f = g*q + r, deg(r) < deg(g).
537RETURN:  poly
538ASSUME:  basering = Q[x] or Q(a)[x]
539"
540{
541  def REMB = basering;
542  def REMR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
543  setring(REMR);
544  export(REMR);
545  poly f = imap(REMB, f);
546  poly g = imap(REMB, g);
547  poly h = RemainderMain(f, g);
548
549  setring(REMB);
550  poly r = imap(REMR, h);
551  kill(REMR);
552  return(r);
553}
554example
555{"EXAMPLE:";  echo = 2;
556 ring R = (0,a), x, lp;
557 minpoly = a2+1;
558 poly f =  x4 - 1;
559 poly g = x3 - 1;
560 Remainder(f, g);
561}
562
563proc RemainderMain(poly f, poly g)
564"USAGE:   RemainderMain(f, g); where f,g are polynomials
565PURPOSE: compute the remainder r s.t. f = g*q + r, deg(r) < deg(g)
566RETURN:  poly
567ASSUME:  basering = Q[x,a] and ideal mpoly is defined (it might be 0),
568         this represents the ring Q(a)[x] together with its minimal polynomial.
569"
570{
571  int dg;
572  intvec wt = 1,0;;
573  poly lc, g1, r;
574
575  if(deg(g, wt) == 0) { return(0); }
576
577  lc = LeadTerm(g, 1)[3];
578  g1 = MultPolys(InvertNumberMain(lc), g);  // make g monic
579
580  return(SimplifyPoly(reduce(f, std(g1))));
581}
582
583///////////////////////////////////////////////////////////////////////////////
584
585proc EGCD(poly f, poly g)
586"USAGE:   EGCD(f, g); where f,g are polynomials
587PURPOSE: compute the polynomial gcd of f and g over Q(a)[x]
588RETURN:  polynomial h s.t. h is a greatest common divisor of f and g (not nec.
589         monic)
590ASSUME:  basering = Q(a)[t]
591EXAMPLE: example  EGCD; shows an example
592"
593{
594  def GCDB = basering;
595  def GCDR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
596  setring GCDR;
597  export(GCDR);
598  poly f = imap(GCDB, f);
599  poly g = imap(GCDB, g);
600  poly h = EGCDMain(f, g);  // squarefree norm of f
601
602  setring(GCDB);
603  poly h = imap(GCDR, h);
604  kill(GCDR);
605  return(h);
606}
607example
608{"EXAMPLE:";  echo = 2;
609 ring R = (0,a), x, lp;
610 minpoly = a2+1;
611 poly f =  x4 - 1;
612 poly g = x2 - 2*a*x - 1;
613 EGCD(f, g);
614}
615
616proc EGCDMain(poly f, poly g)
617"USAGE:   EGCDMain(f, g); where f,g are polynomials
618PURPOSE: compute the polynomial gcd of f and g over Q(a)[x]
619RETURN:  poly
620ASSUME:  basering = Q[x,a] and ideal mpoly is defined (it might be 0),
621         this represents the ring Q(a)[x] together with its minimal polynomial.
622EXAMPLE: example EGCD; shows an example
623"
624// might be extended to return s1, s2 s.t. f*s1 + g*s2 = gcd
625{
626  int i = 1;
627  poly r1, r2, r;
628
629  r1 = f;
630  r2 = g;
631
632  while(r2 != 0) {
633    r  = RemainderMain(r1, r2);
634    r1 = r2;
635    r2 = r;
636  }
637  return(r1);
638}
639
640///////////////////////////////////////////////////////////////////////////////
641
642proc MEGCD(poly f, poly g, int varIndex)
643"USAGE:   MEGCD(f, g, i); poly f, g; int i
644PURPOSE: compute  the polynomial gcd of f and g in the i'th variable
645RETURN:  poly
646ASSUME:  f, g are polynomials in var(i), last variable is the algebraic number
647EXAMPLE: example  MEGCD; shows an example
648"
649// might be extended to return s1, s2 s.t. f*s1 + g*s2 = gc
650// not used !
651{
652  string @str, @sf, @sg, @mp, @parName;
653
654  def @RGCDB = basering;
655
656  @sf = string(f);
657  @sg = string(g);
658  @mp = string(minpoly);
659
660  if(npars(basering) == 0) { @parName = "0";}
661  else { @parName = "(0, " + parstr(basering) + ")"; }
662  @str = "ring @RGCD = " + @parName + ", " + string(var(varIndex)) + ", dp;";
663  execute(@str);
664  if(@mp != "0") { execute ("minpoly = " + @mp + ";"); }
665  execute("poly @f = " + @sf + ";");
666  execute("poly @g = " + @sg + ";");
667  export(@RGCD);
668  poly @h = EGCD(@f, @g);
669  setring(@RGCDB);
670  poly h = imap(@RGCD, @h);
671  kill(@RGCD);
672  return(h);
673}
674
675///////////////////////////////////////////////////////////////////////////////
676
677proc SQFRNorm(poly f)
678"USAGE:   SQFRNorm(f); where f is a polynomial
679PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x].
680RETURN:  list with 3 entries
681  @format
682  _[1] = squarefree norm of g (poly)
683  _[2] = g (= f(x - s*a)) (poly)
684  _[3] = s (int)
685  @end format
686ASSUME:  f must be squarefree, basering = Q(a)[x] and minpoly != 0.
687NOTE:    the norm is an element of Q[x]
688EXAMPLE: example  SQFRNorm; shows an example
689"
690{
691  def SNB = basering;
692  def SNR = TransferRing(SNB);  // new ring with parameter 'a'
693                                // replaced by a variable
694  setring SNR;
695  poly f = imap(SNB, f);
696  list result = SQFRNormMain(f);  // squarefree norm of f
697
698  setring SNB;
699  list result = imap(SNR, result);
700  kill SNR;
701  return(result);
702}
703example
704{"EXAMPLE:";  echo = 2;
705   ring R = (0,a), x, lp;
706   minpoly = a2+1;
707  poly f =  x4 - 2*x + 1;
708  SQFRNorm(f);
709}
710
711proc SQFRNormMain(poly f)
712"USAGE:   SQFRNorm(f); where f is a polynomial
713PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x].
714RETURN:  list with 3 entries
715  @format
716  _[1] = squarefree norm of g (poly)
717  _[2] = g (= f(x - s*a)) (poly)
718  _[3] = s (int)
719  @end format
720ASSUME:  f must be squarefree, basering = Q[x,a] and ideal mpoly is equal to
721         'minpoly',this represents the ring Q(a)[x] together with 'minpoly'.
722NOTE:   the norm is an element of Q[x]
723EXAMPLE: example  SqfrNorm; shows an example
724"
725{
726  def SNRMB = basering;
727  int s = 0;
728  intvec wt = 1,0;
729  ideal mapId;
730  // list result;
731  poly g, N, N1, h;
732  string ringSTR;
733
734  mapId[1] = var(1) - var(2);    // linear transformation
735  mapId[2] = var(2);
736  map Fs = SNRMB, mapId;
737
738  N = resultant(f, mpoly[1], var(2));  // norm of f
739  N1 = diff(N, var(1));
740  g = f;
741
742  ringSTR = "ring SNRM1 = 0, " + string(var(1)) + ", dp;";  // univariate ring
743  execute(ringSTR);
744  poly N, N1, h;        // N, N1 do not contain 'a', use built-in gcd
745  h = gcd(imap(SNRMB, N), imap(SNRMB, N1));
746  setring(SNRMB);
747  h = imap(SNRM1, h);
748  while(deg(h, wt) != 0) {    // while norm is not squarefree
749    s = s + 1;
750    g = reduce(Fs(g), mpoly);
751    N = reduce(resultant(g, mpoly[1], var(2)), mpoly);  // norm of g
752    N1 = reduce(diff(N, var(1)), mpoly);
753    setring(SNRM1);
754    h = gcd(imap(SNRMB, N), imap(SNRMB, N1));
755    setring(SNRMB);
756    h = imap(SNRM1, h);
757  }
758  return(list(N, g, s));
759}
760
761///////////////////////////////////////////////////////////////////////////////
762
763proc Factor(poly f)
764"USAGE:   Factor(f); where f is a polynomial
765PURPOSE: compute the factorization of the squarefree poly f over Q(a)[t]
766RETURN:  list with two entries
767  @format
768  _[1] = factors (monic), first entry is the leading coefficient
769  _[2] = multiplicities (not yet implemented)
770  @end format
771ASSUME:  basering must be the univariate polynomial ring over a field, which
772         is Q or a simple extension of Q given by a minpoly.
773NOTE:    if basering = Q[t] then this is the built-in @code{factorize}
774EXAMPLE: example  Factor; shows an example
775"
776{
777  def FRB = basering;
778  def FRR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
779  setring FRR;
780  export(FRR);
781  poly f = imap(FRB, f);
782  list result = FactorMain(f);  // squarefree norm of f
783
784  setring(FRB);
785  list result = imap(FRR, result);
786  kill(FRR);
787  return(result);
788}
789example
790{"EXAMPLE:";  echo = 2;
791 ring R = (0,a), x, lp;
792 minpoly = a2+1;
793 poly f =  x4 - 1;
794 list fl = Factor(f);
795 fl;
796 fl[1][1]*fl[1][2]*fl[1][3]*fl[1][4]*fl[1][5] - f;
797}
798
799///////////////////////////////////////////////////////////////////////////////
800
801proc FactorMain(poly f)
802"USAGE:   FactorMain(f); where f is a polynomial
803PURPOSE: compute the factorization of the squarefree poly f over Q(a)[t],
804         minpoly  = p(a).
805RETURN:  list with 2 entries
806  @format
807  _[1] = factors, first is a constant
808  _[2] = multiplicities (not yet implemented)
809  @end format
810ASSUME:  basering = Q[x,a], representing Q(a)[x]. An ideal mpoly must
811         be defined, representing the minimal polynomial (it might be 0!).
812EXAMPLE: example  Factor; shows an example
813"
814// extend this by a squarefree factorization !!
815// multiplicities are not valid !!
816{
817  int i, s;
818  list normList, factorList, quo_rem;
819  poly f1, h, h1, H, g, leadCoef, invCoeff;
820  ideal fac1, fac2;
821  map F;
822
823  // if no minimal polynomial is defined then use 'factorize'
824  // FactorOverQ is wrapped around 'factorize'
825
826  if(mpoly[1] == 0) {
827    // print(" factorize : deg = " + string(deg(f, intvec(1,0))));
828    factorList = factorize(f); // FactorOverQ(f);
829    return(factorList);
830  }
831
832  // if mpoly != 0 and f does not contain the algebraic number, a root of
833  // f might be contained in Q(a). Hence one must not use 'factorize'.
834
835  fac1[1] = 1;
836  fac2[1] = 1;
837  normList = SQFRNormMain(f);
838  // print(" factorize : deg = " + string(deg(normList[1], intvec(1,0))));
839  factorList = factorize(normList[1]);     // factor squarefree norm of f over Q[x]
840  g = normList[2];
841  s = normList[3];
842  F[1] = var(1) + s*var(2);      // inverse transformation
843  F[2] = var(2);
844  fac1[1] = factorList[1][1];
845  fac2[1] = factorList[2][1];
846  for(i = 2; i <= size(factorList[1]); i = i + 1) {
847    H = factorList[1][i];
848    h = EGCDMain(H, g);
849    quo_rem = QuotientMain(g, h);
850    g = quo_rem[1];
851    fac1[i] = SimplifyPoly(F(h));
852    fac2[i] = 1;        // to be changed later
853  }
854  return(list(fac1, fac2));
855}
856
857///////////////////////////////////////////////////////////////////////////////
858
859proc ZeroSetMain(ideal I, int primDecQ)
860"USAGE:   ZeroSetMain(ideal I, int opt); ideal I, int opt
861PURPOSE: compute the zero-set of the zero-dim. ideal I, in a simple extension
862         of the ground field.
863RETURN:  list
864         - 'f' is the polynomial f in  Q(a) (a' being substituted by newA)
865         _[1] = zero-set (list), is the list of the zero-set of the ideal I,
866                each entry is an ideal.
867         _[2] = 'newA';  if the ground field is Q(a') and the extension field
868                is Q(a), then 'newA' is the representation of a' in Q(a).
869                If the basering contains a parameter 'a' and the minpoly
870                remains unchanged then 'newA' = 'a'. If the basering does not
871                contain a parameter then 'newA' = 'a' (default).
872         _[3] = 'mpoly' (ideal), the minimal polynomial of the simple extension
873                of the ground field.
874ASSUME:  basering = K[x_1,x_2,...,x_n] where K = Q or a simple extension of Q
875         given by a minpoly; dim(I) = 0.
876NOTE:    opt = 0  no primary decomposition
877         opt > 0  use a primary decomposition
878EXAMPLE: example ZeroSet; shows an example
879"
880{
881  // main work is done in ZeroSetMainWork, here the zero-set of each ideal from the
882  // primary decompostion is coputed by menas of ZeroSetMainWork, and then the
883  // minpoly and the parameter representing the algebraic extension are
884  // transformed according to 'newA', i.e., only bookeeping is done.
885
886  def altring=basering;
887  int i, j, n, noMP, dbPrt;
888  intvec w;
889  list currentSol, result, idealList, primDecList, zeroSet;
890  ideal J;
891  map Fa;
892  poly newA, oldMinPoly;
893
894  dbPrt = printlevel-voice+2;
895  dbprint(dbPrt, "ZeroSet of " + string(I) + ", minpoly = " + string(minpoly));
896
897  n = nvars(basering) - 1;
898  for(i = 1; i <= n; i++) { w[i] = 1;}
899  w[n + 1] = 0;
900
901  if(primDecQ == 0) { return(ZeroSetMainWork(I, w, 0)); }
902
903  newA = var(n + 1);
904  if(mpoly[1] == 0) { noMP = 1;}
905  else {noMP = 0;}
906
907  primDecList = primdecGTZ(I);      // primary decomposition
908  dbprint(dbPrt, "primary decomposition consists of " + string(size(primDecList)) + " primary ideals ");
909  // idealList = PDSort(idealList);    // high degrees first
910
911  for(i = 1; i <= size(primDecList); i = i + 1) {
912    idealList[i] = primDecList[i][2];  // use prime component
913    dbprint(dbPrt, string(i) + "  " + string(idealList[i]));
914  }
915
916  // compute the zero-set of each primary ideal and join them.
917  // If necessary, change the ground field and transform the zero-set
918
919  dbprint(dbPrt, "
920find the zero-set of each primary ideal, form the union
921and keep track of the minimal polynomials ");
922
923  for(i = 1; i <= size(idealList); i = i + 1) {
924    J = idealList[i];
925    idealList[i] = 0;
926    oldMinPoly = mpoly[1];
927    dbprint(dbPrt, " ideal#" + string(i) + " of " + string(size(idealList)) + " = " + string(J));
928    currentSol = ZeroSetMainWork(J, w, 0);
929
930    if(oldMinPoly != currentSol[3]) {   // change minpoly and transform solutions
931      dbprint(dbPrt, " change minpoly to " + string(currentSol[3][1]));
932      dbprint(dbPrt, " new representation of algebraic number = " + string(currentSol[2]));
933      if(!noMP) {      // transform the algebraic number a
934        Fa = basering, maxideal(1);
935        Fa[n + 1] = currentSol[2];
936        newA = SimplifyPoly(Fa(newA));  // new representation of a
937        if(size(zeroSet) > 0) {zeroSet = SimplifyZeroset(Fa(zeroSet)); }
938        if(i < size(idealList)) { idealList = SimplifyZeroset(Fa(idealList)); }
939      }
940      else { noMP = 0;}
941    }
942    zeroSet = zeroSet + currentSol[1];    // add new elements
943  }
944  return(list(zeroSet, newA, mpoly));
945}
946
947///////////////////////////////////////////////////////////////////////////////
948
949proc ZeroSetMainWork(ideal id, intvec wt, int sVars)
950"USAGE:   ZeroSetMainWork(I, wt, sVars);
951PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
952         of the ground field (without multiplicities).
953RETURN:  list, all entries are polynomials
954         _[1] = zeros, each entry is an ideal
955         _[2] = newA; if the ground field is Q(a') this is the rep. of a' w.r.t. a
956         _[3] = minpoly of the algebraic extension of the ground field (ideal)
957         _[4] = name of algebraic number (default = 'a')
958ASSUME:  basering = Q[x_1,x_2,...,x_n,a]
959         ideal mpoly must be defined, it might be 0!
960NOTE:    might change 'mpoly' !!
961EXAMPLE: example  IdealSolve; shows an example
962"
963{
964  def altring=basering;
965  int i, j, k, nrSols, n, noMP;
966  ideal I, generators, gens, solid, partsolid;
967  list linSol, linearSolution, nLinSol, nonlinSolutions, partSol, sol, solutions, result;
968  list linIndex, nlinIndex, index;
969  map Fa, Fsubs;
970  poly oldMinPoly, newA;
971
972  if(mpoly[1] == 0) { noMP = 1;}
973  else { noMP = 0;}
974  n = nvars(basering) - 1;
975  newA = var(n + 1);
976
977  I = std(id);
978
979  // find linear solutions of univariate generators
980
981  linSol = LinearZeroSetMain(I, wt);
982  generators = linSol[3];      // they are a standardbasis
983  linIndex = linSol[2];
984  linearSolution = linSol[1];
985  if(size(linIndex) + sVars == n) {    // all variables solved
986    solid = SubsMapIdeal(linearSolution, linIndex, 0);
987    result[1] = list(solid);
988    result[2] = var(n + 1);
989    result[3] = mpoly;
990    return(result);
991  }
992
993  // find roots of the nonlinear univariate polynomials of generators
994  // if necessary, transform linear solutions w.r.t. newA
995
996  oldMinPoly = mpoly[1];
997  nLinSol =  NonLinearZeroSetMain(generators, wt);    // find solutions of univariate generators
998  nonlinSolutions = nLinSol[1];    // store solutions
999  nlinIndex = nLinSol[4];     // and index of solved variables
1000  generators = nLinSol[5];    // new generators
1001
1002  // change minpoly if necessary and transform the ideal and the partial solutions
1003
1004  if(oldMinPoly != nLinSol[3]) {
1005    newA = nLinSol[2];
1006    if(!noMP && size(linearSolution) > 0) {    // transform the algebraic number a
1007      Fa = basering, maxideal(1);
1008      Fa[n + 1] = newA;
1009      linearSolution = SimplifyData(Fa(linearSolution));  // ...
1010    }
1011  }
1012
1013  // check if all variables are solved.
1014
1015  if(size(linIndex) + size(nlinIndex) == n - sVars) {
1016    solutions = MergeSolutions(linearSolution, linIndex, nonlinSolutions, nlinIndex, list(), n);
1017  }
1018
1019  else {
1020
1021  // some variables are not solved.
1022  // substitute each partial solution in generators and find the
1023  // zero set of the resulting ideal by recursive application
1024  // of ZeroSetMainWork !
1025
1026  index = linIndex + nlinIndex;
1027  nrSols = 0;
1028  for(i = 1; i <=  size(nonlinSolutions); i = i + 1) {
1029    sol = linearSolution + nonlinSolutions[i];
1030    solid = SubsMapIdeal(sol, index, 1);
1031    Fsubs = basering, solid;
1032    gens = std(SimplifyData(Fsubs(generators)));    // substitute partial solution
1033    oldMinPoly = mpoly[1];
1034    partSol = ZeroSetMainWork(gens, wt, size(index) + sVars);
1035
1036    if(oldMinPoly != partSol[3]) {    // minpoly has changed
1037      Fa = basering, maxideal(1);
1038      Fa[n + 1] = partSol[2];    // a -> p(a), representation of a w.r.t. new minpoly
1039      newA = reduce(Fa(newA), mpoly);
1040      generators = std(SimplifyData(Fa(generators)));
1041      if(size(linearSolution) > 0) { linearSolution = SimplifyData(Fa(linearSolution));}
1042      if(size(nonlinSolutions) > 0) {
1043        nonlinSolutions = SimplifyZeroset(Fa(nonlinSolutions));
1044      }
1045      sol = linearSolution + nonlinSolutions[i];
1046    }
1047
1048    for(j = 1; j <= size(partSol[1]); j++) {   // for all partial solutions
1049      partsolid = partSol[1][j];
1050      for(k = 1; k <= size(index); k++) {
1051        partsolid[index[k]] = sol[k];
1052       }
1053      nrSols++;
1054      solutions[nrSols] = partsolid;
1055    }
1056  }
1057
1058  }  // end else
1059  return(list(solutions, newA, mpoly));
1060}
1061
1062///////////////////////////////////////////////////////////////////////////////
1063
1064proc LinearZeroSetMain(ideal I, intvec wt)
1065"USAGE:   LinearZeroSetMain(I, wt)
1066PURPOSE: solve the univariate linear polys in I
1067ASSUME:  basering = Q[x_1,...,x_n,a]
1068RETURN:  list
1069         _[1] = partial solution of I
1070         _[2] = index of solved vars
1071         _[3] = new generators (standardbasis)
1072"
1073{
1074  def altring=basering;
1075  int i, ok, n, found, nrSols;
1076  ideal generators, newGens;
1077  list result, index, totalIndex, vars, sol, temp;
1078  map F;
1079  poly f;
1080
1081  result[1] = index;      // sol[1] should be the empty list
1082  n = nvars(basering) - 1;
1083  generators = I;        // might be wrong, use index !
1084  ok = 1;
1085  nrSols = 0;
1086  while(ok) {
1087    found = 0;
1088    for(i = 1; i <= size(generators); i = i + 1) {
1089      f = generators[i];
1090      vars = Variables(f, n);
1091      if(size(vars) == 1 && deg(f, wt) == 1) {  // univariate,linear
1092        nrSols++; found++;
1093        index[nrSols] = vars[1];
1094        sol[nrSols] = var(vars[1]) - MultPolys(InvertNumberMain(LeadTerm(f, vars[1])[3]), f);
1095      }
1096    }
1097    if(found > 0) {
1098      F = basering, SubsMapIdeal(sol, index, 1);
1099      newGens = std(SimplifyData(F(generators)));    // substitute, simplify alg. number
1100      if(size(newGens) == 0) {ok = 0;}
1101      generators = newGens;
1102    }
1103    else {
1104      ok = 0;
1105    }
1106  }
1107  if(nrSols > 0) { result[1] = sol;}
1108  result[2] = index;
1109  result[3] = generators;
1110  return(result);
1111}
1112
1113///////////////////////////////////////////////////////////////////////////////
1114
1115proc NonLinearZeroSetMain(ideal I, intvec wt)
1116"USAGE:   ZeroSetMainWork(I, wt, sVars);
1117PURPOSE: solves the (nonlinear) univariate polynomials in I
1118         of the ground field (without multiplicities).
1119RETURN:  list, all entries are polynomials
1120         _[1] = list of solutions
1121         _[2] = newA
1122         _[3] = minpoly
1123         _[4] - index of solved variables
1124         _[5] = new representation of I
1125ASSUME:  basering = Q[x_1,x_2,...,x_n,a], ideal 'mpoly' must be defined,
1126         it might be 0 !
1127NOTE:    might change 'mpoly' !!
1128"
1129{
1130  int i, nrSols, ok, n;
1131  ideal generators;
1132  list result, sols, index, vars, partSol;
1133  map F;
1134  poly f, newA;
1135  string ringSTR;
1136
1137  def NLZR = basering;
1138  export(NLZR);
1139
1140  n = nvars(basering) - 1;
1141
1142  generators = I;
1143  newA = var(n + 1);
1144  result[2] = newA;            // default
1145  nrSols = 0;
1146  ok = 1;
1147  i = 1;
1148  while(ok) {
1149
1150    // test if the i-th generator of I is univariate
1151
1152    f = generators[i];
1153    vars = Variables(f, n);
1154    if(size(vars) == 1) {
1155      generators[i] = 0;
1156      generators = simplify(generators, 2);    // remove 0
1157      nrSols++;
1158      index[nrSols] = vars[1];      // store index of solved variable
1159
1160      // create univariate ring
1161
1162      ringSTR = "ring RIS1 = 0, (" + string(var(vars[1])) + ", " + string(var(n+1)) + "), lp;";
1163      execute(ringSTR);
1164      ideal mpoly = std(imap(NLZR, mpoly));
1165      list roots;
1166      poly f = imap(NLZR, f);
1167      export(RIS1);
1168      export(mpoly);
1169      roots = RootsMain(f);
1170
1171      // get "old" basering with new minpoly
1172
1173      setring(NLZR);
1174      partSol = imap(RIS1, roots);
1175      kill(RIS1);
1176      if(mpoly[1] != partSol[3]) {      // change minpoly
1177        mpoly = std(partSol[3]);
1178        F = NLZR, maxideal(1);
1179        F[n + 1] = partSol[2];
1180        if(size(sols) > 0) {sols = SimplifyZeroset(F(sols)); }
1181        newA = reduce(F(newA), mpoly);    // normal form
1182        result[2] = newA;
1183        generators = SimplifyData(F(generators));  // does not remove 0's
1184      }
1185      sols = ExtendSolutions(sols, partSol[1]);
1186    } // end univariate
1187    else {
1188      i = i + 1;
1189    }
1190    if(i > size(generators)) { ok = 0;}
1191  }
1192  result[1] = sols;
1193  result[3] = mpoly;
1194  result[4] = index;
1195  result[5] = std(generators);
1196
1197  kill(NLZR);
1198  return(result);
1199}
1200
1201///////////////////////////////////////////////////////////////////////////////
1202
1203static proc ExtendSolutions(list solutions, list newSolutions)
1204"USAGE:   ExtendSolutions(sols, newSols); list sols, newSols;
1205PURPOSE: extend the entries of 'sols' by the entries of 'newSols',
1206         each entry of 'newSols' is a number.
1207RETURN:  list
1208ASSUME:  basering = Q[x_1,...,x_n,a], ideal 'mpoly' must be defined,
1209         it might be 0 !
1210NOTE:    used by 'NonLinearZeroSetMain'
1211"
1212{
1213  int i, j, k, n, nrSols;
1214  list newSols, temp;
1215
1216  nrSols = size(solutions);
1217  if(nrSols > 0) {n = size(solutions[1]);}
1218  else {
1219    n = 0;
1220    nrSols = 1;
1221  }
1222  k = 0;
1223  for(i = 1; i <= nrSols; i++) {
1224    for(j = 1; j <= size(newSolutions); j++) {
1225      k++;
1226      if(n == 0) { temp[1] = newSolutions[j];}
1227      else {
1228        temp = solutions[i];
1229        temp[n + 1] = newSolutions[j];
1230      }
1231      newSols[k] = temp;
1232    }
1233  }
1234  return(newSols);
1235}
1236
1237///////////////////////////////////////////////////////////////////////////////
1238
1239static proc MergeSolutions(list sol1, list index1, list sol2, list index2)
1240"USAGE:   MergeSolutions(sol1, index1, sol2, index2); all parameters are lists
1241RETURN:  list
1242PURPOSE: create a list of solutions of size n, each entry of 'sol2' must
1243         have size n. 'sol1' is one partial solution (from 'LinearZeroSetMain')
1244         'sol2' is a list of partial solutions (from 'NonLinearZeroSetMain')
1245ASSUME:  'sol2' is not empty
1246NOTE:    used by 'ZeroSetMainWork'
1247{
1248  int i, j, k, m;
1249  ideal sol;
1250  list newSols;
1251
1252  m = 0;
1253  for(i = 1; i <= size(sol2); i++) {
1254    m++;
1255    newSols[m] = SubsMapIdeal(sol1 + sol2[i], index1 + index2, 0);
1256  }
1257  return(newSols);
1258}
1259
1260///////////////////////////////////////////////////////////////////////////////
1261
1262static proc SubsMapIdeal(list sol, list index, int opt)
1263"USAGE:   SubsMapIdeal(sol,index,opt); list sol, index; int opt;
1264PURPOSE: built an ideal I as follows.
1265         if i is contained in 'index' then set I[i] = sol[i]
1266         if i is not contained in 'index' then
1267         - opt = 0: set I[i] = 0
1268         - opt = 1: set I[i] = var(i)
1269         if opt = 1 and n = nvars(basering) then set I[n] = var(n).
1270RETURN:  ideal
1271ASSUME:  size(sol) = size(index) <= nvars(basering)
1272"
1273{
1274  int k = 0;
1275  ideal I;
1276  for(int i = 1; i <= nvars(basering) - 1; i = i + 1) {    // built subs. map
1277    if(ContainedQ(index, i)) {
1278      k++;
1279      I[index[k]] = sol[k];
1280    }
1281    else {
1282      if(opt) { I[i] = var(i); }
1283      else { I[i] = 0; }
1284    }
1285  }
1286  if(opt) {I[nvars(basering)] = var(nvars(basering));}
1287  return(I);
1288}
1289
1290///////////////////////////////////////////////////////////////////////////////
1291
1292proc SimplifyZeroset(data)
1293"USAGE:   SimplifyZeroset(data); list data
1294PURPOSE: reduce the entries of the elements of 'data' w.r.t. the ideal 'mpoly'
1295         'data' is a list of ideals/lists.
1296RETURN:  list
1297ASSUME:  basering = Q[x_1,...,x_n,a], order = lp
1298         'data' is a list of ideals
1299         ideal 'mpoly' must be defined, it might be 0 !
1300"
1301{
1302  int i;
1303  list result;
1304
1305  for(i = 1; i <= size(data); i++) {
1306    result[i] = SimplifyData(data[i]);
1307  }
1308  return(result);
1309}
1310
1311///////////////////////////////////////////////////////////////////////////////
1312
1313proc Variables(poly f, int n)
1314"USAGE:   Variables(f,n); poly f; int n;
1315PURPOSE: list of variables among var(1),...,var(n) which occur in f.
1316RETURN:  list
1317ASSUME:  n <= nvars(basering)
1318"
1319{
1320  int i, nrV;
1321  list index;
1322
1323  nrV = 0;
1324  for(i = 1; i <= n; i = i + 1) {
1325    if(diff(f, var(i)) != 0) { nrV++; index[nrV] = i; }
1326  }
1327  return(index);
1328}
1329
1330///////////////////////////////////////////////////////////////////////////////
1331
1332proc ContainedQ(data, f, list #)
1333"USAGE:    ContainedQ(data, f [, opt]); data=list; f=any type, opt=integer
1334PURPOSE:  test if f is an element of data.
1335RETURN:   int
1336          0 if f not contained in data
1337          1 if f contained in data
1338OPTIONS:  opt = 0 : use '==' for comparing f with elements from data@*
1339          opt = 1 : use @code{SameQ} for comparing f with elements from data
1340"
1341{
1342  int opt, i, found;
1343  if(size(#) > 0) { opt = #[1];}
1344  else { opt = 0; }
1345  i = 1;
1346  found = 0;
1347
1348  while((!found) && (i <= size(data))) {
1349    if(opt == 0) {
1350      if(f == data[i]) { found = 1;}
1351      else {i = i + 1;}
1352    }
1353    else {
1354      if(SameQ(f, data[i])) { found = 1;}
1355      else {i = i + 1;}
1356    }
1357  }
1358  return(found);
1359}
1360
1361//////////////////////////////////////////////////////////////////////////////
1362
1363proc SameQ(a, b)
1364"USAGE:    SameQ(a, b); a,b=list/intvec
1365PURPOSE:  test a == b elementwise, i.e., a[i] = b[i].
1366RETURN:   int
1367          0 if a != b
1368          1 if a == b
1369"
1370{
1371  if(typeof(a) == typeof(b)) {
1372    if(typeof(a) == "list" || typeof(a) == "intvec") {
1373      if(size(a) == size(b)) {
1374        int i = 1;
1375        int ok = 1;
1376        while(ok && (i <= size(a))) {
1377          if(a[i] == b[i]) { i = i + 1;}
1378          else {ok = 0;}
1379        }
1380        return(ok);
1381      }
1382      else { return(0); }
1383    }
1384    else { return(a == b);}
1385  }
1386  else { return(0);}
1387}
1388
1389///////////////////////////////////////////////////////////////////////////////
1390
1391static proc SimplifyPoly(poly f)
1392"USAGE:   SimplifyPoly(f); poly f
1393PURPOSE: reduces the coefficients of f w.r.t. the ideal 'moly' if they contain
1394         the algebraic number 'a'.
1395RETURN:  poly
1396ASSUME:  basering = Q[x_1,...,x_n,a]
1397         ideal mpoly must be defined, it might be 0 !
1398"
1399{
1400  matrix coMx;
1401  poly f1, vp;
1402
1403  vp = 1;
1404  for(int i = 1; i < nvars(basering); i++) { vp = vp * var(i);}
1405
1406  coMx = coef(f, vp);
1407  f1 = 0;
1408  for(i = 1; i <= ncols(coMx); i++) {
1409    f1 = f1 + coMx[1, i] * reduce(coMx[2, i], mpoly);
1410  }
1411  return(f1);
1412}
1413
1414///////////////////////////////////////////////////////////////////////////////
1415
1416static proc SimplifyData(data)
1417"USAGE:   SimplifyData(data); ideal/list data;
1418PURPOSE: reduces the entries of 'data' w.r.t. the ideal 'mpoly' if they contain
1419         the algebraic number 'a'
1420RETURN:  ideal/list
1421ASSUME:  basering = Q[x_1,...,x_n,a]
1422         ideal 'mpoly' must be defined, it might be 0 !
1423"
1424{
1425  def altring=basering;
1426  int n;
1427  poly f;
1428
1429  if(typeof(data) == "ideal") { n = ncols(data); }
1430  else { n = size(data);}
1431
1432  for(int i = 1; i <= n; i++) {
1433    f = data[i];
1434    data[i] = SimplifyPoly(f);
1435  }
1436  return(data);
1437}
1438
1439///////////////////////////////////////////////////////////////////////////////
1440
1441static proc TransferRing(R)
1442"USAGE:   TransferRing(R);
1443PURPOSE: creates a new ring containing the same variables as R, but without
1444         parameters. If R contains a parameter then this parameter is added
1445         as the last variable and 'minpoly' is represented by the ideal 'mpoly'
1446         If the basering does not contain a parameter then 'a' is added and
1447         'mpoly' = 0.
1448RETURN:  ring
1449ASSUME:  R = K[x_1,...,x_n] where K = Q or K = Q(a).
1450NOTE:    Creates the ring needed for all prodecures with name 'proc-name'Main
1451"
1452{
1453  def altring=basering;
1454  string ringSTR, parName, minPoly;
1455
1456  setring(R);
1457
1458  if(npars(basering) == 0) {
1459    parName = "a";
1460    minPoly = "0";
1461  }
1462  else {
1463    parName = parstr(basering);
1464    minPoly = string(minpoly);
1465  }
1466  ringSTR = "ring TR = 0, (" + varstr(basering) + "," + parName + "), lp;";
1467
1468  execute(ringSTR);
1469  execute("ideal mpoly = std(" + minPoly + ");");
1470  export(mpoly);
1471  setring altring;
1472  return(TR);
1473}
1474
1475///////////////////////////////////////////////////////////////////////////////
1476
1477static proc NewBaseRing()
1478"USAGE:   NewBaseRing();
1479PURPOSE: creates a new ring, the last variable is added as a parameter.
1480         minpoly is set to mpoly[1].
1481RETURN:  ring
1482ASSUME:  basering = Q[x_1,...,x_n, a], 'mpoly' must be defined
1483"
1484{
1485  int n = nvars(basering);
1486  int MP;
1487  string ringSTR, parName, varString;
1488
1489  def BR = basering;
1490  if(mpoly[1] != 0) {
1491    parName = "(0, " + string(var(n)) + ")";
1492    MP = 1;
1493  }
1494  else {
1495    parName = "0";
1496    MP = 0;
1497  }
1498
1499
1500  for(int i = 1; i < n - 1; i++) {
1501    varString = varString + string(var(i)) + ",";
1502  }
1503  varString = varString + string(var(n-1));
1504
1505  ringSTR = "ring TR = " + parName + ", (" + varString + "), lp;";
1506  execute(ringSTR);
1507  if(MP) { minpoly = number(imap(BR, mpoly)[1]); }
1508  setring BR;
1509  return(TR);
1510}
1511
1512///////////////////////////////////////////////////////////////////////////////
1513
1514/*
1515                           Examples:
1516
1517
1518// order = 20;
1519ring S1 = 0, (s(1..3)), lp;
1520ideal 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);
1521ideal mpoly = std(0);
1522
1523// order = 10
1524ring S2 = 0, (s(1..5)), lp;
1525ideal 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;
1526ideal mpoly = std(0);
1527
1528//order = 126
1529ring S3 =  0, (s(1..5)), lp;
1530ideal 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;
1531ideal mpoly = std(0);
1532
1533// order = 192
1534ring S4 = 0, (s(1..4)), lp;
1535ideal 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;
1536ideal mpoly = std(0);
1537
1538ring R = (0,a), (x,y,z), lp;
1539minpoly = a2 + 1;
1540ideal I1 = x2 - 1/2, a*z - 1, y - 2;
1541ideal I2 = x3 - 1/2, a*z2 - 3, y - 2*a;
1542
1543*/
Note: See TracBrowser for help on using the repository browser.