source: git/Singular/ipshell.cc @ 2687e2

spielwiese
Last change on this file since 2687e2 was 2687e2, checked in by Hans Schoenemann <hannes@…>, 9 years ago
fix: setting CurrPack for branchTo-calls
  • Property mode set to 100644
File size: 148.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT:
6*/
7
8#include <kernel/mod2.h>
9
10#include <omalloc/omalloc.h>
11
12#include <factory/factory.h>
13
14#include <misc/auxiliary.h>
15#include <misc/options.h>
16#include <misc/mylimits.h>
17#include <misc/intvec.h>
18
19#include <coeffs/numbers.h>
20#include <coeffs/coeffs.h>
21
22#include <coeffs/rmodulon.h>
23#include <coeffs/longrat.h>
24
25#include <polys/monomials/ring.h>
26#include <polys/monomials/maps.h>
27
28#include <polys/prCopy.h>
29#include <polys/matpol.h>
30
31#include <polys/weight.h>
32#include <polys/clapsing.h>
33
34
35#include <polys/ext_fields/algext.h>
36#include <polys/ext_fields/transext.h>
37
38#include <kernel/polys.h>
39#include <kernel/ideals.h>
40
41#include <kernel/numeric/mpr_base.h>
42#include <kernel/numeric/mpr_numeric.h>
43
44#include <kernel/GBEngine/syz.h>
45#include <kernel/GBEngine/kstd1.h>
46
47#include <kernel/combinatorics/stairc.h>
48#include <kernel/combinatorics/hutil.h>
49
50#include <kernel/spectrum/semic.h>
51#include <kernel/spectrum/splist.h>
52#include <kernel/spectrum/spectrum.h>
53
54#include <kernel/oswrapper/feread.h>
55
56#include <Singular/lists.h>
57#include <Singular/attrib.h>
58#include <Singular/ipconv.h>
59#include <Singular/links/silink.h>
60#include <Singular/ipshell.h>
61#include <Singular/maps_ip.h>
62#include <Singular/tok.h>
63#include <Singular/ipid.h>
64#include <Singular/subexpr.h>
65#include <Singular/fevoices.h>
66
67#include <math.h>
68#include <ctype.h>
69
70// define this if you want to use the fast_map routine for mapping ideals
71#define FAST_MAP
72
73#ifdef FAST_MAP
74#include <kernel/maps/fast_maps.h>
75#endif
76
77#ifdef SINGULAR_4_1
78#include <Singular/number2.h>
79#include <coeffs/bigintmat.h>
80#endif
81leftv iiCurrArgs=NULL;
82idhdl iiCurrProc=NULL;
83const char *lastreserved=NULL;
84
85static BOOLEAN iiNoKeepRing=TRUE;
86
87/*0 implementation*/
88
89const char * iiTwoOps(int t)
90{
91  if (t<127)
92  {
93    static char ch[2];
94    switch (t)
95    {
96      case '&':
97        return "and";
98      case '|':
99        return "or";
100      default:
101        ch[0]=t;
102        ch[1]='\0';
103        return ch;
104    }
105  }
106  switch (t)
107  {
108    case COLONCOLON:  return "::";
109    case DOTDOT:      return "..";
110    //case PLUSEQUAL:   return "+=";
111    //case MINUSEQUAL:  return "-=";
112    case MINUSMINUS:  return "--";
113    case PLUSPLUS:    return "++";
114    case EQUAL_EQUAL: return "==";
115    case LE:          return "<=";
116    case GE:          return ">=";
117    case NOTEQUAL:    return "<>";
118    default:          return Tok2Cmdname(t);
119  }
120}
121
122int iiOpsTwoChar(const char *s)
123{
124/* not handling: &&, ||, ** */
125  if (s[1]=='\0') return s[0];
126  else if (s[2]!='\0') return 0;
127  switch(s[0])
128  {
129    case '.': if (s[1]=='.') return DOTDOT;
130              else           return 0;
131    case ':': if (s[1]==':') return COLONCOLON;
132              else           return 0;
133    case '-': if (s[1]=='-') return COLONCOLON;
134              else           return 0;
135    case '+': if (s[1]=='+') return PLUSPLUS;
136              else           return 0;
137    case '=': if (s[1]=='=') return EQUAL_EQUAL;
138              else           return 0;
139    case '<': if (s[1]=='=') return LE;
140              else if (s[1]=='>') return NOTEQUAL;
141              else           return 0;
142    case '>': if (s[1]=='=') return GE;
143              else           return 0;
144    case '!': if (s[1]=='=') return NOTEQUAL;
145              else           return 0;
146  }
147  return 0;
148}
149
150static void list1(const char* s, idhdl h,BOOLEAN c, BOOLEAN fullname)
151{
152  char buffer[22];
153  int l;
154  char buf2[128];
155
156  if(fullname) sprintf(buf2, "%s::%s", "", IDID(h));
157  else sprintf(buf2, "%s", IDID(h));
158
159  Print("%s%-30.30s [%d]  ",s,buf2,IDLEV(h));
160  if (h == currRingHdl) PrintS("*");
161  PrintS(Tok2Cmdname((int)IDTYP(h)));
162
163  ipListFlag(h);
164  switch(IDTYP(h))
165  {
166    case INT_CMD:   Print(" %d",IDINT(h)); break;
167    case INTVEC_CMD:Print(" (%d)",IDINTVEC(h)->length()); break;
168    case INTMAT_CMD:Print(" %d x %d",IDINTVEC(h)->rows(),IDINTVEC(h)->cols());
169                    break;
170    case POLY_CMD:
171    case VECTOR_CMD:if (c)
172                    {
173                      PrintS(" ");wrp(IDPOLY(h));
174                      if(IDPOLY(h) != NULL)
175                      {
176                        Print(", %d monomial(s)",pLength(IDPOLY(h)));
177                      }
178                    }
179                    break;
180    case MODUL_CMD: Print(", rk %d", (int)(IDIDEAL(h)->rank));
181    case IDEAL_CMD: Print(", %u generator(s)",
182                    IDELEMS(IDIDEAL(h))); break;
183    case MAP_CMD:
184                    Print(" from %s",IDMAP(h)->preimage); break;
185    case MATRIX_CMD:Print(" %u x %u"
186                      ,MATROWS(IDMATRIX(h))
187                      ,MATCOLS(IDMATRIX(h))
188                    );
189                    break;
190    case PACKAGE_CMD:
191                    paPrint(IDID(h),IDPACKAGE(h));
192                    break;
193    case PROC_CMD: if((IDPROC(h)->libname!=NULL)
194                   && (strlen(IDPROC(h)->libname)>0))
195                     Print(" from %s",IDPROC(h)->libname);
196                   if(IDPROC(h)->is_static)
197                     PrintS(" (static)");
198                   break;
199    case STRING_CMD:
200                   {
201                     char *s;
202                     l=strlen(IDSTRING(h));
203                     memset(buffer,0,22);
204                     strncpy(buffer,IDSTRING(h),si_min(l,20));
205                     if ((s=strchr(buffer,'\n'))!=NULL)
206                     {
207                       *s='\0';
208                     }
209                     PrintS(" ");
210                     PrintS(buffer);
211                     if((s!=NULL) ||(l>20))
212                     {
213                       Print("..., %d char(s)",l);
214                     }
215                     break;
216                   }
217    case LIST_CMD: Print(", size: %d",IDLIST(h)->nr+1);
218                   break;
219    case QRING_CMD:
220    case RING_CMD:
221                   if ((IDRING(h)==currRing) && (currRingHdl!=h))
222                     PrintS("(*)"); /* this is an alias to currRing */
223#ifdef RDEBUG
224                   if (traceit &TRACE_SHOW_RINGS)
225                     Print(" <%lx>",(long)(IDRING(h)));
226#endif
227                   break;
228#ifdef SINGULAR_4_1
229    case CNUMBER_CMD:
230                   {  number2 n=(number2)IDDATA(h);
231                      Print(" (%s)",nCoeffName(n->cf));
232                      break;
233                   }
234    case CMATRIX_CMD:
235                   {  bigintmat *b=(bigintmat*)IDDATA(h);
236                      Print(" %d x %d (%s)",
237                        b->rows(),b->cols(),
238                        nCoeffName(b->basecoeffs()));
239                      break;
240                   }
241#endif
242    /*default:     break;*/
243  }
244  PrintLn();
245}
246
247void type_cmd(leftv v)
248{
249  BOOLEAN oldShortOut = FALSE;
250
251  if (currRing != NULL)
252  {
253    oldShortOut = currRing->ShortOut;
254    currRing->ShortOut = 1;
255  }
256  int t=v->Typ();
257  Print("// %s %s ",v->Name(),Tok2Cmdname(t));
258  switch (t)
259  {
260    case MAP_CMD:Print(" from %s\n",((map)(v->Data()))->preimage); break;
261    case INTMAT_CMD: Print(" %d x %d\n",((intvec*)(v->Data()))->rows(),
262                                      ((intvec*)(v->Data()))->cols()); break;
263    case MATRIX_CMD:Print(" %u x %u\n" ,
264       MATROWS((matrix)(v->Data())),
265       MATCOLS((matrix)(v->Data())));break;
266    case MODUL_CMD: Print(", rk %d\n", (int)(((ideal)(v->Data()))->rank));break;
267    case LIST_CMD: Print(", size %d\n",((lists)(v->Data()))->nr+1); break;
268
269    case PROC_CMD:
270    case RING_CMD:
271    case IDEAL_CMD:
272    case QRING_CMD: PrintLn(); break;
273
274    //case INT_CMD:
275    //case STRING_CMD:
276    //case INTVEC_CMD:
277    //case POLY_CMD:
278    //case VECTOR_CMD:
279    //case PACKAGE_CMD:
280
281    default:
282      break;
283  }
284  v->Print();
285  if (currRing != NULL)
286    currRing->ShortOut = oldShortOut;
287}
288
289static void killlocals0(int v, idhdl * localhdl, const ring r)
290{
291  idhdl h = *localhdl;
292  while (h!=NULL)
293  {
294    int vv;
295    //Print("consider %s, lev: %d:",IDID(h),IDLEV(h));
296    if ((vv=IDLEV(h))>0)
297    {
298      if (vv < v)
299      {
300        if (iiNoKeepRing)
301        {
302          //PrintS(" break\n");
303          return;
304        }
305        h = IDNEXT(h);
306        //PrintLn();
307      }
308      else //if (vv >= v)
309      {
310        idhdl nexth = IDNEXT(h);
311        killhdl2(h,localhdl,r);
312        h = nexth;
313        //PrintS("kill\n");
314      }
315    }
316    else
317    {
318      h = IDNEXT(h);
319      //PrintLn();
320    }
321  }
322}
323
324void killlocals_rec(idhdl *root,int v, ring r)
325{
326  idhdl h=*root;
327  while (h!=NULL)
328  {
329    if (IDLEV(h)>=v)
330    {
331//      Print("kill %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
332      idhdl n=IDNEXT(h);
333      killhdl2(h,root,r);
334      h=n;
335    }
336    else if (IDTYP(h)==PACKAGE_CMD)
337    {
338 //     Print("into pack %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
339      if (IDPACKAGE(h)!=basePack)
340        killlocals_rec(&(IDRING(h)->idroot),v,r);
341      h=IDNEXT(h);
342    }
343    else if ((IDTYP(h)==RING_CMD)
344    ||(IDTYP(h)==QRING_CMD))
345    {
346      if ((IDRING(h)!=NULL) && (IDRING(h)->idroot!=NULL))
347      // we have to test IDRING(h)!=NULL: qring Q=groebner(...): killlocals
348      {
349  //    Print("into ring %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
350        killlocals_rec(&(IDRING(h)->idroot),v,IDRING(h));
351      }
352      h=IDNEXT(h);
353    }
354    else
355    {
356//      Print("skip %s lev %d for lev %d\n",IDID(h),IDLEV(h),v);
357      h=IDNEXT(h);
358    }
359  }
360}
361BOOLEAN killlocals_list(int v, lists L)
362{
363  if (L==NULL) return FALSE;
364  BOOLEAN changed=FALSE;
365  int n=L->nr;
366  for(;n>=0;n--)
367  {
368    leftv h=&(L->m[n]);
369    void *d=h->data;
370    if (((h->rtyp==RING_CMD) || (h->rtyp==QRING_CMD))
371    && (((ring)d)->idroot!=NULL))
372    {
373      if (d!=currRing) {changed=TRUE;rChangeCurrRing((ring)d);}
374      killlocals0(v,&(((ring)h->data)->idroot),(ring)h->data);
375    }
376    else if (h->rtyp==LIST_CMD)
377      changed|=killlocals_list(v,(lists)d);
378  }
379  return changed;
380}
381void killlocals(int v)
382{
383  BOOLEAN changed=FALSE;
384  idhdl sh=currRingHdl;
385  ring cr=currRing;
386  if (sh!=NULL) changed=((IDLEV(sh)<v) || (IDRING(sh)->ref>0));
387  //if (changed) Print("currRing=%s(%x), lev=%d,ref=%d\n",IDID(sh),IDRING(sh),IDLEV(sh),IDRING(sh)->ref);
388
389  killlocals_rec(&(basePack->idroot),v,currRing);
390
391  if (iiRETURNEXPR_len > myynest)
392  {
393    int t=iiRETURNEXPR.Typ();
394    if ((/*iiRETURNEXPR.Typ()*/ t==RING_CMD)
395    || (/*iiRETURNEXPR.Typ()*/ t==QRING_CMD))
396    {
397      leftv h=&iiRETURNEXPR;
398      if (((ring)h->data)->idroot!=NULL)
399        killlocals0(v,&(((ring)h->data)->idroot),(ring)h->data);
400    }
401    else if (/*iiRETURNEXPR.Typ()*/ t==LIST_CMD)
402    {
403      leftv h=&iiRETURNEXPR;
404      changed |=killlocals_list(v,(lists)h->data);
405    }
406  }
407  if (changed)
408  {
409    currRingHdl=rFindHdl(cr,NULL);
410    if (currRingHdl==NULL)
411      currRing=NULL;
412    else if(cr!=currRing)
413      rChangeCurrRing(cr);
414  }
415
416  if (myynest<=1) iiNoKeepRing=TRUE;
417  //Print("end killlocals  >= %d\n",v);
418  //listall();
419}
420
421void list_cmd(int typ, const char* what, const char *prefix,BOOLEAN iterate, BOOLEAN fullname)
422{
423  package savePack=currPack;
424  idhdl h,start;
425  BOOLEAN all = typ<0;
426  BOOLEAN really_all=FALSE;
427
428  if ( typ==0 )
429  {
430    if (strcmp(what,"all")==0)
431    {
432      if (currPack!=basePack)
433        list_cmd(-1,NULL,prefix,iterate,fullname); // list current package
434      really_all=TRUE;
435      h=basePack->idroot;
436    }
437    else
438    {
439      h = ggetid(what);
440      if (h!=NULL)
441      {
442        if (iterate) list1(prefix,h,TRUE,fullname);
443        if (IDTYP(h)==ALIAS_CMD) PrintS("A");
444        if ((IDTYP(h)==RING_CMD)
445            || (IDTYP(h)==QRING_CMD)
446            //|| (IDTYP(h)==PACKE_CMD)
447        )
448        {
449          h=IDRING(h)->idroot;
450        }
451        else if(IDTYP(h)==PACKAGE_CMD)
452        {
453          currPack=IDPACKAGE(h);
454          //Print("list_cmd:package\n");
455          all=TRUE;typ=PROC_CMD;fullname=TRUE;really_all=TRUE;
456          h=IDPACKAGE(h)->idroot;
457        }
458        else
459        {
460          currPack=savePack;
461          return;
462        }
463      }
464      else
465      {
466        Werror("%s is undefined",what);
467        currPack=savePack;
468        return;
469      }
470    }
471    all=TRUE;
472  }
473  else if (RingDependend(typ))
474  {
475    h = currRing->idroot;
476  }
477  else
478    h = IDROOT;
479  start=h;
480  while (h!=NULL)
481  {
482    if ((all
483      && (IDTYP(h)!=PROC_CMD)
484      &&(IDTYP(h)!=PACKAGE_CMD)
485      && (IDTYP(h)!=CRING_CMD))
486    || (typ == IDTYP(h))
487    || ((typ==RING_CMD) &&(IDTYP(h)==CRING_CMD))
488    || ((IDTYP(h)==QRING_CMD) && (typ==RING_CMD)))
489    {
490      list1(prefix,h,start==currRingHdl, fullname);
491      if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
492        && (really_all || (all && (h==currRingHdl)))
493        && ((IDLEV(h)==0)||(IDLEV(h)==myynest)))
494      {
495        list_cmd(0,IDID(h),"//      ",FALSE);
496      }
497      if (IDTYP(h)==PACKAGE_CMD && really_all)
498      {
499        package save_p=currPack;
500        currPack=IDPACKAGE(h);
501        list_cmd(0,IDID(h),"//      ",FALSE);
502        currPack=save_p;
503      }
504    }
505    h = IDNEXT(h);
506  }
507  currPack=savePack;
508}
509
510void test_cmd(int i)
511{
512  int ii;
513
514  if (i<0)
515  {
516    ii= -i;
517    if (ii < 32)
518    {
519      si_opt_1 &= ~Sy_bit(ii);
520    }
521    else if (ii < 64)
522    {
523      si_opt_2 &= ~Sy_bit(ii-32);
524    }
525    else
526      WerrorS("out of bounds\n");
527  }
528  else if (i<32)
529  {
530    ii=i;
531    if (Sy_bit(ii) & kOptions)
532    {
533      Warn("Gerhard, use the option command");
534      si_opt_1 |= Sy_bit(ii);
535    }
536    else if (Sy_bit(ii) & validOpts)
537      si_opt_1 |= Sy_bit(ii);
538  }
539  else if (i<64)
540  {
541    ii=i-32;
542    si_opt_2 |= Sy_bit(ii);
543  }
544  else
545    WerrorS("out of bounds\n");
546}
547
548int exprlist_length(leftv v)
549{
550  int rc = 0;
551  while (v!=NULL)
552  {
553    switch (v->Typ())
554    {
555      case INT_CMD:
556      case POLY_CMD:
557      case VECTOR_CMD:
558      case NUMBER_CMD:
559        rc++;
560        break;
561      case INTVEC_CMD:
562      case INTMAT_CMD:
563        rc += ((intvec *)(v->Data()))->length();
564        break;
565      case MATRIX_CMD:
566      case IDEAL_CMD:
567      case MODUL_CMD:
568        {
569          matrix mm = (matrix)(v->Data());
570          rc += mm->rows() * mm->cols();
571        }
572        break;
573      case LIST_CMD:
574        rc+=((lists)v->Data())->nr+1;
575        break;
576      default:
577        rc++;
578    }
579    v = v->next;
580  }
581  return rc;
582}
583
584int iiIsPrime0(unsigned p)  /* brute force !!!! */
585{
586  unsigned i,j=0 /*only to avoid compiler warnings*/;
587  if (p<=32749) // max. small prime in factory
588  {
589    int a=0;
590    int e=cf_getNumSmallPrimes()-1;
591    i=e/2;
592    do
593    {
594      j=cf_getSmallPrime(i);
595      if (p==j) return p;
596      if (p<j) e=i-1;
597      else     a=i+1;
598      i=a+(e-a)/2;
599    } while ( a<= e);
600    if (p>j) return j;
601    else     return cf_getSmallPrime(i-1);
602  }
603  unsigned end_i=cf_getNumSmallPrimes()-1;
604  unsigned end_p=(unsigned)sqrt((double)p);
605restart:
606  for (i=0; i<end_i; i++)
607  {
608    j=cf_getSmallPrime(i);
609    if ((p%j) == 0)
610    {
611      if (p<=32751) return iiIsPrime0(p-2);
612      p-=2;
613      goto restart;
614    }
615    if (j > end_p) return p;
616  }
617  if (i>=end_i)
618  {
619    while(j<=end_p)
620    {
621      j+=2;
622      if ((p%j) == 0)
623      {
624        if (p<=32751) return iiIsPrime0(p-2);
625        p-=2;
626        goto restart;
627      }
628    }
629  }
630  return p;
631}
632int IsPrime(int p)  /* brute force !!!! */
633{
634  if      (p == 0)    return 0;
635  else if (p == 1)    return 1/*1*/;
636  else if ((p == 2)||(p==3))    return p;
637  else if (p < 0)     return 2; //(iiIsPrime0((unsigned)(-p)));
638  else if ((p & 1)==0) return iiIsPrime0((unsigned)(p-1));
639  return iiIsPrime0((unsigned)(p));
640}
641
642BOOLEAN iiWRITE(leftv,leftv v)
643{
644  sleftv vf;
645  if (iiConvert(v->Typ(),LINK_CMD,iiTestConvert(v->Typ(),LINK_CMD),v,&vf))
646  {
647    WerrorS("link expected");
648    return TRUE;
649  }
650  si_link l=(si_link)vf.Data();
651  if (vf.next == NULL)
652  {
653    WerrorS("write: need at least two arguments");
654    return TRUE;
655  }
656
657  BOOLEAN b=slWrite(l,vf.next); /* iiConvert preserves next */
658  if (b)
659  {
660    const char *s;
661    if ((l!=NULL)&&(l->name!=NULL)) s=l->name;
662    else                            s=sNoName;
663    Werror("cannot write to %s",s);
664  }
665  vf.CleanUp();
666  return b;
667}
668
669leftv iiMap(map theMap, const char * what)
670{
671  idhdl w,r;
672  leftv v;
673  int i;
674  nMapFunc nMap;
675
676  r=IDROOT->get(theMap->preimage,myynest);
677  if ((currPack!=basePack)
678  &&((r==NULL) || ((r->typ != RING_CMD) && (r->typ != QRING_CMD))))
679    r=basePack->idroot->get(theMap->preimage,myynest);
680  if ((r==NULL) && (currRingHdl!=NULL)
681  && (strcmp(theMap->preimage,IDID(currRingHdl))==0))
682  {
683    r=currRingHdl;
684  }
685  if ((r!=NULL) && ((r->typ == RING_CMD) || (r->typ== QRING_CMD)))
686  {
687    //if ((nMap=nSetMap(rInternalChar(IDRING(r)),
688    //             IDRING(r)->parameter,
689    //             rPar(IDRING(r)),
690    //             IDRING(r)->minpoly)))
691    if ((nMap=n_SetMap(IDRING(r)->cf, currRing->cf))==NULL)
692    {
693////////// WTF?
694//      if (rEqual(IDRING(r),currRing))
695//      {
696//        nMap = n_SetMap(currRing->cf, currRing->cf);
697//      }
698//      else
699//      {
700        Werror("can not map from ground field of %s to current ground field",
701          theMap->preimage);
702        return NULL;
703//      }
704    }
705    if (IDELEMS(theMap)<IDRING(r)->N)
706    {
707      theMap->m=(polyset)omReallocSize((ADDRESS)theMap->m,
708                                 IDELEMS(theMap)*sizeof(poly),
709                                 (IDRING(r)->N)*sizeof(poly));
710      for(i=IDELEMS(theMap);i<IDRING(r)->N;i++)
711        theMap->m[i]=NULL;
712      IDELEMS(theMap)=IDRING(r)->N;
713    }
714    if (what==NULL)
715    {
716      WerrorS("argument of a map must have a name");
717    }
718    else if ((w=IDRING(r)->idroot->get(what,myynest))!=NULL)
719    {
720      char *save_r=NULL;
721      v=(leftv)omAlloc0Bin(sleftv_bin);
722      sleftv tmpW;
723      memset(&tmpW,0,sizeof(sleftv));
724      tmpW.rtyp=IDTYP(w);
725      if (tmpW.rtyp==MAP_CMD)
726      {
727        tmpW.rtyp=IDEAL_CMD;
728        save_r=IDMAP(w)->preimage;
729        IDMAP(w)->preimage=0;
730      }
731      tmpW.data=IDDATA(w);
732#if 0
733      if (((tmpW.rtyp==IDEAL_CMD)||(tmpW.rtyp==MODUL_CMD)) && idIs0(IDIDEAL(w)))
734      {
735        v->rtyp=tmpW.rtyp;
736        v->data=idInit(IDELEMS(IDIDEAL(w)),IDIDEAL(w)->rank);
737      }
738      else
739#endif
740      {
741#ifdef FAST_MAP
742        if ((tmpW.rtyp==IDEAL_CMD) && (nMap == ndCopyMap)
743#ifdef HAVE_PLURAL
744        && (!rIsPluralRing(currRing))
745#endif
746        )
747        {
748          v->rtyp=IDEAL_CMD;
749          v->data=fast_map(IDIDEAL(w), IDRING(r), (ideal)theMap, currRing);
750        }
751        else
752#endif
753        if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,IDRING(r),NULL,NULL,0,nMap))
754        {
755          Werror("cannot map %s(%d)",Tok2Cmdname(w->typ),w->typ);
756          omFreeBin((ADDRESS)v, sleftv_bin);
757          if (save_r!=NULL) IDMAP(w)->preimage=save_r;
758          return NULL;
759        }
760      }
761      if (save_r!=NULL)
762      {
763        IDMAP(w)->preimage=save_r;
764        IDMAP((idhdl)v)->preimage=omStrDup(save_r);
765        v->rtyp=MAP_CMD;
766      }
767      return v;
768    }
769    else
770    {
771      Werror("%s undefined in %s",what,theMap->preimage);
772    }
773  }
774  else
775  {
776    Werror("cannot find preimage %s",theMap->preimage);
777  }
778  return NULL;
779}
780
781#ifdef OLD_RES
782void  iiMakeResolv(resolvente r, int length, int rlen, char * name, int typ0,
783                   intvec ** weights)
784{
785  lists L=liMakeResolv(r,length,rlen,typ0,weights);
786  int i=0;
787  idhdl h;
788  char * s=(char *)omAlloc(strlen(name)+5);
789
790  while (i<=L->nr)
791  {
792    sprintf(s,"%s(%d)",name,i+1);
793    if (i==0)
794      h=enterid(s,myynest,typ0,&(currRing->idroot), FALSE);
795    else
796      h=enterid(s,myynest,MODUL_CMD,&(currRing->idroot), FALSE);
797    if (h!=NULL)
798    {
799      h->data.uideal=(ideal)L->m[i].data;
800      h->attribute=L->m[i].attribute;
801      if (BVERBOSE(V_DEF_RES))
802        Print("//defining: %s as %d-th syzygy module\n",s,i+1);
803    }
804    else
805    {
806      idDelete((ideal *)&(L->m[i].data));
807      Warn("cannot define %s",s);
808    }
809    //L->m[i].data=NULL;
810    //L->m[i].rtyp=0;
811    //L->m[i].attribute=NULL;
812    i++;
813  }
814  omFreeSize((ADDRESS)L->m,(L->nr+1)*sizeof(sleftv));
815  omFreeBin((ADDRESS)L, slists_bin);
816  omFreeSize((ADDRESS)s,strlen(name)+5);
817}
818#endif
819
820//resolvente iiFindRes(char * name, int * len, int *typ0)
821//{
822//  char *s=(char *)omAlloc(strlen(name)+5);
823//  int i=-1;
824//  resolvente r;
825//  idhdl h;
826//
827//  do
828//  {
829//    i++;
830//    sprintf(s,"%s(%d)",name,i+1);
831//    h=currRing->idroot->get(s,myynest);
832//  } while (h!=NULL);
833//  *len=i-1;
834//  if (*len<=0)
835//  {
836//    Werror("no objects %s(1),.. found",name);
837//    omFreeSize((ADDRESS)s,strlen(name)+5);
838//    return NULL;
839//  }
840//  r=(ideal *)omAlloc(/*(len+1)*/ i*sizeof(ideal));
841//  memset(r,0,(*len)*sizeof(ideal));
842//  i=-1;
843//  *typ0=MODUL_CMD;
844//  while (i<(*len))
845//  {
846//    i++;
847//    sprintf(s,"%s(%d)",name,i+1);
848//    h=currRing->idroot->get(s,myynest);
849//    if (h->typ != MODUL_CMD)
850//    {
851//      if ((i!=0) || (h->typ!=IDEAL_CMD))
852//      {
853//        Werror("%s is not of type module",s);
854//        omFreeSize((ADDRESS)r,(*len)*sizeof(ideal));
855//        omFreeSize((ADDRESS)s,strlen(name)+5);
856//        return NULL;
857//      }
858//      *typ0=IDEAL_CMD;
859//    }
860//    if ((i>0) && (idIs0(r[i-1])))
861//    {
862//      *len=i-1;
863//      break;
864//    }
865//    r[i]=IDIDEAL(h);
866//  }
867//  omFreeSize((ADDRESS)s,strlen(name)+5);
868//  return r;
869//}
870
871static resolvente iiCopyRes(resolvente r, int l)
872{
873  int i;
874  resolvente res=(ideal *)omAlloc0((l+1)*sizeof(ideal));
875
876  for (i=0; i<l; i++)
877    res[i]=idCopy(r[i]);
878  return res;
879}
880
881BOOLEAN jjMINRES(leftv res, leftv v)
882{
883  int len=0;
884  int typ0;
885  lists L=(lists)v->Data();
886  intvec *weights=(intvec*)atGet(v,"isHomog",INTVEC_CMD);
887  int add_row_shift = 0;
888  if (weights==NULL)
889    weights=(intvec*)atGet(&(L->m[0]),"isHomog",INTVEC_CMD);
890  if (weights!=NULL)  add_row_shift=weights->min_in();
891  resolvente rr=liFindRes(L,&len,&typ0);
892  if (rr==NULL) return TRUE;
893  resolvente r=iiCopyRes(rr,len);
894
895  syMinimizeResolvente(r,len,0);
896  omFreeSize((ADDRESS)rr,len*sizeof(ideal));
897  len++;
898  res->data=(char *)liMakeResolv(r,len,-1,typ0,NULL,add_row_shift);
899  return FALSE;
900}
901
902BOOLEAN jjBETTI(leftv res, leftv u)
903{
904  sleftv tmp;
905  memset(&tmp,0,sizeof(tmp));
906  tmp.rtyp=INT_CMD;
907  tmp.data=(void *)1;
908  if ((u->Typ()==IDEAL_CMD)
909  || (u->Typ()==MODUL_CMD))
910    return jjBETTI2_ID(res,u,&tmp);
911  else
912    return jjBETTI2(res,u,&tmp);
913}
914
915BOOLEAN jjBETTI2_ID(leftv res, leftv u, leftv v)
916{
917  lists l=(lists) omAllocBin(slists_bin);
918  l->Init(1);
919  l->m[0].rtyp=u->Typ();
920  l->m[0].data=u->Data();
921  attr *a=u->Attribute();
922  if (a!=NULL)
923  l->m[0].attribute=*a;
924  sleftv tmp2;
925  memset(&tmp2,0,sizeof(tmp2));
926  tmp2.rtyp=LIST_CMD;
927  tmp2.data=(void *)l;
928  BOOLEAN r=jjBETTI2(res,&tmp2,v);
929  l->m[0].data=NULL;
930  l->m[0].attribute=NULL;
931  l->m[0].rtyp=DEF_CMD;
932  l->Clean();
933  return r;
934}
935
936BOOLEAN jjBETTI2(leftv res, leftv u, leftv v)
937{
938  resolvente r;
939  int len;
940  int reg,typ0;
941  lists l=(lists)u->Data();
942
943  intvec *weights=NULL;
944  int add_row_shift=0;
945  intvec *ww=(intvec *)atGet(&(l->m[0]),"isHomog",INTVEC_CMD);
946  if (ww!=NULL)
947  {
948     weights=ivCopy(ww);
949     add_row_shift = ww->min_in();
950     (*weights) -= add_row_shift;
951  }
952  //Print("attr:%x\n",weights);
953
954  r=liFindRes(l,&len,&typ0);
955  if (r==NULL) return TRUE;
956  res->data=(char *)syBetti(r,len,&reg,weights,(int)(long)v->Data());
957  omFreeSize((ADDRESS)r,(len)*sizeof(ideal));
958  atSet(res,omStrDup("rowShift"),(void*)(long)add_row_shift,INT_CMD);
959  if (weights!=NULL) delete weights;
960  return FALSE;
961}
962
963int iiRegularity(lists L)
964{
965  int len,reg,typ0;
966
967  resolvente r=liFindRes(L,&len,&typ0);
968
969  if (r==NULL)
970    return -2;
971  intvec *weights=NULL;
972  int add_row_shift=0;
973  intvec *ww=(intvec *)atGet(&(L->m[0]),"isHomog",INTVEC_CMD);
974  if (ww!=NULL)
975  {
976     weights=ivCopy(ww);
977     add_row_shift = ww->min_in();
978     (*weights) -= add_row_shift;
979  }
980  //Print("attr:%x\n",weights);
981
982  intvec *dummy=syBetti(r,len,&reg,weights);
983  if (weights!=NULL) delete weights;
984  delete dummy;
985  omFreeSize((ADDRESS)r,len*sizeof(ideal));
986  return reg+1+add_row_shift;
987}
988
989BOOLEAN iiDebugMarker=TRUE;
990#define BREAK_LINE_LENGTH 80
991void iiDebug()
992{
993  Print("\n-- break point in %s --\n",VoiceName());
994  if (iiDebugMarker) VoiceBackTrack();
995  char * s;
996  iiDebugMarker=FALSE;
997  s = (char *)omAlloc(BREAK_LINE_LENGTH+4);
998  loop
999  {
1000    memset(s,0,80);
1001    fe_fgets_stdin("",s,BREAK_LINE_LENGTH);
1002    if (s[BREAK_LINE_LENGTH-1]!='\0')
1003    {
1004      Print("line too long, max is %d chars\n",BREAK_LINE_LENGTH);
1005    }
1006    else
1007      break;
1008  }
1009  if (*s=='\n')
1010  {
1011    iiDebugMarker=TRUE;
1012  }
1013#if MDEBUG
1014  else if(strncmp(s,"cont;",5)==0)
1015  {
1016    iiDebugMarker=TRUE;
1017  }
1018#endif /* MDEBUG */
1019  else
1020  {
1021    strcat( s, "\n;~\n");
1022    newBuffer(s,BT_execute);
1023  }
1024}
1025
1026lists scIndIndset(ideal S, BOOLEAN all, ideal Q)
1027{
1028  int i;
1029  indset save;
1030  lists res=(lists)omAlloc0Bin(slists_bin);
1031
1032  hexist = hInit(S, Q, &hNexist, currRing);
1033  if (hNexist == 0)
1034  {
1035    intvec *iv=new intvec(rVar(currRing));
1036    for(i=0; i<rVar(currRing); i++) (*iv)[i]=1;
1037    res->Init(1);
1038    res->m[0].rtyp=INTVEC_CMD;
1039    res->m[0].data=(intvec*)iv;
1040    return res;
1041  }
1042  else if (hisModule!=0)
1043  {
1044    res->Init(0);
1045    return res;
1046  }
1047  save = ISet = (indset)omAlloc0Bin(indlist_bin);
1048  hMu = 0;
1049  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1050  hvar = (varset)omAlloc((rVar(currRing) + 1) * sizeof(int));
1051  hpure = (scmon)omAlloc((1 + (rVar(currRing) * rVar(currRing))) * sizeof(long));
1052  hrad = hexist;
1053  hNrad = hNexist;
1054  radmem = hCreate(rVar(currRing) - 1);
1055  hCo = rVar(currRing) + 1;
1056  hNvar = rVar(currRing);
1057  hRadical(hrad, &hNrad, hNvar);
1058  hSupp(hrad, hNrad, hvar, &hNvar);
1059  if (hNvar)
1060  {
1061    hCo = hNvar;
1062    memset(hpure, 0, (rVar(currRing) + 1) * sizeof(long));
1063    hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
1064    hLexR(hrad, hNrad, hvar, hNvar);
1065    hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
1066  }
1067  if (hCo && (hCo < rVar(currRing)))
1068  {
1069    hIndMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
1070  }
1071  if (hMu!=0)
1072  {
1073    ISet = save;
1074    hMu2 = 0;
1075    if (all && (hCo+1 < rVar(currRing)))
1076    {
1077      JSet = (indset)omAlloc0Bin(indlist_bin);
1078      hIndAllMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
1079      i=hMu+hMu2;
1080      res->Init(i);
1081      if (hMu2 == 0)
1082      {
1083        omFreeBin((ADDRESS)JSet, indlist_bin);
1084      }
1085    }
1086    else
1087    {
1088      res->Init(hMu);
1089    }
1090    for (i=0;i<hMu;i++)
1091    {
1092      res->m[i].data = (void *)save->set;
1093      res->m[i].rtyp = INTVEC_CMD;
1094      ISet = save;
1095      save = save->nx;
1096      omFreeBin((ADDRESS)ISet, indlist_bin);
1097    }
1098    omFreeBin((ADDRESS)save, indlist_bin);
1099    if (hMu2 != 0)
1100    {
1101      save = JSet;
1102      for (i=hMu;i<hMu+hMu2;i++)
1103      {
1104        res->m[i].data = (void *)save->set;
1105        res->m[i].rtyp = INTVEC_CMD;
1106        JSet = save;
1107        save = save->nx;
1108        omFreeBin((ADDRESS)JSet, indlist_bin);
1109      }
1110      omFreeBin((ADDRESS)save, indlist_bin);
1111    }
1112  }
1113  else
1114  {
1115    res->Init(0);
1116    omFreeBin((ADDRESS)ISet,  indlist_bin);
1117  }
1118  hKill(radmem, rVar(currRing) - 1);
1119  omFreeSize((ADDRESS)hpure, (1 + (rVar(currRing) * rVar(currRing))) * sizeof(long));
1120  omFreeSize((ADDRESS)hvar, (rVar(currRing) + 1) * sizeof(int));
1121  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1122  hDelete(hexist, hNexist);
1123  return res;
1124}
1125
1126int iiDeclCommand(leftv sy, leftv name, int lev,int t, idhdl* root,BOOLEAN isring, BOOLEAN init_b)
1127{
1128  BOOLEAN res=FALSE;
1129  const char *id = name->name;
1130
1131  memset(sy,0,sizeof(sleftv));
1132  if ((name->name==NULL)||(isdigit(name->name[0])))
1133  {
1134    WerrorS("object to declare is not a name");
1135    res=TRUE;
1136  }
1137  else
1138  {
1139    //if (name->rtyp!=0)
1140    //{
1141    //  Warn("`%s` is already in use",name->name);
1142    //}
1143    {
1144      sy->data = (char *)enterid(id,lev,t,root,init_b);
1145    }
1146    if (sy->data!=NULL)
1147    {
1148      sy->rtyp=IDHDL;
1149      currid=sy->name=IDID((idhdl)sy->data);
1150      // name->name=NULL; /* used in enterid */
1151      //sy->e = NULL;
1152      if (name->next!=NULL)
1153      {
1154        sy->next=(leftv)omAllocBin(sleftv_bin);
1155        res=iiDeclCommand(sy->next,name->next,lev,t,root, isring);
1156      }
1157    }
1158    else res=TRUE;
1159  }
1160  name->CleanUp();
1161  return res;
1162}
1163
1164BOOLEAN iiDefaultParameter(leftv p)
1165{
1166  attr at=NULL;
1167  if (iiCurrProc!=NULL)
1168     at=iiCurrProc->attribute->get("default_arg");
1169  if (at==NULL)
1170    return FALSE;
1171  sleftv tmp;
1172  memset(&tmp,0,sizeof(sleftv));
1173  tmp.rtyp=at->atyp;
1174  tmp.data=at->CopyA();
1175  return iiAssign(p,&tmp);
1176}
1177BOOLEAN iiBranchTo(leftv r, leftv args)
1178{
1179  // <string1...stringN>,<proc>
1180  // known: args!=NULL, l>=1
1181  int l=args->listLength();
1182  int ll=0;
1183  if (iiCurrArgs!=NULL) ll=iiCurrArgs->listLength();
1184  if (ll!=(l-1)) return FALSE;
1185  leftv h=args;
1186  short *t=(short*)omAlloc(l*sizeof(short));
1187  t[0]=l-1;
1188  int b;
1189  int i;
1190  for(i=1;i<l;i++,h=h->next)
1191  {
1192    if (h->Typ()!=STRING_CMD)
1193    {
1194      omFree(t);
1195      Werror("arg %d is not a string",i);
1196      return TRUE;
1197    }
1198    int tt;
1199    b=IsCmd((char *)h->Data(),tt);
1200    if(b) t[i]=tt;
1201    else
1202    {
1203      omFree(t);
1204      Werror("arg %d is not a type name",i);
1205      return TRUE;
1206    }
1207  }
1208  if (h->Typ()!=PROC_CMD)
1209  {
1210    omFree(t);
1211    Werror("last arg is not a proc",i);
1212    return TRUE;
1213  }
1214  b=iiCheckTypes(iiCurrArgs,t,0);
1215  omFree(t);
1216  if (b && (h->rtyp==IDHDL) && (h->e==NULL))
1217  {
1218    BOOLEAN err;
1219    //Print("branchTo: %s\n",h->Name());
1220    iiCurrProc=(idhdl)h->data;
1221    procinfo * pi=IDPROC(iiCurrProc);
1222    if( pi->data.s.body==NULL )
1223    {
1224      iiGetLibProcBuffer(pi);
1225      if (pi->data.s.body==NULL) return TRUE;
1226    }
1227    if ((pi->pack!=NULL)&&(currPack!=pi->pack))
1228    {
1229      currPack=pi->pack;
1230      iiCheckPack(currPack);
1231      currPackHdl=packFindHdl(currPack);
1232      //Print("set pack=%s\n",IDID(currPackHdl));
1233    }
1234    err=iiAllStart(pi,pi->data.s.body,BT_proc,pi->data.s.body_lineno-(iiCurrArgs==NULL));
1235    exitBuffer(BT_proc);
1236    if (iiCurrArgs!=NULL)
1237    {
1238      if (!err) Warn("too many arguments for %s",IDID(iiCurrProc));
1239      iiCurrArgs->CleanUp();
1240      omFreeBin((ADDRESS)iiCurrArgs, sleftv_bin);
1241      iiCurrArgs=NULL;
1242    }
1243    return 2-err;
1244  }
1245  return FALSE;
1246}
1247BOOLEAN iiParameter(leftv p)
1248{
1249  if (iiCurrArgs==NULL)
1250  {
1251    if (strcmp(p->name,"#")==0)
1252      return iiDefaultParameter(p);
1253    Werror("not enough arguments for proc %s",VoiceName());
1254    p->CleanUp();
1255    return TRUE;
1256  }
1257  leftv h=iiCurrArgs;
1258  leftv rest=h->next; /*iiCurrArgs is not NULL here*/
1259  BOOLEAN is_default_list=FALSE;
1260  if (strcmp(p->name,"#")==0)
1261  {
1262    is_default_list=TRUE;
1263    rest=NULL;
1264  }
1265  else
1266  {
1267    h->next=NULL;
1268  }
1269  BOOLEAN res=iiAssign(p,h);
1270  if (is_default_list)
1271  {
1272    iiCurrArgs=NULL;
1273  }
1274  else
1275  {
1276    iiCurrArgs=rest;
1277  }
1278  h->CleanUp();
1279  omFreeBin((ADDRESS)h, sleftv_bin);
1280  return res;
1281}
1282BOOLEAN iiAlias(leftv p)
1283{
1284  if (iiCurrArgs==NULL)
1285  {
1286    Werror("not enough arguments for proc %s",VoiceName());
1287    p->CleanUp();
1288    return TRUE;
1289  }
1290  leftv h=iiCurrArgs;
1291  iiCurrArgs=h->next;
1292  h->next=NULL;
1293  if (h->rtyp!=IDHDL)
1294  {
1295    BOOLEAN res=iiAssign(p,h);
1296    h->CleanUp();
1297    omFreeBin((ADDRESS)h, sleftv_bin);
1298    return res;
1299  }
1300  if (h->Typ()!=p->Typ())
1301  {
1302    WerrorS("type mismatch");
1303    return TRUE;
1304  }
1305  idhdl pp=(idhdl)p->data;
1306  switch(pp->typ)
1307  {
1308#ifdef SINGULAR_4_1
1309      case CRING_CMD:
1310        nKillChar((coeffs)pp);
1311        break;
1312#endif
1313      case INT_CMD:
1314        break;
1315      case INTVEC_CMD:
1316      case INTMAT_CMD:
1317         delete IDINTVEC(pp);
1318         break;
1319      case NUMBER_CMD:
1320         nDelete(&IDNUMBER(pp));
1321         break;
1322      case BIGINT_CMD:
1323         n_Delete(&IDNUMBER(pp),coeffs_BIGINT);
1324         break;
1325      case MAP_CMD:
1326         {
1327           map im = IDMAP(pp);
1328           omFree((ADDRESS)im->preimage);
1329         }
1330         // continue as ideal:
1331      case IDEAL_CMD:
1332      case MODUL_CMD:
1333      case MATRIX_CMD:
1334          idDelete(&IDIDEAL(pp));
1335         break;
1336      case PROC_CMD:
1337      case RESOLUTION_CMD:
1338      case STRING_CMD:
1339         omFree((ADDRESS)IDSTRING(pp));
1340         break;
1341      case LIST_CMD:
1342         IDLIST(pp)->Clean();
1343         break;
1344      case LINK_CMD:
1345         omFreeBin(IDLINK(pp),sip_link_bin);
1346         break;
1347       // case ring: cannot happen
1348       default:
1349         Werror("unknown type %d",p->Typ());
1350         return TRUE;
1351  }
1352  pp->typ=ALIAS_CMD;
1353  IDDATA(pp)=(char*)h->data;
1354  h->CleanUp();
1355  omFreeBin((ADDRESS)h, sleftv_bin);
1356  return FALSE;
1357}
1358
1359static BOOLEAN iiInternalExport (leftv v, int toLev)
1360{
1361  idhdl h=(idhdl)v->data;
1362  //Print("iiInternalExport('%s',%d)%s\n", v->name, toLev,"");
1363  if (IDLEV(h)==0)
1364  {
1365    if (BVERBOSE(V_REDEFINE)) Warn("`%s` is already global",IDID(h));
1366  }
1367  else
1368  {
1369    h=IDROOT->get(v->name,toLev);
1370    idhdl *root=&IDROOT;
1371    if ((h==NULL)&&(currRing!=NULL))
1372    {
1373      h=currRing->idroot->get(v->name,toLev);
1374      root=&currRing->idroot;
1375    }
1376    BOOLEAN keepring=FALSE;
1377    if ((h!=NULL)&&(IDLEV(h)==toLev))
1378    {
1379      if (IDTYP(h)==v->Typ())
1380      {
1381        if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
1382        && (v->Data()==IDDATA(h)))
1383        {
1384          IDRING(h)->ref++;
1385          keepring=TRUE;
1386          IDLEV(h)=toLev;
1387          //WarnS("keepring");
1388          return FALSE;
1389        }
1390        if (BVERBOSE(V_REDEFINE))
1391        {
1392          Warn("redefining %s",IDID(h));
1393        }
1394#ifdef USE_IILOCALRING
1395        if (iiLocalRing[0]==IDRING(h) && (!keepring)) iiLocalRing[0]=NULL;
1396#else
1397        proclevel *p=procstack;
1398        while (p->next!=NULL) p=p->next;
1399        if ((p->cRing==IDRING(h)) && (!keepring))
1400        {
1401          p->cRing=NULL;
1402          p->cRingHdl=NULL;
1403        }
1404#endif
1405        killhdl2(h,root,currRing);
1406      }
1407      else
1408      {
1409        return TRUE;
1410      }
1411    }
1412    h=(idhdl)v->data;
1413    IDLEV(h)=toLev;
1414    if (keepring) IDRING(h)->ref--;
1415    iiNoKeepRing=FALSE;
1416    //Print("export %s\n",IDID(h));
1417  }
1418  return FALSE;
1419}
1420
1421BOOLEAN iiInternalExport (leftv v, int toLev, package rootpack)
1422{
1423  idhdl h=(idhdl)v->data;
1424  if(h==NULL)
1425  {
1426    Warn("'%s': no such identifier\n", v->name);
1427    return FALSE;
1428  }
1429  package frompack=v->req_packhdl;
1430  if (frompack==NULL) frompack=currPack;
1431  if ((RingDependend(IDTYP(h)))
1432  || ((IDTYP(h)==LIST_CMD)
1433     && (lRingDependend(IDLIST(h)))
1434     )
1435  )
1436  {
1437    //Print("// ==> Ringdependent set nesting to 0\n");
1438    return (iiInternalExport(v, toLev));
1439  }
1440  else
1441  {
1442    IDLEV(h)=toLev;
1443    v->req_packhdl=rootpack;
1444    if (h==frompack->idroot)
1445    {
1446      frompack->idroot=h->next;
1447    }
1448    else
1449    {
1450      idhdl hh=frompack->idroot;
1451      while ((hh!=NULL) && (hh->next!=h))
1452        hh=hh->next;
1453      if ((hh!=NULL) && (hh->next==h))
1454        hh->next=h->next;
1455      else
1456      {
1457        Werror("`%s` not found",v->Name());
1458        return TRUE;
1459      }
1460    }
1461    h->next=rootpack->idroot;
1462    rootpack->idroot=h;
1463  }
1464  return FALSE;
1465}
1466
1467BOOLEAN iiExport (leftv v, int toLev)
1468{
1469  BOOLEAN nok=FALSE;
1470  leftv r=v;
1471  while (v!=NULL)
1472  {
1473    if ((v->name==NULL)||(v->rtyp==0)||(v->e!=NULL))
1474    {
1475      Werror("cannot export:%s of internal type %d",v->name,v->rtyp);
1476      nok=TRUE;
1477    }
1478    else
1479    {
1480      if(iiInternalExport(v, toLev))
1481      {
1482        r->CleanUp();
1483        return TRUE;
1484      }
1485    }
1486    v=v->next;
1487  }
1488  r->CleanUp();
1489  return nok;
1490}
1491
1492/*assume root!=idroot*/
1493BOOLEAN iiExport (leftv v, int toLev, package pack)
1494{
1495#ifdef SINGULAR_4_1
1496  if ((pack==basePack)&&(pack!=currPack))
1497  { Warn("'exportto' to Top is depreciated in >>%s<<",my_yylinebuf);}
1498#endif
1499  BOOLEAN nok=FALSE;
1500  leftv rv=v;
1501  while (v!=NULL)
1502  {
1503    if ((v->name==NULL)||(v->rtyp==0)||(v->e!=NULL)
1504    )
1505    {
1506      Werror("cannot export:%s of internal type %d",v->name,v->rtyp);
1507      nok=TRUE;
1508    }
1509    else
1510    {
1511      idhdl old=pack->idroot->get( v->name,toLev);
1512      if (old!=NULL)
1513      {
1514        if ((pack==currPack) && (old==(idhdl)v->data))
1515        {
1516          if (BVERBOSE(V_REDEFINE)) Warn("`%s` is already global",IDID(old));
1517          break;
1518        }
1519        else if (IDTYP(old)==v->Typ())
1520        {
1521          if (BVERBOSE(V_REDEFINE))
1522          {
1523            Warn("redefining %s",IDID(old));
1524          }
1525          v->name=omStrDup(v->name);
1526          killhdl2(old,&(pack->idroot),currRing);
1527        }
1528        else
1529        {
1530          rv->CleanUp();
1531          return TRUE;
1532        }
1533      }
1534      //Print("iiExport: pack=%s\n",IDID(root));
1535      if(iiInternalExport(v, toLev, pack))
1536      {
1537        rv->CleanUp();
1538        return TRUE;
1539      }
1540    }
1541    v=v->next;
1542  }
1543  rv->CleanUp();
1544  return nok;
1545}
1546
1547BOOLEAN iiCheckRing(int i)
1548{
1549  if (currRing==NULL)
1550  {
1551    #ifdef SIQ
1552    if (siq<=0)
1553    {
1554    #endif
1555      if (RingDependend(i))
1556      {
1557        WerrorS("no ring active");
1558        return TRUE;
1559      }
1560    #ifdef SIQ
1561    }
1562    #endif
1563  }
1564  return FALSE;
1565}
1566
1567poly    iiHighCorner(ideal I, int ak)
1568{
1569  int i;
1570  if(!idIsZeroDim(I)) return NULL; // not zero-dim.
1571  poly po=NULL;
1572  if (rHasLocalOrMixedOrdering_currRing())
1573  {
1574    scComputeHC(I,currRing->qideal,ak,po);
1575    if (po!=NULL)
1576    {
1577      pGetCoeff(po)=nInit(1);
1578      for (i=rVar(currRing); i>0; i--)
1579      {
1580        if (pGetExp(po, i) > 0) pDecrExp(po,i);
1581      }
1582      pSetComp(po,ak);
1583      pSetm(po);
1584    }
1585  }
1586  else
1587    po=pOne();
1588  return po;
1589}
1590
1591void iiCheckPack(package &p)
1592{
1593  if (p==basePack) return;
1594
1595  idhdl t=basePack->idroot;
1596
1597  while ((t!=NULL) && (IDTYP(t)!=PACKAGE_CMD) && (IDPACKAGE(t)!=p)) t=t->next;
1598
1599  if (t==NULL)
1600  {
1601    WarnS("package not found\n");
1602    p=basePack;
1603  }
1604  return;
1605}
1606
1607idhdl rDefault(const char *s)
1608{
1609  idhdl tmp=NULL;
1610
1611  if (s!=NULL) tmp = enterid(s, myynest, RING_CMD, &IDROOT);
1612  if (tmp==NULL) return NULL;
1613
1614// if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
1615  if (sLastPrinted.RingDependend())
1616  {
1617    sLastPrinted.CleanUp();
1618    memset(&sLastPrinted,0,sizeof(sleftv));
1619  }
1620
1621  ring r = IDRING(tmp);
1622
1623  r->cf = nInitChar(n_Zp, (void*)32003); //   r->cf->ch = 32003;
1624  r->N      = 3;
1625  /*r->P     = 0; Alloc0 in idhdl::set, ipid.cc*/
1626  /*names*/
1627  r->names = (char **) omAlloc0(3 * sizeof(char_ptr));
1628  r->names[0]  = omStrDup("x");
1629  r->names[1]  = omStrDup("y");
1630  r->names[2]  = omStrDup("z");
1631  /*weights: entries for 3 blocks: NULL*/
1632  r->wvhdl = (int **)omAlloc0(3 * sizeof(int_ptr));
1633  /*order: dp,C,0*/
1634  r->order = (int *) omAlloc(3 * sizeof(int *));
1635  r->block0 = (int *)omAlloc0(3 * sizeof(int *));
1636  r->block1 = (int *)omAlloc0(3 * sizeof(int *));
1637  /* ringorder dp for the first block: var 1..3 */
1638  r->order[0]  = ringorder_dp;
1639  r->block0[0] = 1;
1640  r->block1[0] = 3;
1641  /* ringorder C for the second block: no vars */
1642  r->order[1]  = ringorder_C;
1643  /* the last block: everything is 0 */
1644  r->order[2]  = 0;
1645
1646  /* complete ring intializations */
1647  rComplete(r);
1648  rSetHdl(tmp);
1649  return currRingHdl;
1650}
1651
1652idhdl rFindHdl(ring r, idhdl n)
1653{
1654  idhdl h=rSimpleFindHdl(r,IDROOT,n);
1655  if (h!=NULL)  return h;
1656  if (IDROOT!=basePack->idroot) h=rSimpleFindHdl(r,basePack->idroot,n);
1657  if (h!=NULL)  return h;
1658  proclevel *p=procstack;
1659  while(p!=NULL)
1660  {
1661    if ((p->cPack!=basePack)
1662    && (p->cPack!=currPack))
1663      h=rSimpleFindHdl(r,p->cPack->idroot,n);
1664    if (h!=NULL)  return h;
1665    p=p->next;
1666  }
1667  idhdl tmp=basePack->idroot;
1668  while (tmp!=NULL)
1669  {
1670    if (IDTYP(tmp)==PACKAGE_CMD)
1671      h=rSimpleFindHdl(r,IDPACKAGE(tmp)->idroot,n);
1672    if (h!=NULL)  return h;
1673    tmp=IDNEXT(tmp);
1674  }
1675  return NULL;
1676}
1677
1678void rDecomposeCF(leftv h,const ring r,const ring R)
1679{
1680  lists L=(lists)omAlloc0Bin(slists_bin);
1681  L->Init(4);
1682  h->rtyp=LIST_CMD;
1683  h->data=(void *)L;
1684  // 0: char/ cf - ring
1685  // 1: list (var)
1686  // 2: list (ord)
1687  // 3: qideal
1688  // ----------------------------------------
1689  // 0: char/ cf - ring
1690  L->m[0].rtyp=INT_CMD;
1691  L->m[0].data=(void *)(long)r->cf->ch;
1692  // ----------------------------------------
1693  // 1: list (var)
1694  lists LL=(lists)omAlloc0Bin(slists_bin);
1695  LL->Init(r->N);
1696  int i;
1697  for(i=0; i<r->N; i++)
1698  {
1699    LL->m[i].rtyp=STRING_CMD;
1700    LL->m[i].data=(void *)omStrDup(r->names[i]);
1701  }
1702  L->m[1].rtyp=LIST_CMD;
1703  L->m[1].data=(void *)LL;
1704  // ----------------------------------------
1705  // 2: list (ord)
1706  LL=(lists)omAlloc0Bin(slists_bin);
1707  i=rBlocks(r)-1;
1708  LL->Init(i);
1709  i--;
1710  lists LLL;
1711  for(; i>=0; i--)
1712  {
1713    intvec *iv;
1714    int j;
1715    LL->m[i].rtyp=LIST_CMD;
1716    LLL=(lists)omAlloc0Bin(slists_bin);
1717    LLL->Init(2);
1718    LLL->m[0].rtyp=STRING_CMD;
1719    LLL->m[0].data=(void *)omStrDup(rSimpleOrdStr(r->order[i]));
1720    if (r->block1[i]-r->block0[i] >=0 )
1721    {
1722      j=r->block1[i]-r->block0[i];
1723      if(r->order[i]==ringorder_M) j=(j+1)*(j+1)-1;
1724      iv=new intvec(j+1);
1725      if ((r->wvhdl!=NULL) && (r->wvhdl[i]!=NULL))
1726      {
1727        for(;j>=0; j--) (*iv)[j]=r->wvhdl[i][j];
1728      }
1729      else switch (r->order[i])
1730      {
1731        case ringorder_dp:
1732        case ringorder_Dp:
1733        case ringorder_ds:
1734        case ringorder_Ds:
1735        case ringorder_lp:
1736          for(;j>=0; j--) (*iv)[j]=1;
1737          break;
1738        default: /* do nothing */;
1739      }
1740    }
1741    else
1742    {
1743      iv=new intvec(1);
1744    }
1745    LLL->m[1].rtyp=INTVEC_CMD;
1746    LLL->m[1].data=(void *)iv;
1747    LL->m[i].data=(void *)LLL;
1748  }
1749  L->m[2].rtyp=LIST_CMD;
1750  L->m[2].data=(void *)LL;
1751  // ----------------------------------------
1752  // 3: qideal
1753  L->m[3].rtyp=IDEAL_CMD;
1754  if (nCoeff_is_transExt(R->cf))
1755    L->m[3].data=(void *)idInit(1,1);
1756  else
1757  {
1758    ideal q=idInit(IDELEMS(r->qideal));
1759    q->m[0]=p_Init(R);
1760    pSetCoeff0(q->m[0],(number)(r->qideal->m[0]));
1761    L->m[3].data=(void *)q;
1762//    I->m[0] = pNSet(R->minpoly);
1763  }
1764  // ----------------------------------------
1765}
1766void rDecomposeC(leftv h,const ring R)
1767/* field is R or C */
1768{
1769  lists L=(lists)omAlloc0Bin(slists_bin);
1770  if (rField_is_long_C(R)) L->Init(3);
1771  else                     L->Init(2);
1772  h->rtyp=LIST_CMD;
1773  h->data=(void *)L;
1774  // 0: char/ cf - ring
1775  // 1: list (var)
1776  // 2: list (ord)
1777  // ----------------------------------------
1778  // 0: char/ cf - ring
1779  L->m[0].rtyp=INT_CMD;
1780  L->m[0].data=(void *)0;
1781  // ----------------------------------------
1782  // 1:
1783  lists LL=(lists)omAlloc0Bin(slists_bin);
1784  LL->Init(2);
1785    LL->m[0].rtyp=INT_CMD;
1786    LL->m[0].data=(void *)(long)si_max(R->cf->float_len,SHORT_REAL_LENGTH/2);
1787    LL->m[1].rtyp=INT_CMD;
1788    LL->m[1].data=(void *)(long)si_max(R->cf->float_len2,SHORT_REAL_LENGTH);
1789  L->m[1].rtyp=LIST_CMD;
1790  L->m[1].data=(void *)LL;
1791  // ----------------------------------------
1792  // 2: list (par)
1793  if (rField_is_long_C(R))
1794  {
1795    L->m[2].rtyp=STRING_CMD;
1796    L->m[2].data=(void *)omStrDup(*rParameter(R));
1797  }
1798  // ----------------------------------------
1799}
1800
1801#ifdef HAVE_RINGS
1802void rDecomposeRing(leftv h,const ring R)
1803/* field is R or C */
1804{
1805  lists L=(lists)omAlloc0Bin(slists_bin);
1806  if (rField_is_Ring_Z(R)) L->Init(1);
1807  else                     L->Init(2);
1808  h->rtyp=LIST_CMD;
1809  h->data=(void *)L;
1810  // 0: char/ cf - ring
1811  // 1: list (module)
1812  // ----------------------------------------
1813  // 0: char/ cf - ring
1814  L->m[0].rtyp=STRING_CMD;
1815  L->m[0].data=(void *)omStrDup("integer");
1816  // ----------------------------------------
1817  // 1: module
1818  if (rField_is_Ring_Z(R)) return;
1819  lists LL=(lists)omAlloc0Bin(slists_bin);
1820  LL->Init(2);
1821  LL->m[0].rtyp=BIGINT_CMD;
1822  LL->m[0].data=nlMapGMP((number) R->cf->modBase, R->cf, R->cf); // TODO: what is this?? // extern number nlMapGMP(number from, const coeffs src, const coeffs dst); // FIXME: replace with n_InitMPZ(R->cf->modBase, coeffs_BIGINT); ?
1823  LL->m[1].rtyp=INT_CMD;
1824  LL->m[1].data=(void *) R->cf->modExponent;
1825  L->m[1].rtyp=LIST_CMD;
1826  L->m[1].data=(void *)LL;
1827}
1828#endif
1829
1830
1831lists rDecompose(const ring r)
1832{
1833  assume( r != NULL );
1834  const coeffs C = r->cf;
1835  assume( C != NULL );
1836
1837  // sanity check: require currRing==r for rings with polynomial data
1838  if ( (r!=currRing) && (
1839           (nCoeff_is_algExt(C) && (C != currRing->cf))
1840        || (r->qideal != NULL)
1841#ifdef HAVE_PLURAL
1842        || (rIsPluralRing(r))
1843#endif
1844                        )
1845     )
1846  {
1847    WerrorS("ring with polynomial data must be the base ring or compatible");
1848    return NULL;
1849  }
1850  // 0: char/ cf - ring
1851  // 1: list (var)
1852  // 2: list (ord)
1853  // 3: qideal
1854  // possibly:
1855  // 4: C
1856  // 5: D
1857  lists L=(lists)omAlloc0Bin(slists_bin);
1858  if (rIsPluralRing(r))
1859    L->Init(6);
1860  else
1861    L->Init(4);
1862  // ----------------------------------------
1863  // 0: char/ cf - ring
1864#ifdef SINGULAR_4_1
1865  // 0: char/ cf - ring
1866  L->m[0].rtyp=CRING_CMD;
1867  L->m[0].data=(char*)r->cf; r->cf->ref++;
1868#else
1869  if (rField_is_numeric(r))
1870  {
1871    rDecomposeC(&(L->m[0]),r);
1872  }
1873#ifdef HAVE_RINGS
1874  else if (rField_is_Ring(r))
1875  {
1876    rDecomposeRing(&(L->m[0]),r);
1877  }
1878#endif
1879  else if ( r->cf->extRing!=NULL )// nCoeff_is_algExt(r->cf))
1880  {
1881    rDecomposeCF(&(L->m[0]), r->cf->extRing, r);
1882  }
1883  else if(rField_is_GF(r))
1884  {
1885    lists Lc=(lists)omAlloc0Bin(slists_bin);
1886    Lc->Init(4);
1887    // char:
1888    Lc->m[0].rtyp=INT_CMD;
1889    Lc->m[0].data=(void*)(long)r->cf->m_nfCharQ;
1890    // var:
1891    lists Lv=(lists)omAlloc0Bin(slists_bin);
1892    Lv->Init(1);
1893    Lv->m[0].rtyp=STRING_CMD;
1894    Lv->m[0].data=(void *)omStrDup(*rParameter(r));
1895    Lc->m[1].rtyp=LIST_CMD;
1896    Lc->m[1].data=(void*)Lv;
1897    // ord:
1898    lists Lo=(lists)omAlloc0Bin(slists_bin);
1899    Lo->Init(1);
1900    lists Loo=(lists)omAlloc0Bin(slists_bin);
1901    Loo->Init(2);
1902    Loo->m[0].rtyp=STRING_CMD;
1903    Loo->m[0].data=(void *)omStrDup(rSimpleOrdStr(ringorder_lp));
1904
1905    intvec *iv=new intvec(1); (*iv)[0]=1;
1906    Loo->m[1].rtyp=INTVEC_CMD;
1907    Loo->m[1].data=(void *)iv;
1908
1909    Lo->m[0].rtyp=LIST_CMD;
1910    Lo->m[0].data=(void*)Loo;
1911
1912    Lc->m[2].rtyp=LIST_CMD;
1913    Lc->m[2].data=(void*)Lo;
1914    // q-ideal:
1915    Lc->m[3].rtyp=IDEAL_CMD;
1916    Lc->m[3].data=(void *)idInit(1,1);
1917    // ----------------------
1918    L->m[0].rtyp=LIST_CMD;
1919    L->m[0].data=(void*)Lc;
1920  }
1921  else
1922  {
1923    L->m[0].rtyp=INT_CMD;
1924    L->m[0].data=(void *)(long)r->cf->ch;
1925  }
1926#endif
1927  // ----------------------------------------
1928  // 1: list (var)
1929  lists LL=(lists)omAlloc0Bin(slists_bin);
1930  LL->Init(r->N);
1931  int i;
1932  for(i=0; i<r->N; i++)
1933  {
1934    LL->m[i].rtyp=STRING_CMD;
1935    LL->m[i].data=(void *)omStrDup(r->names[i]);
1936  }
1937  L->m[1].rtyp=LIST_CMD;
1938  L->m[1].data=(void *)LL;
1939  // ----------------------------------------
1940  // 2: list (ord)
1941  LL=(lists)omAlloc0Bin(slists_bin);
1942  i=rBlocks(r)-1;
1943  LL->Init(i);
1944  i--;
1945  lists LLL;
1946  for(; i>=0; i--)
1947  {
1948    intvec *iv;
1949    int j;
1950    LL->m[i].rtyp=LIST_CMD;
1951    LLL=(lists)omAlloc0Bin(slists_bin);
1952    LLL->Init(2);
1953    LLL->m[0].rtyp=STRING_CMD;
1954    LLL->m[0].data=(void *)omStrDup(rSimpleOrdStr(r->order[i]));
1955
1956    if(r->order[i] == ringorder_IS) //  || r->order[i] == ringorder_s || r->order[i] == ringorder_S)
1957    {
1958      assume( r->block0[i] == r->block1[i] );
1959      const int s = r->block0[i];
1960      assume( -2 < s && s < 2);
1961
1962      iv=new intvec(1);
1963      (*iv)[0] = s;
1964    }
1965    else if (r->block1[i]-r->block0[i] >=0 )
1966    {
1967      int bl=j=r->block1[i]-r->block0[i];
1968      if (r->order[i]==ringorder_M)
1969      {
1970        j=(j+1)*(j+1)-1;
1971        bl=j+1;
1972      }
1973      else if (r->order[i]==ringorder_am)
1974      {
1975        j+=r->wvhdl[i][bl+1];
1976      }
1977      iv=new intvec(j+1);
1978      if ((r->wvhdl!=NULL) && (r->wvhdl[i]!=NULL))
1979      {
1980        for(;j>=0; j--) (*iv)[j]=r->wvhdl[i][j+(j>bl)];
1981      }
1982      else switch (r->order[i])
1983      {
1984        case ringorder_dp:
1985        case ringorder_Dp:
1986        case ringorder_ds:
1987        case ringorder_Ds:
1988        case ringorder_lp:
1989          for(;j>=0; j--) (*iv)[j]=1;
1990          break;
1991        default: /* do nothing */;
1992      }
1993    }
1994    else
1995    {
1996      iv=new intvec(1);
1997    }
1998    LLL->m[1].rtyp=INTVEC_CMD;
1999    LLL->m[1].data=(void *)iv;
2000    LL->m[i].data=(void *)LLL;
2001  }
2002  L->m[2].rtyp=LIST_CMD;
2003  L->m[2].data=(void *)LL;
2004  // ----------------------------------------
2005  // 3: qideal
2006  L->m[3].rtyp=IDEAL_CMD;
2007  if (r->qideal==NULL)
2008    L->m[3].data=(void *)idInit(1,1);
2009  else
2010    L->m[3].data=(void *)idCopy(r->qideal);
2011  // ----------------------------------------
2012#ifdef HAVE_PLURAL // NC! in rDecompose
2013  if (rIsPluralRing(r))
2014  {
2015    L->m[4].rtyp=MATRIX_CMD;
2016    L->m[4].data=(void *)mp_Copy(r->GetNC()->C, r, r);
2017    L->m[5].rtyp=MATRIX_CMD;
2018    L->m[5].data=(void *)mp_Copy(r->GetNC()->D, r, r);
2019  }
2020#endif
2021  return L;
2022}
2023
2024void rComposeC(lists L, ring R)
2025/* field is R or C */
2026{
2027  // ----------------------------------------
2028  // 0: char/ cf - ring
2029  if ((L->m[0].rtyp!=INT_CMD) || (L->m[0].data!=(char *)0))
2030  {
2031    Werror("invald coeff. field description, expecting 0");
2032    return;
2033  }
2034//  R->cf->ch=0;
2035  // ----------------------------------------
2036  // 1:
2037  if (L->m[1].rtyp!=LIST_CMD)
2038    Werror("invald coeff. field description, expecting precision list");
2039  lists LL=(lists)L->m[1].data;
2040  int r1=(int)(long)LL->m[0].data;
2041  int r2=(int)(long)LL->m[1].data;
2042  if (L->nr==2) // complex
2043    R->cf = nInitChar(n_long_C, NULL);
2044  else if ((r1<=SHORT_REAL_LENGTH)
2045  && (r2=SHORT_REAL_LENGTH))
2046    R->cf = nInitChar(n_R, NULL);
2047  else
2048  {
2049    LongComplexInfo* p = (LongComplexInfo *)omAlloc0(sizeof(LongComplexInfo));
2050    p->float_len=r1;
2051    p->float_len2=r2;
2052    R->cf = nInitChar(n_long_R, NULL);
2053  }
2054
2055  if ((r1<=SHORT_REAL_LENGTH)   // should go into nInitChar
2056  && (r2=SHORT_REAL_LENGTH))
2057  {
2058    R->cf->float_len=SHORT_REAL_LENGTH/2;
2059    R->cf->float_len2=SHORT_REAL_LENGTH;
2060  }
2061  else
2062  {
2063    R->cf->float_len=si_min(r1,32767);
2064    R->cf->float_len2=si_min(r2,32767);
2065  }
2066  // ----------------------------------------
2067  // 2: list (par)
2068  if (L->nr==2)
2069  {
2070    //R->cf->extRing->N=1;
2071    if (L->m[2].rtyp!=STRING_CMD)
2072    {
2073      Werror("invald coeff. field description, expecting parameter name");
2074      return;
2075    }
2076    //(rParameter(R))=(char**)omAlloc0(rPar(R)*sizeof(char_ptr));
2077    rParameter(R)[0]=omStrDup((char *)L->m[2].data);
2078  }
2079  // ----------------------------------------
2080}
2081
2082#ifdef HAVE_RINGS
2083void rComposeRing(lists L, ring R)
2084/* field is R or C */
2085{
2086  // ----------------------------------------
2087  // 0: string: integer
2088  // no further entries --> Z
2089  mpz_ptr modBase = NULL;
2090  unsigned int modExponent = 1;
2091
2092  modBase = (mpz_ptr) omAlloc(sizeof(mpz_t));
2093  if (L->nr == 0)
2094  {
2095    mpz_init_set_ui(modBase,0);
2096    modExponent = 1;
2097  }
2098  // ----------------------------------------
2099  // 1:
2100  else
2101  {
2102    if (L->m[1].rtyp!=LIST_CMD) Werror("invald data, expecting list of numbers");
2103    lists LL=(lists)L->m[1].data;
2104    if ((LL->nr >= 0) && LL->m[0].rtyp == BIGINT_CMD)
2105    {
2106      number tmp= (number) LL->m[0].data; // never use CopyD() on list elements
2107                                    // assume that tmp is integer, not rational
2108      n_MPZ (modBase, tmp, coeffs_BIGINT);
2109    }
2110    else if (LL->nr >= 0 && LL->m[0].rtyp == INT_CMD)
2111    {
2112      mpz_init_set_ui(modBase,(unsigned long) LL->m[0].data);
2113    }
2114    else
2115    {
2116      mpz_init_set_ui(modBase,0);
2117    }
2118    if (LL->nr >= 1)
2119    {
2120      modExponent = (unsigned long) LL->m[1].data;
2121    }
2122    else
2123    {
2124      modExponent = 1;
2125    }
2126  }
2127  // ----------------------------------------
2128  if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_cmp_ui(modBase, 0) < 0))
2129  {
2130    Werror("Wrong ground ring specification (module is 1)");
2131    return;
2132  }
2133  if (modExponent < 1)
2134  {
2135    Werror("Wrong ground ring specification (exponent smaller than 1");
2136    return;
2137  }
2138  // module is 0 ---> integers
2139  if (mpz_cmp_ui(modBase, 0) == 0)
2140  {
2141    R->cf=nInitChar(n_Z,NULL);
2142  }
2143  // we have an exponent
2144  else if (modExponent > 1)
2145  {
2146    //R->cf->ch = R->cf->modExponent;
2147    if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
2148    {
2149      /* this branch should be active for modExponent = 2..32 resp. 2..64,
2150           depending on the size of a long on the respective platform */
2151      R->cf=nInitChar(n_Z2m,(void*)(long)modExponent);       // Use Z/2^ch
2152      omFreeSize (modBase, sizeof(mpz_t));
2153    }
2154    else
2155    {
2156      //ringtype 3
2157      ZnmInfo info;
2158      info.base= modBase;
2159      info.exp= modExponent;
2160      R->cf=nInitChar(n_Znm,(void*) &info);
2161    }
2162  }
2163  // just a module m > 1
2164  else
2165  {
2166    //ringtype = 2;
2167    //const int ch = mpz_get_ui(modBase);
2168    ZnmInfo info;
2169    info.base= modBase;
2170    info.exp= modExponent;
2171    R->cf=nInitChar(n_Zn,(void*) &info);
2172  }
2173}
2174#endif
2175
2176static void rRenameVars(ring R)
2177{
2178  int i,j;
2179  BOOLEAN ch;
2180  do
2181  {
2182    ch=0;
2183    for(i=0;i<R->N-1;i++)
2184    {
2185      for(j=i+1;j<R->N;j++)
2186      {
2187        if (strcmp(R->names[i],R->names[j])==0)
2188        {
2189          ch=TRUE;
2190          Warn("name conflict var(%d) and var(%d): `%s`, rename to `@%s`",i+1,j+1,R->names[i],R->names[i]);
2191          omFree(R->names[j]);
2192          R->names[j]=(char *)omAlloc(2+strlen(R->names[i]));
2193          sprintf(R->names[j],"@%s",R->names[i]);
2194        }
2195      }
2196    }
2197  }
2198  while (ch);
2199  for(i=0;i<rPar(R); i++)
2200  {
2201    for(j=0;j<R->N;j++)
2202    {
2203      if (strcmp(rParameter(R)[i],R->names[j])==0)
2204      {
2205        Warn("name conflict par(%d) and var(%d): `%s`, renaming the VARIABLE to `@@(%d)`",i+1,j+1,R->names[j],i+1);
2206//        omFree(rParameter(R)[i]);
2207//        rParameter(R)[i]=(char *)omAlloc(10);
2208//        sprintf(rParameter(R)[i],"@@(%d)",i+1);
2209        omFree(R->names[j]);
2210        R->names[j]=(char *)omAlloc(10);
2211        sprintf(R->names[j],"@@(%d)",i+1);
2212      }
2213    }
2214  }
2215}
2216
2217ring rCompose(const lists  L, const BOOLEAN check_comp)
2218{
2219  if ((L->nr!=3)
2220#ifdef HAVE_PLURAL
2221  &&(L->nr!=5)
2222#endif
2223  )
2224    return NULL;
2225  int is_gf_char=0;
2226  // 0: char/ cf - ring
2227  // 1: list (var)
2228  // 2: list (ord)
2229  // 3: qideal
2230  // possibly:
2231  // 4: C
2232  // 5: D
2233
2234  ring R = (ring) omAlloc0Bin(sip_sring_bin);
2235
2236
2237  // ------------------------------------------------------------------
2238  // 0: char:
2239#ifdef SINGULAR_4_1
2240  if (L->m[0].Typ()==CRING_CMD)
2241  {
2242    R->cf=(coeffs)L->m[0].Data();
2243    R->cf->ref++;
2244  }
2245  else
2246#endif
2247  if (L->m[0].Typ()==INT_CMD)
2248  {
2249    int ch = (int)(long)L->m[0].Data();
2250    assume( ch >= 0 );
2251
2252    if (ch == 0) // Q?
2253      R->cf = nInitChar(n_Q, NULL);
2254    else
2255    {
2256      int l = IsPrime(ch); // Zp?
2257      if( l != ch )
2258      {
2259        Warn("%d is invalid characteristic of ground field. %d is used.", ch, l);
2260        ch = l;
2261      }
2262      R->cf = nInitChar(n_Zp, (void*)(long)ch);
2263    }
2264  }
2265  else if (L->m[0].Typ()==LIST_CMD) // something complicated...
2266  {
2267    lists LL=(lists)L->m[0].Data();
2268
2269#ifdef HAVE_RINGS
2270    if (LL->m[0].Typ() == STRING_CMD) // 1st comes a string?
2271    {
2272      rComposeRing(LL, R); // Ring!?
2273    }
2274    else
2275#endif
2276    if (LL->nr < 3)
2277      rComposeC(LL,R); // R, long_R, long_C
2278    else
2279    {
2280      if (LL->m[0].Typ()==INT_CMD)
2281      {
2282        int ch = (int)(long)LL->m[0].Data();
2283        while ((ch!=fftable[is_gf_char]) && (fftable[is_gf_char])) is_gf_char++;
2284        if (fftable[is_gf_char]==0) is_gf_char=-1;
2285
2286        if(is_gf_char!= -1)
2287        {
2288          GFInfo param;
2289
2290          param.GFChar = ch;
2291          param.GFDegree = 1;
2292          param.GFPar_name = (const char*)(((lists)(LL->m[1].Data()))->m[0].Data());
2293
2294          // nfInitChar should be able to handle the case when ch is in fftables!
2295          R->cf = nInitChar(n_GF, (void*)&param);
2296        }
2297      }
2298
2299      if( R->cf == NULL )
2300      {
2301        ring extRing = rCompose((lists)L->m[0].Data(),FALSE);
2302
2303        if (extRing==NULL)
2304        {
2305          WerrorS("could not create the specified coefficient field");
2306          goto rCompose_err;
2307        }
2308
2309        if( extRing->qideal != NULL ) // Algebraic extension
2310        {
2311          AlgExtInfo extParam;
2312
2313          extParam.r = extRing;
2314
2315          R->cf = nInitChar(n_algExt, (void*)&extParam);
2316        }
2317        else // Transcendental extension
2318        {
2319          TransExtInfo extParam;
2320          extParam.r = extRing;
2321          assume( extRing->qideal == NULL );
2322
2323          R->cf = nInitChar(n_transExt, &extParam);
2324        }
2325      }
2326    }
2327  }
2328  else
2329  {
2330    WerrorS("coefficient field must be described by `int` or `list`");
2331    goto rCompose_err;
2332  }
2333
2334  if( R->cf == NULL )
2335  {
2336    WerrorS("could not create coefficient field described by the input!");
2337    goto rCompose_err;
2338  }
2339
2340  // ------------------------- VARS ---------------------------
2341  if (L->m[1].Typ()==LIST_CMD)
2342  {
2343    lists v=(lists)L->m[1].Data();
2344    R->N = v->nr+1;
2345    if (R->N<=0)
2346    {
2347      WerrorS("no ring variables");
2348      goto rCompose_err;
2349    }
2350    R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
2351    int i;
2352    for(i=0;i<R->N;i++)
2353    {
2354      if (v->m[i].Typ()==STRING_CMD)
2355        R->names[i]=omStrDup((char *)v->m[i].Data());
2356      else if (v->m[i].Typ()==POLY_CMD)
2357      {
2358        poly p=(poly)v->m[i].Data();
2359        int nr=pIsPurePower(p);
2360        if (nr>0)
2361          R->names[i]=omStrDup(currRing->names[nr-1]);
2362        else
2363        {
2364          Werror("var name %d must be a string or a ring variable",i+1);
2365          goto rCompose_err;
2366        }
2367      }
2368      else
2369      {
2370        Werror("var name %d must be `string`",i+1);
2371        goto rCompose_err;
2372      }
2373    }
2374  }
2375  else
2376  {
2377    WerrorS("variable must be given as `list`");
2378    goto rCompose_err;
2379  }
2380  // ------------------------ ORDER ------------------------------
2381  if (L->m[2].Typ()==LIST_CMD)
2382  {
2383    lists v=(lists)L->m[2].Data();
2384    int n= v->nr+2;
2385    int j;
2386    // initialize fields of R
2387    R->order=(int *)omAlloc0(n*sizeof(int));
2388    R->block0=(int *)omAlloc0(n*sizeof(int));
2389    R->block1=(int *)omAlloc0(n*sizeof(int));
2390    R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
2391    // init order, so that rBlocks works correctly
2392    for (j=0; j < n-1; j++)
2393      R->order[j] = (int) ringorder_unspec;
2394    // orderings
2395    for(j=0;j<n-1;j++)
2396    {
2397    // todo: a(..), M
2398      if (v->m[j].Typ()!=LIST_CMD)
2399      {
2400        WerrorS("ordering must be list of lists");
2401        goto rCompose_err;
2402      }
2403      lists vv=(lists)v->m[j].Data();
2404      if ((vv->nr!=1)
2405      || (vv->m[0].Typ()!=STRING_CMD)
2406      || ((vv->m[1].Typ()!=INTVEC_CMD) && (vv->m[1].Typ()!=INT_CMD)))
2407      {
2408        WerrorS("ordering name must be a (string,intvec)");
2409        goto rCompose_err;
2410      }
2411      R->order[j]=rOrderName(omStrDup((char*)vv->m[0].Data())); // assume STRING
2412
2413      if (j==0) R->block0[0]=1;
2414      else
2415      {
2416         int jj=j-1;
2417         while((jj>=0)
2418         && ((R->order[jj]== ringorder_a)
2419            || (R->order[jj]== ringorder_aa)
2420            || (R->order[jj]== ringorder_am)
2421            || (R->order[jj]== ringorder_c)
2422            || (R->order[jj]== ringorder_C)
2423            || (R->order[jj]== ringorder_s)
2424            || (R->order[jj]== ringorder_S)
2425         ))
2426         {
2427           //Print("jj=%, skip %s\n",rSimpleOrdStr(R->order[jj]));
2428           jj--;
2429         }
2430         if (jj<0) R->block0[j]=1;
2431         else       R->block0[j]=R->block1[jj]+1;
2432      }
2433      intvec *iv;
2434      if (vv->m[1].Typ()==INT_CMD)
2435        iv=new intvec((int)(long)vv->m[1].Data(),(int)(long)vv->m[1].Data());
2436      else
2437        iv=ivCopy((intvec*)vv->m[1].Data()); //assume INTVEC
2438      int iv_len=iv->length();
2439      R->block1[j]=si_max(R->block0[j],R->block0[j]+iv_len-1);
2440      if (R->block1[j]>R->N)
2441      {
2442        R->block1[j]=R->N;
2443        iv_len=R->block1[j]-R->block0[j]+1;
2444      }
2445      //Print("block %d from %d to %d\n",j,R->block0[j], R->block1[j]);
2446      int i;
2447      switch (R->order[j])
2448      {
2449         case ringorder_ws:
2450         case ringorder_Ws:
2451            R->OrdSgn=-1;
2452         case ringorder_aa:
2453         case ringorder_a:
2454         case ringorder_wp:
2455         case ringorder_Wp:
2456           R->wvhdl[j] =( int *)omAlloc(iv_len*sizeof(int));
2457           for (i=0; i<iv_len;i++)
2458           {
2459             R->wvhdl[j][i]=(*iv)[i];
2460           }
2461           break;
2462         case ringorder_am:
2463           R->wvhdl[j] =( int *)omAlloc((iv->length()+1)*sizeof(int));
2464           for (i=0; i<iv_len;i++)
2465           {
2466             R->wvhdl[j][i]=(*iv)[i];
2467           }
2468           R->wvhdl[j][i]=iv->length() - iv_len;
2469           //printf("ivlen:%d,iv->len:%d,mod:%d\n",iv_len,iv->length(),R->wvhdl[j][i]);
2470           for (; i<iv->length(); i++)
2471           {
2472              R->wvhdl[j][i+1]=(*iv)[i];
2473           }
2474           break;
2475         case ringorder_M:
2476           R->wvhdl[j] =( int *)omAlloc((iv->length())*sizeof(int));
2477           for (i=0; i<iv->length();i++) R->wvhdl[j][i]=(*iv)[i];
2478           R->block1[j]=si_max(R->block0[j],R->block0[j]+(int)sqrt((double)(iv->length()-1)));
2479           if (R->block1[j]>R->N)
2480           {
2481             WerrorS("ordering matrix too big");
2482             goto rCompose_err;
2483           }
2484           break;
2485         case ringorder_ls:
2486         case ringorder_ds:
2487         case ringorder_Ds:
2488         case ringorder_rs:
2489           R->OrdSgn=-1;
2490         case ringorder_lp:
2491         case ringorder_dp:
2492         case ringorder_Dp:
2493         case ringorder_rp:
2494           break;
2495         case ringorder_S:
2496           break;
2497         case ringorder_c:
2498         case ringorder_C:
2499           R->block1[j]=R->block0[j]=0;
2500           break;
2501
2502         case ringorder_s:
2503           break;
2504
2505         case ringorder_IS:
2506         {
2507           R->block1[j] = R->block0[j] = 0;
2508           if( iv->length() > 0 )
2509           {
2510             const int s = (*iv)[0];
2511             assume( -2 < s && s < 2 );
2512             R->block1[j] = R->block0[j] = s;
2513           }
2514           break;
2515         }
2516         case 0:
2517         case ringorder_unspec:
2518           break;
2519      }
2520      delete iv;
2521    }
2522    // sanity check
2523    j=n-2;
2524    if ((R->order[j]==ringorder_c)
2525    || (R->order[j]==ringorder_C)
2526    || (R->order[j]==ringorder_unspec)) j--;
2527    if (R->block1[j] != R->N)
2528    {
2529      if (((R->order[j]==ringorder_dp) ||
2530           (R->order[j]==ringorder_ds) ||
2531           (R->order[j]==ringorder_Dp) ||
2532           (R->order[j]==ringorder_Ds) ||
2533           (R->order[j]==ringorder_rp) ||
2534           (R->order[j]==ringorder_rs) ||
2535           (R->order[j]==ringorder_lp) ||
2536           (R->order[j]==ringorder_ls))
2537          &&
2538            R->block0[j] <= R->N)
2539      {
2540        R->block1[j] = R->N;
2541      }
2542      else
2543      {
2544        Werror("ordering incomplete: size (%d) should be %d",R->block1[j],R->N);
2545        goto rCompose_err;
2546      }
2547    }
2548    if (R->block0[j]>R->N)
2549    {
2550      Werror("not enough variables (%d) for ordering block %d, scanned so far:",R->N,j+1);
2551      for(int ii=0;ii<=j;ii++)
2552        Werror("ord[%d]: %s from v%d to v%d",ii+1,rSimpleOrdStr(R->order[ii]),R->block0[ii],R->block1[ii]);
2553      goto rCompose_err;
2554    }
2555    if (check_comp)
2556    {
2557      BOOLEAN comp_order=FALSE;
2558      int jj;
2559      for(jj=0;jj<n;jj++)
2560      {
2561        if ((R->order[jj]==ringorder_c) ||
2562            (R->order[jj]==ringorder_C)) { comp_order=TRUE; break; }
2563      }
2564      if (!comp_order)
2565      {
2566        R->order=(int*)omRealloc0Size(R->order,n*sizeof(int),(n+1)*sizeof(int));
2567        R->block0=(int*)omRealloc0Size(R->block0,n*sizeof(int),(n+1)*sizeof(int));
2568        R->block1=(int*)omRealloc0Size(R->block1,n*sizeof(int),(n+1)*sizeof(int));
2569        R->wvhdl=(int**)omRealloc0Size(R->wvhdl,n*sizeof(int_ptr),(n+1)*sizeof(int_ptr));
2570        R->order[n-1]=ringorder_C;
2571        R->block0[n-1]=0;
2572        R->block1[n-1]=0;
2573        R->wvhdl[n-1]=NULL;
2574        n++;
2575      }
2576    }
2577  }
2578  else
2579  {
2580    WerrorS("ordering must be given as `list`");
2581    goto rCompose_err;
2582  }
2583
2584  // ------------------------ ??????? --------------------
2585
2586  rRenameVars(R);
2587  rComplete(R);
2588
2589/*#ifdef HAVE_RINGS
2590// currently, coefficients which are ring elements require a global ordering:
2591  if (rField_is_Ring(R) && (R->OrdSgn==-1))
2592  {
2593    WerrorS("global ordering required for these coefficients");
2594    goto rCompose_err;
2595  }
2596#endif*/
2597
2598
2599  // ------------------------ Q-IDEAL ------------------------
2600
2601  if (L->m[3].Typ()==IDEAL_CMD)
2602  {
2603    ideal q=(ideal)L->m[3].Data();
2604    if (q->m[0]!=NULL)
2605    {
2606      if (R->cf != currRing->cf) //->cf->ch!=currRing->cf->ch)
2607      {
2608      #if 0
2609            WerrorS("coefficient fields must be equal if q-ideal !=0");
2610            goto rCompose_err;
2611      #else
2612        ring orig_ring=currRing;
2613        rChangeCurrRing(R);
2614        int *perm=NULL;
2615        int *par_perm=NULL;
2616        int par_perm_size=0;
2617        nMapFunc nMap;
2618
2619        if ((nMap=nSetMap(orig_ring->cf))==NULL)
2620        {
2621          if (rEqual(orig_ring,currRing))
2622          {
2623            nMap=n_SetMap(currRing->cf, currRing->cf);
2624          }
2625          else
2626          // Allow imap/fetch to be make an exception only for:
2627          if ( (rField_is_Q_a(orig_ring) &&  // Q(a..) -> Q(a..) || Q || Zp || Zp(a)
2628            (rField_is_Q(currRing) || rField_is_Q_a(currRing) ||
2629             (rField_is_Zp(currRing) || rField_is_Zp_a(currRing))))
2630           ||
2631           (rField_is_Zp_a(orig_ring) &&  // Zp(a..) -> Zp(a..) || Zp
2632            (rField_is_Zp(currRing, rInternalChar(orig_ring)) ||
2633             rField_is_Zp_a(currRing, rInternalChar(orig_ring)))) )
2634          {
2635            par_perm_size=rPar(orig_ring);
2636
2637//            if ((orig_ring->minpoly != NULL) || (orig_ring->qideal != NULL))
2638//              naSetChar(rInternalChar(orig_ring),orig_ring);
2639//            else ntSetChar(rInternalChar(orig_ring),orig_ring);
2640
2641            nSetChar(currRing->cf);
2642          }
2643          else
2644          {
2645            WerrorS("coefficient fields must be equal if q-ideal !=0");
2646            goto rCompose_err;
2647          }
2648        }
2649        perm=(int *)omAlloc0((orig_ring->N+1)*sizeof(int));
2650        if (par_perm_size!=0)
2651          par_perm=(int *)omAlloc0(par_perm_size*sizeof(int));
2652        int i;
2653        #if 0
2654        // use imap:
2655        maFindPerm(orig_ring->names,orig_ring->N,orig_ring->parameter,orig_ring->P,
2656          currRing->names,currRing->N,currRing->parameter, currRing->P,
2657          perm,par_perm, currRing->ch);
2658        #else
2659        // use fetch
2660        if ((rPar(orig_ring)>0) && (rPar(currRing)==0))
2661        {
2662          for(i=si_min(rPar(orig_ring),rVar(currRing))-1;i>=0;i--) par_perm[i]=i+1;
2663        }
2664        else if (par_perm_size!=0)
2665          for(i=si_min(rPar(orig_ring),rPar(currRing))-1;i>=0;i--) par_perm[i]=-(i+1);
2666        for(i=si_min(orig_ring->N,rVar(currRing));i>0;i--) perm[i]=i;
2667        #endif
2668        ideal dest_id=idInit(IDELEMS(q),1);
2669        for(i=IDELEMS(q)-1; i>=0; i--)
2670        {
2671          dest_id->m[i]=p_PermPoly(q->m[i],perm,orig_ring, currRing,nMap,
2672                                  par_perm,par_perm_size);
2673          //  PrintS("map:");pWrite(dest_id->m[i]);PrintLn();
2674          pTest(dest_id->m[i]);
2675        }
2676        R->qideal=dest_id;
2677        if (perm!=NULL)
2678          omFreeSize((ADDRESS)perm,(orig_ring->N+1)*sizeof(int));
2679        if (par_perm!=NULL)
2680          omFreeSize((ADDRESS)par_perm,par_perm_size*sizeof(int));
2681        rChangeCurrRing(orig_ring);
2682      #endif
2683      }
2684      else
2685        R->qideal=idrCopyR(q,currRing,R);
2686    }
2687  }
2688  else
2689  {
2690    WerrorS("q-ideal must be given as `ideal`");
2691    goto rCompose_err;
2692  }
2693
2694
2695  // ---------------------------------------------------------------
2696  #ifdef HAVE_PLURAL
2697  if (L->nr==5)
2698  {
2699    if (nc_CallPlural((matrix)L->m[4].Data(),
2700                      (matrix)L->m[5].Data(),
2701                      NULL,NULL,
2702                      R,
2703                      true, // !!!
2704                      true, false,
2705                      currRing, FALSE)) goto rCompose_err;
2706    // takes care about non-comm. quotient! i.e. calls "nc_SetupQuotient" due to last true
2707  }
2708  #endif
2709  return R;
2710
2711rCompose_err:
2712  if (R->N>0)
2713  {
2714    int i;
2715    if (R->names!=NULL)
2716    {
2717      i=R->N-1;
2718      while (i>=0) { if (R->names[i]!=NULL) omFree(R->names[i]); i--; }
2719      omFree(R->names);
2720    }
2721  }
2722  if (R->order!=NULL) omFree(R->order);
2723  if (R->block0!=NULL) omFree(R->block0);
2724  if (R->block1!=NULL) omFree(R->block1);
2725  if (R->wvhdl!=NULL) omFree(R->wvhdl);
2726  omFree(R);
2727  return NULL;
2728}
2729
2730// from matpol.cc
2731
2732/*2
2733* compute the jacobi matrix of an ideal
2734*/
2735BOOLEAN mpJacobi(leftv res,leftv a)
2736{
2737  int     i,j;
2738  matrix result;
2739  ideal id=(ideal)a->Data();
2740
2741  result =mpNew(IDELEMS(id),rVar(currRing));
2742  for (i=1; i<=IDELEMS(id); i++)
2743  {
2744    for (j=1; j<=rVar(currRing); j++)
2745    {
2746      MATELEM(result,i,j) = pDiff(id->m[i-1],j);
2747    }
2748  }
2749  res->data=(char *)result;
2750  return FALSE;
2751}
2752
2753/*2
2754* returns the Koszul-matrix of degree d of a vectorspace with dimension n
2755* uses the first n entrees of id, if id <> NULL
2756*/
2757BOOLEAN mpKoszul(leftv res,leftv c/*ip*/, leftv b/*in*/, leftv id)
2758{
2759  int n=(int)(long)b->Data();
2760  int d=(int)(long)c->Data();
2761  int     k,l,sign,row,col;
2762  matrix  result;
2763  ideal temp;
2764  BOOLEAN bo;
2765  poly    p;
2766
2767  if ((d>n) || (d<1) || (n<1))
2768  {
2769    res->data=(char *)mpNew(1,1);
2770    return FALSE;
2771  }
2772  int *choise = (int*)omAlloc(d*sizeof(int));
2773  if (id==NULL)
2774    temp=idMaxIdeal(1);
2775  else
2776    temp=(ideal)id->Data();
2777
2778  k = binom(n,d);
2779  l = k*d;
2780  l /= n-d+1;
2781  result =mpNew(l,k);
2782  col = 1;
2783  idInitChoise(d,1,n,&bo,choise);
2784  while (!bo)
2785  {
2786    sign = 1;
2787    for (l=1;l<=d;l++)
2788    {
2789      if (choise[l-1]<=IDELEMS(temp))
2790      {
2791        p = pCopy(temp->m[choise[l-1]-1]);
2792        if (sign == -1) p = pNeg(p);
2793        sign *= -1;
2794        row = idGetNumberOfChoise(l-1,d,1,n,choise);
2795        MATELEM(result,row,col) = p;
2796      }
2797    }
2798    col++;
2799    idGetNextChoise(d,n,&bo,choise);
2800  }
2801  if (id==NULL) idDelete(&temp);
2802
2803  res->data=(char *)result;
2804  return FALSE;
2805}
2806
2807// from syz1.cc
2808/*2
2809* read out the Betti numbers from resolution
2810* (interpreter interface)
2811*/
2812BOOLEAN syBetti2(leftv res, leftv u, leftv w)
2813{
2814  syStrategy syzstr=(syStrategy)u->Data();
2815
2816  BOOLEAN minim=(int)(long)w->Data();
2817  int row_shift=0;
2818  int add_row_shift=0;
2819  intvec *weights=NULL;
2820  intvec *ww=(intvec *)atGet(u,"isHomog",INTVEC_CMD);
2821  if (ww!=NULL)
2822  {
2823     weights=ivCopy(ww);
2824     add_row_shift = ww->min_in();
2825     (*weights) -= add_row_shift;
2826  }
2827
2828  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift,weights);
2829  //row_shift += add_row_shift;
2830  //Print("row_shift=%d, add_row_shift=%d\n",row_shift,add_row_shift);
2831  atSet(res,omStrDup("rowShift"),(void*)(long)add_row_shift,INT_CMD);
2832
2833  return FALSE;
2834}
2835BOOLEAN syBetti1(leftv res, leftv u)
2836{
2837  sleftv tmp;
2838  memset(&tmp,0,sizeof(tmp));
2839  tmp.rtyp=INT_CMD;
2840  tmp.data=(void *)1;
2841  return syBetti2(res,u,&tmp);
2842}
2843
2844/*3
2845* converts a resolution into a list of modules
2846*/
2847lists syConvRes(syStrategy syzstr,BOOLEAN toDel,int add_row_shift)
2848{
2849  resolvente fullres = syzstr->fullres;
2850  resolvente minres = syzstr->minres;
2851
2852  const int length = syzstr->length;
2853
2854  if ((fullres==NULL) && (minres==NULL))
2855  {
2856    if (syzstr->hilb_coeffs==NULL)
2857    { // La Scala
2858      fullres = syReorder(syzstr->res, length, syzstr);
2859    }
2860    else
2861    { // HRES
2862      minres = syReorder(syzstr->orderedRes, length, syzstr);
2863      syKillEmptyEntres(minres, length);
2864    }
2865  }
2866
2867  resolvente tr;
2868  int typ0=IDEAL_CMD;
2869
2870  if (minres!=NULL)
2871    tr = minres;
2872  else
2873    tr = fullres;
2874
2875  resolvente trueres=NULL; intvec ** w=NULL;
2876
2877  if (length>0)
2878  {
2879    trueres = (resolvente)omAlloc0((length)*sizeof(ideal));
2880    for (int i=(length)-1;i>=0;i--)
2881    {
2882      if (tr[i]!=NULL)
2883      {
2884        trueres[i] = idCopy(tr[i]);
2885      }
2886    }
2887    if ( id_RankFreeModule(trueres[0], currRing) > 0)
2888      typ0 = MODUL_CMD;
2889    if (syzstr->weights!=NULL)
2890    {
2891      w = (intvec**)omAlloc0(length*sizeof(intvec*));
2892      for (int i=length-1;i>=0;i--)
2893      {
2894        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
2895      }
2896    }
2897  }
2898
2899  lists li = liMakeResolv(trueres, length, syzstr->list_length,typ0,
2900                          w, add_row_shift);
2901
2902  if (w != NULL) omFreeSize(w, length*sizeof(intvec*));
2903
2904  if (toDel)
2905    syKillComputation(syzstr);
2906  else
2907  {
2908    if( fullres != NULL && syzstr->fullres == NULL )
2909      syzstr->fullres = fullres;
2910
2911    if( minres != NULL && syzstr->minres == NULL )
2912      syzstr->minres = minres;
2913  }
2914
2915  return li;
2916
2917
2918}
2919
2920/*3
2921* converts a list of modules into a resolution
2922*/
2923syStrategy syConvList(lists li,BOOLEAN toDel)
2924{
2925  int typ0;
2926  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2927
2928  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
2929  if (fr != NULL)
2930  {
2931
2932    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2933    for (int i=result->length-1;i>=0;i--)
2934    {
2935      if (fr[i]!=NULL)
2936        result->fullres[i] = idCopy(fr[i]);
2937    }
2938    result->list_length=result->length;
2939    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2940  }
2941  else
2942  {
2943    omFreeSize(result, sizeof(ssyStrategy));
2944    result = NULL;
2945  }
2946  if (toDel) li->Clean();
2947  return result;
2948}
2949
2950/*3
2951* converts a list of modules into a minimal resolution
2952*/
2953syStrategy syForceMin(lists li)
2954{
2955  int typ0;
2956  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2957
2958  resolvente fr = liFindRes(li,&(result->length),&typ0);
2959  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2960  for (int i=result->length-1;i>=0;i--)
2961  {
2962    if (fr[i]!=NULL)
2963      result->minres[i] = idCopy(fr[i]);
2964  }
2965  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2966  return result;
2967}
2968// from weight.cc
2969BOOLEAN kWeight(leftv res,leftv id)
2970{
2971  ideal F=(ideal)id->Data();
2972  intvec * iv = new intvec(rVar(currRing));
2973  polyset s;
2974  int  sl, n, i;
2975  int  *x;
2976
2977  res->data=(char *)iv;
2978  s = F->m;
2979  sl = IDELEMS(F) - 1;
2980  n = rVar(currRing);
2981  double wNsqr = (double)2.0 / (double)n;
2982  wFunctional = wFunctionalBuch;
2983  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
2984  wCall(s, sl, x, wNsqr, currRing);
2985  for (i = n; i!=0; i--)
2986    (*iv)[i-1] = x[i + n + 1];
2987  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
2988  return FALSE;
2989}
2990
2991BOOLEAN kQHWeight(leftv res,leftv v)
2992{
2993  res->data=(char *)id_QHomWeight((ideal)v->Data(), currRing);
2994  if (res->data==NULL)
2995    res->data=(char *)new intvec(rVar(currRing));
2996  return FALSE;
2997}
2998/*==============================================================*/
2999// from clapsing.cc
3000#if 0
3001BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
3002{
3003  BOOLEAN b=singclap_factorize((poly)(u->CopyD()), &v, 0);
3004  res->data=(void *)b;
3005}
3006#endif
3007
3008BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
3009{
3010  res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
3011                  (poly)w->CopyD(), currRing);
3012  return errorreported;
3013}
3014
3015BOOLEAN jjCHARSERIES(leftv res, leftv u)
3016{
3017  res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
3018  return (res->data==NULL);
3019}
3020
3021// from semic.cc
3022#ifdef HAVE_SPECTRUM
3023
3024// ----------------------------------------------------------------------------
3025//  Initialize a  spectrum  deep from a  singular  lists
3026// ----------------------------------------------------------------------------
3027
3028void copy_deep( spectrum& spec, lists l )
3029{
3030    spec.mu = (int)(long)(l->m[0].Data( ));
3031    spec.pg = (int)(long)(l->m[1].Data( ));
3032    spec.= (int)(long)(l->m[2].Data( ));
3033
3034    spec.copy_new( spec.n );
3035
3036    intvec  *num = (intvec*)l->m[3].Data( );
3037    intvec  *den = (intvec*)l->m[4].Data( );
3038    intvec  *mul = (intvec*)l->m[5].Data( );
3039
3040    for( int i=0; i<spec.n; i++ )
3041    {
3042        spec.s[i] = (Rational)((*num)[i])/(Rational)((*den)[i]);
3043        spec.w[i] = (*mul)[i];
3044    }
3045}
3046
3047// ----------------------------------------------------------------------------
3048//  singular lists  constructor for  spectrum
3049// ----------------------------------------------------------------------------
3050
3051spectrum /*former spectrum::spectrum ( lists l )*/
3052spectrumFromList( lists l )
3053{
3054    spectrum result;
3055    copy_deep( result, l );
3056    return result;
3057}
3058
3059// ----------------------------------------------------------------------------
3060//  generate a Singular  lists  from a spectrum
3061// ----------------------------------------------------------------------------
3062
3063/* former spectrum::thelist ( void )*/
3064lists   getList( spectrum& spec )
3065{
3066    lists   L  = (lists)omAllocBin( slists_bin);
3067
3068    L->Init( 6 );
3069
3070    intvec            *num  = new intvec( spec.n );
3071    intvec            *den  = new intvec( spec.n );
3072    intvec            *mult = new intvec( spec.n );
3073
3074    for( int i=0; i<spec.n; i++ )
3075    {
3076        (*num) [i] = spec.s[i].get_num_si( );
3077        (*den) [i] = spec.s[i].get_den_si( );
3078        (*mult)[i] = spec.w[i];
3079    }
3080
3081    L->m[0].rtyp = INT_CMD;    //  milnor number
3082    L->m[1].rtyp = INT_CMD;    //  geometrical genus
3083    L->m[2].rtyp = INT_CMD;    //  # of spectrum numbers
3084    L->m[3].rtyp = INTVEC_CMD; //  numerators
3085    L->m[4].rtyp = INTVEC_CMD; //  denomiantors
3086    L->m[5].rtyp = INTVEC_CMD; //  multiplicities
3087
3088    L->m[0].data = (void*)(long)spec.mu;
3089    L->m[1].data = (void*)(long)spec.pg;
3090    L->m[2].data = (void*)(long)spec.n;
3091    L->m[3].data = (void*)num;
3092    L->m[4].data = (void*)den;
3093    L->m[5].data = (void*)mult;
3094
3095    return  L;
3096}
3097// from spectrum.cc
3098// ----------------------------------------------------------------------------
3099//  print out an error message for a spectrum list
3100// ----------------------------------------------------------------------------
3101
3102typedef enum
3103{
3104    semicOK,
3105    semicMulNegative,
3106
3107    semicListTooShort,
3108    semicListTooLong,
3109
3110    semicListFirstElementWrongType,
3111    semicListSecondElementWrongType,
3112    semicListThirdElementWrongType,
3113    semicListFourthElementWrongType,
3114    semicListFifthElementWrongType,
3115    semicListSixthElementWrongType,
3116
3117    semicListNNegative,
3118    semicListWrongNumberOfNumerators,
3119    semicListWrongNumberOfDenominators,
3120    semicListWrongNumberOfMultiplicities,
3121
3122    semicListMuNegative,
3123    semicListPgNegative,
3124    semicListNumNegative,
3125    semicListDenNegative,
3126    semicListMulNegative,
3127
3128    semicListNotSymmetric,
3129    semicListNotMonotonous,
3130
3131    semicListMilnorWrong,
3132    semicListPGWrong
3133
3134} semicState;
3135
3136void    list_error( semicState state )
3137{
3138    switch( state )
3139    {
3140        case semicListTooShort:
3141            WerrorS( "the list is too short" );
3142            break;
3143        case semicListTooLong:
3144            WerrorS( "the list is too long" );
3145            break;
3146
3147        case semicListFirstElementWrongType:
3148            WerrorS( "first element of the list should be int" );
3149            break;
3150        case semicListSecondElementWrongType:
3151            WerrorS( "second element of the list should be int" );
3152            break;
3153        case semicListThirdElementWrongType:
3154            WerrorS( "third element of the list should be int" );
3155            break;
3156        case semicListFourthElementWrongType:
3157            WerrorS( "fourth element of the list should be intvec" );
3158            break;
3159        case semicListFifthElementWrongType:
3160            WerrorS( "fifth element of the list should be intvec" );
3161            break;
3162        case semicListSixthElementWrongType:
3163            WerrorS( "sixth element of the list should be intvec" );
3164            break;
3165
3166        case semicListNNegative:
3167            WerrorS( "first element of the list should be positive" );
3168            break;
3169        case semicListWrongNumberOfNumerators:
3170            WerrorS( "wrong number of numerators" );
3171            break;
3172        case semicListWrongNumberOfDenominators:
3173            WerrorS( "wrong number of denominators" );
3174            break;
3175        case semicListWrongNumberOfMultiplicities:
3176            WerrorS( "wrong number of multiplicities" );
3177            break;
3178
3179        case semicListMuNegative:
3180            WerrorS( "the Milnor number should be positive" );
3181            break;
3182        case semicListPgNegative:
3183            WerrorS( "the geometrical genus should be nonnegative" );
3184            break;
3185        case semicListNumNegative:
3186            WerrorS( "all numerators should be positive" );
3187            break;
3188        case semicListDenNegative:
3189            WerrorS( "all denominators should be positive" );
3190            break;
3191        case semicListMulNegative:
3192            WerrorS( "all multiplicities should be positive" );
3193            break;
3194
3195        case semicListNotSymmetric:
3196            WerrorS( "it is not symmetric" );
3197            break;
3198        case semicListNotMonotonous:
3199            WerrorS( "it is not monotonous" );
3200            break;
3201
3202        case semicListMilnorWrong:
3203            WerrorS( "the Milnor number is wrong" );
3204            break;
3205        case semicListPGWrong:
3206            WerrorS( "the geometrical genus is wrong" );
3207            break;
3208
3209        default:
3210            WerrorS( "unspecific error" );
3211            break;
3212    }
3213}
3214// ----------------------------------------------------------------------------
3215//  this is the main spectrum computation function
3216// ----------------------------------------------------------------------------
3217
3218enum    spectrumState
3219{
3220    spectrumOK,
3221    spectrumZero,
3222    spectrumBadPoly,
3223    spectrumNoSingularity,
3224    spectrumNotIsolated,
3225    spectrumDegenerate,
3226    spectrumWrongRing,
3227    spectrumNoHC,
3228    spectrumUnspecErr
3229};
3230
3231// from splist.cc
3232// ----------------------------------------------------------------------------
3233//  Compute the spectrum of a  spectrumPolyList
3234// ----------------------------------------------------------------------------
3235
3236/* former spectrumPolyList::spectrum ( lists*, int) */
3237spectrumState   spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3238{
3239  spectrumPolyNode  **node = &speclist.root;
3240  spectrumPolyNode  *search;
3241
3242  poly              f,tmp;
3243  int               found,cmp;
3244
3245  Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3246                 ( fast==2 ? 2 : 1 ) );
3247
3248  Rational weight_prev( 0,1 );
3249
3250  int     mu = 0;          // the milnor number
3251  int     pg = 0;          // the geometrical genus
3252  int     n  = 0;          // number of different spectral numbers
3253  int     z  = 0;          // number of spectral number equal to smax
3254
3255  while( (*node)!=(spectrumPolyNode*)NULL &&
3256         ( fast==0 || (*node)->weight<=smax ) )
3257  {
3258        // ---------------------------------------
3259        //  determine the first normal form which
3260        //  contains the monomial  node->mon
3261        // ---------------------------------------
3262
3263    found  = FALSE;
3264    search = *node;
3265
3266    while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3267    {
3268      if( search->nf!=(poly)NULL )
3269      {
3270        f = search->nf;
3271
3272        do
3273        {
3274                    // --------------------------------
3275                    //  look for  (*node)->mon  in   f
3276                    // --------------------------------
3277
3278          cmp = pCmp( (*node)->mon,f );
3279
3280          if( cmp<0 )
3281          {
3282            f = pNext( f );
3283          }
3284          else if( cmp==0 )
3285          {
3286                        // -----------------------------
3287                        //  we have found a normal form
3288                        // -----------------------------
3289
3290            found = TRUE;
3291
3292                        //  normalize coefficient
3293
3294            number inv = nInvers( pGetCoeff( f ) );
3295            pMult_nn( search->nf,inv );
3296            nDelete( &inv );
3297
3298                        //  exchange  normal forms
3299
3300            tmp         = (*node)->nf;
3301            (*node)->nf = search->nf;
3302            search->nf  = tmp;
3303          }
3304        }
3305        while( cmp<0 && f!=(poly)NULL );
3306      }
3307      search = search->next;
3308    }
3309
3310    if( found==FALSE )
3311    {
3312            // ------------------------------------------------
3313            //  the weight of  node->mon  is a spectrum number
3314            // ------------------------------------------------
3315
3316      mu++;
3317
3318      if( (*node)->weight<=(Rational)1 )              pg++;
3319      if( (*node)->weight==smax )           z++;
3320      if( (*node)->weight>weight_prev )     n++;
3321
3322      weight_prev = (*node)->weight;
3323      node = &((*node)->next);
3324    }
3325    else
3326    {
3327            // -----------------------------------------------
3328            //  determine all other normal form which contain
3329            //  the monomial  node->mon
3330            //  replace for  node->mon  its normal form
3331            // -----------------------------------------------
3332
3333      while( search!=(spectrumPolyNode*)NULL )
3334      {
3335        if( search->nf!=(poly)NULL )
3336        {
3337          f = search->nf;
3338
3339          do
3340          {
3341                        // --------------------------------
3342                        //  look for  (*node)->mon  in   f
3343                        // --------------------------------
3344
3345            cmp = pCmp( (*node)->mon,f );
3346
3347            if( cmp<0 )
3348            {
3349              f = pNext( f );
3350            }
3351            else if( cmp==0 )
3352            {
3353              search->nf = pSub( search->nf,
3354                                 ppMult_nn( (*node)->nf,pGetCoeff( f ) ) );
3355              pNorm( search->nf );
3356            }
3357          }
3358          while( cmp<0 && f!=(poly)NULL );
3359        }
3360        search = search->next;
3361      }
3362      speclist.delete_node( node );
3363    }
3364
3365  }
3366
3367    // --------------------------------------------------------
3368    //  fast computation exploits the symmetry of the spectrum
3369    // --------------------------------------------------------
3370
3371  if( fast==2 )
3372  {
3373    mu = 2*mu - z;
3374    n  = ( z > 0 ? 2*n - 1 : 2*n );
3375  }
3376
3377    // --------------------------------------------------------
3378    //  compute the spectrum numbers with their multiplicities
3379    // --------------------------------------------------------
3380
3381  intvec            *nom  = new intvec( n );
3382  intvec            *den  = new intvec( n );
3383  intvec            *mult = new intvec( n );
3384
3385  int count         = 0;
3386  int multiplicity  = 1;
3387
3388  for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3389              ( fast==0 || search->weight<=smax );
3390       search=search->next )
3391  {
3392    if( search->next==(spectrumPolyNode*)NULL ||
3393        search->weight<search->next->weight )
3394    {
3395      (*nom) [count] = search->weight.get_num_si( );
3396      (*den) [count] = search->weight.get_den_si( );
3397      (*mult)[count] = multiplicity;
3398
3399      multiplicity=1;
3400      count++;
3401    }
3402    else
3403    {
3404      multiplicity++;
3405    }
3406  }
3407
3408    // --------------------------------------------------------
3409    //  fast computation exploits the symmetry of the spectrum
3410    // --------------------------------------------------------
3411
3412  if( fast==2 )
3413  {
3414    int n1,n2;
3415    for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3416    {
3417      (*nom) [n2] = rVar(currRing)*(*den)[n1]-(*nom)[n1];
3418      (*den) [n2] = (*den)[n1];
3419      (*mult)[n2] = (*mult)[n1];
3420    }
3421  }
3422
3423    // -----------------------------------
3424    //  test if the spectrum is symmetric
3425    // -----------------------------------
3426
3427  if( fast==0 || fast==1 )
3428  {
3429    int symmetric=TRUE;
3430
3431    for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3432    {
3433      if( (*mult)[n1]!=(*mult)[n2] ||
3434          (*den) [n1]!= (*den)[n2] ||
3435          (*nom)[n1]+(*nom)[n2]!=rVar(currRing)*(*den) [n1] )
3436      {
3437        symmetric = FALSE;
3438      }
3439    }
3440
3441    if( symmetric==FALSE )
3442    {
3443            // ---------------------------------------------
3444            //  the spectrum is not symmetric => degenerate
3445            //  principal part
3446            // ---------------------------------------------
3447
3448      *L = (lists)omAllocBin( slists_bin);
3449      (*L)->Init( 1 );
3450      (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3451      (*L)->m[0].data = (void*)(long)mu;
3452
3453      return spectrumDegenerate;
3454    }
3455  }
3456
3457  *L = (lists)omAllocBin( slists_bin);
3458
3459  (*L)->Init( 6 );
3460
3461  (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3462  (*L)->m[1].rtyp = INT_CMD;    //  geometrical genus
3463  (*L)->m[2].rtyp = INT_CMD;    //  number of spectrum values
3464  (*L)->m[3].rtyp = INTVEC_CMD; //  nominators
3465  (*L)->m[4].rtyp = INTVEC_CMD; //  denomiantors
3466  (*L)->m[5].rtyp = INTVEC_CMD; //  multiplicities
3467
3468  (*L)->m[0].data = (void*)(long)mu;
3469  (*L)->m[1].data = (void*)(long)pg;
3470  (*L)->m[2].data = (void*)(long)n;
3471  (*L)->m[3].data = (void*)nom;
3472  (*L)->m[4].data = (void*)den;
3473  (*L)->m[5].data = (void*)mult;
3474
3475  return  spectrumOK;
3476}
3477
3478spectrumState   spectrumCompute( poly h,lists *L,int fast )
3479{
3480  int i;
3481
3482  #ifdef SPECTRUM_DEBUG
3483  #ifdef SPECTRUM_PRINT
3484  #ifdef SPECTRUM_IOSTREAM
3485    cout << "spectrumCompute\n";
3486    if( fast==0 ) cout << "    no optimization" << endl;
3487    if( fast==1 ) cout << "    weight optimization" << endl;
3488    if( fast==2 ) cout << "    symmetry optimization" << endl;
3489  #else
3490    fprintf( stdout,"spectrumCompute\n" );
3491    if( fast==0 ) fprintf( stdout,"    no optimization\n" );
3492    if( fast==1 ) fprintf( stdout,"    weight optimization\n" );
3493    if( fast==2 ) fprintf( stdout,"    symmetry optimization\n" );
3494  #endif
3495  #endif
3496  #endif
3497
3498  // ----------------------
3499  //  check if  h  is zero
3500  // ----------------------
3501
3502  if( h==(poly)NULL )
3503  {
3504    return  spectrumZero;
3505  }
3506
3507  // ----------------------------------
3508  //  check if  h  has a constant term
3509  // ----------------------------------
3510
3511  if( hasConstTerm( h, currRing ) )
3512  {
3513    return  spectrumBadPoly;
3514  }
3515
3516  // --------------------------------
3517  //  check if  h  has a linear term
3518  // --------------------------------
3519
3520  if( hasLinearTerm( h, currRing ) )
3521  {
3522    *L = (lists)omAllocBin( slists_bin);
3523    (*L)->Init( 1 );
3524    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3525    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3526
3527    return  spectrumNoSingularity;
3528  }
3529
3530  // ----------------------------------
3531  //  compute the jacobi ideal of  (h)
3532  // ----------------------------------
3533
3534  ideal J = NULL;
3535  J = idInit( rVar(currRing),1 );
3536
3537  #ifdef SPECTRUM_DEBUG
3538  #ifdef SPECTRUM_PRINT
3539  #ifdef SPECTRUM_IOSTREAM
3540    cout << "\n   computing the Jacobi ideal...\n";
3541  #else
3542    fprintf( stdout,"\n   computing the Jacobi ideal...\n" );
3543  #endif
3544  #endif
3545  #endif
3546
3547  for( i=0; i<rVar(currRing); i++ )
3548  {
3549    J->m[i] = pDiff( h,i+1); //j );
3550
3551    #ifdef SPECTRUM_DEBUG
3552    #ifdef SPECTRUM_PRINT
3553    #ifdef SPECTRUM_IOSTREAM
3554      cout << "        ";
3555    #else
3556      fprintf( stdout,"        " );
3557    #endif
3558      pWrite( J->m[i] );
3559    #endif
3560    #endif
3561  }
3562
3563  // --------------------------------------------
3564  //  compute a standard basis  stdJ  of  jac(h)
3565  // --------------------------------------------
3566
3567  #ifdef SPECTRUM_DEBUG
3568  #ifdef SPECTRUM_PRINT
3569  #ifdef SPECTRUM_IOSTREAM
3570    cout << endl;
3571    cout << "    computing a standard basis..." << endl;
3572  #else
3573    fprintf( stdout,"\n" );
3574    fprintf( stdout,"    computing a standard basis...\n" );
3575  #endif
3576  #endif
3577  #endif
3578
3579  ideal stdJ = kStd(J,currRing->qideal,isNotHomog,NULL);
3580  idSkipZeroes( stdJ );
3581
3582  #ifdef SPECTRUM_DEBUG
3583  #ifdef SPECTRUM_PRINT
3584    for( i=0; i<IDELEMS(stdJ); i++ )
3585    {
3586      #ifdef SPECTRUM_IOSTREAM
3587        cout << "        ";
3588      #else
3589        fprintf( stdout,"        " );
3590      #endif
3591
3592      pWrite( stdJ->m[i] );
3593    }
3594  #endif
3595  #endif
3596
3597  idDelete( &J );
3598
3599  // ------------------------------------------
3600  //  check if the  h  has a singularity
3601  // ------------------------------------------
3602
3603  if( hasOne( stdJ, currRing ) )
3604  {
3605    // -------------------------------
3606    //  h is smooth in the origin
3607    //  return only the Milnor number
3608    // -------------------------------
3609
3610    *L = (lists)omAllocBin( slists_bin);
3611    (*L)->Init( 1 );
3612    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3613    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3614
3615    return  spectrumNoSingularity;
3616  }
3617
3618  // ------------------------------------------
3619  //  check if the singularity  h  is isolated
3620  // ------------------------------------------
3621
3622  for( i=rVar(currRing); i>0; i-- )
3623  {
3624    if( hasAxis( stdJ,i, currRing )==FALSE )
3625    {
3626      return  spectrumNotIsolated;
3627    }
3628  }
3629
3630  // ------------------------------------------
3631  //  compute the highest corner  hc  of  stdJ
3632  // ------------------------------------------
3633
3634  #ifdef SPECTRUM_DEBUG
3635  #ifdef SPECTRUM_PRINT
3636  #ifdef SPECTRUM_IOSTREAM
3637    cout << "\n    computing the highest corner...\n";
3638  #else
3639    fprintf( stdout,"\n    computing the highest corner...\n" );
3640  #endif
3641  #endif
3642  #endif
3643
3644  poly hc = (poly)NULL;
3645
3646  scComputeHC( stdJ,currRing->qideal, 0,hc );
3647
3648  if( hc!=(poly)NULL )
3649  {
3650    pGetCoeff(hc) = nInit(1);
3651
3652    for( i=rVar(currRing); i>0; i-- )
3653    {
3654      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3655    }
3656    pSetm( hc );
3657  }
3658  else
3659  {
3660    return  spectrumNoHC;
3661  }
3662
3663  #ifdef SPECTRUM_DEBUG
3664  #ifdef SPECTRUM_PRINT
3665  #ifdef SPECTRUM_IOSTREAM
3666    cout << "       ";
3667  #else
3668    fprintf( stdout,"       " );
3669  #endif
3670    pWrite( hc );
3671  #endif
3672  #endif
3673
3674  // ----------------------------------------
3675  //  compute the Newton polygon  nph  of  h
3676  // ----------------------------------------
3677
3678  #ifdef SPECTRUM_DEBUG
3679  #ifdef SPECTRUM_PRINT
3680  #ifdef SPECTRUM_IOSTREAM
3681    cout << "\n    computing the newton polygon...\n";
3682  #else
3683    fprintf( stdout,"\n    computing the newton polygon...\n" );
3684  #endif
3685  #endif
3686  #endif
3687
3688  newtonPolygon nph( h, currRing );
3689
3690  #ifdef SPECTRUM_DEBUG
3691  #ifdef SPECTRUM_PRINT
3692    cout << nph;
3693  #endif
3694  #endif
3695
3696  // -----------------------------------------------
3697  //  compute the weight corner  wc  of  (stdj,nph)
3698  // -----------------------------------------------
3699
3700  #ifdef SPECTRUM_DEBUG
3701  #ifdef SPECTRUM_PRINT
3702  #ifdef SPECTRUM_IOSTREAM
3703    cout << "\n    computing the weight corner...\n";
3704  #else
3705    fprintf( stdout,"\n    computing the weight corner...\n" );
3706  #endif
3707  #endif
3708  #endif
3709
3710  poly    wc = ( fast==0 ? pCopy( hc ) :
3711               ( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
3712              /* fast==2 */computeWC( nph,
3713                      ((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
3714
3715  #ifdef SPECTRUM_DEBUG
3716  #ifdef SPECTRUM_PRINT
3717  #ifdef SPECTRUM_IOSTREAM
3718    cout << "        ";
3719  #else
3720    fprintf( stdout,"        " );
3721  #endif
3722    pWrite( wc );
3723  #endif
3724  #endif
3725
3726  // -------------
3727  //  compute  NF
3728  // -------------
3729
3730  #ifdef SPECTRUM_DEBUG
3731  #ifdef SPECTRUM_PRINT
3732  #ifdef SPECTRUM_IOSTREAM
3733    cout << "\n    computing NF...\n" << endl;
3734  #else
3735    fprintf( stdout,"\n    computing NF...\n" );
3736  #endif
3737  #endif
3738  #endif
3739
3740  spectrumPolyList NF( &nph );
3741
3742  computeNF( stdJ,hc,wc,&NF, currRing );
3743
3744  #ifdef SPECTRUM_DEBUG
3745  #ifdef SPECTRUM_PRINT
3746    cout << NF;
3747  #ifdef SPECTRUM_IOSTREAM
3748    cout << endl;
3749  #else
3750    fprintf( stdout,"\n" );
3751  #endif
3752  #endif
3753  #endif
3754
3755  // ----------------------------
3756  //  compute the spectrum of  h
3757  // ----------------------------
3758//  spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
3759
3760  return spectrumStateFromList(NF, L, fast );
3761}
3762
3763// ----------------------------------------------------------------------------
3764//  this procedure is called from the interpreter
3765// ----------------------------------------------------------------------------
3766//  first  = polynomial
3767//  result = list of spectrum numbers
3768// ----------------------------------------------------------------------------
3769
3770void spectrumPrintError(spectrumState state)
3771{
3772  switch( state )
3773  {
3774    case spectrumZero:
3775      WerrorS( "polynomial is zero" );
3776      break;
3777    case spectrumBadPoly:
3778      WerrorS( "polynomial has constant term" );
3779      break;
3780    case spectrumNoSingularity:
3781      WerrorS( "not a singularity" );
3782      break;
3783    case spectrumNotIsolated:
3784      WerrorS( "the singularity is not isolated" );
3785      break;
3786    case spectrumNoHC:
3787      WerrorS( "highest corner cannot be computed" );
3788      break;
3789    case spectrumDegenerate:
3790      WerrorS( "principal part is degenerate" );
3791      break;
3792    case spectrumOK:
3793      break;
3794
3795    default:
3796      WerrorS( "unknown error occurred" );
3797      break;
3798  }
3799}
3800
3801BOOLEAN spectrumProc( leftv result,leftv first )
3802{
3803  spectrumState state = spectrumOK;
3804
3805  // -------------------
3806  //  check consistency
3807  // -------------------
3808
3809  //  check for a local ring
3810
3811  if( !ringIsLocal(currRing ) )
3812  {
3813    WerrorS( "only works for local orderings" );
3814    state = spectrumWrongRing;
3815  }
3816
3817  //  no quotient rings are allowed
3818
3819  else if( currRing->qideal != NULL )
3820  {
3821    WerrorS( "does not work in quotient rings" );
3822    state = spectrumWrongRing;
3823  }
3824  else
3825  {
3826    lists   L    = (lists)NULL;
3827    int     flag = 1; // weight corner optimization is safe
3828
3829    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3830
3831    if( state==spectrumOK )
3832    {
3833      result->rtyp = LIST_CMD;
3834      result->data = (char*)L;
3835    }
3836    else
3837    {
3838      spectrumPrintError(state);
3839    }
3840  }
3841
3842  return  (state!=spectrumOK);
3843}
3844
3845// ----------------------------------------------------------------------------
3846//  this procedure is called from the interpreter
3847// ----------------------------------------------------------------------------
3848//  first  = polynomial
3849//  result = list of spectrum numbers
3850// ----------------------------------------------------------------------------
3851
3852BOOLEAN spectrumfProc( leftv result,leftv first )
3853{
3854  spectrumState state = spectrumOK;
3855
3856  // -------------------
3857  //  check consistency
3858  // -------------------
3859
3860  //  check for a local polynomial ring
3861
3862  if( currRing->OrdSgn != -1 )
3863  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
3864  // or should we use:
3865  //if( !ringIsLocal( ) )
3866  {
3867    WerrorS( "only works for local orderings" );
3868    state = spectrumWrongRing;
3869  }
3870  else if( currRing->qideal != NULL )
3871  {
3872    WerrorS( "does not work in quotient rings" );
3873    state = spectrumWrongRing;
3874  }
3875  else
3876  {
3877    lists   L    = (lists)NULL;
3878    int     flag = 2; // symmetric optimization
3879
3880    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3881
3882    if( state==spectrumOK )
3883    {
3884      result->rtyp = LIST_CMD;
3885      result->data = (char*)L;
3886    }
3887    else
3888    {
3889      spectrumPrintError(state);
3890    }
3891  }
3892
3893  return  (state!=spectrumOK);
3894}
3895
3896// ----------------------------------------------------------------------------
3897//  check if a list is a spectrum
3898//  check for:
3899//      list has 6 elements
3900//      1st element is int (mu=Milnor number)
3901//      2nd element is int (pg=geometrical genus)
3902//      3rd element is int (n =number of different spectrum numbers)
3903//      4th element is intvec (num=numerators)
3904//      5th element is intvec (den=denomiantors)
3905//      6th element is intvec (mul=multiplicities)
3906//      exactly n numerators
3907//      exactly n denominators
3908//      exactly n multiplicities
3909//      mu>0
3910//      pg>=0
3911//      n>0
3912//      num>0
3913//      den>0
3914//      mul>0
3915//      symmetriy with respect to numberofvariables/2
3916//      monotony
3917//      mu = sum of all multiplicities
3918//      pg = sum of all multiplicities where num/den<=1
3919// ----------------------------------------------------------------------------
3920
3921semicState  list_is_spectrum( lists l )
3922{
3923    // -------------------
3924    //  check list length
3925    // -------------------
3926
3927    if( l->nr < 5 )
3928    {
3929        return  semicListTooShort;
3930    }
3931    else if( l->nr > 5 )
3932    {
3933        return  semicListTooLong;
3934    }
3935
3936    // -------------
3937    //  check types
3938    // -------------
3939
3940    if( l->m[0].rtyp != INT_CMD )
3941    {
3942        return  semicListFirstElementWrongType;
3943    }
3944    else if( l->m[1].rtyp != INT_CMD )
3945    {
3946        return  semicListSecondElementWrongType;
3947    }
3948    else if( l->m[2].rtyp != INT_CMD )
3949    {
3950        return  semicListThirdElementWrongType;
3951    }
3952    else if( l->m[3].rtyp != INTVEC_CMD )
3953    {
3954        return  semicListFourthElementWrongType;
3955    }
3956    else if( l->m[4].rtyp != INTVEC_CMD )
3957    {
3958        return  semicListFifthElementWrongType;
3959    }
3960    else if( l->m[5].rtyp != INTVEC_CMD )
3961    {
3962        return  semicListSixthElementWrongType;
3963    }
3964
3965    // -------------------------
3966    //  check number of entries
3967    // -------------------------
3968
3969    int     mu = (int)(long)(l->m[0].Data( ));
3970    int     pg = (int)(long)(l->m[1].Data( ));
3971    int     n  = (int)(long)(l->m[2].Data( ));
3972
3973    if( n <= 0 )
3974    {
3975        return  semicListNNegative;
3976    }
3977
3978    intvec  *num = (intvec*)l->m[3].Data( );
3979    intvec  *den = (intvec*)l->m[4].Data( );
3980    intvec  *mul = (intvec*)l->m[5].Data( );
3981
3982    if( n != num->length( ) )
3983    {
3984        return  semicListWrongNumberOfNumerators;
3985    }
3986    else if( n != den->length( ) )
3987    {
3988        return  semicListWrongNumberOfDenominators;
3989    }
3990    else if( n != mul->length( ) )
3991    {
3992        return  semicListWrongNumberOfMultiplicities;
3993    }
3994
3995    // --------
3996    //  values
3997    // --------
3998
3999    if( mu <= 0 )
4000    {
4001        return  semicListMuNegative;
4002    }
4003    if( pg < 0 )
4004    {
4005        return  semicListPgNegative;
4006    }
4007
4008    int i;
4009
4010    for( i=0; i<n; i++ )
4011    {
4012        if( (*num)[i] <= 0 )
4013        {
4014            return  semicListNumNegative;
4015        }
4016        if( (*den)[i] <= 0 )
4017        {
4018            return  semicListDenNegative;
4019        }
4020        if( (*mul)[i] <= 0 )
4021        {
4022            return  semicListMulNegative;
4023        }
4024    }
4025
4026    // ----------------
4027    //  check symmetry
4028    // ----------------
4029
4030    int     j;
4031
4032    for( i=0, j=n-1; i<=j; i++,j-- )
4033    {
4034        if( (*num)[i] != rVar(currRing)*((*den)[i]) - (*num)[j] ||
4035            (*den)[i] != (*den)[j] ||
4036            (*mul)[i] != (*mul)[j] )
4037        {
4038            return  semicListNotSymmetric;
4039        }
4040    }
4041
4042    // ----------------
4043    //  check monotony
4044    // ----------------
4045
4046    for( i=0, j=1; i<n/2; i++,j++ )
4047    {
4048        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
4049        {
4050            return  semicListNotMonotonous;
4051        }
4052    }
4053
4054    // ---------------------
4055    //  check Milnor number
4056    // ---------------------
4057
4058    for( mu=0, i=0; i<n; i++ )
4059    {
4060        mu += (*mul)[i];
4061    }
4062
4063    if( mu != (int)(long)(l->m[0].Data( )) )
4064    {
4065        return  semicListMilnorWrong;
4066    }
4067
4068    // -------------------------
4069    //  check geometrical genus
4070    // -------------------------
4071
4072    for( pg=0, i=0; i<n; i++ )
4073    {
4074        if( (*num)[i]<=(*den)[i] )
4075        {
4076            pg += (*mul)[i];
4077        }
4078    }
4079
4080    if( pg != (int)(long)(l->m[1].Data( )) )
4081    {
4082        return  semicListPGWrong;
4083    }
4084
4085    return  semicOK;
4086}
4087
4088// ----------------------------------------------------------------------------
4089//  this procedure is called from the interpreter
4090// ----------------------------------------------------------------------------
4091//  first  = list of spectrum numbers
4092//  second = list of spectrum numbers
4093//  result = sum of the two lists
4094// ----------------------------------------------------------------------------
4095
4096BOOLEAN spaddProc( leftv result,leftv first,leftv second )
4097{
4098    semicState  state;
4099
4100    // -----------------
4101    //  check arguments
4102    // -----------------
4103
4104    lists l1 = (lists)first->Data( );
4105    lists l2 = (lists)second->Data( );
4106
4107    if( (state=list_is_spectrum( l1 )) != semicOK )
4108    {
4109        WerrorS( "first argument is not a spectrum:" );
4110        list_error( state );
4111    }
4112    else if( (state=list_is_spectrum( l2 )) != semicOK )
4113    {
4114        WerrorS( "second argument is not a spectrum:" );
4115        list_error( state );
4116    }
4117    else
4118    {
4119        spectrum s1= spectrumFromList ( l1 );
4120        spectrum s2= spectrumFromList ( l2 );
4121        spectrum sum( s1+s2 );
4122
4123        result->rtyp = LIST_CMD;
4124        result->data = (char*)(getList(sum));
4125    }
4126
4127    return  (state!=semicOK);
4128}
4129
4130// ----------------------------------------------------------------------------
4131//  this procedure is called from the interpreter
4132// ----------------------------------------------------------------------------
4133//  first  = list of spectrum numbers
4134//  second = integer
4135//  result = the multiple of the first list by the second factor
4136// ----------------------------------------------------------------------------
4137
4138BOOLEAN spmulProc( leftv result,leftv first,leftv second )
4139{
4140    semicState  state;
4141
4142    // -----------------
4143    //  check arguments
4144    // -----------------
4145
4146    lists   l = (lists)first->Data( );
4147    int     k = (int)(long)second->Data( );
4148
4149    if( (state=list_is_spectrum( l ))!=semicOK )
4150    {
4151        WerrorS( "first argument is not a spectrum" );
4152        list_error( state );
4153    }
4154    else if( k < 0 )
4155    {
4156        WerrorS( "second argument should be positive" );
4157        state = semicMulNegative;
4158    }
4159    else
4160    {
4161        spectrum s= spectrumFromList( l );
4162        spectrum product( k*s );
4163
4164        result->rtyp = LIST_CMD;
4165        result->data = (char*)getList(product);
4166    }
4167
4168    return  (state!=semicOK);
4169}
4170
4171// ----------------------------------------------------------------------------
4172//  this procedure is called from the interpreter
4173// ----------------------------------------------------------------------------
4174//  first  = list of spectrum numbers
4175//  second = list of spectrum numbers
4176//  result = semicontinuity index
4177// ----------------------------------------------------------------------------
4178
4179BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
4180{
4181  semicState  state;
4182  BOOLEAN qh=(((int)(long)w->Data())==1);
4183
4184  // -----------------
4185  //  check arguments
4186  // -----------------
4187
4188  lists l1 = (lists)u->Data( );
4189  lists l2 = (lists)v->Data( );
4190
4191  if( (state=list_is_spectrum( l1 ))!=semicOK )
4192  {
4193    WerrorS( "first argument is not a spectrum" );
4194    list_error( state );
4195  }
4196  else if( (state=list_is_spectrum( l2 ))!=semicOK )
4197  {
4198    WerrorS( "second argument is not a spectrum" );
4199    list_error( state );
4200  }
4201  else
4202  {
4203    spectrum s1= spectrumFromList( l1 );
4204    spectrum s2= spectrumFromList( l2 );
4205
4206    res->rtyp = INT_CMD;
4207    if (qh)
4208      res->data = (void*)(long)(s1.mult_spectrumh( s2 ));
4209    else
4210      res->data = (void*)(long)(s1.mult_spectrum( s2 ));
4211  }
4212
4213  // -----------------
4214  //  check status
4215  // -----------------
4216
4217  return  (state!=semicOK);
4218}
4219BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
4220{
4221  sleftv tmp;
4222  memset(&tmp,0,sizeof(tmp));
4223  tmp.rtyp=INT_CMD;
4224  /* tmp.data = (void *)0;  -- done by memset */
4225
4226  return  semicProc3(res,u,v,&tmp);
4227}
4228
4229#endif
4230
4231BOOLEAN loNewtonP( leftv res, leftv arg1 )
4232{
4233  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4234  return FALSE;
4235}
4236
4237BOOLEAN loSimplex( leftv res, leftv args )
4238{
4239  if ( !(rField_is_long_R(currRing)) )
4240  {
4241    WerrorS("Ground field not implemented!");
4242    return TRUE;
4243  }
4244
4245  simplex * LP;
4246  matrix m;
4247
4248  leftv v= args;
4249  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4250    return TRUE;
4251  else
4252    m= (matrix)(v->CopyD());
4253
4254  LP = new simplex(MATROWS(m),MATCOLS(m));
4255  LP->mapFromMatrix(m);
4256
4257  v= v->next;
4258  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
4259    return TRUE;
4260  else
4261    LP->m= (int)(long)(v->Data());
4262
4263  v= v->next;
4264  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
4265    return TRUE;
4266  else
4267    LP->n= (int)(long)(v->Data());
4268
4269  v= v->next;
4270  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
4271    return TRUE;
4272  else
4273    LP->m1= (int)(long)(v->Data());
4274
4275  v= v->next;
4276  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
4277    return TRUE;
4278  else
4279    LP->m2= (int)(long)(v->Data());
4280
4281  v= v->next;
4282  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
4283    return TRUE;
4284  else
4285    LP->m3= (int)(long)(v->Data());
4286
4287#ifdef mprDEBUG_PROT
4288  Print("m (constraints) %d\n",LP->m);
4289  Print("n (columns) %d\n",LP->n);
4290  Print("m1 (<=) %d\n",LP->m1);
4291  Print("m2 (>=) %d\n",LP->m2);
4292  Print("m3 (==) %d\n",LP->m3);
4293#endif
4294
4295  LP->compute();
4296
4297  lists lres= (lists)omAlloc( sizeof(slists) );
4298  lres->Init( 6 );
4299
4300  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4301  lres->m[0].data=(void*)LP->mapToMatrix(m);
4302
4303  lres->m[1].rtyp= INT_CMD;   // found a solution?
4304  lres->m[1].data=(void*)(long)LP->icase;
4305
4306  lres->m[2].rtyp= INTVEC_CMD;
4307  lres->m[2].data=(void*)LP->posvToIV();
4308
4309  lres->m[3].rtyp= INTVEC_CMD;
4310  lres->m[3].data=(void*)LP->zrovToIV();
4311
4312  lres->m[4].rtyp= INT_CMD;
4313  lres->m[4].data=(void*)(long)LP->m;
4314
4315  lres->m[5].rtyp= INT_CMD;
4316  lres->m[5].data=(void*)(long)LP->n;
4317
4318  res->data= (void*)lres;
4319
4320  return FALSE;
4321}
4322
4323BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4324{
4325  ideal gls = (ideal)(arg1->Data());
4326  int imtype= (int)(long)arg2->Data();
4327
4328  uResultant::resMatType mtype= determineMType( imtype );
4329
4330  // check input ideal ( = polynomial system )
4331  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4332  {
4333    return TRUE;
4334  }
4335
4336  uResultant *resMat= new uResultant( gls, mtype, false );
4337  if (resMat!=NULL)
4338  {
4339    res->rtyp = MODUL_CMD;
4340    res->data= (void*)resMat->accessResMat()->getMatrix();
4341    if (!errorreported) delete resMat;
4342  }
4343  return errorreported;
4344}
4345
4346BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4347{
4348
4349  poly gls;
4350  gls= (poly)(arg1->Data());
4351  int howclean= (int)(long)arg3->Data();
4352
4353  if ( !(rField_is_R(currRing) ||
4354         rField_is_Q(currRing) ||
4355         rField_is_long_R(currRing) ||
4356         rField_is_long_C(currRing)) )
4357  {
4358    WerrorS("Ground field not implemented!");
4359    return TRUE;
4360  }
4361
4362  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4363                          rField_is_long_C(currRing)) )
4364  {
4365    unsigned long int ii = (unsigned long int)arg2->Data();
4366    setGMPFloatDigits( ii, ii );
4367  }
4368
4369  if ( gls == NULL || pIsConstant( gls ) )
4370  {
4371    WerrorS("Input polynomial is constant!");
4372    return TRUE;
4373  }
4374
4375  int ldummy;
4376  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4377  //  int deg= pDeg( gls );
4378  //  int len= pLength( gls );
4379  int i,vpos=0;
4380  poly piter;
4381  lists elist;
4382  lists rlist;
4383
4384  elist= (lists)omAlloc( sizeof(slists) );
4385  elist->Init( 0 );
4386
4387  if ( rVar(currRing) > 1 )
4388  {
4389    piter= gls;
4390    for ( i= 1; i <= rVar(currRing); i++ )
4391      if ( pGetExp( piter, i ) )
4392      {
4393        vpos= i;
4394        break;
4395      }
4396    while ( piter )
4397    {
4398      for ( i= 1; i <= rVar(currRing); i++ )
4399        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4400        {
4401          WerrorS("The input polynomial must be univariate!");
4402          return TRUE;
4403        }
4404      pIter( piter );
4405    }
4406  }
4407
4408  rootContainer * roots= new rootContainer();
4409  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4410  piter= gls;
4411  for ( i= deg; i >= 0; i-- )
4412  {
4413    //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4414    if ( piter && pTotaldegree(piter) == i )
4415    {
4416      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4417      //nPrint( pcoeffs[i] );PrintS("  ");
4418      pIter( piter );
4419    }
4420    else
4421    {
4422      pcoeffs[i]= nInit(0);
4423    }
4424  }
4425
4426#ifdef mprDEBUG_PROT
4427  for (i=deg; i >= 0; i--)
4428  {
4429    nPrint( pcoeffs[i] );PrintS("  ");
4430  }
4431  PrintLn();
4432#endif
4433
4434  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4435  roots->solver( howclean );
4436
4437  int elem= roots->getAnzRoots();
4438  char *dummy;
4439  int j;
4440
4441  rlist= (lists)omAlloc( sizeof(slists) );
4442  rlist->Init( elem );
4443
4444  if (rField_is_long_C(currRing))
4445  {
4446    for ( j= 0; j < elem; j++ )
4447    {
4448      rlist->m[j].rtyp=NUMBER_CMD;
4449      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4450      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4451    }
4452  }
4453  else
4454  {
4455    for ( j= 0; j < elem; j++ )
4456    {
4457      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4458      rlist->m[j].rtyp=STRING_CMD;
4459      rlist->m[j].data=(void *)dummy;
4460    }
4461  }
4462
4463  elist->Clean();
4464  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4465
4466  // this is (via fillContainer) the same data as in root
4467  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4468  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4469
4470  delete roots;
4471
4472  res->rtyp= LIST_CMD;
4473  res->data= (void*)rlist;
4474
4475  return FALSE;
4476}
4477
4478BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4479{
4480  int i;
4481  ideal p,w;
4482  p= (ideal)arg1->Data();
4483  w= (ideal)arg2->Data();
4484
4485  // w[0] = f(p^0)
4486  // w[1] = f(p^1)
4487  // ...
4488  // p can be a vector of numbers (multivariate polynom)
4489  //   or one number (univariate polynom)
4490  // tdg = deg(f)
4491
4492  int n= IDELEMS( p );
4493  int m= IDELEMS( w );
4494  int tdg= (int)(long)arg3->Data();
4495
4496  res->data= (void*)NULL;
4497
4498  // check the input
4499  if ( tdg < 1 )
4500  {
4501    WerrorS("Last input parameter must be > 0!");
4502    return TRUE;
4503  }
4504  if ( n != rVar(currRing) )
4505  {
4506    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4507    return TRUE;
4508  }
4509  if ( m != (int)pow((double)tdg+1,(double)n) )
4510  {
4511    Werror("Size of second input ideal must be equal to %d!",
4512      (int)pow((double)tdg+1,(double)n));
4513    return TRUE;
4514  }
4515  if ( !(rField_is_Q(currRing) /* ||
4516         rField_is_R() || rField_is_long_R() ||
4517         rField_is_long_C()*/ ) )
4518         {
4519    WerrorS("Ground field not implemented!");
4520    return TRUE;
4521  }
4522
4523  number tmp;
4524  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4525  for ( i= 0; i < n; i++ )
4526  {
4527    pevpoint[i]=nInit(0);
4528    if (  (p->m)[i] )
4529    {
4530      tmp = pGetCoeff( (p->m)[i] );
4531      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4532      {
4533        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4534        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4535        return TRUE;
4536      }
4537    } else tmp= NULL;
4538    if ( !nIsZero(tmp) )
4539    {
4540      if ( !pIsConstant((p->m)[i]))
4541      {
4542        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4543        WerrorS("Elements of first input ideal must be numbers!");
4544        return TRUE;
4545      }
4546      pevpoint[i]= nCopy( tmp );
4547    }
4548  }
4549
4550  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4551  for ( i= 0; i < m; i++ )
4552  {
4553    wresults[i]= nInit(0);
4554    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4555    {
4556      if ( !pIsConstant((w->m)[i]))
4557      {
4558        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4559        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4560        WerrorS("Elements of second input ideal must be numbers!");
4561        return TRUE;
4562      }
4563      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4564    }
4565  }
4566
4567  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4568  number *ncpoly= vm.interpolateDense( wresults );
4569  // do not free ncpoly[]!!
4570  poly rpoly= vm.numvec2poly( ncpoly );
4571
4572  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4573  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4574
4575  res->data= (void*)rpoly;
4576  return FALSE;
4577}
4578
4579BOOLEAN nuUResSolve( leftv res, leftv args )
4580{
4581  leftv v= args;
4582
4583  ideal gls;
4584  int imtype;
4585  int howclean;
4586
4587  // get ideal
4588  if ( v->Typ() != IDEAL_CMD )
4589    return TRUE;
4590  else gls= (ideal)(v->Data());
4591  v= v->next;
4592
4593  // get resultant matrix type to use (0,1)
4594  if ( v->Typ() != INT_CMD )
4595    return TRUE;
4596  else imtype= (int)(long)v->Data();
4597  v= v->next;
4598
4599  if (imtype==0)
4600  {
4601    ideal test_id=idInit(1,1);
4602    int j;
4603    for(j=IDELEMS(gls)-1;j>=0;j--)
4604    {
4605      if (gls->m[j]!=NULL)
4606      {
4607        test_id->m[0]=gls->m[j];
4608        intvec *dummy_w=id_QHomWeight(test_id, currRing);
4609        if (dummy_w!=NULL)
4610        {
4611          WerrorS("Newton polytope not of expected dimension");
4612          delete dummy_w;
4613          return TRUE;
4614        }
4615      }
4616    }
4617  }
4618
4619  // get and set precision in digits ( > 0 )
4620  if ( v->Typ() != INT_CMD )
4621    return TRUE;
4622  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4623                          rField_is_long_C(currRing)) )
4624  {
4625    unsigned long int ii=(unsigned long int)v->Data();
4626    setGMPFloatDigits( ii, ii );
4627  }
4628  v= v->next;
4629
4630  // get interpolation steps (0,1,2)
4631  if ( v->Typ() != INT_CMD )
4632    return TRUE;
4633  else howclean= (int)(long)v->Data();
4634
4635  uResultant::resMatType mtype= determineMType( imtype );
4636  int i,count;
4637  lists listofroots= NULL;
4638  number smv= NULL;
4639  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4640
4641  //emptylist= (lists)omAlloc( sizeof(slists) );
4642  //emptylist->Init( 0 );
4643
4644  //res->rtyp = LIST_CMD;
4645  //res->data= (void *)emptylist;
4646
4647  // check input ideal ( = polynomial system )
4648  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4649  {
4650    return TRUE;
4651  }
4652
4653  uResultant * ures;
4654  rootContainer ** iproots;
4655  rootContainer ** muiproots;
4656  rootArranger * arranger;
4657
4658  // main task 1: setup of resultant matrix
4659  ures= new uResultant( gls, mtype );
4660  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4661  {
4662    WerrorS("Error occurred during matrix setup!");
4663    return TRUE;
4664  }
4665
4666  // if dense resultant, check if minor nonsingular
4667  if ( mtype == uResultant::denseResMat )
4668  {
4669    smv= ures->accessResMat()->getSubDet();
4670#ifdef mprDEBUG_PROT
4671    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4672#endif
4673    if ( nIsZero(smv) )
4674    {
4675      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4676      return TRUE;
4677    }
4678  }
4679
4680  // main task 2: Interpolate specialized resultant polynomials
4681  if ( interpolate_det )
4682    iproots= ures->interpolateDenseSP( false, smv );
4683  else
4684    iproots= ures->specializeInU( false, smv );
4685
4686  // main task 3: Interpolate specialized resultant polynomials
4687  if ( interpolate_det )
4688    muiproots= ures->interpolateDenseSP( true, smv );
4689  else
4690    muiproots= ures->specializeInU( true, smv );
4691
4692#ifdef mprDEBUG_PROT
4693  int c= iproots[0]->getAnzElems();
4694  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4695  c= muiproots[0]->getAnzElems();
4696  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4697#endif
4698
4699  // main task 4: Compute roots of specialized polys and match them up
4700  arranger= new rootArranger( iproots, muiproots, howclean );
4701  arranger->solve_all();
4702
4703  // get list of roots
4704  if ( arranger->success() )
4705  {
4706    arranger->arrange();
4707    listofroots= listOfRoots(arranger, gmp_output_digits );
4708  }
4709  else
4710  {
4711    WerrorS("Solver was unable to find any roots!");
4712    return TRUE;
4713  }
4714
4715  // free everything
4716  count= iproots[0]->getAnzElems();
4717  for (i=0; i < count; i++) delete iproots[i];
4718  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4719  count= muiproots[0]->getAnzElems();
4720  for (i=0; i < count; i++) delete muiproots[i];
4721  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4722
4723  delete ures;
4724  delete arranger;
4725  nDelete( &smv );
4726
4727  res->data= (void *)listofroots;
4728
4729  //emptylist->Clean();
4730  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4731
4732  return FALSE;
4733}
4734
4735// from mpr_numeric.cc
4736lists listOfRoots( rootArranger* self, const unsigned int oprec )
4737{
4738  int i,j;
4739  int count= self->roots[0]->getAnzRoots(); // number of roots
4740  int elem= self->roots[0]->getAnzElems();  // number of koordinates per root
4741
4742  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
4743
4744  if ( self->found_roots )
4745  {
4746    listofroots->Init( count );
4747
4748    for (i=0; i < count; i++)
4749    {
4750      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
4751      onepoint->Init(elem);
4752      for ( j= 0; j < elem; j++ )
4753      {
4754        if ( !rField_is_long_C(currRing) )
4755        {
4756          onepoint->m[j].rtyp=STRING_CMD;
4757          onepoint->m[j].data=(void *)complexToStr((*self->roots[j])[i],oprec, currRing->cf);
4758        }
4759        else
4760        {
4761          onepoint->m[j].rtyp=NUMBER_CMD;
4762          onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
4763        }
4764        onepoint->m[j].next= NULL;
4765        onepoint->m[j].name= NULL;
4766      }
4767      listofroots->m[i].rtyp=LIST_CMD;
4768      listofroots->m[i].data=(void *)onepoint;
4769      listofroots->m[j].next= NULL;
4770      listofroots->m[j].name= NULL;
4771    }
4772
4773  }
4774  else
4775  {
4776    listofroots->Init( 0 );
4777  }
4778
4779  return listofroots;
4780}
4781
4782// from ring.cc
4783void rSetHdl(idhdl h)
4784{
4785  ring rg = NULL;
4786  if (h!=NULL)
4787  {
4788//   Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
4789    rg = IDRING(h);
4790    if (rg==NULL) return; //id <>NULL, ring==NULL
4791    omCheckAddrSize((ADDRESS)h,sizeof(idrec));
4792    if (IDID(h))  // OB: ????
4793      omCheckAddr((ADDRESS)IDID(h));
4794    rTest(rg);
4795  }
4796
4797  // clean up history
4798  if (sLastPrinted.RingDependend())
4799  {
4800    sLastPrinted.CleanUp();
4801    memset(&sLastPrinted,0,sizeof(sleftv));
4802  }
4803
4804  // test for valid "currRing":
4805  if ((rg!=NULL) && (rg->idroot==NULL))
4806  {
4807    ring old=rg;
4808    rg=rAssure_HasComp(rg);
4809    if (old!=rg)
4810    {
4811      rKill(old);
4812      IDRING(h)=rg;
4813    }
4814  }
4815   /*------------ change the global ring -----------------------*/
4816  rChangeCurrRing(rg);
4817  currRingHdl = h;
4818}
4819
4820BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
4821{
4822  int last = 0, o=0, n = 1, i=0, typ = 1, j;
4823  sleftv *sl = ord;
4824
4825  // determine nBlocks
4826  while (sl!=NULL)
4827  {
4828    intvec *iv = (intvec *)(sl->data);
4829    if (((*iv)[1]==ringorder_c)||((*iv)[1]==ringorder_C))
4830      i++;
4831    else if ((*iv)[1]==ringorder_L)
4832    {
4833      R->bitmask=(*iv)[2];
4834      n--;
4835    }
4836    else if (((*iv)[1]!=ringorder_a)
4837    && ((*iv)[1]!=ringorder_a64)
4838    && ((*iv)[1]!=ringorder_am))
4839      o++;
4840    n++;
4841    sl=sl->next;
4842  }
4843  // check whether at least one real ordering
4844  if (o==0)
4845  {
4846    WerrorS("invalid combination of orderings");
4847    return TRUE;
4848  }
4849  // if no c/C ordering is given, increment n
4850  if (i==0) n++;
4851  else if (i != 1)
4852  {
4853    // throw error if more than one is given
4854    WerrorS("more than one ordering c/C specified");
4855    return TRUE;
4856  }
4857
4858  // initialize fields of R
4859  R->order=(int *)omAlloc0(n*sizeof(int));
4860  R->block0=(int *)omAlloc0(n*sizeof(int));
4861  R->block1=(int *)omAlloc0(n*sizeof(int));
4862  R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
4863
4864  int *weights=(int*)omAlloc0((R->N+1)*sizeof(int));
4865
4866  // init order, so that rBlocks works correctly
4867  for (j=0; j < n-1; j++)
4868    R->order[j] = (int) ringorder_unspec;
4869  // set last _C order, if no c/C order was given
4870  if (i == 0) R->order[n-2] = ringorder_C;
4871
4872  /* init orders */
4873  sl=ord;
4874  n=-1;
4875  while (sl!=NULL)
4876  {
4877    intvec *iv;
4878    iv = (intvec *)(sl->data);
4879    if ((*iv)[1]!=ringorder_L)
4880    {
4881      n++;
4882
4883      /* the format of an ordering:
4884       *  iv[0]: factor
4885       *  iv[1]: ordering
4886       *  iv[2..end]: weights
4887       */
4888      R->order[n] = (*iv)[1];
4889      typ=1;
4890      switch ((*iv)[1])
4891      {
4892          case ringorder_ws:
4893          case ringorder_Ws:
4894            typ=-1;
4895          case ringorder_wp:
4896          case ringorder_Wp:
4897            R->wvhdl[n]=(int*)omAlloc((iv->length()-1)*sizeof(int));
4898            R->block0[n] = last+1;
4899            for (i=2; i<iv->length(); i++)
4900            {
4901              R->wvhdl[n][i-2] = (*iv)[i];
4902              last++;
4903              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4904            }
4905            R->block1[n] = last;
4906            break;
4907          case ringorder_ls:
4908          case ringorder_ds:
4909          case ringorder_Ds:
4910          case ringorder_rs:
4911            typ=-1;
4912          case ringorder_lp:
4913          case ringorder_dp:
4914          case ringorder_Dp:
4915          case ringorder_rp:
4916            R->block0[n] = last+1;
4917            if (iv->length() == 3) last+=(*iv)[2];
4918            else last += (*iv)[0];
4919            R->block1[n] = last;
4920            //if ((R->block0[n]>R->block1[n])
4921            //|| (R->block1[n]>rVar(R)))
4922            //{
4923            //  R->block1[n]=rVar(R);
4924            //  //WerrorS("ordering larger than number of variables");
4925            //  break;
4926            //}
4927            if (rCheckIV(iv)) return TRUE;
4928            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4929            {
4930              if (weights[i]==0) weights[i]=typ;
4931            }
4932            break;
4933
4934          case ringorder_s: // no 'rank' params!
4935          {
4936
4937            if(iv->length() > 3)
4938              return TRUE;
4939
4940            if(iv->length() == 3)
4941            {
4942              const int s = (*iv)[2];
4943              R->block0[n] = s;
4944              R->block1[n] = s;
4945            }
4946            break;
4947          }
4948          case ringorder_IS:
4949          {
4950            if(iv->length() != 3) return TRUE;
4951
4952            const int s = (*iv)[2];
4953
4954            if( 1 < s || s < -1 ) return TRUE;
4955
4956            R->block0[n] = s;
4957            R->block1[n] = s;
4958            break;
4959          }
4960          case ringorder_S:
4961          case ringorder_c:
4962          case ringorder_C:
4963          {
4964            if (rCheckIV(iv)) return TRUE;
4965            break;
4966          }
4967          case ringorder_aa:
4968          case ringorder_a:
4969          {
4970            R->block0[n] = last+1;
4971            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4972            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int));
4973            for (i=2; i<iv->length(); i++)
4974            {
4975              R->wvhdl[n][i-2]=(*iv)[i];
4976              last++;
4977              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4978            }
4979            last=R->block0[n]-1;
4980            break;
4981          }
4982          case ringorder_am:
4983          {
4984            R->block0[n] = last+1;
4985            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4986            R->wvhdl[n] = (int*)omAlloc(iv->length()*sizeof(int));
4987            if (R->block1[n]- R->block0[n]+2>=iv->length())
4988               WarnS("missing module weights");
4989            for (i=2; i<=(R->block1[n]-R->block0[n]+2); i++)
4990            {
4991              R->wvhdl[n][i-2]=(*iv)[i];
4992              last++;
4993              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4994            }
4995            R->wvhdl[n][i-2]=iv->length() -3 -(R->block1[n]- R->block0[n]);
4996            for (; i<iv->length(); i++)
4997            {
4998              R->wvhdl[n][i-1]=(*iv)[i];
4999            }
5000            last=R->block0[n]-1;
5001            break;
5002          }
5003          case ringorder_a64:
5004          {
5005            R->block0[n] = last+1;
5006            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
5007            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int64));
5008            int64 *w=(int64 *)R->wvhdl[n];
5009            for (i=2; i<iv->length(); i++)
5010            {
5011              w[i-2]=(*iv)[i];
5012              last++;
5013              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5014            }
5015            last=R->block0[n]-1;
5016            break;
5017          }
5018          case ringorder_M:
5019          {
5020            int Mtyp=rTypeOfMatrixOrder(iv);
5021            if (Mtyp==0) return TRUE;
5022            if (Mtyp==-1) typ = -1;
5023
5024            R->wvhdl[n] =( int *)omAlloc((iv->length()-1)*sizeof(int));
5025            for (i=2; i<iv->length();i++)
5026              R->wvhdl[n][i-2]=(*iv)[i];
5027
5028            R->block0[n] = last+1;
5029            last += (int)sqrt((double)(iv->length()-2));
5030            R->block1[n] = last;
5031            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
5032            {
5033              if (weights[i]==0) weights[i]=typ;
5034            }
5035            break;
5036          }
5037
5038          case ringorder_no:
5039            R->order[n] = ringorder_unspec;
5040            return TRUE;
5041
5042          default:
5043            Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
5044            R->order[n] = ringorder_unspec;
5045            return TRUE;
5046      }
5047    }
5048    sl=sl->next;
5049  }
5050
5051  // check for complete coverage
5052  while ( n >= 0 && (
5053          (R->order[n]==ringorder_c)
5054      ||  (R->order[n]==ringorder_C)
5055      ||  (R->order[n]==ringorder_s)
5056      ||  (R->order[n]==ringorder_S)
5057      ||  (R->order[n]==ringorder_IS)
5058                    )) n--;
5059
5060  assume( n >= 0 );
5061
5062  if (R->block1[n] != R->N)
5063  {
5064    if (((R->order[n]==ringorder_dp) ||
5065         (R->order[n]==ringorder_ds) ||
5066         (R->order[n]==ringorder_Dp) ||
5067         (R->order[n]==ringorder_Ds) ||
5068         (R->order[n]==ringorder_rp) ||
5069         (R->order[n]==ringorder_rs) ||
5070         (R->order[n]==ringorder_lp) ||
5071         (R->order[n]==ringorder_ls))
5072        &&
5073        R->block0[n] <= R->N)
5074    {
5075      R->block1[n] = R->N;
5076    }
5077    else
5078    {
5079      Werror("mismatch of number of vars (%d) and ordering (%d vars)",
5080             R->N,R->block1[n]);
5081      return TRUE;
5082    }
5083  }
5084  // find OrdSgn:
5085  R->OrdSgn = 1;
5086  for(i=1;i<=R->N;i++)
5087  { if (weights[i]<0) { R->OrdSgn=-1;break; }}
5088  omFree(weights);
5089  return FALSE;
5090}
5091
5092BOOLEAN rSleftvList2StringArray(sleftv* sl, char** p)
5093{
5094
5095  while(sl!=NULL)
5096  {
5097    if (sl->Name() == sNoName)
5098    {
5099      if (sl->Typ()==POLY_CMD)
5100      {
5101        sleftv s_sl;
5102        iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
5103        if (s_sl.Name() != sNoName)
5104          *p = omStrDup(s_sl.Name());
5105        else
5106          *p = NULL;
5107        sl->next = s_sl.next;
5108        s_sl.next = NULL;
5109        s_sl.CleanUp();
5110        if (*p == NULL) return TRUE;
5111      }
5112      else
5113        return TRUE;
5114    }
5115    else
5116      *p = omStrDup(sl->Name());
5117    p++;
5118    sl=sl->next;
5119  }
5120  return FALSE;
5121}
5122
5123const short MAX_SHORT = 32767; // (1 << (sizeof(short)*8)) - 1;
5124
5125////////////////////
5126//
5127// rInit itself:
5128//
5129// INPUT:  s: name, pn: ch & parameter (names), rv: variable (names)
5130//         ord: ordering
5131// RETURN: currRingHdl on success
5132//         NULL        on error
5133// NOTE:   * makes new ring to current ring, on success
5134//         * considers input sleftv's as read-only
5135//idhdl rInit(char *s, sleftv* pn, sleftv* rv, sleftv* ord)
5136ring rInit(sleftv* pn, sleftv* rv, sleftv* ord)
5137{
5138#ifdef HAVE_RINGS
5139  //unsigned int ringtype = 0;
5140  mpz_ptr modBase = NULL;
5141  unsigned int modExponent = 1;
5142#endif
5143  int float_len=0;
5144  int float_len2=0;
5145  ring R = NULL;
5146  //BOOLEAN ffChar=FALSE;
5147
5148  /* ch -------------------------------------------------------*/
5149  // get ch of ground field
5150
5151  // allocated ring
5152  R = (ring) omAlloc0Bin(sip_sring_bin);
5153
5154  coeffs cf = NULL;
5155
5156  assume( pn != NULL );
5157  const int P = pn->listLength();
5158
5159  if ((pn->Typ()==CRING_CMD)&&(P==1))
5160  {
5161    cf=(coeffs)pn->CopyD();
5162    assume( cf != NULL );
5163  }
5164  else if (pn->Typ()==INT_CMD)
5165  {
5166    int ch = (int)(long)pn->Data();
5167
5168    /* parameter? -------------------------------------------------------*/
5169    pn = pn->next;
5170
5171    if (pn == NULL) // no params!?
5172    {
5173      if (ch!=0)
5174      {
5175        int ch2=IsPrime(ch);
5176        if ((ch<2)||(ch!=ch2))
5177        {
5178          Warn("%d is invalid as characteristic of the ground field. 32003 is used.", ch);
5179          ch=32003;
5180        }
5181        cf = nInitChar(n_Zp, (void*)(long)ch);
5182      }
5183      else
5184        cf = nInitChar(n_Q, (void*)(long)ch);
5185    }
5186    else
5187    {
5188      const int pars = pn->listLength();
5189
5190      assume( pars > 0 );
5191
5192      // predefined finite field: (p^k, a)
5193      if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
5194      {
5195        GFInfo param;
5196
5197        param.GFChar = ch;
5198        param.GFDegree = 1;
5199        param.GFPar_name = pn->name;
5200
5201        cf = nInitChar(n_GF, &param);
5202      }
5203      else // (0/p, a, b, ..., z)
5204      {
5205        if ((ch!=0) && (ch!=IsPrime(ch)))
5206        {
5207          WerrorS("too many parameters");
5208          goto rInitError;
5209        }
5210
5211        char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5212
5213        if (rSleftvList2StringArray(pn, names))
5214        {
5215          WerrorS("parameter expected");
5216          goto rInitError;
5217        }
5218
5219        TransExtInfo extParam;
5220
5221        extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5222        for(int i=pars-1; i>=0;i--)
5223        {
5224          omFree(names[i]);
5225        }
5226        omFree(names);
5227
5228        cf = nInitChar(n_transExt, &extParam);
5229      }
5230    }
5231
5232//    if (cf==NULL) goto rInitError;
5233    assume( cf != NULL );
5234  }
5235  else if ((pn->name != NULL)
5236  && ((strcmp(pn->name,"real")==0) || (strcmp(pn->name,"complex")==0)))
5237  {
5238    BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5239    if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5240    {
5241      float_len=(int)(long)pn->next->Data();
5242      float_len2=float_len;
5243      pn=pn->next;
5244      if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5245      {
5246        float_len2=(int)(long)pn->next->Data();
5247        pn=pn->next;
5248      }
5249    }
5250
5251    if (!complex_flag)
5252      complex_flag= pn->next != NULL;
5253    if( !complex_flag && (float_len2 <= (short)SHORT_REAL_LENGTH))
5254       cf=nInitChar(n_R, NULL);
5255    else // longR or longC?
5256    {
5257       LongComplexInfo param;
5258
5259       param.float_len = si_min (float_len, 32767);
5260       param.float_len2 = si_min (float_len2, 32767);
5261
5262       // set the parameter name
5263       if (complex_flag)
5264       {
5265         if (param.float_len < SHORT_REAL_LENGTH)
5266         {
5267           param.float_len= SHORT_REAL_LENGTH;
5268           param.float_len2= SHORT_REAL_LENGTH;
5269         }
5270         if (pn->next == NULL)
5271           param.par_name=(const char*)"i"; //default to i
5272         else
5273           param.par_name = (const char*)pn->next->name;
5274       }
5275
5276       cf = nInitChar(complex_flag ? n_long_C: n_long_R, (void*)&param);
5277    }
5278    assume( cf != NULL );
5279  }
5280#ifdef HAVE_RINGS
5281  else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5282  {
5283    // TODO: change to use coeffs_BIGINT!?
5284    modBase = (mpz_ptr) omAlloc(sizeof(mpz_t));
5285    mpz_init_set_si(modBase, 0);
5286    if (pn->next!=NULL)
5287    {
5288      if (pn->next->Typ()==INT_CMD)
5289      {
5290        mpz_set_ui(modBase, (int)(long) pn->next->Data());
5291        pn=pn->next;
5292        if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5293        {
5294          modExponent = (long) pn->next->Data();
5295          pn=pn->next;
5296        }
5297        while ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5298        {
5299          mpz_mul_ui(modBase, modBase, (int)(long) pn->next->Data());
5300          pn=pn->next;
5301        }
5302      }
5303      else if (pn->next->Typ()==BIGINT_CMD)
5304      {
5305        number p=(number)pn->next->CopyD(); // FIXME: why CopyD() here if nlGMP should not overtake p!?
5306        nlGMP(p,(number)modBase,coeffs_BIGINT); // TODO? // extern void   nlGMP(number &i, number n, const coeffs r); // FIXME: n_MPZ( modBase, p, coeffs_BIGINT); ?
5307        n_Delete(&p,coeffs_BIGINT);
5308      }
5309    }
5310    else
5311      cf=nInitChar(n_Z,NULL);
5312
5313    if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_cmp_ui(modBase, 0) < 0))
5314    {
5315      Werror("Wrong ground ring specification (module is 1)");
5316      goto rInitError;
5317    }
5318    if (modExponent < 1)
5319    {
5320      Werror("Wrong ground ring specification (exponent smaller than 1");
5321      goto rInitError;
5322    }
5323    // module is 0 ---> integers ringtype = 4;
5324    // we have an exponent
5325    if (modExponent > 1 && cf == NULL)
5326    {
5327      if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
5328      {
5329        /* this branch should be active for modExponent = 2..32 resp. 2..64,
5330           depending on the size of a long on the respective platform */
5331        //ringtype = 1;       // Use Z/2^ch
5332        cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5333        mpz_clear(modBase);
5334        omFreeSize (modBase, sizeof (mpz_t));
5335      }
5336      else
5337      {
5338        if (mpz_cmp_ui(modBase,0)==0)
5339        {
5340          WerrorS("modulus must not be 0 or parameter not allowed");
5341          goto rInitError;
5342        }
5343        //ringtype = 3;
5344        ZnmInfo info;
5345        info.base= modBase;
5346        info.exp= modExponent;
5347        cf=nInitChar(n_Znm,(void*) &info); //exponent is missing
5348      }
5349    }
5350    // just a module m > 1
5351    else if (cf == NULL)
5352    {
5353      if (mpz_cmp_ui(modBase,0)==0)
5354      {
5355        WerrorS("modulus must not be 0 or parameter not allowed");
5356        goto rInitError;
5357      }
5358      //ringtype = 2;
5359      ZnmInfo info;
5360      info.base= modBase;
5361      info.exp= modExponent;
5362      cf=nInitChar(n_Zn,(void*) &info);
5363    }
5364    assume( cf != NULL );
5365  }
5366#endif
5367  // ring NEW = OLD, (), (); where OLD is a polynomial ring...
5368  else if ((pn->Typ()==RING_CMD) && (P == 1))
5369  {
5370    TransExtInfo extParam;
5371    extParam.r = (ring)pn->Data();
5372    cf = nInitChar(n_transExt, &extParam);
5373  }
5374  else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5375  {
5376    AlgExtInfo extParam;
5377    extParam.r = (ring)pn->Data();
5378
5379    cf = nInitChar(n_algExt, &extParam);   // Q[a]/<minideal>
5380  }
5381  else
5382  {
5383    Werror("Wrong or unknown ground field specification");
5384#ifndef SING_NDEBUG
5385    sleftv* p = pn;
5386    while (p != NULL)
5387    {
5388      Print( "pn[%p]: type: %d [%s]: %p, name: %s", (void*)p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5389      PrintLn();
5390      p = p->next;
5391    }
5392#endif
5393    goto rInitError;
5394  }
5395//  pn=pn->next;
5396
5397  /*every entry in the new ring is initialized to 0*/
5398
5399  /* characteristic -----------------------------------------------*/
5400  /* input: 0 ch=0 : Q     parameter=NULL    ffChar=FALSE   float_len
5401   *         0    1 : Q(a,...)        *names         FALSE
5402   *         0   -1 : R               NULL           FALSE  0
5403   *         0   -1 : R               NULL           FALSE  prec. >6
5404   *         0   -1 : C               *names         FALSE  prec. 0..?
5405   *         p    p : Fp              NULL           FALSE
5406   *         p   -p : Fp(a)           *names         FALSE
5407   *         q    q : GF(q=p^n)       *names         TRUE
5408  */
5409  if (cf==NULL)
5410  {
5411    Werror("Invalid ground field specification");
5412    goto rInitError;
5413//    const int ch=32003;
5414//    cf=nInitChar(n_Zp, (void*)(long)ch);
5415  }
5416
5417  assume( R != NULL );
5418
5419  R->cf = cf;
5420
5421  /* names and number of variables-------------------------------------*/
5422  {
5423    int l=rv->listLength();
5424
5425    if (l>MAX_SHORT)
5426    {
5427      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5428       goto rInitError;
5429    }
5430    R->N = l; /*rv->listLength();*/
5431  }
5432  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5433  if (rSleftvList2StringArray(rv, R->names))
5434  {
5435    WerrorS("name of ring variable expected");
5436    goto rInitError;
5437  }
5438
5439  /* check names and parameters for conflicts ------------------------- */
5440  rRenameVars(R); // conflicting variables will be renamed
5441  /* ordering -------------------------------------------------------------*/
5442  if (rSleftvOrdering2Ordering(ord, R))
5443    goto rInitError;
5444
5445  // Complete the initialization
5446  if (rComplete(R,1))
5447    goto rInitError;
5448
5449/*#ifdef HAVE_RINGS
5450// currently, coefficients which are ring elements require a global ordering:
5451  if (rField_is_Ring(R) && (R->OrdSgn==-1))
5452  {
5453    WerrorS("global ordering required for these coefficients");
5454    goto rInitError;
5455  }
5456#endif*/
5457
5458  rTest(R);
5459
5460  // try to enter the ring into the name list
5461  // need to clean up sleftv here, before this ring can be set to
5462  // new currRing or currRing can be killed beacuse new ring has
5463  // same name
5464  if (pn != NULL) pn->CleanUp();
5465  if (rv != NULL) rv->CleanUp();
5466  if (ord != NULL) ord->CleanUp();
5467  //if ((tmp = enterid(s, myynest, RING_CMD, &IDROOT))==NULL)
5468  //  goto rInitError;
5469
5470  //memcpy(IDRING(tmp),R,sizeof(*R));
5471  // set current ring
5472  //omFreeBin(R,  ip_sring_bin);
5473  //return tmp;
5474  return R;
5475
5476  // error case:
5477  rInitError:
5478  if  ((R != NULL)&&(R->cf!=NULL)) rDelete(R);
5479  if (pn != NULL) pn->CleanUp();
5480  if (rv != NULL) rv->CleanUp();
5481  if (ord != NULL) ord->CleanUp();
5482  return NULL;
5483}
5484
5485ring rSubring(ring org_ring, sleftv* rv)
5486{
5487  ring R = rCopy0(org_ring);
5488  int *perm=(int *)omAlloc0((org_ring->N+1)*sizeof(int));
5489  int n = rBlocks(org_ring), i=0, j;
5490
5491  /* names and number of variables-------------------------------------*/
5492  {
5493    int l=rv->listLength();
5494    if (l>MAX_SHORT)
5495    {
5496      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5497       goto rInitError;
5498    }
5499    R->N = l; /*rv->listLength();*/
5500  }
5501  omFree(R->names);
5502  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5503  if (rSleftvList2StringArray(rv, R->names))
5504  {
5505    WerrorS("name of ring variable expected");
5506    goto rInitError;
5507  }
5508
5509  /* check names for subring in org_ring ------------------------- */
5510  {
5511    i=0;
5512
5513    for(j=0;j<R->N;j++)
5514    {
5515      for(;i<org_ring->N;i++)
5516      {
5517        if (strcmp(org_ring->names[i],R->names[j])==0)
5518        {
5519          perm[i+1]=j+1;
5520          break;
5521        }
5522      }
5523      if (i>org_ring->N)
5524      {
5525        Werror("variable %d (%s) not in basering",j+1,R->names[j]);
5526        break;
5527      }
5528    }
5529  }
5530  //Print("perm=");
5531  //for(i=1;i<org_ring->N;i++) Print("v%d -> v%d\n",i,perm[i]);
5532  /* ordering -------------------------------------------------------------*/
5533
5534  for(i=0;i<n;i++)
5535  {
5536    int min_var=-1;
5537    int max_var=-1;
5538    for(j=R->block0[i];j<=R->block1[i];j++)
5539    {
5540      if (perm[j]>0)
5541      {
5542        if (min_var==-1) min_var=perm[j];
5543        max_var=perm[j];
5544      }
5545    }
5546    if (min_var!=-1)
5547    {
5548      //Print("block %d: old %d..%d, now:%d..%d\n",
5549      //      i,R->block0[i],R->block1[i],min_var,max_var);
5550      R->block0[i]=min_var;
5551      R->block1[i]=max_var;
5552      if (R->wvhdl[i]!=NULL)
5553      {
5554        omFree(R->wvhdl[i]);
5555        R->wvhdl[i]=(int*)omAlloc0((max_var-min_var+1)*sizeof(int));
5556        for(j=org_ring->block0[i];j<=org_ring->block1[i];j++)
5557        {
5558          if (perm[j]>0)
5559          {
5560            R->wvhdl[i][perm[j]-R->block0[i]]=
5561                org_ring->wvhdl[i][j-org_ring->block0[i]];
5562            //Print("w%d=%d (orig_w%d)\n",perm[j],R->wvhdl[i][perm[j]-R->block0[i]],j);
5563          }
5564        }
5565      }
5566    }
5567    else
5568    {
5569      if(R->block0[i]>0)
5570      {
5571        //Print("skip block %d\n",i);
5572        R->order[i]=ringorder_unspec;
5573        if (R->wvhdl[i] !=NULL) omFree(R->wvhdl[i]);
5574        R->wvhdl[i]=NULL;
5575      }
5576      //else Print("keep block %d\n",i);
5577    }
5578  }
5579  i=n-1;
5580  while(i>0)
5581  {
5582    // removed unneded blocks
5583    if(R->order[i-1]==ringorder_unspec)
5584    {
5585      for(j=i;j<=n;j++)
5586      {
5587        R->order[j-1]=R->order[j];
5588        R->block0[j-1]=R->block0[j];
5589        R->block1[j-1]=R->block1[j];
5590        if (R->wvhdl[j-1] !=NULL) omFree(R->wvhdl[j-1]);
5591        R->wvhdl[j-1]=R->wvhdl[j];
5592      }
5593      R->order[n]=ringorder_unspec;
5594      n--;
5595    }
5596    i--;
5597  }
5598  n=rBlocks(org_ring)-1;
5599  while (R->order[n]==0)  n--;
5600  while (R->order[n]==ringorder_unspec)  n--;
5601  if ((R->order[n]==ringorder_c) ||  (R->order[n]==ringorder_C)) n--;
5602  if (R->block1[n] != R->N)
5603  {
5604    if (((R->order[n]==ringorder_dp) ||
5605         (R->order[n]==ringorder_ds) ||
5606         (R->order[n]==ringorder_Dp) ||
5607         (R->order[n]==ringorder_Ds) ||
5608         (R->order[n]==ringorder_rp) ||
5609         (R->order[n]==ringorder_rs) ||
5610         (R->order[n]==ringorder_lp) ||
5611         (R->order[n]==ringorder_ls))
5612        &&
5613        R->block0[n] <= R->N)
5614    {
5615      R->block1[n] = R->N;
5616    }
5617    else
5618    {
5619      Werror("mismatch of number of vars (%d) and ordering (%d vars) in block %d",
5620             R->N,R->block1[n],n);
5621      return NULL;
5622    }
5623  }
5624  omFree(perm);
5625  // find OrdSgn:
5626  R->OrdSgn = org_ring->OrdSgn; // IMPROVE!
5627  //for(i=1;i<=R->N;i++)
5628  //{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
5629  //omFree(weights);
5630  // Complete the initialization
5631  if (rComplete(R,1))
5632    goto rInitError;
5633
5634  rTest(R);
5635
5636  if (rv != NULL) rv->CleanUp();
5637
5638  return R;
5639
5640  // error case:
5641  rInitError:
5642  if  (R != NULL) rDelete(R);
5643  if (rv != NULL) rv->CleanUp();
5644  return NULL;
5645}
5646
5647void rKill(ring r)
5648{
5649  if ((r->ref<=0)&&(r->order!=NULL))
5650  {
5651#ifdef RDEBUG
5652    if (traceit &TRACE_SHOW_RINGS) Print("kill ring %lx\n",(long)r);
5653#endif
5654    if (r->qideal!=NULL)
5655    {
5656      id_Delete(&r->qideal, r);
5657      r->qideal = NULL;
5658    }
5659    int j;
5660#ifdef USE_IILOCALRING
5661    for (j=0;j<myynest;j++)
5662    {
5663      if (iiLocalRing[j]==r)
5664      {
5665        if (j+1==myynest) Warn("killing the basering for level %d",j);
5666        iiLocalRing[j]=NULL;
5667      }
5668    }
5669#else /* USE_IILOCALRING */
5670//#endif /* USE_IILOCALRING */
5671    {
5672      proclevel * nshdl = procstack;
5673      int lev=myynest-1;
5674
5675      for(; nshdl != NULL; nshdl = nshdl->next)
5676      {
5677        if (nshdl->cRing==r)
5678        {
5679          Warn("killing the basering for level %d",lev);
5680          nshdl->cRing=NULL;
5681          nshdl->cRingHdl=NULL;
5682        }
5683      }
5684    }
5685#endif /* USE_IILOCALRING */
5686// any variables depending on r ?
5687    while (r->idroot!=NULL)
5688    {
5689      r->idroot->lev=myynest; // avoid warning about kill global objects
5690      killhdl2(r->idroot,&(r->idroot),r);
5691    }
5692    if (r==currRing)
5693    {
5694      // all dependend stuff is done, clean global vars:
5695      if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
5696      if (sLastPrinted.RingDependend())
5697      {
5698        sLastPrinted.CleanUp();
5699      }
5700      //if ((myynest>0) && (iiRETURNEXPR.RingDependend()))
5701      //{
5702      //  WerrorS("return value depends on local ring variable (export missing ?)");
5703      //  iiRETURNEXPR.CleanUp();
5704      //}
5705      currRing=NULL;
5706      currRingHdl=NULL;
5707    }
5708
5709    /* nKillChar(r); will be called from inside of rDelete */
5710    rDelete(r);
5711    return;
5712  }
5713  r->ref--;
5714}
5715
5716void rKill(idhdl h)
5717{
5718  ring r = IDRING(h);
5719  int ref=0;
5720  if (r!=NULL)
5721  {
5722    ref=r->ref;
5723    rKill(r);
5724  }
5725  if (h==currRingHdl)
5726  {
5727    if (ref<=0) { currRing=NULL; currRingHdl=NULL;}
5728    else
5729    {
5730      currRingHdl=rFindHdl(r,currRingHdl);
5731    }
5732  }
5733}
5734
5735idhdl rSimpleFindHdl(ring r, idhdl root, idhdl n)
5736{
5737  //idhdl next_best=NULL;
5738  idhdl h=root;
5739  while (h!=NULL)
5740  {
5741    if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
5742    && (h!=n)
5743    && (IDRING(h)==r)
5744    )
5745    {
5746   //   if (IDLEV(h)==myynest)
5747   //     return h;
5748   //   if ((IDLEV(h)==0) || (next_best==NULL))
5749   //     next_best=h;
5750   //   else if (IDLEV(next_best)<IDLEV(h))
5751   //     next_best=h;
5752      return h;
5753    }
5754    h=IDNEXT(h);
5755  }
5756  //return next_best;
5757  return NULL;
5758}
5759
5760extern BOOLEAN jjPROC(leftv res, leftv u, leftv v);
5761ideal kGroebner(ideal F, ideal Q)
5762{
5763  //test|=Sy_bit(OPT_PROT);
5764  idhdl save_ringhdl=currRingHdl;
5765  ideal resid;
5766  idhdl new_ring=NULL;
5767  if ((currRingHdl==NULL) || (IDRING(currRingHdl)!=currRing))
5768  {
5769    currRingHdl=enterid(omStrDup(" GROEBNERring"),0,RING_CMD,&IDROOT,FALSE);
5770    new_ring=currRingHdl;
5771    IDRING(currRingHdl)=currRing;
5772  }
5773  sleftv v; memset(&v,0,sizeof(v)); v.rtyp=IDEAL_CMD; v.data=(char *) F;
5774  idhdl h=ggetid("groebner");
5775  sleftv u; memset(&u,0,sizeof(u)); u.rtyp=IDHDL; u.data=(char *) h;
5776            u.name=IDID(h);
5777
5778  sleftv res; memset(&res,0,sizeof(res));
5779  if(jjPROC(&res,&u,&v))
5780  {
5781    resid=kStd(F,Q,testHomog,NULL);
5782  }
5783  else
5784  {
5785    //printf("typ:%d\n",res.rtyp);
5786    resid=(ideal)(res.data);
5787  }
5788  // cleanup GROEBNERring, save_ringhdl, u,v,(res )
5789  if (new_ring!=NULL)
5790  {
5791    idhdl h=IDROOT;
5792    if (h==new_ring) IDROOT=h->next;
5793    else
5794    {
5795      while ((h!=NULL) &&(h->next!=new_ring)) h=h->next;
5796      if (h!=NULL) h->next=h->next->next;
5797    }
5798    if (h!=NULL) omFreeSize(h,sizeof(*h));
5799  }
5800  currRingHdl=save_ringhdl;
5801  u.CleanUp();
5802  v.CleanUp();
5803  return resid;
5804}
5805
5806static void jjINT_S_TO_ID(int n,int *e, leftv res)
5807{
5808  if (n==0) n=1;
5809  ideal l=idInit(n,1);
5810  int i;
5811  poly p;
5812  for(i=rVar(currRing);i>0;i--)
5813  {
5814    if (e[i]>0)
5815    {
5816      n--;
5817      p=pOne();
5818      pSetExp(p,i,1);
5819      pSetm(p);
5820      l->m[n]=p;
5821      if (n==0) break;
5822    }
5823  }
5824  res->data=(char*)l;
5825  setFlag(res,FLAG_STD);
5826  omFreeSize((ADDRESS)e,(rVar(currRing)+1)*sizeof(int));
5827}
5828BOOLEAN jjVARIABLES_P(leftv res, leftv u)
5829{
5830  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5831  int n=pGetVariables((poly)u->Data(),e);
5832  jjINT_S_TO_ID(n,e,res);
5833  return FALSE;
5834}
5835
5836BOOLEAN jjVARIABLES_ID(leftv res, leftv u)
5837{
5838  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5839  ideal I=(ideal)u->Data();
5840  int i;
5841  int n=0;
5842  for(i=I->nrows*I->ncols-1;i>=0;i--)
5843  {
5844    int n0=pGetVariables(I->m[i],e);
5845    if (n0>n) n=n0;
5846  }
5847  jjINT_S_TO_ID(n,e,res);
5848  return FALSE;
5849}
5850
5851void paPrint(const char *n,package p)
5852{
5853  Print(" %s (",n);
5854  switch (p->language)
5855  {
5856    case LANG_SINGULAR: PrintS("S"); break;
5857    case LANG_C:        PrintS("C"); break;
5858    case LANG_TOP:      PrintS("T"); break;
5859    case LANG_NONE:     PrintS("N"); break;
5860    default:            PrintS("U");
5861  }
5862  if(p->libname!=NULL)
5863  Print(",%s", p->libname);
5864  PrintS(")");
5865}
5866
5867BOOLEAN iiApplyINTVEC(leftv res, leftv a, int op, leftv proc)
5868{
5869  intvec *aa=(intvec*)a->Data();
5870  intvec *r=ivCopy(aa);
5871  sleftv tmp_out;
5872  sleftv tmp_in;
5873  BOOLEAN bo=FALSE;
5874  for(int i=0;i<aa->length(); i++)
5875  {
5876    memset(&tmp_in,0,sizeof(tmp_in));
5877    tmp_in.rtyp=INT_CMD;
5878    tmp_in.data=(void*)(long)(*aa)[i];
5879    if (proc==NULL)
5880      bo=iiExprArith1(&tmp_out,&tmp_in,op);
5881    else
5882      bo=jjPROC(&tmp_out,proc,&tmp_in);
5883    if (bo || (tmp_out.rtyp!=INT_CMD))
5884    {
5885      if (r!=NULL) delete r;
5886      Werror("apply fails at index %d",i+1);
5887      return TRUE;
5888    }
5889    (*r)[i]=(int)(long)tmp_out.data;
5890  }
5891  res->data=(void*)r;
5892  return FALSE;
5893}
5894BOOLEAN iiApplyBIGINTMAT(leftv res, leftv a, int op, leftv proc)
5895{
5896  WerrorS("not implemented");
5897  return TRUE;
5898}
5899BOOLEAN iiApplyIDEAL(leftv res, leftv a, int op, leftv proc)
5900{
5901  WerrorS("not implemented");
5902  return TRUE;
5903}
5904BOOLEAN iiApplyLIST(leftv res, leftv a, int op, leftv proc)
5905{
5906  lists aa=(lists)a->Data();
5907  lists r=(lists)omAlloc0Bin(slists_bin); r->Init(aa->nr+1);
5908  sleftv tmp_out;
5909  sleftv tmp_in;
5910  BOOLEAN bo=FALSE;
5911  for(int i=0;i<=aa->nr; i++)
5912  {
5913    memset(&tmp_in,0,sizeof(tmp_in));
5914    tmp_in.Copy(&(aa->m[i]));
5915    if (proc==NULL)
5916      bo=iiExprArith1(&tmp_out,&tmp_in,op);
5917    else
5918      bo=jjPROC(&tmp_out,proc,&tmp_in);
5919    tmp_in.CleanUp();
5920    if (bo)
5921    {
5922      if (r!=NULL) r->Clean();
5923      Werror("apply fails at index %d",i+1);
5924      return TRUE;
5925    }
5926    memcpy(&(r->m[i]),&tmp_out,sizeof(sleftv));
5927  }
5928  res->data=(void*)r;
5929  return FALSE;
5930}
5931BOOLEAN iiApply(leftv res, leftv a, int op, leftv proc)
5932{
5933  memset(res,0,sizeof(sleftv));
5934  res->rtyp=a->Typ();
5935  switch (res->rtyp /*a->Typ()*/)
5936  {
5937    case INTVEC_CMD:
5938    case INTMAT_CMD:
5939        return iiApplyINTVEC(res,a,op,proc);
5940    case BIGINTMAT_CMD:
5941        return iiApplyBIGINTMAT(res,a,op,proc);
5942    case IDEAL_CMD:
5943    case MODUL_CMD:
5944    case MATRIX_CMD:
5945        return iiApplyIDEAL(res,a,op,proc);
5946    case LIST_CMD:
5947        return iiApplyLIST(res,a,op,proc);
5948  }
5949  WerrorS("first argument to `apply` must allow an index");
5950  return TRUE;
5951}
5952
5953BOOLEAN iiTestAssume(leftv a, leftv b)
5954{
5955  // assume a: level
5956  if ((a->Typ()==INT_CMD)&&((long)a->Data()>=0))
5957  {
5958    if ((TEST_V_ALLWARN) && (myynest==0)) WarnS("ASSUME at top level is of no use: see documentation");
5959    char       assume_yylinebuf[80];
5960    strncpy(assume_yylinebuf,my_yylinebuf,79);
5961    int lev=(long)a->Data();
5962    int startlev=0;
5963    idhdl h=ggetid("assumeLevel");
5964    if ((h!=NULL)&&(IDTYP(h)==INT_CMD)) startlev=(long)IDINT(h);
5965    if(lev <=startlev)
5966    {
5967      BOOLEAN bo=b->Eval();
5968      if (bo) { WerrorS("syntax error in ASSUME");return TRUE;}
5969      if (b->Typ()!=INT_CMD) { WerrorS("ASUMME(<level>,<int expr>)");return TRUE; }
5970      if (b->Data()==NULL) { Werror("ASSUME failed:%s",assume_yylinebuf);return TRUE;}
5971    }
5972  }
5973  b->CleanUp();
5974  a->CleanUp();
5975  return FALSE;
5976}
5977
5978#include "libparse.h"
5979
5980BOOLEAN iiARROW(leftv r, char* a,