source: git/Singular/ipshell.cc @ 9051ae

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