source: git/Singular/LIB/solve.lib @ 0bc582c

spielwiese
Last change on this file since 0bc582c was 0bc582c, checked in by Frank Seelisch <seelisch@…>, 15 years ago
removed some docu bugs prior to release 3-1-0 git-svn-id: file:///usr/local/Singular/svn/trunk@11624 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 57.8 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: solve.lib,v 1.38 2009-04-06 09:17:01 seelisch Exp $";
3category="Symbolic-numerical solving";
4info="
5LIBRARY: solve.lib     Complex Solving of Polynomial Systems
6AUTHOR:  Moritz Wenk,  email: wenk@mathematik.uni-kl.de
7         Wilfred Pohl, email: pohl@mathematik.uni-kl.de
8
9PROCEDURES:
10laguerre_solve(p,[..]); find all roots of univariate polynomial p
11solve(i,[..]);          all roots of 0-dim. ideal i using triangular sets
12ures_solve(i,[..]);     find all roots of 0-dimensional ideal i with resultants
13mp_res_mat(i,[..]);     multipolynomial resultant matrix of ideal i
14interpolate(p,v,d);     interpolate poly from evaluation points i and results j
15fglm_solve(i,[..]);     find roots of 0-dim. ideal using FGLM and lex_solve
16lex_solve(i,p,[..]);    find roots of reduced lexicographic standard basis
17simplexOut(l);          prints solution of simplex in nice format
18triangLf_solve(l,[..]); find roots using triangular sys. (factorizing Lazard)
19triangM_solve(l,[..]);  find roots of given triangular system (Moeller)
20triangL_solve(l,[..]);  find roots using triangular system (Lazard)
21triang_solve(l,p,[..]); find roots of given triangular system
22";
23
24LIB "triang.lib";    // needed for triang_solve
25
26///////////////////////////////////////////////////////////////////////////////
27
28proc laguerre_solve( poly f, list # )
29"USAGE:   laguerre_solve(f [, m, l, n, s] ); f = polynomial,@*
30         m, l, n, s = integers (control parameters of the method)@*
31 m: precision of output in digits ( 4 <= m), if basering is not ring of
32      complex numbers;
33 l: precision of internal computation in decimal digits ( l >=8 )
34      only if the basering is not complex or complex with smaller precision;@*
35 n: control of multiplicity of roots or of splitting of f into
36      squarefree factors
37      n < 0, no split of f (good, if all roots are simple)
38      n >= 0, try to split
39      n = 0, return only different roots
40      n > 0, find all roots (with multiplicity)
41 s: s != 0, returns ERROR if  | f(root) | > 0.1^m (when computing in the
42      current ring)
43 ( default: m, l, n, s = 8, 30, 1, 0 )
44ASSUME:  f is a univariate polynomial;@*
45         basering has characteristic 0 and is either complex or without
46         parameters.
47RETURN:  list of (complex) roots of the polynomial f, depending on n. The
48         entries of the result are of type@*
49          string: if the basering is not complex,@*
50          number: otherwise.
51NOTE:    If printlevel >0: displays comments ( default = 0 ).
52         If s != 0 and if the procedure stops with ERROR, try a higher
53         internal precision m.
54EXAMPLE: example laguerre_solve; shows an example
55"
56{
57  if (char(basering)!=0){ERROR("characteristic of basering not 0");}
58  if ((charstr(basering)[1]=="0") and (npars(basering)!=0))
59    {ERROR("basering has parameters");}
60  int OLD_COMPLEX=0;
61  int iv=checkv(f);  // check for variable appearing in f
62  if(iv==0){ERROR("Wrong polynomial!");}
63  poly v=var(iv);    // f univariate in v
64
65  int solutionprec=8;// set the control
66  int numberprec=30;
67  int splitcontrol=1;
68  int rootcheck=0;
69  if(size(#)>0){solutionprec=#[1];if(solutionprec<4){solutionprec=4;}}
70  if(size(#)>1){numberprec=#[2];if(numberprec<8){numberprec=8;}}
71  if(solutionprec>numberprec){numberprec=solutionprec;}
72  if(size(#)>2){splitcontrol=#[3];}
73  if(size(#)>3){rootcheck=#[4];}
74  int prot=printlevel-voice+2;
75  int ringprec=0;
76
77  poly p=divzero(f,iv); // divide out zeros as solution
78  int iz=deg(f)-deg(p); // multiplicity of zero solution
79  if(prot!=0)
80  {
81    string pout;
82    string nl=newline;
83    pout="//BEGIN laguerre_solve";
84    if(iz!=0){pout=pout+nl+"//zeros: divide out "+string(iz);}
85    dbprint(prot,pout);
86  }
87  string ss,tt,oo;
88  ss="";oo=ss;
89  if(npars(basering)==1)
90  {
91    if(OLD_COMPLEX)
92    {
93      tt="1,"+string(par(1));
94      if(tt==charstr(basering))
95      {ss=tt;ringprec=system("getPrecDigits");}
96    }
97    else
98    {
99      tt=charstr(basering);
100      if(size(tt)>7)
101      {
102        if(string(tt[1..7])=="complex")
103        {
104          ss=tt;
105          ringprec=system("getPrecDigits");
106        }
107      }
108    }
109  }
110
111  list roots,simple;
112  if(deg(p)==0) // only zero was root
113  {
114    roots=addzero(roots,ss,iz,splitcontrol);
115    if(prot!=0){dbprint(prot,"//END laguerre_solve");}
116    return(roots);
117  }
118
119  if(prot!=0)// more informations
120  {
121    pout="//control: complex ring with precision "+string(numberprec);
122    if(size(ss)==0){pout=pout+nl+
123        "//         basering not complex, hence solutiontype string";
124    if(solutionprec<numberprec){pout=pout+nl+
125        "//         with precision "+string(solutionprec);}}
126    if(splitcontrol<0){pout=pout+nl+ "//       no spliting";}
127    if(splitcontrol==0){pout=pout+nl+"//       output without multiple roots";}
128    if(rootcheck){pout=pout+nl+
129        "//         check roots with precision "+string(solutionprec);}
130    dbprint(prot,pout);
131  }
132
133  def rn = basering;// set the complex ground field
134  if (ringprec<numberprec)
135  {
136    tt="ring lagc=(complex,"+string(numberprec)+","+string(numberprec)+
137       "),"+string(var(iv))+",lp;";
138    execute(tt);
139    poly p=imap(rn,p);
140    poly v=var(1);
141  }
142  int ima=0;
143  if(size(ss)!=0){ima=checkim(p);}
144  number prc=0.1;// set precision of the solution
145  prc=prc^solutionprec;
146  if(prot!=0)
147  {
148    if(ringprec<numberprec){pout="//working in:  "+tt;}
149    if((size(ss)!=0)&&(ima!=0)){pout=pout+nl+
150        "//         polynomial has complex coefficients";}
151    dbprint(prot,pout);
152  }
153
154  int i1=1;
155  int i2=1;
156  ideal SPLIT=p;
157  if(splitcontrol>=0)// splitting
158  {
159    if(prot!=0){dbprint(prot,"//split in working ring:");}
160    SPLIT=splitsqrfree(p,v);
161    i1=size(SPLIT);
162    if((i1==1)&&(charstr(rn)=="0"))
163    {
164      if(prot!=0){dbprint(prot,"//split exact in basering:");}
165      setring rn;
166      if(v>1)
167      {
168        ideal SQQQQ=splitsqrfree(p,v);
169        setring lagc;
170        SPLIT=imap(rn,SQQQQ);
171      }
172      else
173      {
174        oo="ring exa=0,"+string(var(1))+",lp;";
175        execute(oo);
176        ideal SQQQQ=splitsqrfree(imap(rn,p),var(1));
177        setring lagc;
178        SPLIT=imap(exa,SQQQQ);
179        kill exa;
180      }
181      i1=size(SPLIT);
182    }
183    if(prot!=0)
184    {
185      if(i1>1)
186      {
187        int i3=deg(SPLIT[1]);
188        pout="//results of split(the squarefree factors):";
189        if(i3>0){pout=pout+nl+
190           "//  multiplicity "+string(i2)+", degree "+string(i3);}
191        while(i2<i1)
192        {
193          i2++;
194          i3=deg(SPLIT[i2]);
195          if(i3>0){pout=pout+nl+
196             "//  multiplicity "+string(i2)+", degree "+string(i3);}
197        }
198        dbprint(prot,pout);
199        i2=1;
200      }
201      else
202      {
203        if(charstr(rn)=="0"){dbprint(prot,"// polynomial is squarefree");}
204        else{dbprint(prot,"//split without result");}
205      }
206    }
207  }
208
209  p=SPLIT[1];// the first part
210  if(deg(p)>0)
211  {
212    roots=laguerre(p,numberprec,1);// the ring is already complex, hence numberprec is dummy
213    if((size(roots)==0)||(string(roots[1])=="0")){ERROR("laguerre: no roots found");}
214    if(rootcheck){checkroots(p,v,roots,ima,prc);}
215  }
216  while(i2<i1)
217  {
218    i2++;
219    p=SPLIT[i2];// the part with multiplicity i2
220    if(deg(p)>0)
221    {
222      simple=laguerre(p,numberprec,1);
223      if((size(simple)==0)||(string(simple[1])=="0")){ERROR("laguerre: no roots found");}
224      if(rootcheck){checkroots(p,v,simple,ima,prc);}
225      if(splitcontrol==0)// no multiple roots
226      {
227        roots=roots+simple;
228      }
229      else// multiple roots
230      {
231        roots=roots+makemult(simple,i2);
232      }
233    }
234  }
235
236  if((solutionprec<numberprec)&&(size(ss)==0))// shorter output
237  {
238    oo="ring lout=(complex,"+string(solutionprec)+",1),"
239    +string(var(1))+",lp;";
240    execute(oo);
241    list roots=imap(lagc,roots);
242    roots=transroots(roots);
243    if(iz>0){roots=addzero(roots,ss,iz,splitcontrol);}
244    if(prot!=0){dbprint(prot,"//END laguerre_solve");}
245    return(roots);
246  }
247  if(size(ss)==0){roots=transroots(roots);}// transform to string
248  else         // or map in basering
249  {
250    if(ringprec<numberprec)
251    {
252      setring rn;
253      list roots=imap(lagc,roots);
254    }
255  }
256  if(iz>0){roots=addzero(roots,ss,iz,splitcontrol);}
257  if(prot!=0){dbprint(prot,"//END laguerre_solve");}
258  return(roots);
259}
260example
261{
262    "EXAMPLE:";echo=2;
263    // Find all roots of an univariate polynomial using Laguerre's method:
264    ring rs1= 0,(x,y),lp;
265    poly f = 15x5 + x3 + x2 - 10;
266    // 10 digits precision
267    laguerre_solve(f,10);
268
269    // Now with complex coefficients,
270    // internal precision is 30 digits (default)
271    printlevel=2;
272    ring rsc= (real,10,i),x,lp;
273    poly f = (15.4+i*5)*x^5 + (25.0e-2+i*2)*x^3 + x2 - 10*i;
274    list l = laguerre_solve(f);
275    l;
276    // check result, value of substituted poly should be near to zero
277    // remember that l contains a list of strings
278    // in the case of a different ring
279    subst(f,x,l[1]);
280    subst(f,x,l[2]);
281}
282//////////////////////////////////////////////////////////////////////////////
283//                subprocedures for laguerre_solve
284/*
285*  if p depends only on var(i)
286*       returns i
287*  otherwise    0
288*/
289static proc checkv(poly p)
290{
291  int n=nvars(basering);
292  int i=0;
293  int v;
294
295  while (n>0)
296  {
297    if ((p-subst(p,var(n),0))!=0)
298    {
299      i++;
300      if (i>1){return(0);}
301      v=n;
302    }
303    n--;
304  }
305  return(v);
306}
307/*
308*  if p has only real coefficients
309*       returns 0
310*  otherwise    1
311*/
312static proc checkim(poly p)
313{
314  poly q=p;
315
316  while(q!=0)
317  {
318    if(impart(leadcoef(q))!=0){return(1);}
319    q=q-lead(q);
320  }
321  return(0);
322}
323/*
324*  make multiplicity m
325*/
326static proc makemult(list si,int m)
327{
328  int k0=0;
329  int k1=size(si);
330  int k2,k3;
331  number ro;
332  list msi;
333
334  for(k2=1;k2<=k1;k2++)
335  {
336    ro=si[k2];
337    for(k3=m;k3>0;k3--){k0++;msi[k0]=ro;}
338  }
339  return(msi);
340}
341/*
342*  returns 1 for n<1
343*/
344static proc cmp1(number n)
345{
346  number r=repart(n);
347  number i=impart(n);
348  number c=r*r+i*i;
349  if(c>1){return(1);}
350  else{return(0);}
351}
352/*
353*  exact division of polys f/g
354*  (should be internal)
355*/
356static proc exdiv(poly f,poly g,poly v)
357{
358  int d1=deg(f);
359  int d2=deg(g);
360  poly r0=f;
361  poly rf=0;
362  poly h;
363  number n,m;
364
365  m=leadcoef(g);
366  while ((r0!=0)&&(d1>=d2))
367  {
368    n=leadcoef(r0)/m;
369    h=n*v^(d1-d2);
370    rf=rf+h;
371    r0=r0-h*g;
372    d1=deg(r0);
373  }
374  return(cleardenom(rf));
375}
376/*
377*  p is univariant in x
378*  perform a split of p into squarefree factors
379*  such that the returned ideal 'split' consists of
380*  the faktors, i.e.
381*    p = n * product ( split[i]^i ) , n a number
382*/
383static proc splitsqrfree(poly p, poly x)
384{
385  int dd=deg(p);
386  if(dd==1){return(p);}
387  int i=1;
388  int j;
389  ideal h,split;
390  poly high;
391
392  h=interred(ideal(p,diff(p,x)));
393  if(deg(h[1])==0){return(p);}
394  high=h[1];
395  split[1]=exdiv(p,high,x);
396  while(1)
397  {
398    h=interred(ideal(split[i],high));
399    j=deg(h[1]);
400    if(j==0){return(p);}
401    if(deg(h[1])==deg(split[i]))
402    {
403      split=split,split[i];
404      split[i]=1;
405    }
406    else
407    {
408      split[i]=exdiv(split[i],h[1],x);
409      split=split,h[1];
410      dd=dd-deg(split[i])*i;
411    }
412    j=j*(i+1);
413    if(j==dd){break;}
414    if(j>dd){return(p);}
415    high=exdiv(high,h[1],x);
416    if(deg(high)==0){return(p);};
417    i++;
418  }
419  return(split);
420}
421/*
422*  see checkroots
423*/
424static proc nerr(number n,number m)
425{
426  int r;
427  number z=0;
428  number nr=repart(n);
429  number ni=impart(n);
430  if(nr<z){nr=z-nr;}
431  if(ni<z){ni=nr-ni;}
432  else{ni=nr+ni;}
433  if(ni<m){r=0;}
434  else{r=1;}
435  return(r);
436}
437/*
438*  returns ERROR for nerr(p(r[i]))>=pr
439*/
440static proc checkroots(poly p,poly v,list r,int ima,number pr)
441{
442  int i=0;
443  int j;
444  number n,m;
445  ideal li;
446
447  while(i<size(r))
448  {
449    i++;
450    n=r[i];
451    j=cmp1(n);
452    if(j!=0){li[1]=v/n-1;m=1;}
453    else{li[1]=v-n;m=n;}
454    if((ima==0)&&(impart(n)!=0))
455    {
456      i++;
457      n=r[i];
458      if(j!=0){li[1]=li[1]*(v/n-1);}
459      else{li[1]=li[1]*(v-n);m=m*n;}
460    }
461    attrib(li,"isSB",1);
462    n=leadcoef(reduce(p,li));n=n/m;
463    if(n!=0)
464    {if(nerr(n,pr)!=0){ERROR("Unsufficient accuracy!");}}
465  }
466}
467/*
468*  transforms thr result to string
469*/
470static proc transroots(list r)
471{
472  int i=size(r);
473  while (i>0)
474  {
475    r[i]=string(r[i]);
476    i--;
477  }
478  return(r);
479}
480/*
481* returns a poly without zeroroots
482*/
483static proc divzero(poly f,int iv)
484{
485  poly p=f;
486  poly q=p;
487  poly r;
488  while(p==q)
489  {
490    q=p/var(iv);
491    r=q*var(iv);
492    if(r==p){p=q;}
493  }
494  return(p);
495}
496/*
497*  add zeros to solution
498*/
499static proc addzero(list zz,string ss,int iz,int a)
500{
501    int i=1;
502    int j=size(zz);
503
504    if(size(ss)==0){zz[j+1]="0";}
505    else{zz[j+1]=number(0);}
506    if(a==0){return(zz);}
507    while(i<iz)
508    {
509      i++;
510      if(size(ss)==0){zz[j+i]="0";}
511      else{zz[j+i]=number(0);}
512    }
513    return(zz);
514}
515///////////////////////////////////////////////////////////////////////////////
516
517proc solve( ideal G, list # )
518"USAGE:   solve(G [, m, n [, l]] [,\"oldring\"] [,\"nodisplay\"] ); G = ideal,
519         m, n, l = integers (control parameters of the method), outR ring,@*
520         m: precision of output in digits ( 4 <= m) and of the generated ring
521            of complex numbers;
522         n: control of multiplicity
523            n = 0, return all different roots
524            n != 0, find all roots (with multiplicity)
525         l: precision of internal computation in decimal digits ( l >=8 )
526            only if the basering is not complex or complex with smaller
527            precision, @*
528         [default: (m,n,l) = (8,0,30), or if only (m,n) are set explicitly
529          with n!=0, then (m,n,l) = (m,n,60) ]
530ASSUME:  the ideal is 0-dimensional;@*
531         basering has characteristic 0 and is either complex or without
532         parameters;
533RETURN:  (1) If called without the additional parameter @code{\"oldring\"}: @*
534         ring @code{R} with the same number of variables but with complex
535         coefficients (and precision m). @code{R} comes with a list
536         @code{SOL} of numbers, in which complex roots of G are stored: @*
537         * If n  = 0, @code{SOL} is the list of all different solutions, each
538         of them being represented by a list of numbers. @*
539         * If n != 0, @code{SOL} is a list of two list: SOL[i][1] is the list
540         of all different solutions with the multiplicity SOL[i][2].@*
541         SOL is ordered w.r.t. multiplicity (the smallest first). @*
542         (2) If called with the additional parameter @code{\"oldring\"}, the
543         procedure looks for an appropriate ring (at top level) in which
544         the solutions can be stored (interactive). @*
545         The user may then select an appropriate ring and choose a name for
546         the output list in this ring. The list is exported directly to the
547         selected ring and the return value is a string \"result exported to\"
548         + name of the selected ring.
549NOTE:    If the problem is not 0-dim. the procedure stops with ERROR. If the
550         ideal G is not a lexicographic Groebner basis, the lexicographic
551         Groebner basis is computed internally (Hilbert driven).  @*
552         The computed solutions are displayed, unless @code{solve} is called
553         with the additional parameter @code{\"nodisplay\"}.
554EXAMPLE: example solve; shows an example
555"
556{
557// test if basering admissible
558  if (char(basering)!=0){ERROR("characteristic of basering not 0");}
559  if ((charstr(basering)[1]=="0") and (npars(basering)!=0))
560  { ERROR("basering has parameters"); }
561
562// some global settings and control
563  int oldr, nodisp, ii, jj;
564  list LL;
565  int outprec = 8;
566  int mu = 0;
567  int prec = 30;
568  // check additional parameters...
569  if (size(#)>0)
570  {
571    int sofar=1;
572    if (typeof(#[1])=="int")
573    {
574      outprec = #[1];
575      if (outprec<4){outprec = 4;}
576      if (size(#)>1)
577      {
578        if (typeof(#[2])=="int")
579        {
580          mu = #[2];
581          if (size(#)>2)
582          {
583            if (typeof(#[3])=="int")
584            {
585              prec = #[3];
586              if (prec<8){prec = 8;}
587            }
588            else
589            {
590              if(mu!=0){prec = 60;}
591              if (#[3]=="oldring"){ oldr=1; }
592              if (#[3]=="nodisplay"){ nodisp=1; }
593            }
594            sofar=3;
595          }
596        }
597        else
598        {
599          if (#[2]=="oldring"){ oldr=1; }
600          if (#[2]=="nodisplay"){ nodisp=1; }
601        }
602        sofar=2;
603      }
604    }
605    else
606    {
607      if (#[1]=="oldring"){ oldr=1; }
608      if (#[1]=="nodisplay"){ nodisp=1; }
609    }
610    for (ii=sofar+1;ii<=size(#);ii++)
611    { // check for additional strings
612      if (typeof(#[ii])=="string")
613      {
614        if (#[ii]=="oldring"){ oldr=1; }
615        if (#[ii]=="nodisplay"){ nodisp=1; }
616      }
617    }
618  }
619  if (outprec>prec){prec = outprec;}
620  // if interaktive version is chosen -- choice of basering (Top::`outR`)
621  // and name for list of solutions (outL):
622  if (oldr==1)
623  {
624    list Out;
625    LL=names(Top);
626    for (ii=1;ii<=size(LL);ii++)
627    {
628      if (typeof(`LL[ii]`)=="ring")
629      {
630        if (find(charstr(`LL[ii]`),"complex,"+string(outprec)))
631        {
632          jj++;
633          Out[jj]=LL[ii];
634        }
635      }
636    }
637    if (size(Out)>0)
638    {
639      print("// *** You may select between the following rings for storing "+
640            "the list of");
641      print("// *** complex solutions:");
642      Out;
643      print("// *** Enter the number of the chosen ring");
644      print("// ***  (0: none of them => new ring created and returned)");
645      string chosen;
646      while (chosen=="") { chosen=read(""); }
647      execute("def tchosen = "+chosen);
648      if (typeof(tchosen)=="int")
649      {
650        if ((tchosen>0) and (tchosen<=size(Out)))
651        {
652          string outR = Out[tchosen];
653          print("// *** You have chosen the ring "+ outR +". In this ring"
654                +" the following objects");
655          print("//*** are defined:");
656          listvar(Top::`outR`);
657          print("// *** Enter a name for the list of solutions (different "+
658                "from existing names):");
659          string outL;
660          while (outL==""){ outL=read(""); }
661        }
662      }
663    }
664    else
665    {
666      print("No appropriate ring for storing the list of solutions found " +
667             "=> new ring created and returned");
668    }
669    if (not(defined(outR))) { oldr=0; }
670  }
671
672//  string rinC = nameof(basering)+"C";
673  string sord = ordstr(basering);
674  int nv = nvars(basering);
675  def rin = basering;
676  intvec ovec = option(get);
677  option(redSB);
678  option(returnSB);
679  int sb = attrib(G,"isSB");
680  int lp = 0;
681  if (size(sord)==size("C,lp()"+string(nv)))
682  {
683    lp = find(sord,"lp");
684  }
685
686// ERROR
687  if (sb){if (dim(G)!=0){ERROR("ideal not zero-dimensional");}}
688
689// the trivial homogeneous case (unique solution: (0,...0))
690  if (homog(G))
691  {
692    if (sb==0)
693    {
694      execute("ring dphom=("+charstr(rin)+"),("+
695      varstr(rin)+"),dp;");
696      ideal G = std(imap(rin,G));
697      if (dim(G)!=0){ERROR("ideal not zero-dimensional");}
698      int vdG=vdim(G);
699    }
700    if (oldr!=1)
701    {
702      execute("ring rinC =(complex,"+string(outprec)+
703                 "),("+varstr(basering)+"),lp;");
704      list SOL;
705      if (mu==0){SOL[1] = zerolist(nv);}
706      else{SOL[1] = list(zerolist(nv),list(vdG));}
707      export SOL;
708      if (nodisp==0) { print(SOL); }
709      option(set,ovec);
710      dbprint( printlevel-voice+3,"
711// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
712// is stored.
713// To access the list of complex solutions, type (if the name R was assigned
714// to the return value):
715        setring R; SOL; ");
716      return(rinC);
717    }
718    else
719    {
720      setring (Top::`outR`);
721      list SOL;
722      if (mu==0){SOL[1] = zerolist(nv);}
723      else{SOL[1] = list(zerolist(nv),list(vdG));}
724      execute("def "+outL + "=SOL;");
725      execute("export "+outL+";");
726      if (nodisp==0) { print(SOL); }
727      option(set,ovec);
728      kill SOL;
729      return("result exported to "+outR+" as list "+outL);
730    }
731  }
732
733// look for reduced standard basis in lex
734  if (sb*lp==0)
735  {
736    if (sb==0)
737    {
738      execute("ring dphilb=("+charstr(rin)+"),("+ varstr(rin)+"),dp;");
739      ideal G = imap(rin,G);
740      G = std(G);
741      if (dim(G)!=0){ERROR("ideal not zero-dimensional");}
742    }
743    else
744    {
745      def dphilb = basering;
746      G=interred(G);
747      attrib(G,"isSB",1);
748    }
749    execute("ring lexhilb=("+charstr(rin)+"),("+ varstr(rin)+"),lp;");
750    option(redTail);
751    ideal H = fglm(dphilb,G);
752    kill dphilb;
753    H = simplify(H,2);
754    if (lp){setring rin;}
755    else
756    {
757      execute("ring lplex=("+charstr(rin)+"),("+varstr(rin)+"),lp;");
758    }
759    ideal H = imap(lexhilb,H);
760    kill lexhilb;
761  }
762  else{ideal H = interred(G);}
763
764// only 1 variable
765  def hr = basering;
766  if (nv==1)
767  {
768    if ((mu==0) and (charstr(basering)[1]=="0"))
769    { // special case
770      list L = laguerre_solve(H[1],prec,prec,mu,0); // list of strings
771      if (oldr!=1)
772      {
773        execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
774        list SOL;
775        for (ii=1; ii<=size(L); ii++ ) { execute("SOL[ii]=number("+L[ii]+");"); }
776        export SOL;
777        if (nodisp==0) { print(SOL); }
778        option(set,ovec);
779        dbprint( printlevel-voice+3,"
780// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
781// is stored.
782// To access the list of complex solutions, type (if the name R was assigned
783// to the return value):
784        setring R; SOL; ");
785        return(rinC);
786      }
787      else
788      {
789        setring (Top::`outR`);
790        list SOL;
791        for (ii=1; ii<=size(L); ii++ ) { execute("SOL[ii]="+L[ii]+";"); }
792        execute("def "+outL + "=SOL;");
793        execute("export "+outL+";");
794        if (nodisp==0) { print(SOL); }
795        option(set,ovec);
796        kill SOL;
797        return("result exported to "+outR+" as list "+outL);
798      }
799    }
800    else
801    {
802      execute("ring internC=(complex,"+string(prec)+"),("+varstr(basering)+"),lp;");
803      ideal H = imap(hr,H);
804      list sp = splittolist(splitsqrfree(H[1],var(1)));
805      jj = size(sp);
806      while(jj>0)
807      {
808        sp[jj][1] = laguerre(sp[jj][1],prec,1);
809        jj--;
810      }
811      setring hr;
812      if (oldr!=1)
813      {
814        execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
815        list SOL;
816        list sp=imap(internC,sp);
817        if(mu!=0){ SOL=sp; }
818        else
819        {
820          jj = size(sp);
821          SOL=sp[jj][1];
822          while(jj>1)
823          {
824            jj--;
825            SOL = sp[jj][1]+SOL;
826          }
827        }
828        export SOL;
829        if (nodisp==0) { print(SOL); }
830        option(set,ovec);
831        dbprint( printlevel-voice+3,"
832// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
833// is stored.
834// To access the list of complex solutions, type (if the name R was assigned
835// to the return value):
836        setring R; SOL; ");
837        return(rinC);
838      }
839      else
840      {
841        setring (Top::`outR`);
842        list SOL;
843        list sp=imap(internC,sp);
844        if(mu!=0){ SOL=sp; }
845        else
846        {
847          jj = size(sp);
848          SOL=sp[jj][1];
849          while(jj>1)
850          {
851            jj--;
852            SOL = sp[jj][1]+SOL;
853          }
854        }
855        kill sp;
856        execute("def "+outL + "=SOL;");
857        execute("export "+outL+";");
858        if (nodisp==0) { print(SOL); }
859        option(set,ovec);
860        kill SOL;
861        return("result exported to "+outR+" as list "+outL);
862      }
863    }
864  }
865
866// the triangular sets (not univariate case)
867  attrib(H,"isSB",1);
868  if (mu==0)
869  {
870    list sp = triangMH(H); // faster, but destroy multiplicity
871  }
872  else
873  {
874    list sp = triangM(H);
875  }
876
877//   create the complex ring and map the result
878  if (outprec<prec)
879  {
880    execute("ring internC=(complex,"+string(prec)+"),("+varstr(hr)+"),lp;");
881  }
882  else
883  {
884    execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
885  }
886  list triC = imap(hr,sp);
887
888// solve the tridiagonal systems
889  int js = size(triC);
890  list ret1;
891  if (mu==0)
892  {
893    ret1 = trisolve(list(),triC[1],prec);
894    while (js>1)
895    {
896      ret1 = trisolve(list(),triC[js],prec)+ret1;
897      js--;
898    }
899  }
900  else
901  {
902    ret1 = mutrisolve(list(),triC[1],prec);
903    while (js>1)
904    {
905      ret1 = addlist(mutrisolve(list(),triC[js],prec),ret1,1);
906      js--;
907    }
908    ret1 = finalclear(ret1);
909  }
910
911// final computations
912  option(set,ovec);
913  if (outprec==prec)
914  { // we are in ring rinC
915    if (oldr!=1)
916    {
917      list SOL=ret1;
918      export SOL;
919      if (nodisp==0) { print(SOL); }
920      dbprint( printlevel-voice+3,"
921// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
922// is stored.
923// To access the list of complex solutions, type (if the name R was assigned
924// to the return value):
925        setring R; SOL; ");
926      return(rinC);
927    }
928    else
929    {
930      setring (Top::`outR`);
931      list SOL=imap(rinC,ret1);
932      execute("def "+outL + "=SOL;");
933      execute("export "+outL+";");
934      if (nodisp==0) { print(SOL); }
935      kill SOL;
936      return("result exported to "+outR+" as list "+outL);
937    }
938  }
939  else
940  {
941    if (oldr!=1)
942    {
943      execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
944      list SOL=imap(internC,ret1);
945      export SOL;
946      if (nodisp==0) { print(SOL); }
947      dbprint( printlevel-voice+3,"
948// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
949// is stored.
950// To access the list of complex solutions, type (if the name R was assigned
951// to the return value):
952        setring R; SOL; ");
953      return(rinC);
954    }
955    else
956    {
957      setring (Top::`outR`);
958      list SOL=imap(internC,ret1);
959      execute("def "+outL + "=SOL;");
960      execute("export "+outL+";");
961      if (nodisp==0) { print(SOL); }
962      kill SOL;
963      return("result exported to "+outR+" as list "+outL);
964    }
965  }
966}
967example
968{
969    "EXAMPLE:";echo=2;
970    // Find all roots of a multivariate ideal using triangular sets:
971    int d,t,s = 4,3,2 ;
972    int i;
973    ring A=0,x(1..d),dp;
974    poly p=-1;
975    for (i=d; i>0; i--) { p=p+x(i)^s; }
976    ideal I = x(d)^t-x(d)^s+p;
977    for (i=d-1; i>0; i--) { I=x(i)^t-x(i)^s+p,I; }
978    I;
979    // the multiplicity is
980    vdim(std(I));
981    def AC=solve(I,6,0,"nodisplay");  // solutions should not be displayed
982    // list of solutions is stored in AC as the list SOL (default name)
983    setring AC;
984    size(SOL);               // number of different solutions
985    SOL[5];                  // the 5th solution
986    // you must start with char. 0
987    setring A;
988    def AC1=solve(I,6,1,"nodisplay");
989    setring AC1;
990    size(SOL);               // number of different multiplicities
991    SOL[1][1][1];            // a solution with
992    SOL[1][2];               // multiplicity 1
993    SOL[2][1][1];            // a solution with
994    SOL[2][2];               // multiplicity 12
995    // the number of different solutions is equal to
996    size(SOL[1][1])+size(SOL[2][1]);
997    // the number of complex solutions (counted with multiplicities) is
998    size(SOL[1][1])*SOL[1][2]+size(SOL[2][1])*SOL[2][2];
999}
1000//////////////////////////////////////////////////////////////////////////////
1001//                subprocedures for solve
1002
1003
1004/*
1005* return one zero-solution
1006*/
1007static proc zerolist(int nv)
1008{
1009  list ret;
1010  int i;
1011  number o=0;
1012
1013  for (i=nv;i>0;i--){ret[i] = o;}
1014  return(ret);
1015}
1016
1017/* ----------------------- check solution ----------------------- */
1018static proc multsol(list ff, int c)
1019{
1020  int i,j;
1021
1022  i = 0;
1023  j = size(ff);
1024  while (j>0)
1025  {
1026    if(c){i = i+ff[j][2]*size(ff[j][1]);}
1027    else{i = i+size(ff[j][1]);}
1028    j--;
1029  }
1030  return(i);
1031}
1032
1033/*
1034* the inputideal A => zero ?
1035*/
1036static proc checksol(ideal A, list lr)
1037{
1038  int d = nvars(basering);
1039  list ro;
1040  ideal re,h;
1041  int i,j,k;
1042
1043  for (i=size(lr);i>0;i--)
1044  {
1045    ro = lr[i];
1046    for (j=d;j>0;j--)
1047    {
1048      re[j] = var(j)-ro[j];
1049    }
1050    attrib(re,"isSB",1);
1051    k = size(reduce(A,re,1));
1052    if (k){return(i);}
1053  }
1054  return(0);
1055}
1056
1057/*
1058*  compare 2 solutions: returns 0 for equal
1059*/
1060static proc cmpn(list a,list b)
1061{
1062  int ii;
1063
1064  for(ii=size(a);ii>0;ii--){if(a[ii]!=b[ii]) break;}
1065  return(ii);
1066}
1067
1068/*
1069*  delete equal solutions in the list
1070*/
1071static proc delequal(list r, int w)
1072{
1073  list h;
1074  int i,j,k,c;
1075
1076  if (w)
1077  {
1078    k = size(r);
1079    h = r[k][1];
1080    k--;
1081    while (k>0)
1082    {
1083      h = r[k][1]+h;
1084      k--;
1085    }
1086  }
1087  else{h = r;}
1088  k=size(h);
1089  i=1;
1090  while(i<k)
1091  {
1092    j=k;
1093    while(j>i)
1094    {
1095      c=cmpn(h[i],h[j]);
1096      if(c==0)
1097      {
1098        h=delete(h,j);
1099        k--;
1100      }
1101      j--;
1102    }
1103    i++;
1104  }
1105  return(h);
1106}
1107
1108/* ----------------------- substitution ----------------------- */
1109/*
1110* instead of subst(T,var(v),n), much faster
1111*   need option(redSB) !
1112*/
1113static proc linreduce(ideal T, int v, number n)
1114{
1115  ideal re = var(v)-n;
1116  attrib (re,"isSB",1);
1117  return (reduce(T,re));
1118}
1119
1120/* ----------------------- triangular solution ----------------------- */
1121/*
1122* solution of one tridiagonal system T
1123*   with precision prec
1124*   T[1] is univariant in var(1)
1125*   list o is empty for the first call
1126*/
1127static proc trisolve(list o, ideal T, int prec)
1128{
1129  list lroots,ll;
1130  ideal S;
1131  int i,d;
1132
1133  d = size(T);
1134  S = interred(ideal(T[1],diff(T[1],var(d))));
1135  if (deg(S[1]))
1136  {
1137    T[1] = exdiv(T[1],S[1],var(d));
1138  }
1139  ll = laguerre(T[1],prec,1);
1140  for (i=size(ll);i>0;i--){ll[i] = list(ll[i])+o;}
1141  if (d==1){return(ll);}
1142  for (i=size(ll);i>0;i--)
1143  {
1144    S = linreduce(ideal(T[2..d]),d,ll[i][1]);
1145    lroots = trisolve(ll[i],S,prec)+lroots;
1146  }
1147  return(lroots);
1148}
1149
1150/* ------------------- triangular solution (mult) ------------------- */
1151/*
1152*  recompute equal solutions w.r.t. multiplicity
1153*/
1154static proc finalclear(list b)
1155{
1156  list a = b;
1157  list r;
1158  int i,l,ju,j,k,ku,mu,c;
1159
1160// a[i] only
1161  i = 1;
1162  while (i<=size(a))
1163  {
1164    ju = size(a[i][1]);
1165    j = 1;
1166    while (j<=ju)
1167    {
1168      mu = 1;
1169      k = j+1;
1170      while (k<=ju)
1171      {
1172        c = cmpn(a[i][1][j],a[i][1][k]);
1173        if (c==0)
1174        {
1175          a[i][1] = delete(a[i][1],k);
1176          ju--;
1177          mu++;
1178        }
1179        else{k++;}
1180      }
1181      if (mu>1)
1182      {
1183        r[1] = a[i];
1184        r[1][1] = list(a[i][1][j]);
1185        a[i][1] = delete(a[i][1],j);
1186        a = addlist(r,a,mu);
1187        ju--;
1188      }
1189      else{j++;}
1190    }
1191    if (ju==0){a = delete(a,i);}
1192    else{i++;}
1193  }
1194// a[i], a[l]
1195  i = 1;
1196  while (i<size(a))
1197  {
1198    ju = size(a[i][1]);
1199    l = i+1;
1200    while (l<=size(a))
1201    {
1202      ku = size(a[l][1]);
1203      j = 1;
1204      while (j<=ju)
1205      {
1206        mu = 0;
1207        k = 1;
1208        while (k<=ku)
1209        {
1210          c = cmpn(a[i][1][j],a[l][1][k]);
1211          if (c==0)
1212          {
1213            mu = a[i][2]+a[l][2];
1214            r[1] = a[l];
1215            r[1][1] = list(a[l][1][k]);
1216            r[1][2] = mu;
1217            a[l][1] = delete(a[l][1],k);
1218            a = addlist(r,a,1);
1219            ku--;
1220            break;
1221          }
1222          else{k++;}
1223        }
1224        if (mu)
1225        {
1226          a[i][1] = delete(a[i][1],j);
1227          ju--;
1228        }
1229        else{j++;}
1230      }
1231      if (ku){l++;}
1232      else
1233      {
1234        a = delete(a,l);
1235      }
1236    }
1237    if (ju){i++;}
1238    else
1239    {
1240      a = delete(a,i);
1241    }
1242  }
1243  return(a);
1244}
1245
1246/*
1247* convert to list
1248*/
1249static proc splittolist(ideal sp)
1250{
1251  int j = size(sp);
1252  list spl = list(list(sp[j],j));
1253
1254  j--;
1255  while (j>0)
1256  {
1257    if (deg(sp[j]))
1258    {
1259      spl = list(list(sp[j],j))+spl;
1260    }
1261    j--;
1262  }
1263  return(spl);
1264}
1265
1266/*
1267*  multiply the multiplicity
1268*/
1269static proc multlist(list a, int m)
1270{
1271  int i;
1272  for (i=size(a);i>0;i--){a[i][2] = a[i][2]*m;}
1273  return(a);
1274}
1275
1276/*
1277*  a+b w.r.t. to multiplicity as ordering
1278*    (programming like spolys)
1279*/
1280static proc addlist(list a, list b, int m)
1281{
1282  int i,j,k,l,s;
1283  list r = list();
1284
1285  if (m>1){a = multlist(a,m);}
1286  k = size(a);
1287  l = size(b);
1288  i = 1;
1289  j = 1;
1290  while ((i<=k)&&(j<=l))
1291  {
1292    s = a[i][2]-b[j][2];
1293    if (s>=0)
1294    {
1295      r = r+list(b[j]);
1296      j++;
1297      if (s==0)
1298      {
1299        s = size(r);
1300        r[s][1] = r[s][1]+a[i][1];
1301        i++;
1302      }
1303    }
1304    else
1305    {
1306      r = r+list(a[i]);
1307      i++;
1308    }
1309  }
1310  if (i>k)
1311  {
1312    if (j<=l){r = r+list(b[j..l]);}
1313  }
1314  else{r = r+list(a[i..k]);}
1315  return(r);
1316}
1317
1318/*
1319* solution of one tridiagonal system T with multiplicity
1320*   with precision prec
1321*   T[1] is univariant in var(1)
1322*   list o is empty for the first call
1323*/
1324static proc mutrisolve(list o, ideal T, int prec)
1325{
1326  list lroots,ll,sp;
1327  ideal S,h;
1328  int i,d,m,z;
1329
1330  d = size(T);
1331  sp = splittolist(splitsqrfree(T[1],var(d)));
1332  if (d==1){return(l_mutrisolve(sp,o,prec));}
1333  z = size(sp);
1334  while (z>0)
1335  {
1336    m = sp[z][2];
1337    ll = laguerre(sp[z][1],prec,1);
1338    i = size(ll);
1339    while(i>0)
1340    {
1341      h = linreduce(ideal(T[2..d]),d,ll[i]);
1342      if (size(lroots))
1343      {
1344        lroots = addlist(mutrisolve(list(ll[i])+o,h,prec),lroots,m);
1345      }
1346      else
1347      {
1348        lroots = mutrisolve(list(ll[i])+o,h,prec);
1349        if (m>1){lroots=multlist(lroots,m);}
1350      }
1351      i--;
1352    }
1353    z--;
1354  }
1355  return(lroots);
1356}
1357
1358/*
1359*  the last call, we are ready
1360*/
1361static proc l_mutrisolve(list sp, list o, int prec)
1362{
1363  list lroots,ll;
1364  int z,m,i;
1365
1366  z = size(sp);
1367  while (z>0)
1368  {
1369    m = sp[z][2];
1370    ll = laguerre(sp[z][1],prec,1);
1371    for (i=size(ll);i>0;i--){ll[i] = list(ll[i])+o;}
1372    if (size(lroots))
1373    {
1374      lroots = addlist(list(list(ll,m)),lroots,1);
1375    }
1376    else
1377    {
1378      lroots = list(list(ll,m));
1379    }
1380    z--;
1381  }
1382  return(lroots);
1383}
1384///////////////////////////////////////////////////////////////////////////////
1385
1386proc ures_solve( ideal gls, list # )
1387"USAGE:   ures_solve(i [, k, p] ); i = ideal, k, p = integers
1388   k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky, @*
1389   k=1: use resultant matrix of Macaulay which works only for
1390          homogeneous ideals,@*
1391   p>0: defines precision of the long floats for internal computation
1392          if the basering is not complex (in decimal digits),
1393   (default: k=0, p=30)
1394ASSUME:  i is a zerodimensional ideal given by a quadratic system, that is,@*
1395         nvars(basering) = ncols(i) = number of vars actually occurring in i,
1396RETURN:  If the ground field is the field of complex numbers: list of numbers
1397         (the complex roots of the polynomial system i=0). @*
1398         Otherwise: ring @code{R} with the same number of variables but with
1399         complex coefficients (and precision p). @code{R} comes with a list
1400         @code{SOL} of numbers, in which complex roots of the polynomial
1401         system i are stored: @*
1402EXAMPLE: example ures_solve; shows an example
1403"
1404{
1405  int typ=0;// defaults
1406  int prec=30;
1407
1408  if ( size(#) > 0 )
1409  {
1410    typ= #[1];
1411    if ( typ < 0 || typ > 1 )
1412    {
1413      ERROR("Valid values for second parameter k are:
1414                   0: use sparse Resultant (default)
1415                   1: use Macaulay Resultant");
1416    }
1417  }
1418  if ( size(#) > 1 )
1419  {
1420    prec= #[2];
1421    if ( prec < 8 )
1422    {
1423      prec = 8;
1424    }
1425  }
1426
1427  list LL=uressolve(gls,typ,prec,1);
1428  int sizeLL=size(LL);
1429  if (sizeLL==0)
1430  {
1431    dbprint(printlevel-voice+3,"No solution found!");
1432    return(list());
1433  }
1434  if (typeof(LL[1][1])=="string")
1435  {
1436    int ii,jj;
1437    int nv=size(LL[1]);
1438    execute("ring rinC =(complex,"+string(prec)+",I),("
1439                           +varstr(basering)+"),lp;");
1440    list SOL,SOLnew;
1441    for (ii=1; ii<=sizeLL; ii++)
1442    {
1443      SOLnew=list();
1444      for (jj=1; jj<=nv; jj++)
1445      {
1446        execute("SOLnew["+string(jj)+"]="+LL[ii][jj]+";");
1447      }
1448      SOL[ii]=SOLnew;
1449    }
1450    kill SOLnew;
1451    export SOL;
1452    dbprint( printlevel-voice+3,"
1453// 'ures_solve' created a ring, in which a list SOL of numbers (the complex
1454// solutions) is stored.
1455// To access the list of complex solutions, type (if the name R was assigned
1456// to the return value):
1457        setring R; SOL; ");
1458    return(rinC);
1459  }
1460  else
1461  {
1462    return(LL);
1463  }
1464}
1465example
1466{
1467    "EXAMPLE:";echo=2;
1468    // compute the intersection points of two curves
1469    ring rsq = 0,(x,y),lp;
1470    ideal gls=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
1471    def R=ures_solve(gls,0,16);
1472    setring R; SOL;
1473}
1474///////////////////////////////////////////////////////////////////////////////
1475
1476proc mp_res_mat( ideal i, list # )
1477"USAGE:   mp_res_mat(i [, k] ); i ideal, k integer,
1478    k=0: sparse resultant matrix of Gelfand, Kapranov and Zelevinsky,@*
1479    k=1: resultant matrix of Macaulay (k=0 is default)
1480ASSUME:  The number of elements in the input system must be the number of
1481         variables in the basering plus one;
1482         if k=1 then i must be homogeneous.
1483RETURN:  module representing the multipolynomial resultant matrix
1484EXAMPLE: example mp_res_mat; shows an example
1485"
1486{
1487  int typ=0;
1488
1489  if ( size(#) > 0 )
1490  {
1491    typ= #[1];
1492    if ( typ < 0 || typ > 1 )
1493    {
1494      ERROR("Valid values for third parameter are:
1495                   0: sparse resultant (default)
1496                   1: Macaulay resultant");
1497    }
1498  }
1499  return(mpresmat(i,typ));
1500}
1501example
1502{
1503    "EXAMPLE:";echo=2;
1504    // compute resultant matrix in ring with parameters (sparse resultant matrix)
1505    ring rsq= (0,u0,u1,u2),(x1,x2),lp;
1506    ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
1507    module m = mp_res_mat(i);
1508    print(m);
1509    // computing sparse resultant
1510    det(m);
1511
1512    // compute resultant matrix (Macaulay resultant matrix)
1513    ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
1514    ideal h=  homog(imap(rsq,i),x0);
1515    h;
1516
1517    module m = mp_res_mat(h,1);
1518    print(m);
1519    // computing Macaulay resultant (should be the same as above!)
1520    det(m);
1521
1522    // compute numerical sparse resultant matrix
1523    setring rsq;
1524    ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
1525    module mn = mp_res_mat(ir);
1526    print(mn);
1527    // computing sparse resultant
1528    det(mn);
1529}
1530///////////////////////////////////////////////////////////////////////////////
1531
1532proc interpolate( ideal p, ideal w, int d )
1533"USAGE:   interpolate(p,v,d); p,v=ideals of numbers, d=integer
1534ASSUME:  Ground field K is the field of rational numbers, p and v are lists
1535         of elements of the ground field K with p[j] != -1,0,1, size(p) = n
1536         (= number of vars) and size(v)=N=(d+1)^n.
1537RETURN:  poly f, the unique polynomial f of degree n*d with prescribed values
1538         v[i] at the points p(i)=(p[1]^(i-1),..,p[n]^(i-1)), i=1,..,N.
1539NOTE:    mainly useful when n=1, i.e. f is satisfying f(p^(i-1)) = v[i],
1540         i=1..d+1.
1541SEE ALSO: vandermonde.
1542EXAMPLE: example interpolate; shows an example
1543"
1544{
1545    return(vandermonde(p,w,d));
1546}
1547example
1548{
1549    "EXAMPLE:";  echo=2;
1550    ring r1 = 0,(x),lp;
1551    // determine f with deg(f) = 4 and
1552    // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
1553    ideal v=16,0,11376,1046880,85949136;
1554    interpolate( 3, v, 4 );
1555}
1556
1557///////////////////////////////////////////////////////////////////////////////
1558// changed for Singular 3
1559// Return value is now a list: (rlist, rn@)
1560static proc psubst( int d, int dd, int n, list resl,
1561                    ideal fi, int elem, int nv, int prec, int rn@, list rlist)
1562{
1563    //   nv: number of ring variables         (fixed value)
1564    // elem: number of elements in ideal fi   (fixed value)
1565    //   fi: input ideal                      (fixed value)
1566    //   rl: output list of roots
1567    // resl: actual list of roots
1568    //    n:
1569    //   dd: actual element of fi
1570    //    d: actual variable
1571
1572    list LL;
1573    int pdebug;
1574    int olddd=dd;
1575
1576    dbprint(printlevel-voice+2, "// 0 step "+string(dd)+" of "+string(elem) );
1577
1578    if ( dd <= elem )
1579    {
1580        int loop = 1;
1581        int k;
1582        list lsr,lh;
1583        poly ps;
1584        int thedd;
1585
1586        dbprint( printlevel-voice+1,"// 1 dd = "+string(dd) );
1587
1588        thedd=0;
1589        while ( (dd+1 <= elem) && loop )
1590        {
1591            ps= fi[dd+1];
1592
1593            if ( n-1 > 0 )
1594            {
1595                dbprint( printlevel-voice,
1596                    "// 2 ps=fi["+string(dd+1)+"]"+" size="
1597                        +string(size(coeffs(ps,var(n-1))))
1598                        +"  leadexp(ps)="+string(leadexp(ps)) );
1599
1600                if ( size(coeffs(ps,var(n-1))) == 1 )
1601                {
1602                    dd++;
1603                    // hier Leading-Exponent pruefen???
1604                    // oder ist das Poly immer als letztes in der Liste?!?
1605                    // leadexp(ps)
1606                }
1607                else
1608                {
1609                    loop=0;
1610                }
1611            }
1612            else
1613            {
1614                dbprint( printlevel-voice,
1615                    "// 2 ps=fi["+string(dd+1)+"]"+"  leadexp(ps)="
1616                        +string(leadexp(ps)) );
1617                dd++;
1618            }
1619        }
1620        thedd=dd;
1621        ps= fi[thedd];
1622
1623        dbprint( printlevel-voice+1,
1624            "// 3    fi["+string(thedd-1)+"]"+"  leadexp(fi[thedd-1])="
1625                +string(leadexp(fi[thedd-1])) );
1626        dbprint( printlevel-voice+1,
1627            "// 3 ps=fi["+string(thedd)+"]"+"  leadexp(ps)="
1628                +string(leadexp(ps)) );
1629
1630        for ( k= nv; k > nv-d; k-- )
1631        {
1632            dbprint( printlevel-voice,
1633                "// 4 subst(fi["+string(thedd)+"],"
1634                    +string(var(k))+","+string(resl[k])+");" );
1635            ps = subst(ps,var(k),resl[k]);
1636        }
1637
1638        dbprint( printlevel-voice, "// 5 substituted ps="+string(ps) );
1639
1640        if ( ps != 0 )
1641        {
1642            lsr= laguerre_solve( ps, prec, prec, 0 );
1643        }
1644        else
1645        {
1646            dbprint( printlevel-voice+1,"// 30 ps == 0, thats not cool...");
1647            lsr=list(number(0));
1648        }
1649
1650        dbprint( printlevel-voice+1,
1651         "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]" );
1652
1653        if ( size(lsr) > 1 )
1654        {
1655            dbprint( printlevel-voice+1,
1656                "// 10 checking roots found before, range "
1657                    +string(dd-olddd)+" -- "+string(dd) );
1658            dbprint( printlevel-voice+1,
1659                "// 10 thedd = "+string(thedd) );
1660
1661            int i,j,l;
1662            int ls=size(lsr);
1663            int lss;
1664            poly pss;
1665            list nares;
1666            int rroot;
1667            int nares_size;
1668
1669
1670            for ( i = 1; i <= ls; i++ ) // lsr[1..ls]
1671            {
1672                rroot=1;
1673
1674                if ( pdebug>=2 )
1675                {"// 13 root lsr["+string(i)+"] = "+string(lsr[i]);}
1676                for ( l = 0; l <= dd-olddd; l++ )
1677                {
1678                    if ( l+olddd != thedd )
1679                    {
1680                        if ( pdebug>=2 )
1681                        {"// 11 checking ideal element "+string(l+olddd);}
1682                        ps=fi[l+olddd];
1683                        if ( pdebug>=3 )
1684                        {"// 14 ps=fi["+string(l+olddd)+"]";}
1685                        for ( k= nv; k > nv-d; k-- )
1686                        {
1687                            if ( pdebug>=3 )
1688                            {
1689                                "// 11 subst(fi["+string(olddd+l)+"],"
1690                                    +string(var(k))+","+string(resl[k])+");";
1691                            }
1692                            ps = subst(ps,var(k),resl[k]);
1693
1694                        }
1695
1696                        pss=subst(ps,var(k),lsr[i]); // k=nv-d
1697                        if ( pdebug>=3 )
1698                        { "// 15 0 == "+string(pss); }
1699                        if ( pss != 0 )
1700                        {
1701                            if ( system("complexNearZero",
1702                                        leadcoef(pss),
1703                                        prec) )
1704                            {
1705                                if ( pdebug>=2 )
1706                                { "// 16 root "+string(i)+" is a real root"; }
1707                            }
1708                            else
1709                            {
1710                                if ( pdebug>=2 )
1711                                { "// 17 0 == "+string(pss); }
1712                                rroot=0;
1713                            }
1714                        }
1715
1716                    }
1717                }
1718
1719                if ( rroot == 1 ) // add root to list ?
1720                {
1721                    if ( size(nares) > 0 )
1722                    {
1723                        nares=nares[1..size(nares)],lsr[i];
1724                    }
1725                    else
1726                    {
1727                        nares=lsr[i];
1728                    }
1729                    if ( pdebug>=2 )
1730                    { "// 18 added root to list nares"; }
1731                }
1732            }
1733
1734            nares_size=size(nares);
1735            if ( nares_size == 0 )
1736            {
1737                "Numerical problem: No root found...";
1738                "Output may be incorrect!";
1739                nares=list(number(0));
1740            }
1741
1742            if ( pdebug>=1 )
1743            { "// 20 found <"+string(size(nares))+"> roots"; }
1744
1745            for ( i= 1; i <= nares_size; i++ )
1746            {
1747                resl[nv-d]= nares[i];
1748
1749                if ( dd < elem )
1750                {
1751                    if ( i > 1 )
1752                    {
1753                        rn@++;
1754                    }
1755                    LL = psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec,
1756                                 rn@, rlist );
1757                    rlist = LL[1];
1758                    rn@ = LL[2];
1759                }
1760                else
1761                {
1762                   if ( i > 1 ) {  rn@++; }  //bug found by O.Labs
1763                   if ( pdebug>=1 )
1764                   {"// 30_1 <"+string(rn@)+"> "+string(size(resl))+" <-----";}
1765                   if ( pdebug>=2 ){ resl; }
1766                   rlist[rn@]=resl;
1767                }
1768            }
1769        }
1770        else
1771        {
1772            if ( pdebug>=2 )
1773            { "// 21 found root to be: "+string(lsr[1]); }
1774            resl[nv-d]= lsr[1];
1775
1776            if ( dd < elem )
1777            {
1778                LL= psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec,
1779                            rn@, rlist );
1780                rlist = LL[1];
1781                rn@ = LL[2];
1782            }
1783            else
1784            {
1785                if ( pdebug>=1 )
1786                { "// 30_2 <"+string(rn@)+"> "+string(size(resl))+" <-----";}
1787                if ( pdebug>=2 )
1788                { resl; }
1789                rlist[rn@]=resl;
1790            }
1791        }
1792    }
1793    return(list(rlist,rn@));
1794}
1795
1796///////////////////////////////////////////////////////////////////////////////
1797
1798proc fglm_solve( ideal fi, list # )
1799"USAGE:   fglm_solve(i [, p] ); i ideal, p integer
1800ASSUME:  the ground field has char 0.
1801RETURN:  ring @code{R} with the same number of variables but with complex
1802         coefficients (and precision p). @code{R} comes with a list
1803         @code{rlist} of numbers, in which the complex roots of i are stored.@*
1804         p>0: gives precision of complex numbers in decimal digits [default:
1805         p=30].
1806NOTE:    The procedure uses a standard basis of i to determine all complex
1807         roots of i.
1808EXAMPLE: example fglm_solve; shows an example
1809"
1810{
1811    int prec=30;
1812
1813    if ( size(#)>=1  && typeof(#[1])=="int")
1814    {
1815        prec=#[1];
1816    }
1817
1818    def R = lex_solve(stdfglm(fi),prec);
1819    dbprint( printlevel-voice+3,"
1820// 'fglm_solve' created a ring, in which a list rlist of numbers (the
1821// complex solutions) is stored.
1822// To access the list of complex solutions, type (if the name R was assigned
1823// to the return value):
1824        setring R; rlist; ");
1825    return(R);
1826}
1827example
1828{
1829    "EXAMPLE:";echo=2;
1830    ring r = 0,(x,y),lp;
1831    // compute the intersection points of two curves
1832    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
1833    def R = fglm_solve(s,10);
1834    setring R; rlist;
1835}
1836
1837///////////////////////////////////////////////////////////////////////////////
1838
1839proc lex_solve( ideal fi, list # )
1840"USAGE:   lex_solve( i[,p] ); i=ideal, p=integer,
1841  p>0: gives precision of complex numbers in decimal digits (default: p=30).
1842ASSUME:  i is a reduced lexicographical Groebner bases of a zero-dimensional
1843         ideal, sorted by increasing leading terms.
1844RETURN:  ring @code{R} with the same number of variables but with complex
1845         coefficients (and precision p). @code{R} comes with a list
1846         @code{rlist} of numbers, in which the complex roots of i are stored.
1847EXAMPLE: example lex_solve; shows an example
1848"
1849{
1850    int prec=30;
1851    list LL;
1852
1853    if ( size(#)>=1  && typeof(#[1])=="int")
1854    {
1855        prec=#[1];
1856    }
1857
1858    if ( !defined(pdebug) ) { int pdebug; }
1859    def oring= basering;
1860
1861    // change the ground field to complex numbers
1862    string nrings= "ring RC =(complex,"+string(prec)
1863        +"),("+varstr(basering)+"),lp;";
1864    execute(nrings);
1865
1866    // map fi from old to new ring
1867    ideal fi= imap(oring,fi);
1868
1869    int idelem= size(fi);
1870    int nv= nvars(basering);
1871    int i,j,k,lis;
1872    list resl,li;
1873
1874    if ( !defined(rlist) )
1875    {
1876        list rlist;
1877        export rlist;
1878    }
1879
1880    li= laguerre_solve(fi[1],prec,prec,0);
1881    lis= size(li);
1882
1883    dbprint(printlevel-voice+2,"// laguerre found roots: "+string(size(li)));
1884    int rn@;
1885
1886    for ( j= 1; j <= lis; j++ )
1887    {
1888        dbprint(printlevel-voice+1,"// root "+string(j) );
1889        rn@++;
1890        resl[nv]= li[j];
1891        LL = psubst( 1, 2, nv-1, resl, fi, idelem, nv, prec, rn@, rlist );
1892        rlist=LL[1];
1893        rn@=LL[2];
1894    }
1895
1896    dbprint( printlevel-voice+3,"
1897// 'lex_solve' created a ring, in which a list rlist of numbers (the
1898// complex solutions) is stored.
1899// To access the list of complex solutions, type (if the name R was assigned
1900// to the return value):
1901        setring R; rlist; ");
1902
1903    return(RC);
1904}
1905example
1906{
1907    "EXAMPLE:";echo=2;
1908    ring r = 0,(x,y),lp;
1909    // compute the intersection points of two curves
1910    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
1911    def R = lex_solve(stdfglm(s),10);
1912    setring R; rlist;
1913}
1914
1915///////////////////////////////////////////////////////////////////////////////
1916
1917proc triangLf_solve( ideal fi, list # )
1918"USAGE:   triangLf_solve(i [, p] ); i ideal, p integer,
1919         p>0: gives precision of complex numbers in digits (default: p=30).
1920ASSUME:  the ground field has char 0;  i is a zero-dimensional ideal
1921RETURN:  ring @code{R} with the same number of variables but with complex
1922         coefficients (and precision p). @code{R} comes with a list
1923         @code{rlist} of numbers, in which the complex roots of i are stored.
1924NOTE:    The procedure uses a triangular system (Lazard's Algorithm with
1925         factorization) computed from a standard basis to determine
1926         recursively all complex roots of the input ideal i with Laguerre's
1927         algorithm.
1928EXAMPLE: example triangLf_solve; shows an example
1929"
1930{
1931    int prec=30;
1932
1933    if ( size(#)>=1  && typeof(#[1])=="int")
1934    {
1935        prec=#[1];
1936    }
1937
1938    def R=triang_solve(triangLfak(stdfglm(fi)),prec);
1939    dbprint( printlevel-voice+3,"
1940// 'triangLf_solve' created a ring, in which a list rlist of numbers (the
1941// complex solutions) is stored.
1942// To access the list of complex solutions, type (if the name R was assigned
1943// to the return value):
1944        setring R; rlist; ");
1945    return(R);
1946}
1947example
1948{
1949    "EXAMPLE:";echo=2;
1950    ring r = 0,(x,y),lp;
1951    // compute the intersection points of two curves
1952    ideal s = x2 + y2 - 10, x2 + xy + 2y2 - 16;
1953    def R = triangLf_solve(s,10);
1954    setring R; rlist;
1955}
1956
1957///////////////////////////////////////////////////////////////////////////////
1958
1959proc triangM_solve( ideal fi, list # )
1960"USAGE:   triangM_solve(i [, p ] ); i=ideal, p=integer,
1961         p>0: gives precision of complex numbers in digits (default: p=30).
1962ASSUME:  the ground field has char 0;@*
1963         i zero-dimensional ideal
1964RETURN:  ring @code{R} with the same number of variables but with complex
1965         coefficients (and precision p). @code{R} comes with a list
1966         @code{rlist} of numbers, in which the complex roots of i are stored.
1967NOTE:    The procedure uses a triangular system (Moellers Algorithm) computed
1968         from a standard basis  of input ideal i to determine recursively all
1969         complex roots with Laguerre's algorithm.
1970EXAMPLE: example triangM_solve; shows an example
1971"
1972{
1973    int prec=30;
1974
1975    if ( size(#)>=1  && typeof(#[1])=="int")
1976    {
1977        prec=#[1];
1978    }
1979
1980    def R = triang_solve(triangM(stdfglm(fi)),prec);
1981    dbprint( printlevel-voice+3,"
1982// 'triangM_solve' created a ring, in which a list rlist of numbers (the
1983// complex solutions) is stored.
1984// To access the list of complex solutions, type (if the name R was assigned
1985// to the return value):
1986        setring R; rlist; ");
1987    return(R);
1988}
1989example
1990{
1991    "EXAMPLE:";echo=2;
1992    ring r = 0,(x,y),lp;
1993    // compute the intersection points of two curves
1994    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
1995    def R = triangM_solve(s,10);
1996    setring R; rlist;
1997}
1998
1999///////////////////////////////////////////////////////////////////////////////
2000
2001proc triangL_solve( ideal fi, list # )
2002"USAGE:   triangL_solve(i [, p] ); i=ideal, p=integer,@*
2003         p>0: gives precision of complex numbers in digits (default: p=30).
2004ASSUME:  the ground field has char 0; i is a zero-dimensional ideal.
2005RETURN:  ring @code{R} with the same number of variables, but with complex
2006         coefficients (and precision p). @code{R} comes with a list
2007         @code{rlist} of numbers, in which the complex roots of i are stored.
2008NOTE:    The procedure uses a triangular system (Lazard's Algorithm) computed
2009         from a standard basis of input ideal i to determine recursively all
2010         complex roots with Laguerre's algorithm.
2011EXAMPLE: example triangL_solve; shows an example
2012"
2013{
2014    int prec=30;
2015
2016    if ( size(#)>=1  && typeof(#[1])=="int")
2017    {
2018        prec=#[1];
2019    }
2020
2021    def R=triang_solve(triangL(stdfglm(fi)),prec);
2022    dbprint( printlevel-voice+3,"
2023// 'triangL_solve' created a ring, in which a list rlist of numbers (the
2024// complex solutions) is stored.
2025// To access the list of complex solutions, type (if the name R was assigned
2026// to the return value):
2027        setring R; rlist; ");
2028    return(R);
2029}
2030example
2031{
2032    "EXAMPLE:";echo=2;
2033    ring r = 0,(x,y),lp;
2034    // compute the intersection points of two curves
2035    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
2036    def R = triangL_solve(s,10);
2037    setring R; rlist;
2038}
2039
2040
2041///////////////////////////////////////////////////////////////////////////////
2042
2043proc triang_solve( list lfi, int prec, list # )
2044"USAGE:   triang_solve(l,p [,d] ); l=list, p,d=integers@*
2045         l is a list of finitely many triangular systems, such that the union of
2046         their varieties equals the variety of the initial ideal.@*
2047         p>0: gives precision of complex numbers in digits,@*
2048         d>0: gives precision (1<d<p) for near-zero-determination,@*
2049         (default: d=1/2*p).
2050ASSUME:  the ground field has char 0;@*
2051         l was computed using the algorithm of Lazard or the algorithm of
2052         Moeller (see triang.lib).
2053RETURN:  ring @code{R} with the same number of variables, but with complex
2054         coefficients (and precision p). @code{R} comes with a list
2055         @code{rlist} of numbers, in which the complex roots of l are stored.@*
2056EXAMPLE: example triang_solve; shows an example
2057"
2058{
2059    def oring= basering;
2060    list LL;
2061
2062    // change the ground field to complex numbers
2063    string nrings= "ring RC =(complex,"+string(prec)
2064        +",I),("+varstr(basering)+"),lp;";
2065    execute(nrings);
2066
2067    // list with entry 0 (number)
2068    number nn=0;
2069
2070    // set number of digits for zero-comparison of roots
2071    if ( !defined(myCompDigits) )
2072    {
2073        int myCompDigits;
2074    }
2075    if ( size(#)>=1  && typeof(#[1])=="int" )
2076    {
2077        myCompDigits=#[1];
2078    }
2079    else
2080    {
2081        myCompDigits=(system("getPrecDigits"));
2082    }
2083
2084    dbprint( printlevel-voice+2,"// myCompDigits="+string(myCompDigits) );
2085
2086    int idelem;
2087    int nv= nvars(basering);
2088    int i,j,lis;
2089    list resu,li;
2090
2091    if ( !defined(rlist) )
2092    {
2093        list rlist;
2094        export rlist;
2095    }
2096
2097    int rn@=0;
2098
2099    // map the list
2100    list lfi= imap(oring,lfi);
2101    int slfi= size(lfi);
2102
2103    ideal fi;
2104    for ( i= 1; i <= slfi; i++ )
2105    {
2106        // map fi from old to new ring
2107        fi= lfi[i]; //imap(oring,lfi[i]);
2108
2109        idelem= size(fi);
2110
2111        // solve fi[1]
2112        li= laguerre_solve(fi[1],myCompDigits,myCompDigits,0);
2113        lis= size(li);
2114
2115        dbprint( printlevel-voice+2,"// laguerre found roots: "+string(lis) );
2116
2117        for ( j= 1; j <= lis; j++ )
2118        {
2119            dbprint( printlevel-voice+2,"// root "+string(j) );
2120            rn@++;
2121            resu[nv]= li[j];
2122            LL = psubst( 1, 2, nv-1, resu, fi, idelem, nv, myCompDigits,
2123                         rn@, rlist );
2124            rlist = LL[1];
2125            rn@ = LL[2];
2126        }
2127    }
2128
2129    dbprint( printlevel-voice+3,"
2130// 'triang_solve' created a ring, in which a list rlist of numbers (the
2131// complex solutions) is stored.
2132// To access the list of complex solutions, type (if the name R was assigned
2133// to the return value):
2134        setring R; rlist; ");
2135
2136    return(RC);
2137}
2138example
2139{
2140    "EXAMPLE:";echo=2;
2141    ring r = 0,(x,y),lp;
2142    // compute the intersection points of two curves
2143    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
2144    def R=triang_solve(triangLfak(stdfglm(s)),10);
2145    setring R; rlist;
2146}
2147
2148///////////////////////////////////////////////////////////////////////////////
2149
2150proc simplexOut(list l)
2151"USAGE:   simplexOut(l); l list
2152ASSUME:  l is the output of simplex.
2153RETURN:  Nothing. The procedure prints the computed solution of simplex
2154         (as strings) in a nice format.
2155SEE ALSO: simplex
2156EXAMPLE: example simplexOut; shows an example
2157"
2158{
2159  int i,j;
2160  matrix m= l[1];
2161  intvec iposv= l[3];
2162  int icase= l[2];
2163
2164  int cols= ncols(m);
2165  int rows= nrows(m);
2166
2167  int N= l[6];
2168
2169  if ( 1 == icase )  // objective function is unbound
2170  {
2171    "objective function is unbound";
2172    return();
2173  }
2174  if ( -1 == icase )  // no solution satisfies the given constraints
2175  {
2176    "no solution satisfies the given constraints";
2177    return();
2178  }
2179  if ( -2 == icase )  // other error
2180  {
2181    "an error occurred during simplex computation!";
2182    return();
2183  }
2184
2185  for ( i = 1; i <= rows; i++ )
2186  {
2187    if (i == 1)
2188    {
2189      "z = "+string(m[1][1]);
2190    }
2191    else
2192    {
2193      if ( iposv[i-1] <= N )
2194      {
2195        "x"+string(iposv[i-1])+" = "+string(m[i,1]);
2196      }
2197//        else
2198//        {
2199//               "Y"; iposv[i-1]-N+1;
2200//        }
2201    }
2202  }
2203}
2204example
2205{
2206  "EXAMPLE:";echo=2;
2207  ring r = (real,10),(x),lp;
2208
2209  // consider the max. problem:
2210  //
2211  //    maximize  x(1) + x(2) + 3*x(3) - 0.5*x(4)
2212  //
2213  //  with constraints:   x(1) +          2*x(3)          <= 740
2214  //                             2*x(2)          - 7*x(4) <=   0
2215  //                               x(2) -   x(3) + 2*x(4) >=   0.5
2216  //                      x(1) +   x(2) +   x(3) +   x(4)  =   9
2217  //
2218  matrix sm[5][5]=   0, 1, 1, 3,-0.5,
2219                   740,-1, 0,-2, 0,
2220                     0, 0,-2, 0, 7,
2221                   0.5, 0,-1, 1,-2,
2222                     9,-1,-1,-1,-1;
2223
2224  int n = 4;  // number of constraints
2225  int m = 4;  // number of variables
2226  int m1= 2;  // number of <= constraints
2227  int m2= 1;  // number of >= constraints
2228  int m3= 1;  // number of == constraints
2229
2230  list sol=simplex(sm, n, m, m1, m2, m3);
2231  simplexOut(sol);
2232}
2233
2234
2235// local Variables: ***
2236// c-set-style: bsd ***
2237// End: ***
Note: See TracBrowser for help on using the repository browser.