source: git/Singular/ipshell.cc @ a2e447

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