source: git/Singular/ipshell.cc

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