source: git/Singular/LIB/zeroset.lib @ fba86f8

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