source: git/Singular/ipshell.cc @ f13c85

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