source: git/Singular/LIB/realrad.lib @ 224420b

spielwiese
Last change on this file since 224420b was 224420b, checked in by Frank Seelisch <seelisch@…>, 14 years ago
removed duplicate proc declarations when doing LIB "all.lib"; i.e., no more "redefine" prompts git-svn-id: file:///usr/local/Singular/svn/trunk@12529 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.0 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="real algebra";
4info="
5LIBRARY:  realrad.lib   Computation of real radicals
6AUTHOR :  Silke Spang
7
8OVERVIEW:
9   Algorithms about the computation of the real
10   radical of an arbitary ideal over the rational numbers
11   and transcendetal extensions thereof
12
13PROCEDURES:
14 realpoly(f);     Computes the real part of the univariate polynomial f
15 realzero(j);     Computes the real radical of the zerodimensional ideal j
16 realrad(j);      Computes the real radical of an arbitary ideal over
17                  transcendental extension of the rational numbers
18";
19
20LIB "inout.lib";
21LIB "poly.lib";
22LIB "matrix.lib";
23LIB "general.lib";
24LIB "rootsur.lib";
25LIB "algebra.lib";
26LIB "standard.lib";
27LIB "primdec.lib";
28LIB "elim.lib";
29
30///////////////////////////////////////////////////////////////////////////////
31
32///////////////////////////////////////////////////////////////////////////////
33///////////////////////////////////////////////////////////////////////////////
34//// the main procedure //////////////////////////////////////////////////////
35//////////////////////////////////////////////////////////////////////////////
36proc realrad(ideal id)
37"USAGE: realrad(id), id an ideal of arbitary dimension
38RETURN: the real radical of id
39EXAMPE: example realrad; shows an example"
40{
41
42  def r=basering;
43  int n=nvars(basering);
44  // for faster Groebner basis and dimension compuations
45  string neuring ="ring schnell=("+charstr(r)+"),("+varstr(r)+"),dp;";
46  execute(neuring);
47  def ri=basering;
48
49  list reddim;//reduct dimension to 0
50  list lpar,lvar,sub;//for the ringchange
51  string pari,vari;
52  int i,siz,l,j;
53  string less="list lessvar="+varstr(r)+";";
54  execute(less);
55  ideal id=imap(r,id);
56  l=size(id);
57  for (i=1;i<=l;i++)
58  {
59    id[i]=simplify_gen(id[i]);
60  }
61  id=groebner(id);
62  if (dim(id)<=0)
63  {
64    id=realzero(id);
65    setring r;
66    id=imap(ri,id);
67    return(id);
68  }
69  //sub are the subsets of {x_1,...,x_n}
70  sub=subsets(n);
71  siz=size(sub)-1;//we dont want to localize on all variables
72
73  //for the empty set
74  reddim[1]=zeroreduct(id);
75  reddim[1]=realzero(reddim[1]);
76  for (i=1;i<=siz;i++)
77  {
78
79    lvar=lessvar;
80    lpar=list();
81    l=size(sub[i]);
82    for (j=1;j<=l;j++)
83    {
84      lpar=lpar+list(lvar[sub[i][j]-j+1]);
85      lvar=delete(lvar,sub[i][j]-j+1);
86    }
87    for(j=1;j<=l;j++)//there are l entries in lpar
88    {
89      pari=pari+","+string(lpar[j]);
90    }
91    l=n-l;//there are the remaining n-l entries in lvar
92    for(j=1;j<=l;j++)//there are l entries in lpar
93    {
94      vari=vari+","+string(lvar[j]);
95    }
96    vari=vari[2..size(vari)];
97    neuring="ring neu=("+charstr(r)+pari+"),("+vari+"),dp;";
98    execute(neuring);
99    ideal id=imap(r,id);
100    ideal buffer=zeroreduct(id);
101    buffer=realzero(buffer);
102    setring ri;
103    reddim[i+1]=imap(neu,buffer);
104    kill neu;
105    //compute the intersection of buffer with r
106    reddim[i+1]=contnonloc(reddim[i+1],pari,vari);
107    vari="";
108    pari="";
109  }
110  id=intersect(reddim[1..(siz+1)]);
111  id=timeStd(id,301);//simplify the output
112  setring r;
113  id=imap(ri,id);
114  return(id);
115
116}
117example
118{ "EXAMPLE:"; echo = 2;
119  ring r1=0,(x,y,z),lp;
120  //dimension 0
121  ideal i0=(x2+1)*(x3-2),(y3-2)*(y2+y+1),z3+2;
122  //dimension 1
123  ideal i1=(y3+3y2+y+1)*(y2+4y+4)*(x2+1),(x2+y)*(x2-y2)*(x2+2xy+y2)*(y2+y+1);
124  ideal i=intersect(i0,i1);
125  realrad(i);
126}
127
128
129/*static*/ proc zeroreduct(ideal i)
130"USAGE:zeroreduct(i), i an arbitary ideal
131RETURN: an ideal j of dimension <=0 s.th. i is contained in
132        j and j is contained in i_{Iso} which is the zariski closure
133        of all real isolated points of i
134"
135{
136 list equi;
137 int d,n,di;
138 n=nvars(basering);
139 def r=basering;
140
141 //chance ring to get faster groebner bases computation for dimensions
142
143 string rneu="ring neu=("+charstr(r)+"),("+varstr(r)+"),dp;";
144 execute(rneu);
145 ideal i=imap(r,i);
146
147 i=groebner(i);
148 while (dim(i)> 0)
149 {
150   equi=equidim(i);
151   d=size(equi);
152   equi[d]=radical(equi[d]);
153   di=dim(std(equi[d]));
154   equi[d]=equi[d],minor(jacob(equi[d]),n-di);
155   equi[d]=radical(equi[d]);
156   i=intersect(equi[1..d]);
157   i=groebner(i);
158 }
159
160 setring r;
161 i=imap(neu,i);
162 i=timeStd(i,301);
163 return(i);
164}
165//////////////////////////////////////////////////////////////////////////////
166///////the zero-dimensional case /////////////////////////////////////////////
167//////////////////////////////////////////////////////////////////////////////
168proc realzero(ideal j)
169"USAGE:    realzero(j); a zero-dimensional ideal j
170RETURN:    j: a zero dimensional ideal, which is the real radical
171              of i, if dim(i)=0
172           0: otherwise
173           this acts via
174           primary decomposition (i=1)
175           listdecomp (i=2) or facstd (i=3)
176EXAMPLE:   example realzero; shows an example"
177
178
179{
180 list prim,prepared,nonshape,realu;
181 int r;//counter
182 int l;//number of first polynomial with degree >1 or even
183 l=size(j);
184 for (r=1;r<=l;r++)
185 {
186   j[r]=simplify_gen(j[r]);
187   if (j[r]==1)
188   {
189     return(ideal(1));
190   }
191 }
192 option(redSB);
193 //j=groebner(j);
194 //special case
195 //if (j==1)
196 //{
197 //  return(j);
198 //}
199 if (nvars(basering)==1)
200 {
201   j=groebner(j);
202   j=realpoly(j[1]);
203   return(j);
204 }
205
206
207 //if (dim(j)>0) {return(0);}
208
209 def r_alt=basering;
210 //store the ring
211 //for a ring chance to the ordering lp;
212 execute("ring r_neu =("+charstr(basering)+"),("+varstr(basering)+"),lp;");
213 setring r_neu;
214 ideal boeser,max;
215 prepared[1]=ideal(1);
216 ideal j=imap(r_alt,j);
217 //ideal j=fglm(r_alt,j);
218 prim=primdecGTZ(j);
219 for (r=1;r<=size(prim);r++)
220 {
221   max=prim[r][2];
222   max=groebner(max);
223   realu=prepare_max(max);
224   max=realu[1];
225   if (max!=1)
226   {
227     if (realu[2]==1)
228     {
229       prepared=insert(prepared,max);
230     }
231     else
232     {
233       nonshape=insert(nonshape,max);
234     }
235   }
236 }
237 j=intersect(prepared[1..size(prepared)]);
238
239 //use a variable change into general position to obtain
240 //the shape via radzero
241 if (size(nonshape)>0)
242 {
243   boeser=GeneralPos(nonshape);
244   j=intersect(j,boeser);
245 }
246 j=timeStd(j,301);
247 setring r_alt;
248 j=fetch(r_neu,j);
249 return(j);
250}
251example
252{ "EXAMPLE:"; echo = 2;
253  //in non parametric fields
254  ring r=0,(x,y),dp;
255  ideal i=(y3+3y2+y+1)*(y2+4y+4)*(x2+1),(x2+y)*(x2-y2)*(x2+2xy+y2)*(y2+y+1);
256  realzero(i);
257  ideal j=(y3+3y2+y+1)*(y2-2y+1),(x2+y)*(x2-y2);
258  realzero(j);
259
260  //to get every path
261  ring r1=(0,t),(x,y),lp;
262  ideal m1=x2+1-t,y3+t2;
263  ideal m2=x2+t2+1,y2+t;
264  ideal m3=x2+1-t,y2-t;
265  ideal m4=x^2+1+t,y2-t;
266  ideal i=intersect(m1,m2,m3,m4);
267  realzero(i);
268
269}
270
271static proc GeneralPos(list buffer)
272"USAGE:    GeneralPos(buffer);
273           buffer a list of maximal ideals which failed the prepare_max-test
274RETURN:    j: the intersection of their realradicals
275EXAMPLE:   example radzero; shows no example"
276{
277 def r=basering;
278 int n,ll;
279 //for the mapping in general position
280 map phi,psi;
281 ideal j;
282 ideal jmap=randomLast(20);
283 string ri;
284 intvec @hilb;
285 ideal trans,transprep;// the transformation ideals
286 int nva=nvars(r);
287 int zz,k,l;//counter
288 poly randp;
289 for (zz=1;zz<nva;zz++)
290 {
291   if (npars(basering)>0)
292   {
293     randp=randp+(random(0,5)*par(1)+random(0,5)*par(1)^2+random(0,5))*var(zz);
294   }
295   else
296   {
297     randp=randp+random(0,5)*var(zz);
298   }
299 }
300 randp=randp+var(nva);
301
302 //now they are all irreducible in the non univariate case and
303 //real in the univariate case
304
305 int m=size(buffer);
306 for (l=1;l<=m;l++)
307 {
308   //searching first non univariate polynomial with an even degree
309   //for odd degree we could use the fundamental theorem of algebra and
310   //get real zeros
311
312   //this will act via a coordinate chance into general position
313   //denote that this random chance doesn't work allways
314   //the ideas for the transformation into general position are
315   //used from the primdec.lib
316   transprep=buffer[l];
317   if (voice>=10)
318   {
319     jmap[size(jmap)]=randp;
320   }
321
322
323   for (k=2;k<=n;k++)
324   {
325     if (ord(buffer[l][k])==1)
326     {
327       for (zz=1;zz<=nva;zz++)
328       {
329         if (lead(buffer[l][k])/var(zz)!=0)
330         {
331           transprep[k]=var(zz);
332         }
333       }
334       jmap[nva]=subst(jmap[nva],lead(buffer[l][k]),0);
335     }
336   }
337   phi =r,jmap;
338   for (k=1;k<=nva;k++)
339   {
340     jmap[k]=-(jmap[k]-2*var(k));
341   }
342   psi =r,jmap;
343
344   //coordinate chance
345   trans=phi(transprep);
346
347   //acting with the chanced ideal
348
349   trans=groebner(trans);
350   trans[1]=realpoly(trans[1]);
351
352   //special case
353   if (trans==1)
354   {
355     buffer[l]=trans;
356   }
357   else
358   {
359     ri="ring rhelp=("+charstr(r)+ "),(" +varstr(r)+ ",@t),dp;";
360     execute(ri);
361     ideal trans=homog(imap(r,trans),@t);
362
363     ideal trans1=std(trans);
364     @hilb=hilb(trans1,1);
365     ri= "ring rhelp1=("
366                   +charstr(r)+ "),(" +varstr(rhelp)+ "),lp;";
367     execute(ri);
368     ideal trans=homog(imap(r,trans),@t);
369     kill rhelp;
370     trans=std(trans,@hilb);
371     trans=subst(trans,@t,1);//dehomogenising
372     setring r;
373     trans=imap(rhelp1,trans);
374     kill rhelp1;
375     trans=std(trans);
376     attrib(trans,"isSB",1);
377
378     trans=realzero(trans);
379
380     //going back
381     buffer[l]=psi(trans);
382     buffer[l]=timeStd(buffer[l],301);//timelimit for std computation
383   }
384 }
385 //option(returnSB);
386 j=intersect(buffer[1..m]);
387 return(j);
388
389}
390
391/*proc minAssReal(ideal i, int erg)
392{
393  int l,m,d,e,r,fac;
394  ideal buffer,factor;
395  list minreal;
396  l=size(i);
397  for (r=1;r<=l;r++)
398  {
399    i[r]=simplify_gen(i[r]);
400
401  }
402
403  list pr=primdecGTZ(i);
404  m=size(pr);
405  for (l=1;l<=m;l++)
406  {
407     d=dim(std(pr[l][2]));
408     buffer=realrad(pr[l][2]);
409     buffer=std(buffer);
410     e=dim(buffer);
411     if (d==e)
412     {
413         minreal=minreal+list(pr[l]);
414     }
415  }
416  if (erg==0)
417  {
418   return(minreal);
419  }
420  else
421  {
422     pr=list();
423     m=size(minreal);
424     for (l=1;l<=m;l++)
425     {
426       pr=insert(pr,minreal[l][2]);
427     }
428     i=intersect(pr[1..m]);
429     i=timeStd(i,301);
430     list realmin=minreal+list(i);
431     return(realmin);
432  }
433}*/
434//////////////////////////////////////////////////////////////////////////////
435///////the univariate case ///////////////////////////////////////////////////
436//////////////////////////////////////////////////////////////////////////////
437proc realpoly(poly f)
438"USAGE:    realpoly(f); a univariate polynomial f;
439RETURN:    poly f, where f is the real part of the input f
440EXAMPLE:   example realpoly; shows an example"
441{
442 def r=basering;
443 int tester;
444 if (size(parstr(r))!=0)
445 {
446   string changering="ring rneu=0,("+parstr(r)+","+varstr(r)+"),lp";
447   execute(changering);
448   poly f=imap(r,f);
449   tester=1;
450 }
451 f=simplify(f,1);//wlog f is monic
452 if (f==1)
453 {
454   setring r;
455   return(f);
456 }
457 ideal j=factorize(f,1);//for getting the squarefree factorization
458 poly erg=1;
459 for (int i=1;i<=size(j);i=i+1)
460 {
461  if (is_real(j[i])==1) {erg=erg*j[i];}
462  //we only need real primes
463 }
464 if (tester==1)
465 {
466   setring(r);
467   poly erg=imap(rneu,erg);
468 }
469 return(erg);
470}
471example
472{ "EXAMPLE:"; echo = 2;
473   ring r1 = 0,x,dp;
474   poly f=x5+16x2+x+1;
475   realpoly(f);
476   realpoly(f*(x4+2));
477   ring r2=0,(x,y),dp;
478   poly f=x6-3x4y2 + y6 + x2y2 -6y+5;
479   realpoly(f);
480   ring r3=0,(x,y,z),dp;
481   poly f=x4y4-2x5y3z2+x6y2z4+2x2y3z-4x3y2z3+2x4yz5+z2y2-2z4yx+z6x2;
482   realpoly(f);
483   realpoly(f*(x2+y2+1));
484}
485
486
487
488
489///////////////////////////////////////////////////////////////////////////////
490//// for semi-definiteness/////////////////////////////////////////////////////
491///////////////////////////////////////////////////////////////////////////////
492proc decision(poly f)
493" USAGE: decission(f); a multivariate polynomial f in Q[x_1,..,x_n] and lc f=0
494 RETURN: assume that the basering has a lexicographical ordering,
495         1 if f is positive semidefinite  0 if f is indefinite
496EXAMPLE: decision shows an example
497{
498   string ri,lessvar,parvar,perm;
499   ideal jac;
500   list varlist,buffer,isol,@s,lhelp,lhelp1,lfac,worklist;
501   poly p,g;
502   def rbuffer;
503   def r=basering;
504   //diverse zaehler
505   int @z,zz,count,tester;
506   int n=nvars(r);
507   //specialcases
508
509   if (leadcoef(f)<0)
510   {
511     return(0);
512   }
513   lfac=factorize(f,2);
514   ideal factor=lfac[1];
515   intvec @ex=lfac[2];
516   factor=factor[1];
517   zz=size(factor);
518   f=1;
519   for (@z=1;@z<=zz;@z++)
520   {
521     if ((@ex[@z] mod 2)==1)
522     {
523        f=f*factor[@z];
524     }
525   }
526   if (deg(f)<=0)
527   {
528     if (leadcoef(f)>=0)
529     {
530       return(1);
531     }
532     return(0);
533   }
534   //for recursion
535   if (n==1)
536   {
537    if (sturm(f,-length(f),length(f))==0)
538    {
539      return(1);
540    }
541    return(0);
542   }
543   //search for a p in Q[x_n] such that f is pos. sem. definite
544   //if and only if for every isolating setting S={a_1,...,a_r} holds that
545   //every f(x_1,..,x_n-1, a_i) is positiv semidefinite
546   //recursion of variables
547   ///////////////////////////////////////////////////////////////////////////
548   ///////////////////////////////////////////////////////////////////////////
549   perm="varlist="+varstr(r)+";";
550   execute(perm);
551   varlist=delete(varlist,n);
552   for (@z=1;@z<n;@z++)
553   {
554     lessvar=lessvar+","+string(varlist[@z]);
555   }
556   lessvar=lessvar[2..size(lessvar)];
557   parvar=string(var(n));
558   ri="ring r_neu="+charstr(r)+",(@t,"+parvar+","+lessvar+"),dp;";
559   execute(ri);
560   poly f=imap(r,f);
561   list varlist=imap(r,varlist);
562   ideal jac=jacob(@t+f);
563   jac=jac[3..(n+1)];
564   ideal eins=std(jac);
565   ideal i=@t+f,jac;
566   //use Wu method
567   if (eins==1)
568   {
569     zz=0;
570   }
571   else
572   {
573     matrix m=char_series(i);
574     zz=nrows(m);//number of rows
575   }
576   poly p=1;
577   for (@z=1;@z<=zz;@z++)
578   {
579     p=p*m[@z,1];
580   }
581   //trailing coefficient of p
582   p=subst(p,@t,0);
583   p=realpoly(p);
584   @s=subsets(n-1);
585   ideal jacs;
586   for (@z=1;@z<=size(@s);@z++)
587   {
588     perm="";
589     lhelp=list();
590
591     worklist=varlist;
592     buffer=jac[1..(n-1)];
593     //vorbereitungen fuer den Ringwechsel
594     //setze worklist=x_1,..,x_(n-1)
595
596     for (zz=1;zz<=size(@s[@z]);zz++)
597     {
598      buffer =delete(buffer ,@s[@z][zz]-zz+1);
599      worklist=delete(worklist,@s[@z][zz]-zz+1);
600      lhelp=lhelp+list(string(var(@s[@z][zz]+2)));
601      lhelp1=insert(lhelp,string(var(@s[@z][zz]+2)));
602     }
603     //worklist=(x_1,...,x_n-1)\(x_i1,...,x_ik)
604     //lhelp =(x_i1,...,x_ik)
605     //buffer=diff(f,x_i) i not in (i1,..,ik);
606
607     worklist=list("@t",string(var(2)))+lhelp+worklist;
608     for (zz=1;zz<=n+1;zz++)
609     {
610       perm=perm+","+string(worklist[zz]);
611     }
612     perm=perm[2..size(perm)];
613     if (size(buffer)!=0)
614     {
615      jacs=buffer[1..size(buffer)];
616      jacs=@t+f,jacs;
617     }
618     else
619     {
620      jacs=@t+f;
621     }
622     rbuffer=basering;
623     //perm=@t,x_n,x_1,..,x_ik,x\(x_i1,..,x_ik)
624     ri="ring rh=0,("+perm+"),dp;";
625     execute(ri);
626     ideal jacs=imap(rbuffer,jacs);
627     poly p=imap(rbuffer,p);
628     matrix m=char_series(jacs);
629     poly e=1;
630     for (count=1;count<=nrows(m);count++)
631     {
632       e=e*m[count,1];
633     }
634     //search for the leading coefficient of e in
635     //Q(@t,x_n)[x_@s[@z][1],..,x_@s[@z][size(@s[@z])]
636     intmat l[n-1][n-1];
637     for (zz=1;zz<n;zz++)
638     {
639       l[zz,n-zz]=1;
640     }
641     ri="ring rcoef="+"(0,@t,"+parvar+"),
642     ("+lessvar+"),M(l);";
643     execute(ri);
644     kill l;
645     poly e=imap(rh,e);
646     e=leadcoef(e);
647     setring rh;
648     e=imap(rcoef,e);
649     e=subst(e,@t,0);
650     e=realpoly(e);
651     p=p*e;
652     setring r_neu;
653     p=imap(rh,p);
654     kill rh,rcoef;
655   }
656   setring r;
657   p=imap(r_neu,p);
658   ///////////////////////////////////////////////////////////////////////////
659   ///////////found polynomial p /////////////////////////////////////////////
660   ///////////////////////////////////////////////////////////////////////////
661   //Compute an isolating set for p
662   ri="ring iso="+charstr(r)+","+parvar+",lp;";
663   execute(ri);
664   poly p=imap(r,p);
665   isol=isolset(p);
666   setring r;
667   list isol=imap(iso,isol);
668   tester=1;
669   for (@z=1;@z<=size(isol);@z++)
670   {
671     ri="ring rless="+charstr(r)+",("+lessvar+"),lp;";
672     g=subst(f,var(n),isol[@z]);
673     execute(ri);
674     poly g=imap(r,g);
675     tester=tester*decision(g);
676     setring r;
677     kill rless;
678   }
679   return(tester);
680}
681
682
683proc isolset(poly f)
684"USAGE:  isolset(f); f a univariate polynomial over the rational numbers
685RETURN: An isolating set of f
686NOTE: algorithm can be found in M-F. Roy,R: Pollack, S. Basu page 373
687EXAMPLE: example isolset; shows an example"
688{
689 int i,case;
690 number m;
691 list buffer;
692 //only real roots count
693 f=realpoly(f);
694 poly seppart=f;
695 seppart=simplify(seppart,1);
696 //int N=binlog(length(seppart));
697 //number zweihochN=exp(2,N+1);
698 number zweihochN=length(f);
699 //a special case
700 if (deg(seppart)==0)
701 {
702   return(list(number(0)));
703 }
704 if (sturm(seppart,-zweihochN,zweihochN)==1)
705 {
706  return(list(-zweihochN,zweihochN));
707 }
708 //getting bernstein coeffs
709 ideal id=isuni(f)-zweihochN;
710 map jmap=basering,id;
711 seppart=jmap(seppart);
712
713 id=2*zweihochN*var(1);
714 jmap=basering,id;
715 seppart=jmap(seppart);
716
717 matrix c=coeffs(seppart,var(1));
718 int s=size(c);
719 poly recproc;
720 //Reciprocal polynomial
721 for (i=1;i<=s;i++)
722 {
723  recproc=recproc+c[s+1-i,1]*(var(1)^(i-1));
724 }
725 jmap=basering,var(1)+1;
726 seppart=jmap(recproc);
727 list bernsteincoeffs,bern;
728 c=coeffs(seppart,var(1));
729 for (i=1;i<=s;i++)
730 {
731   bern[i]=number(c[s+1-i,1])/binomial(s-1,i-1);
732 }
733 bernsteincoeffs=bern,list(-zweihochN,zweihochN);
734 list POS;
735 POS[1]=bernsteincoeffs;
736 list L;
737 while (size(POS)!=0)
738 {
739  if (varsigns(POS[1][1])<2)
740  {
741    case=varsigns(POS[1][1]);
742  }
743  else
744  {
745    case=2;
746  }
747  //case Anweisung
748  buffer=POS[1];
749  POS=delete(POS,1);
750  while(1)
751  {
752   if (case==1)
753   {
754    L=L+buffer[2];
755    break;
756   }
757
758   if (case==2)
759   {
760    m=number(buffer[2][1]+buffer[2][2])/2;
761    bern=BernsteinCoefficients(buffer[1],buffer[2],m);
762    POS=bern+POS;
763    if (leadcoef(sign(leadcoef(subst(f,isuni(f),m))))==0)
764    {
765      number epsilon=1/10;
766      while (sturm(f,m-epsilon,m+epsilon)!=1)
767      {
768        epsilon=epsilon/10;
769      }
770      L=L+list(m-epsilon,m+epsilon);
771    }
772    break;
773   }
774  break;
775  }
776 }
777 i=1;
778 while (i<size(L))
779 {
780   if (L[i]==L[i+1])
781   {
782     L=delete(L,i);
783   }
784   else
785   {
786     i=i+1;
787   }
788 }
789 return(L);
790}
791
792static proc BernsteinCoefficients(list bern,list lr,number m)
793"USAGE :BernsteinCoefficients(bern,lr,m);
794        a list bern=b_0,...,b_p representing a polynomial P of degree <=p
795        in the Bernstein basis pf lr=(l,r) an a number m in Q
796 RETURN:a list erg=erg1,erg2 s.th. erg1=erg1[1],erg[2] and erg1[1] are
797        the bernstein coefficients of P w.r.t. to erg1[2]=(l,m) and erg2[1]
798        is one for erg2[2]=(m,r)
799 EXAMPLE: Bernsteincoefficients shows no example
800"
801{
802 //Zaehler
803 int i,j;
804 list erg,erg1,erg2;
805 number a=(lr[2]-m)/(lr[2]-lr[1]);
806 number b=(m-lr[1])/(lr[2]-lr[1]);
807 int p=size(bern);
808 list berns,buffer,buffer2;
809 berns[1]=bern;
810 for (i=2;i<=p;i++)
811 {
812  for (j=1;j<=p+1-i;j++)
813  {
814   buffer[j]=a*berns[i-1][j]+b*berns[i-1][j+1];
815  }
816  berns[i]=buffer;
817  buffer=list();
818 }
819
820 for (i=1;i<=p;i++)
821 {
822  buffer[i]=berns[i][1];
823  buffer2[i]=berns[p+1-i][i];
824 }
825 erg1=buffer,list(lr[1],m);
826 erg2=buffer2,list(m,lr[2]);
827 erg=erg1,erg2;
828 return(erg);
829}
830
831static proc binlog(number i)
832{
833 int erg;
834 if (i<2) {return(0);}
835 else
836 {
837  erg=1+binlog(i/2);
838  return(erg);
839 }
840}
841static proc exp(int a, int b)
842{
843   number res = 1;
844   while (b > 0) {
845   if (b mod 2 == 1) {
846       res = res * a;
847       b = b-1;
848   }
849   a = a * a;
850   b = b div 2;
851   }
852   return(res);
853}
854
855
856
857//////////////////////////////////////////////////////////////////////////////
858///////diverse Hilfsprozeduren ///////////////////////////////////////////////
859//////////////////////////////////////////////////////////////////////////////
860
861/////////////////////////////////////////////////////////////////////////////
862/////wichtig fuers Verstaendnis//////////////////////////////////////////////
863/////////////////////////////////////////////////////////////////////////////
864static proc is_real(poly f)
865"USAGE:     is_real(f);a univariate irreducible polynomial f;
866RETURN:    1: if f is real
867           0: is f is not real
868EXAMPLE:   example is_real; shows an example"
869
870{
871  int d,anz,i;
872  def r=basering;
873
874  if (f==1) {return(1);}
875  if (isuniv(f)==0)
876  {
877   for (i=1;i<=nvars(r);i++)
878   {
879     d=size(coeffs(f,var(i)))+1;
880     if ((d mod 2)==1)
881     {
882       return(1);
883     }
884   }
885   d=1-decision(f);
886   return(d);
887  }
888  d=deg(f) mod 2;
889  if (d==1)
890  {
891     return(1);//because of fundamental theorem of algebra
892  }
893  else
894  {
895   f=simplify(f,1);//wlog we can assume that f is monic
896   number a=leadcoef(sign(leadcoef(subst(f,isuni(f),-length(f)))));
897   number b=leadcoef(sign(leadcoef(subst(f,isuni(f),length(f)))));
898   if
899   (a*b!=1)
900   //polynomials are contineous so the image is an interval
901   //referres to analysis
902   {
903      return(1);
904   }
905   else
906   {
907      anz=sturm(f,-length(f),length(f));
908      if (anz==0) {return(0);}
909      else {return(1);}
910   }
911  }
912}
913example
914{ "EXAMPLE:"; echo = 2;
915   ring r1 = 0,x,dp;
916   poly f=x2+1;
917   is_real(f);
918
919}
920
921
922static proc prepare_max(ideal m)
923"USAGE: prepare_max(m); m a maximal ideal in Q(y_1,...,y_m)[x_1,...,x_n]
924RETURN: a list erg=(id,j); where id is the real radical of m if j=1 (i.e. m
925        satisfies the shape lemma in one variable x_i) else id=m and j=0;
926EXAMPLE: is_in_shape shows an exmaple;
927"
928
929{
930  int j,k,i,l,fakul;
931  def r=basering;
932  int n=nvars(r);
933  list erg,varlist,perm;
934  string wechsler,vari;
935  //option(redSB);
936
937  for (i=1;i<=n;i++)
938  {
939   varlist=varlist+list(var(i));
940  }
941  perm=permutation(varlist);
942  fakul=size(perm);
943  for (i=1;i<=fakul;i++)
944  {
945    for (j=1;j<=n;j++)
946    {
947      vari=vari+","+string(perm[i][j]);
948    }
949    vari=vari[2..size(vari)];
950    wechsler="ring r_neu=("+charstr(r)+"),("+vari+"),lp;";
951    execute(wechsler);
952    ideal id=imap(r,m);
953    id=groebner(id);
954    k=search_first(id,2,2);
955    setring r;
956    m=imap(r_neu,id);
957    m[1]=realpoly(m[1]);
958    if (m[1]==1)
959    {
960      erg[1]=ideal(1);
961      erg[2]=1;
962      return(erg);
963    }
964    if (k>n)
965    {
966      erg[1]=m;
967      erg[2]=1;
968      return(erg);
969    }
970    else
971    {
972      for (l=k;l<=n;l++)
973      {
974        if (realpoly(m[l])==1)
975        {
976          erg[1]=ideal(1);
977          erg[2]=1;
978          return(erg);
979        }
980      }
981    }
982    vari="";
983    kill r_neu;
984  }
985  if (size(parstr(r))==0)
986  {
987    erg[1]=m;
988    j=1;
989    for (i=1;i<=n;i++)
990    {
991      j=j*isuniv(m[i]);
992    }
993    erg[2]=j;
994    return(erg);
995  }
996  erg[1]=m;
997  erg[2]=0;
998  return(erg);
999}
1000
1001static proc length(poly f)
1002"USAGE:    length(f); poly f;
1003RETURN:    sum of the absolute Value of all coeffients of an irreducible
1004           poly nomial f
1005EXAMPLE:   example length; shows an example"
1006
1007{
1008 number erg,buffer;
1009 f=simplify(f,1);//wlog f is monic
1010 int n=size(f);
1011 for (int i=1;i<=n;i=i+1)
1012 {
1013   buffer= leadcoef(f[i]);
1014   erg=erg + absValue(buffer);
1015 }
1016
1017 return(erg);
1018}
1019example
1020{ "EXAMPLE:"; echo = 2;
1021   ring r1 = 0,x,dp;
1022   poly f=x4-6x3+x2+1;
1023   norm(f);
1024
1025   ring r2=0,(x,y),dp;
1026   poly g=x2-y3;
1027   length(g);
1028
1029}
1030//////////////////////////////////////////////////////////////////////////////
1031//////////////weniger wichtig fuers Verstaendnis//////////////////////////////
1032//////////////////////////////////////////////////////////////////////////////
1033static proc isuniv(poly f)
1034{
1035  int erg;
1036  if (f==0)
1037  {
1038    erg=1;
1039  }
1040  else
1041  {
1042  erg=(isuni(f)!=0);
1043  }
1044  return(erg);
1045}
1046static proc search_first(ideal j,int start, int i)
1047"USAGE:    searchfirst(j, start, i);
1048           id a reduced groebner basis w.r.t. lex
1049RETURN:    if i=1 then turns the number of the first non univariate entry
1050           with order >1 in its leading term after start
1051           else the first non univariate of even order
1052EXAMPLE:   example norm; shows no example"
1053{
1054  int n=size(j);
1055  int k=start;//counter
1056  j=j,0;
1057  if (i==1)
1058  {
1059    while
1060    ((k<=n)&&(ord(j[k])==1))
1061    {
1062      k=k+1;
1063    }
1064  }
1065  else
1066  {
1067    while
1068    ((k<=n)&&(ord(j[k]) mod 2==1))
1069    {
1070      k=k+1;
1071    }
1072
1073  }
1074 return(k);
1075}
1076
1077static proc subsets(int n)
1078"USAGE :subsets(n); n>=0 in Z
1079RETURN :l a list of all non-empty subsets of {1,..,n}
1080EXAMPLE:subsets(n) shows an example;
1081"
1082{
1083 list l,buffer;
1084 int i,j,binzahl;
1085 if (n<=0)
1086 {
1087   return(l);
1088 }
1089 int grenze=int(exp(2,n))-1;
1090 for (i=1;i<=grenze;i++)
1091 {
1092  binzahl=i;
1093  for (j=1;j<=n;j++)
1094  {
1095   if ((binzahl mod 2)==1)
1096   {
1097     buffer=buffer+list(j);
1098   }
1099   binzahl=binzahl div 2;
1100  }
1101  l[i]=buffer;
1102  buffer=list();
1103 }
1104 return(l);
1105}
1106example
1107{ "EXAMPLE:"; echo = 2;
1108  subsets(3);
1109  subsets(4);
1110}
1111
1112proc permutation(list L)
1113" USAGE: permutation(L); L a list
1114 OUTPUT: a list of all permutation lists of L
1115EXAMPLE: permutation(L) gives an example"
1116{
1117  list erg,buffer,permi,einfueger;
1118  int i,j,l;
1119  int n=size(L);
1120  if (n==0)
1121  {
1122    return(erg);
1123  }
1124  if (n==1)
1125  {
1126    erg=list(L);
1127    return(erg);
1128  }
1129  for (i=1;i<=n;i++)
1130  {
1131    buffer=delete(L,i);
1132    einfueger=permutation(buffer);
1133    l=size(einfueger);
1134    for (j=1;j<=l;j++)
1135    {
1136      permi=list(L[i])+einfueger[j];
1137      erg=insert(erg,permi);
1138    }
1139  }
1140  return(erg);
1141}
1142example
1143{ "EXAMPLE:"; echo = 2;
1144  list L1="Just","an","example";
1145  permutation(L1);
1146  list L2=1,2,3,4;
1147  permutation(L2);
1148}
1149static proc simplify_gen(poly f)
1150"USAGE : simplify_gen(f); f a polymimial in Q(y_1,..,y_m)[x_1,..,x_n]
1151RETURN : a polynomial g such that g is the square-free part of f  and
1152        every real univariate factor of f is cancelled out
1153EXAMPLE:simplify_gen gives no example"
1154{
1155  int i,l;
1156  ideal factor;
1157  poly g=1;
1158  factor=factorize(f,2)[1];
1159  l=size(factor);
1160  for (i=1;i<=l;i++)
1161  {
1162    if (isuniv(factor[i]))
1163    {
1164        g=g*realpoly(factor[i]);
1165    }
1166    else
1167    {
1168      g=g*factor[i];
1169    }
1170  }
1171  return(g);
1172}
1173static proc contnonloc(ideal id,string pari, string vari)
1174"INPUT : a radical ideal id in in F[pari+vari] which is radical in
1175         F(pari)[vari), pari and vari strings of variables
1176OUTPUT : the contraction ideal of id, i.e. idF(pari)[vari]\cap F[pari+vari]
1177EXAMPLE: contnonloc shows an example
1178"
1179{
1180  list pr;
1181  list contractpr;
1182  int i,l,tester;
1183  ideal primcomp;
1184  def r=basering;
1185  string neu="ring r_neu=("+charstr(r)+pari+"),("+vari+"),dp;";
1186  execute(neu);
1187  def r1=basering;
1188  ideal buffer;
1189  setring r;
1190  pr=primdecGTZ(id);
1191  l=size(pr);
1192  contractpr[1]=ideal(1);
1193  for (i=1;i<=l;i++)
1194  {
1195    primcomp=pr[i][2];
1196    setring r1;
1197    buffer=imap(r,primcomp);
1198    buffer=groebner(buffer);
1199    if (buffer==1)
1200    {
1201     tester=0;
1202    }
1203    else
1204    {
1205     tester=1;
1206    }
1207    setring r;
1208
1209    //id only consits of non units in F(pari)
1210    if (tester==1)
1211    {
1212     contractpr=insert(contractpr,primcomp);
1213    }
1214  }
1215  l=size(contractpr);
1216  id=intersect(contractpr[1..l]);
1217  return(id);
1218}
1219example
1220{ "EXAMPLE:"; echo = 2;
1221   ring  r = 0,(a,b,c),lp;
1222   ideal i=b3+c5,ab2+c3;
1223   ideal j=contnonloc(i,",b","a,c");
1224   j;
1225}
Note: See TracBrowser for help on using the repository browser.