source: git/Singular/LIB/realrad.lib @ 33b509

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