source: git/Singular/ipshell.cc @ 40b65a3

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