source: git/Singular/ipshell.cc @ 97181f8

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