source: git/Singular/LIB/realrad.lib

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