source: git/Singular/ipshell.cc @ d770e6

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