source: git/Singular/ipshell.cc @ bdc9282

spielwiese
Last change on this file since bdc9282 was bdc9282, checked in by Hans Schoenemann <hannes@…>, 4 years ago
fix: ringlist->ring for real/complex cf.
  • Property mode set to 100644
File size: 161.5 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          for(;j>=0; j--) (*iv)[j]=1;
1695          break;
1696        default: /* do nothing */;
1697      }
1698    }
1699    else
1700    {
1701      iv=new intvec(1);
1702    }
1703    LLL->m[1].rtyp=INTVEC_CMD;
1704    LLL->m[1].data=(void *)iv;
1705    LL->m[i].data=(void *)LLL;
1706  }
1707  L->m[2].rtyp=LIST_CMD;
1708  L->m[2].data=(void *)LL;
1709  // ----------------------------------------
1710  // 3: qideal
1711  L->m[3].rtyp=IDEAL_CMD;
1712  if (nCoeff_is_transExt(R->cf))
1713    L->m[3].data=(void *)idInit(1,1);
1714  else
1715  {
1716    ideal q=idInit(IDELEMS(r->qideal));
1717    q->m[0]=p_Init(R);
1718    pSetCoeff0(q->m[0],(number)(r->qideal->m[0]));
1719    L->m[3].data=(void *)q;
1720//    I->m[0] = pNSet(R->minpoly);
1721  }
1722  // ----------------------------------------
1723}
1724static void rDecomposeC_41(leftv h,const coeffs C)
1725/* field is R or C */
1726{
1727  lists L=(lists)omAlloc0Bin(slists_bin);
1728  if (nCoeff_is_long_C(C)) L->Init(3);
1729  else                     L->Init(2);
1730  h->rtyp=LIST_CMD;
1731  h->data=(void *)L;
1732  // 0: char/ cf - ring
1733  // 1: list (var)
1734  // 2: list (ord)
1735  // ----------------------------------------
1736  // 0: char/ cf - ring
1737  L->m[0].rtyp=INT_CMD;
1738  L->m[0].data=(void *)0;
1739  // ----------------------------------------
1740  // 1:
1741  lists LL=(lists)omAlloc0Bin(slists_bin);
1742  LL->Init(2);
1743    LL->m[0].rtyp=INT_CMD;
1744    LL->m[0].data=(void *)(long)si_max(C->float_len,SHORT_REAL_LENGTH/2);
1745    LL->m[1].rtyp=INT_CMD;
1746    LL->m[1].data=(void *)(long)si_max(C->float_len2,SHORT_REAL_LENGTH);
1747  L->m[1].rtyp=LIST_CMD;
1748  L->m[1].data=(void *)LL;
1749  // ----------------------------------------
1750  // 2: list (par)
1751  if (nCoeff_is_long_C(C))
1752  {
1753    L->m[2].rtyp=STRING_CMD;
1754    L->m[2].data=(void *)omStrDup(*n_ParameterNames(C));
1755  }
1756  // ----------------------------------------
1757}
1758static void rDecomposeC(leftv h,const ring R)
1759/* field is R or C */
1760{
1761  lists L=(lists)omAlloc0Bin(slists_bin);
1762  if (rField_is_long_C(R)) L->Init(3);
1763  else                     L->Init(2);
1764  h->rtyp=LIST_CMD;
1765  h->data=(void *)L;
1766  // 0: char/ cf - ring
1767  // 1: list (var)
1768  // 2: list (ord)
1769  // ----------------------------------------
1770  // 0: char/ cf - ring
1771  L->m[0].rtyp=INT_CMD;
1772  L->m[0].data=(void *)0;
1773  // ----------------------------------------
1774  // 1:
1775  lists LL=(lists)omAlloc0Bin(slists_bin);
1776  LL->Init(2);
1777    LL->m[0].rtyp=INT_CMD;
1778    LL->m[0].data=(void *)(long)si_max(R->cf->float_len,SHORT_REAL_LENGTH/2);
1779    LL->m[1].rtyp=INT_CMD;
1780    LL->m[1].data=(void *)(long)si_max(R->cf->float_len2,SHORT_REAL_LENGTH);
1781  L->m[1].rtyp=LIST_CMD;
1782  L->m[1].data=(void *)LL;
1783  // ----------------------------------------
1784  // 2: list (par)
1785  if (rField_is_long_C(R))
1786  {
1787    L->m[2].rtyp=STRING_CMD;
1788    L->m[2].data=(void *)omStrDup(*rParameter(R));
1789  }
1790  // ----------------------------------------
1791}
1792
1793#ifdef HAVE_RINGS
1794void rDecomposeRing_41(leftv h,const coeffs C)
1795/* field is R or C */
1796{
1797  lists L=(lists)omAlloc0Bin(slists_bin);
1798  if (nCoeff_is_Ring(C)) L->Init(1);
1799  else                   L->Init(2);
1800  h->rtyp=LIST_CMD;
1801  h->data=(void *)L;
1802  // 0: char/ cf - ring
1803  // 1: list (module)
1804  // ----------------------------------------
1805  // 0: char/ cf - ring
1806  L->m[0].rtyp=STRING_CMD;
1807  L->m[0].data=(void *)omStrDup("integer");
1808  // ----------------------------------------
1809  // 1: modulo
1810  if (nCoeff_is_Z(C)) return;
1811  lists LL=(lists)omAlloc0Bin(slists_bin);
1812  LL->Init(2);
1813  LL->m[0].rtyp=BIGINT_CMD;
1814  LL->m[0].data=nlMapGMP((number) C->modBase, C, coeffs_BIGINT);
1815  LL->m[1].rtyp=INT_CMD;
1816  LL->m[1].data=(void *) C->modExponent;
1817  L->m[1].rtyp=LIST_CMD;
1818  L->m[1].data=(void *)LL;
1819}
1820#endif
1821
1822void rDecomposeRing(leftv h,const ring R)
1823/* field is R or C */
1824{
1825#ifdef HAVE_RINGS
1826  lists L=(lists)omAlloc0Bin(slists_bin);
1827  if (rField_is_Z(R)) L->Init(1);
1828  else                     L->Init(2);
1829  h->rtyp=LIST_CMD;
1830  h->data=(void *)L;
1831  // 0: char/ cf - ring
1832  // 1: list (module)
1833  // ----------------------------------------
1834  // 0: char/ cf - ring
1835  L->m[0].rtyp=STRING_CMD;
1836  L->m[0].data=(void *)omStrDup("integer");
1837  // ----------------------------------------
1838  // 1: module
1839  if (rField_is_Z(R)) return;
1840  lists LL=(lists)omAlloc0Bin(slists_bin);
1841  LL->Init(2);
1842  LL->m[0].rtyp=BIGINT_CMD;
1843  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); ?
1844  LL->m[1].rtyp=INT_CMD;
1845  LL->m[1].data=(void *) R->cf->modExponent;
1846  L->m[1].rtyp=LIST_CMD;
1847  L->m[1].data=(void *)LL;
1848#else
1849  WerrorS("rDecomposeRing");
1850#endif
1851}
1852
1853
1854BOOLEAN rDecompose_CF(leftv res,const coeffs C)
1855{
1856  assume( C != NULL );
1857
1858  // sanity check: require currRing==r for rings with polynomial data
1859  if ( nCoeff_is_algExt(C) && (C != currRing->cf))
1860  {
1861    WerrorS("ring with polynomial data must be the base ring or compatible");
1862    return TRUE;
1863  }
1864  if (nCoeff_is_numeric(C))
1865  {
1866    rDecomposeC_41(res,C);
1867  }
1868#ifdef HAVE_RINGS
1869  else if (nCoeff_is_Ring(C))
1870  {
1871    rDecomposeRing_41(res,C);
1872  }
1873#endif
1874  else if ( C->extRing!=NULL )// nCoeff_is_algExt(r->cf))
1875  {
1876    rDecomposeCF(res, C->extRing, currRing);
1877  }
1878  else if(nCoeff_is_GF(C))
1879  {
1880    lists Lc=(lists)omAlloc0Bin(slists_bin);
1881    Lc->Init(4);
1882    // char:
1883    Lc->m[0].rtyp=INT_CMD;
1884    Lc->m[0].data=(void*)(long)C->m_nfCharQ;
1885    // var:
1886    lists Lv=(lists)omAlloc0Bin(slists_bin);
1887    Lv->Init(1);
1888    Lv->m[0].rtyp=STRING_CMD;
1889    Lv->m[0].data=(void *)omStrDup(*n_ParameterNames(C));
1890    Lc->m[1].rtyp=LIST_CMD;
1891    Lc->m[1].data=(void*)Lv;
1892    // ord:
1893    lists Lo=(lists)omAlloc0Bin(slists_bin);
1894    Lo->Init(1);
1895    lists Loo=(lists)omAlloc0Bin(slists_bin);
1896    Loo->Init(2);
1897    Loo->m[0].rtyp=STRING_CMD;
1898    Loo->m[0].data=(void *)omStrDup(rSimpleOrdStr(ringorder_lp));
1899
1900    intvec *iv=new intvec(1); (*iv)[0]=1;
1901    Loo->m[1].rtyp=INTVEC_CMD;
1902    Loo->m[1].data=(void *)iv;
1903
1904    Lo->m[0].rtyp=LIST_CMD;
1905    Lo->m[0].data=(void*)Loo;
1906
1907    Lc->m[2].rtyp=LIST_CMD;
1908    Lc->m[2].data=(void*)Lo;
1909    // q-ideal:
1910    Lc->m[3].rtyp=IDEAL_CMD;
1911    Lc->m[3].data=(void *)idInit(1,1);
1912    // ----------------------
1913    res->rtyp=LIST_CMD;
1914    res->data=(void*)Lc;
1915  }
1916  else
1917  {
1918    res->rtyp=INT_CMD;
1919    res->data=(void *)(long)C->ch;
1920  }
1921  // ----------------------------------------
1922  return FALSE;
1923}
1924
1925lists rDecompose_list_cf(const ring r)
1926{
1927  assume( r != NULL );
1928  const coeffs C = r->cf;
1929  assume( C != NULL );
1930
1931  // sanity check: require currRing==r for rings with polynomial data
1932  if ( (r!=currRing) && (
1933        (r->qideal != NULL)
1934#ifdef HAVE_PLURAL
1935        || (rIsPluralRing(r))
1936#endif
1937                        )
1938     )
1939  {
1940    WerrorS("ring with polynomial data must be the base ring or compatible");
1941    return NULL;
1942  }
1943  // 0: char/ cf - ring
1944  // 1: list (var)
1945  // 2: list (ord)
1946  // 3: qideal
1947  // possibly:
1948  // 4: C
1949  // 5: D
1950  lists L=(lists)omAlloc0Bin(slists_bin);
1951  if (rIsPluralRing(r))
1952    L->Init(6);
1953  else
1954    L->Init(4);
1955  // ----------------------------------------
1956  // 0: char/ cf - ring
1957  L->m[0].rtyp=CRING_CMD;
1958  L->m[0].data=(char*)r->cf; r->cf->ref++;
1959  // ----------------------------------------
1960  // 1: list (var)
1961  lists LL=(lists)omAlloc0Bin(slists_bin);
1962  LL->Init(r->N);
1963  int i;
1964  for(i=0; i<r->N; i++)
1965  {
1966    LL->m[i].rtyp=STRING_CMD;
1967    LL->m[i].data=(void *)omStrDup(r->names[i]);
1968  }
1969  L->m[1].rtyp=LIST_CMD;
1970  L->m[1].data=(void *)LL;
1971  // ----------------------------------------
1972  // 2: list (ord)
1973  LL=(lists)omAlloc0Bin(slists_bin);
1974  i=rBlocks(r)-1;
1975  LL->Init(i);
1976  i--;
1977  lists LLL;
1978  for(; i>=0; i--)
1979  {
1980    intvec *iv;
1981    int j;
1982    LL->m[i].rtyp=LIST_CMD;
1983    LLL=(lists)omAlloc0Bin(slists_bin);
1984    LLL->Init(2);
1985    LLL->m[0].rtyp=STRING_CMD;
1986    LLL->m[0].data=(void *)omStrDup(rSimpleOrdStr(r->order[i]));
1987
1988    if(r->order[i] == ringorder_IS) //  || r->order[i] == ringorder_s || r->order[i] == ringorder_S)
1989    {
1990      assume( r->block0[i] == r->block1[i] );
1991      const int s = r->block0[i];
1992      assume( -2 < s && s < 2);
1993
1994      iv=new intvec(1);
1995      (*iv)[0] = s;
1996    }
1997    else if (r->block1[i]-r->block0[i] >=0 )
1998    {
1999      int bl=j=r->block1[i]-r->block0[i];
2000      if (r->order[i]==ringorder_M)
2001      {
2002        j=(j+1)*(j+1)-1;
2003        bl=j+1;
2004      }
2005      else if (r->order[i]==ringorder_am)
2006      {
2007        j+=r->wvhdl[i][bl+1];
2008      }
2009      iv=new intvec(j+1);
2010      if ((r->wvhdl!=NULL) && (r->wvhdl[i]!=NULL))
2011      {
2012        for(;j>=0; j--) (*iv)[j]=r->wvhdl[i][j+(j>bl)];
2013      }
2014      else switch (r->order[i])
2015      {
2016        case ringorder_dp:
2017        case ringorder_Dp:
2018        case ringorder_ds:
2019        case ringorder_Ds:
2020        case ringorder_lp:
2021          for(;j>=0; j--) (*iv)[j]=1;
2022          break;
2023        default: /* do nothing */;
2024      }
2025    }
2026    else
2027    {
2028      iv=new intvec(1);
2029    }
2030    LLL->m[1].rtyp=INTVEC_CMD;
2031    LLL->m[1].data=(void *)iv;
2032    LL->m[i].data=(void *)LLL;
2033  }
2034  L->m[2].rtyp=LIST_CMD;
2035  L->m[2].data=(void *)LL;
2036  // ----------------------------------------
2037  // 3: qideal
2038  L->m[3].rtyp=IDEAL_CMD;
2039  if (r->qideal==NULL)
2040    L->m[3].data=(void *)idInit(1,1);
2041  else
2042    L->m[3].data=(void *)idCopy(r->qideal);
2043  // ----------------------------------------
2044#ifdef HAVE_PLURAL // NC! in rDecompose
2045  if (rIsPluralRing(r))
2046  {
2047    L->m[4].rtyp=MATRIX_CMD;
2048    L->m[4].data=(void *)mp_Copy(r->GetNC()->C, r, r);
2049    L->m[5].rtyp=MATRIX_CMD;
2050    L->m[5].data=(void *)mp_Copy(r->GetNC()->D, r, r);
2051  }
2052#endif
2053  return L;
2054}
2055
2056lists rDecompose(const ring r)
2057{
2058  assume( r != NULL );
2059  const coeffs C = r->cf;
2060  assume( C != NULL );
2061
2062  // sanity check: require currRing==r for rings with polynomial data
2063  if ( (r!=currRing) && (
2064           (nCoeff_is_algExt(C) && (C != currRing->cf))
2065        || (r->qideal != NULL)
2066#ifdef HAVE_PLURAL
2067        || (rIsPluralRing(r))
2068#endif
2069                        )
2070     )
2071  {
2072    WerrorS("ring with polynomial data must be the base ring or compatible");
2073    return NULL;
2074  }
2075  // 0: char/ cf - ring
2076  // 1: list (var)
2077  // 2: list (ord)
2078  // 3: qideal
2079  // possibly:
2080  // 4: C
2081  // 5: D
2082  lists L=(lists)omAlloc0Bin(slists_bin);
2083  if (rIsPluralRing(r))
2084    L->Init(6);
2085  else
2086    L->Init(4);
2087  // ----------------------------------------
2088  // 0: char/ cf - ring
2089  if (rField_is_numeric(r))
2090  {
2091    rDecomposeC(&(L->m[0]),r);
2092  }
2093  else if (rField_is_Ring(r))
2094  {
2095    rDecomposeRing(&(L->m[0]),r);
2096  }
2097  else if ( r->cf->extRing!=NULL )// nCoeff_is_algExt(r->cf))
2098  {
2099    rDecomposeCF(&(L->m[0]), r->cf->extRing, r);
2100  }
2101  else if(rField_is_GF(r))
2102  {
2103    lists Lc=(lists)omAlloc0Bin(slists_bin);
2104    Lc->Init(4);
2105    // char:
2106    Lc->m[0].rtyp=INT_CMD;
2107    Lc->m[0].data=(void*)(long)r->cf->m_nfCharQ;
2108    // var:
2109    lists Lv=(lists)omAlloc0Bin(slists_bin);
2110    Lv->Init(1);
2111    Lv->m[0].rtyp=STRING_CMD;
2112    Lv->m[0].data=(void *)omStrDup(*rParameter(r));
2113    Lc->m[1].rtyp=LIST_CMD;
2114    Lc->m[1].data=(void*)Lv;
2115    // ord:
2116    lists Lo=(lists)omAlloc0Bin(slists_bin);
2117    Lo->Init(1);
2118    lists Loo=(lists)omAlloc0Bin(slists_bin);
2119    Loo->Init(2);
2120    Loo->m[0].rtyp=STRING_CMD;
2121    Loo->m[0].data=(void *)omStrDup(rSimpleOrdStr(ringorder_lp));
2122
2123    intvec *iv=new intvec(1); (*iv)[0]=1;
2124    Loo->m[1].rtyp=INTVEC_CMD;
2125    Loo->m[1].data=(void *)iv;
2126
2127    Lo->m[0].rtyp=LIST_CMD;
2128    Lo->m[0].data=(void*)Loo;
2129
2130    Lc->m[2].rtyp=LIST_CMD;
2131    Lc->m[2].data=(void*)Lo;
2132    // q-ideal:
2133    Lc->m[3].rtyp=IDEAL_CMD;
2134    Lc->m[3].data=(void *)idInit(1,1);
2135    // ----------------------
2136    L->m[0].rtyp=LIST_CMD;
2137    L->m[0].data=(void*)Lc;
2138  }
2139  else
2140  {
2141    L->m[0].rtyp=INT_CMD;
2142    L->m[0].data=(void *)(long)r->cf->ch;
2143  }
2144  // ----------------------------------------
2145  // 1: list (var)
2146  lists LL=(lists)omAlloc0Bin(slists_bin);
2147  LL->Init(r->N);
2148  int i;
2149  for(i=0; i<r->N; i++)
2150  {
2151    LL->m[i].rtyp=STRING_CMD;
2152    LL->m[i].data=(void *)omStrDup(r->names[i]);
2153  }
2154  L->m[1].rtyp=LIST_CMD;
2155  L->m[1].data=(void *)LL;
2156  // ----------------------------------------
2157  // 2: list (ord)
2158  LL=(lists)omAlloc0Bin(slists_bin);
2159  i=rBlocks(r)-1;
2160  LL->Init(i);
2161  i--;
2162  lists LLL;
2163  for(; i>=0; i--)
2164  {
2165    intvec *iv;
2166    int j;
2167    LL->m[i].rtyp=LIST_CMD;
2168    LLL=(lists)omAlloc0Bin(slists_bin);
2169    LLL->Init(2);
2170    LLL->m[0].rtyp=STRING_CMD;
2171    LLL->m[0].data=(void *)omStrDup(rSimpleOrdStr(r->order[i]));
2172
2173    if((r->order[i] == ringorder_IS)
2174    || (r->order[i] == ringorder_s)) //|| r->order[i] == ringorder_S)
2175    {
2176      assume( r->block0[i] == r->block1[i] );
2177      const int s = r->block0[i];
2178      assume( (-2 < s && s < 2)||(r->order[i] != ringorder_IS));
2179
2180      iv=new intvec(1);
2181      (*iv)[0] = s;
2182    }
2183    else if (r->block1[i]-r->block0[i] >=0 )
2184    {
2185      int bl=j=r->block1[i]-r->block0[i];
2186      if (r->order[i]==ringorder_M)
2187      {
2188        j=(j+1)*(j+1)-1;
2189        bl=j+1;
2190      }
2191      else if (r->order[i]==ringorder_am)
2192      {
2193        j+=r->wvhdl[i][bl+1];
2194      }
2195      iv=new intvec(j+1);
2196      if ((r->wvhdl!=NULL) && (r->wvhdl[i]!=NULL))
2197      {
2198        for(;j>=0; j--) (*iv)[j]=r->wvhdl[i][j+(j>bl)];
2199      }
2200      else switch (r->order[i])
2201      {
2202        case ringorder_dp:
2203        case ringorder_Dp:
2204        case ringorder_ds:
2205        case ringorder_Ds:
2206        case ringorder_lp:
2207          for(;j>=0; j--) (*iv)[j]=1;
2208          break;
2209        default: /* do nothing */;
2210      }
2211    }
2212    else
2213    {
2214      iv=new intvec(1);
2215    }
2216    LLL->m[1].rtyp=INTVEC_CMD;
2217    LLL->m[1].data=(void *)iv;
2218    LL->m[i].data=(void *)LLL;
2219  }
2220  L->m[2].rtyp=LIST_CMD;
2221  L->m[2].data=(void *)LL;
2222  // ----------------------------------------
2223  // 3: qideal
2224  L->m[3].rtyp=IDEAL_CMD;
2225  if (r->qideal==NULL)
2226    L->m[3].data=(void *)idInit(1,1);
2227  else
2228    L->m[3].data=(void *)idCopy(r->qideal);
2229  // ----------------------------------------
2230#ifdef HAVE_PLURAL // NC! in rDecompose
2231  if (rIsPluralRing(r))
2232  {
2233    L->m[4].rtyp=MATRIX_CMD;
2234    L->m[4].data=(void *)mp_Copy(r->GetNC()->C, r, r);
2235    L->m[5].rtyp=MATRIX_CMD;
2236    L->m[5].data=(void *)mp_Copy(r->GetNC()->D, r, r);
2237  }
2238#endif
2239  return L;
2240}
2241
2242void rComposeC(lists L, ring R)
2243/* field is R or C */
2244{
2245  // ----------------------------------------
2246  // 0: char/ cf - ring
2247  if ((L->m[0].rtyp!=INT_CMD) || (L->m[0].data!=(char *)0))
2248  {
2249    WerrorS("invalid coeff. field description, expecting 0");
2250    return;
2251  }
2252//  R->cf->ch=0;
2253  // ----------------------------------------
2254  // 0, (r1,r2) [, "i" ]
2255  if (L->m[1].rtyp!=LIST_CMD)
2256  {
2257    WerrorS("invalid coeff. field description, expecting precision list");
2258    return;
2259  }
2260  lists LL=(lists)L->m[1].data;
2261  if ((LL->nr!=1)
2262    || (LL->m[0].rtyp!=INT_CMD)
2263    || (LL->m[1].rtyp!=INT_CMD))
2264  {
2265    WerrorS("invalid coeff. field description list, expected list(`int`,`int`)");
2266    return;
2267  }
2268  int r1=(int)(long)LL->m[0].data;
2269  int r2=(int)(long)LL->m[1].data;
2270  r1=si_min(r1,32767);
2271  r2=si_min(r2,32767);
2272  LongComplexInfo par; memset(&par, 0, sizeof(par));
2273  par.float_len=r1;
2274  par.float_len2=r2;
2275  if (L->nr==2) // complex
2276  {
2277    if (L->m[2].rtyp!=STRING_CMD)
2278    {
2279      WerrorS("invalid coeff. field description, expecting parameter name");
2280      return;
2281    }
2282    par.par_name=(char*)L->m[2].data;
2283    R->cf = nInitChar(n_long_C, &par);
2284  }
2285  else if ((r1<=SHORT_REAL_LENGTH) && (r2<=SHORT_REAL_LENGTH)) /* && L->nr==1*/
2286    R->cf = nInitChar(n_R, NULL);
2287  else /* && L->nr==1*/
2288  {
2289    R->cf = nInitChar(n_long_R, &par);
2290  }
2291}
2292
2293#ifdef HAVE_RINGS
2294void rComposeRing(lists L, ring R)
2295/* field is R or C */
2296{
2297  // ----------------------------------------
2298  // 0: string: integer
2299  // no further entries --> Z
2300  mpz_t modBase;
2301  unsigned int modExponent = 1;
2302
2303  if (L->nr == 0)
2304  {
2305    mpz_init_set_ui(modBase,0);
2306    modExponent = 1;
2307  }
2308  // ----------------------------------------
2309  // 1:
2310  else
2311  {
2312    if (L->m[1].rtyp!=LIST_CMD) WerrorS("invalid data, expecting list of numbers");
2313    lists LL=(lists)L->m[1].data;
2314    if ((LL->nr >= 0) && LL->m[0].rtyp == BIGINT_CMD)
2315    {
2316      number tmp= (number) LL->m[0].data; // never use CopyD() on list elements
2317                                    // assume that tmp is integer, not rational
2318      mpz_init(modBase);
2319      n_MPZ (modBase, tmp, coeffs_BIGINT);
2320    }
2321    else if (LL->nr >= 0 && LL->m[0].rtyp == INT_CMD)
2322    {
2323      mpz_init_set_ui(modBase,(unsigned long) LL->m[0].data);
2324    }
2325    else
2326    {
2327      mpz_init_set_ui(modBase,0);
2328    }
2329    if (LL->nr >= 1)
2330    {
2331      modExponent = (unsigned long) LL->m[1].data;
2332    }
2333    else
2334    {
2335      modExponent = 1;
2336    }
2337  }
2338  // ----------------------------------------
2339  if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_sgn1(modBase) < 0))
2340  {
2341    WerrorS("Wrong ground ring specification (module is 1)");
2342    return;
2343  }
2344  if (modExponent < 1)
2345  {
2346    WerrorS("Wrong ground ring specification (exponent smaller than 1)");
2347    return;
2348  }
2349  // module is 0 ---> integers
2350  if (mpz_sgn1(modBase) == 0)
2351  {
2352    R->cf=nInitChar(n_Z,NULL);
2353  }
2354  // we have an exponent
2355  else if (modExponent > 1)
2356  {
2357    //R->cf->ch = R->cf->modExponent;
2358    if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
2359    {
2360      /* this branch should be active for modExponent = 2..32 resp. 2..64,
2361           depending on the size of a long on the respective platform */
2362      R->cf=nInitChar(n_Z2m,(void*)(long)modExponent);       // Use Z/2^ch
2363    }
2364    else
2365    {
2366      //ringtype 3
2367      ZnmInfo info;
2368      info.base= modBase;
2369      info.exp= modExponent;
2370      R->cf=nInitChar(n_Znm,(void*) &info);
2371    }
2372  }
2373  // just a module m > 1
2374  else
2375  {
2376    //ringtype = 2;
2377    //const int ch = mpz_get_ui(modBase);
2378    ZnmInfo info;
2379    info.base= modBase;
2380    info.exp= modExponent;
2381    R->cf=nInitChar(n_Zn,(void*) &info);
2382  }
2383  mpz_clear(modBase);
2384}
2385#endif
2386
2387static void rRenameVars(ring R)
2388{
2389  int i,j;
2390  BOOLEAN ch;
2391  do
2392  {
2393    ch=0;
2394    for(i=0;i<R->N-1;i++)
2395    {
2396      for(j=i+1;j<R->N;j++)
2397      {
2398        if (strcmp(R->names[i],R->names[j])==0)
2399        {
2400          ch=TRUE;
2401          Warn("name conflict var(%d) and var(%d): `%s`, rename to `@%s`",i+1,j+1,R->names[i],R->names[i]);
2402          omFree(R->names[j]);
2403          R->names[j]=(char *)omAlloc(2+strlen(R->names[i]));
2404          sprintf(R->names[j],"@%s",R->names[i]);
2405        }
2406      }
2407    }
2408  }
2409  while (ch);
2410  for(i=0;i<rPar(R); i++)
2411  {
2412    for(j=0;j<R->N;j++)
2413    {
2414      if (strcmp(rParameter(R)[i],R->names[j])==0)
2415      {
2416        Warn("name conflict par(%d) and var(%d): `%s`, renaming the VARIABLE to `@@(%d)`",i+1,j+1,R->names[j],i+1);
2417//        omFree(rParameter(R)[i]);
2418//        rParameter(R)[i]=(char *)omAlloc(10);
2419//        sprintf(rParameter(R)[i],"@@(%d)",i+1);
2420        omFree(R->names[j]);
2421        R->names[j]=(char *)omAlloc(10);
2422        sprintf(R->names[j],"@@(%d)",i+1);
2423      }
2424    }
2425  }
2426}
2427
2428static inline BOOLEAN rComposeVar(const lists  L, ring R)
2429{
2430  assume(R!=NULL);
2431  if (L->m[1].Typ()==LIST_CMD)
2432  {
2433    lists v=(lists)L->m[1].Data();
2434    R->N = v->nr+1;
2435    if (R->N<=0)
2436    {
2437      WerrorS("no ring variables");
2438      return TRUE;
2439    }
2440    R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
2441    int i;
2442    for(i=0;i<R->N;i++)
2443    {
2444      if (v->m[i].Typ()==STRING_CMD)
2445        R->names[i]=omStrDup((char *)v->m[i].Data());
2446      else if (v->m[i].Typ()==POLY_CMD)
2447      {
2448        poly p=(poly)v->m[i].Data();
2449        int nr=pIsPurePower(p);
2450        if (nr>0)
2451          R->names[i]=omStrDup(currRing->names[nr-1]);
2452        else
2453        {
2454          Werror("var name %d must be a string or a ring variable",i+1);
2455          return TRUE;
2456        }
2457      }
2458      else
2459      {
2460        Werror("var name %d must be `string` (not %d)",i+1, v->m[i].Typ());
2461        return TRUE;
2462      }
2463    }
2464  }
2465  else
2466  {
2467    WerrorS("variable must be given as `list`");
2468    return TRUE;
2469  }
2470  return FALSE;
2471}
2472
2473static inline BOOLEAN rComposeOrder(const lists  L, const BOOLEAN check_comp, ring R)
2474{
2475  assume(R!=NULL);
2476  long bitmask=0L;
2477  if (L->m[2].Typ()==LIST_CMD)
2478  {
2479    lists v=(lists)L->m[2].Data();
2480    int n= v->nr+2;
2481    int j_in_R,j_in_L;
2482    // do we have an entry "L",... ?: set bitmask
2483    for (int j=0; j < n-1; j++)
2484    {
2485      if (v->m[j].Typ()==LIST_CMD)
2486      {
2487        lists vv=(lists)v->m[j].Data();
2488        if ((vv->nr==1)
2489        &&(vv->m[0].Typ()==STRING_CMD)
2490        &&(strcmp((char*)vv->m[0].Data(),"L")==0))
2491        {
2492          number nn=(number)vv->m[1].Data();
2493          if (vv->m[1].Typ()==BIGINT_CMD)
2494            bitmask=n_Int(nn,coeffs_BIGINT);
2495          else if (vv->m[1].Typ()==INT_CMD)
2496            bitmask=(long)nn;
2497          else
2498          {
2499            Werror("illegal argument for pseudo ordering L: %d",vv->m[1].Typ());
2500            return TRUE;
2501          }
2502          break;
2503        }
2504      }
2505    }
2506    if (bitmask!=0) n--;
2507
2508    // initialize fields of R
2509    R->order=(rRingOrder_t *)omAlloc0(n*sizeof(rRingOrder_t));
2510    R->block0=(int *)omAlloc0(n*sizeof(int));
2511    R->block1=(int *)omAlloc0(n*sizeof(int));
2512    R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
2513    // init order, so that rBlocks works correctly
2514    for (j_in_R= n-2; j_in_R>=0; j_in_R--)
2515      R->order[j_in_R] = ringorder_unspec;
2516    // orderings
2517    for(j_in_R=0,j_in_L=0;j_in_R<n-1;j_in_R++,j_in_L++)
2518    {
2519    // todo: a(..), M
2520      if (v->m[j_in_L].Typ()!=LIST_CMD)
2521      {
2522        WerrorS("ordering must be list of lists");
2523        return TRUE;
2524      }
2525      lists vv=(lists)v->m[j_in_L].Data();
2526      if ((vv->nr==1)
2527      && (vv->m[0].Typ()==STRING_CMD))
2528      {
2529        if (strcmp((char*)vv->m[0].Data(),"L")==0)
2530        {
2531          j_in_R--;
2532          continue;
2533        }
2534        if ((vv->m[1].Typ()!=INTVEC_CMD) && (vv->m[1].Typ()!=INT_CMD))
2535        {
2536          PrintS(lString(vv));
2537          WerrorS("ordering name must be a (string,intvec)(1)");
2538          return TRUE;
2539        }
2540        R->order[j_in_R]=rOrderName(omStrDup((char*)vv->m[0].Data())); // assume STRING
2541
2542        if (j_in_R==0) R->block0[0]=1;
2543        else
2544        {
2545          int jj=j_in_R-1;
2546          while((jj>=0)
2547          && ((R->order[jj]== ringorder_a)
2548             || (R->order[jj]== ringorder_aa)
2549             || (R->order[jj]== ringorder_am)
2550             || (R->order[jj]== ringorder_c)
2551             || (R->order[jj]== ringorder_C)
2552             || (R->order[jj]== ringorder_s)
2553             || (R->order[jj]== ringorder_S)
2554          ))
2555          {
2556            //Print("jj=%, skip %s\n",rSimpleOrdStr(R->order[jj]));
2557            jj--;
2558          }
2559          if (jj<0) R->block0[j_in_R]=1;
2560          else      R->block0[j_in_R]=R->block1[jj]+1;
2561        }
2562        intvec *iv;
2563        if (vv->m[1].Typ()==INT_CMD)
2564          iv=new intvec((int)(long)vv->m[1].Data(),(int)(long)vv->m[1].Data());
2565        else
2566          iv=ivCopy((intvec*)vv->m[1].Data()); //assume INTVEC
2567        int iv_len=iv->length();
2568        if (iv_len==0)
2569        {
2570          Werror("empty intvec for ordering %d (%s)",j_in_R,rSimpleOrdStr(R->order[j_in_R]));
2571          return TRUE;
2572        }
2573        if ((R->order[j_in_R]!=ringorder_s)
2574        &&(R->order[j_in_R]!=ringorder_c)
2575        &&(R->order[j_in_R]!=ringorder_C))
2576        {
2577          R->block1[j_in_R]=si_max(R->block0[j_in_R],R->block0[j_in_R]+iv_len-1);
2578          if (R->block1[j_in_R]>R->N)
2579          {
2580            if (R->block0[j_in_R]>R->N)
2581            {
2582              Werror("not enough variables for ordering %d (%s)",j_in_R,rSimpleOrdStr(R->order[j_in_R]));
2583              return TRUE;
2584            }
2585            R->block1[j_in_R]=R->N;
2586            iv_len=R->block1[j_in_R]-R->block0[j_in_R]+1;
2587          }
2588          //Print("block %d from %d to %d\n",j,R->block0[j], R->block1[j]);
2589        }
2590        int i;
2591        switch (R->order[j_in_R])
2592        {
2593           case ringorder_ws:
2594           case ringorder_Ws:
2595              R->OrdSgn=-1;
2596           case ringorder_aa:
2597           case ringorder_a:
2598           case ringorder_wp:
2599           case ringorder_Wp:
2600             R->wvhdl[j_in_R] =( int *)omAlloc(iv_len*sizeof(int));
2601             for (i=0; i<iv_len;i++)
2602             {
2603               R->wvhdl[j_in_R][i]=(*iv)[i];
2604             }
2605             break;
2606           case ringorder_am:
2607             R->wvhdl[j_in_R] =( int *)omAlloc((iv->length()+1)*sizeof(int));
2608             for (i=0; i<iv_len;i++)
2609             {
2610               R->wvhdl[j_in_R][i]=(*iv)[i];
2611             }
2612             R->wvhdl[j_in_R][i]=iv->length() - iv_len;
2613             //printf("ivlen:%d,iv->len:%d,mod:%d\n",iv_len,iv->length(),R->wvhdl[j][i]);
2614             for (; i<iv->length(); i++)
2615             {
2616                R->wvhdl[j_in_R][i+1]=(*iv)[i];
2617             }
2618             break;
2619           case ringorder_M:
2620             R->wvhdl[j_in_R] =( int *)omAlloc((iv->length())*sizeof(int));
2621             for (i=0; i<iv->length();i++) R->wvhdl[j_in_R][i]=(*iv)[i];
2622             R->block1[j_in_R]=si_max(R->block0[j_in_R],R->block0[j_in_R]+(int)sqrt((double)(iv->length()-1)));
2623             if (R->block1[j_in_R]>R->N)
2624             {
2625               WerrorS("ordering matrix too big");
2626               return TRUE;
2627             }
2628             break;
2629           case ringorder_ls:
2630           case ringorder_ds:
2631           case ringorder_Ds:
2632           case ringorder_rs:
2633             R->OrdSgn=-1;
2634           case ringorder_lp:
2635           case ringorder_dp:
2636           case ringorder_Dp:
2637           case ringorder_rp:
2638             break;
2639           case ringorder_S:
2640             break;
2641           case ringorder_c:
2642           case ringorder_C:
2643             R->block1[j_in_R]=R->block0[j_in_R]=0;
2644             break;
2645
2646           case ringorder_s:
2647             R->block1[j_in_R]=R->block0[j_in_R]=(*iv)[0];
2648             rSetSyzComp(R->block0[j_in_R],R);
2649             break;
2650
2651           case ringorder_IS:
2652           {
2653             R->block1[j_in_R] = R->block0[j_in_R] = 0;
2654             if( iv->length() > 0 )
2655             {
2656               const int s = (*iv)[0];
2657               assume( -2 < s && s < 2 );
2658               R->block1[j_in_R] = R->block0[j_in_R] = s;
2659             }
2660             break;
2661           }
2662           case 0:
2663           case ringorder_unspec:
2664             break;
2665           case ringorder_L: /* cannot happen */
2666           case ringorder_a64: /*not implemented */
2667             WerrorS("ring order not implemented");
2668             return TRUE;
2669        }
2670        delete iv;
2671      }
2672      else
2673      {
2674        PrintS(lString(vv));
2675        WerrorS("ordering name must be a (string,intvec)");
2676        return TRUE;
2677      }
2678    }
2679    // sanity check
2680    j_in_R=n-2;
2681    if ((R->order[j_in_R]==ringorder_c)
2682    || (R->order[j_in_R]==ringorder_C)
2683    || (R->order[j_in_R]==ringorder_unspec)) j_in_R--;
2684    if (R->block1[j_in_R] != R->N)
2685    {
2686      if (((R->order[j_in_R]==ringorder_dp) ||
2687           (R->order[j_in_R]==ringorder_ds) ||
2688           (R->order[j_in_R]==ringorder_Dp) ||
2689           (R->order[j_in_R]==ringorder_Ds) ||
2690           (R->order[j_in_R]==ringorder_rp) ||
2691           (R->order[j_in_R]==ringorder_rs) ||
2692           (R->order[j_in_R]==ringorder_lp) ||
2693           (R->order[j_in_R]==ringorder_ls))
2694          &&
2695            R->block0[j_in_R] <= R->N)
2696      {
2697        R->block1[j_in_R] = R->N;
2698      }
2699      else
2700      {
2701        Werror("ordering incomplete: size (%d) should be %d",R->block1[j_in_R],R->N);
2702        return TRUE;
2703      }
2704    }
2705    if (R->block0[j_in_R]>R->N)
2706    {
2707      Werror("not enough variables (%d) for ordering block %d, scanned so far:",R->N,j_in_R+1);
2708      for(int ii=0;ii<=j_in_R;ii++)
2709        Werror("ord[%d]: %s from v%d to v%d",ii+1,rSimpleOrdStr(R->order[ii]),R->block0[ii],R->block1[ii]);
2710      return TRUE;
2711    }
2712    if (check_comp)
2713    {
2714      BOOLEAN comp_order=FALSE;
2715      int jj;
2716      for(jj=0;jj<n;jj++)
2717      {
2718        if ((R->order[jj]==ringorder_c) ||
2719            (R->order[jj]==ringorder_C)) { comp_order=TRUE; break; }
2720      }
2721      if (!comp_order)
2722      {
2723        R->order=(rRingOrder_t*)omRealloc0Size(R->order,n*sizeof(rRingOrder_t),(n+1)*sizeof(rRingOrder_t));
2724        R->block0=(int*)omRealloc0Size(R->block0,n*sizeof(int),(n+1)*sizeof(int));
2725        R->block1=(int*)omRealloc0Size(R->block1,n*sizeof(int),(n+1)*sizeof(int));
2726        R->wvhdl=(int**)omRealloc0Size(R->wvhdl,n*sizeof(int_ptr),(n+1)*sizeof(int_ptr));
2727        R->order[n-1]=ringorder_C;
2728        R->block0[n-1]=0;
2729        R->block1[n-1]=0;
2730        R->wvhdl[n-1]=NULL;
2731        n++;
2732      }
2733    }
2734  }
2735  else
2736  {
2737    WerrorS("ordering must be given as `list`");
2738    return TRUE;
2739  }
2740  if (bitmask!=0) R->bitmask=bitmask*2;
2741  return FALSE;
2742}
2743
2744ring rCompose(const lists  L, const BOOLEAN check_comp, const long bitmask,const int isLetterplace)
2745{
2746  if ((L->nr!=3)
2747#ifdef HAVE_PLURAL
2748  &&(L->nr!=5)
2749#endif
2750  )
2751    return NULL;
2752  int is_gf_char=0;
2753  // 0: char/ cf - ring
2754  // 1: list (var)
2755  // 2: list (ord)
2756  // 3: qideal
2757  // possibly:
2758  // 4: C
2759  // 5: D
2760
2761  ring R = (ring) omAlloc0Bin(sip_sring_bin);
2762
2763  // ------------------------------------------------------------------
2764  // 0: char:
2765  if (L->m[0].Typ()==CRING_CMD)
2766  {
2767    R->cf=(coeffs)L->m[0].Data();
2768    R->cf->ref++;
2769  }
2770  else if (L->m[0].Typ()==INT_CMD)
2771  {
2772    int ch = (int)(long)L->m[0].Data();
2773    assume( ch >= 0 );
2774
2775    if (ch == 0) // Q?
2776      R->cf = nInitChar(n_Q, NULL);
2777    else
2778    {
2779      int l = IsPrime(ch); // Zp?
2780      if( l != ch )
2781      {
2782        Warn("%d is invalid characteristic of ground field. %d is used.", ch, l);
2783        ch = l;
2784      }
2785      #ifndef TEST_ZN_AS_ZP
2786      R->cf = nInitChar(n_Zp, (void*)(long)ch);
2787      #else
2788      mpz_t modBase;
2789      mpz_init_set_ui(modBase,(long) ch);
2790      ZnmInfo info;
2791      info.base= modBase;
2792      info.exp= 1;
2793      R->cf=nInitChar(n_Zn,(void*) &info); //exponent is missing
2794      R->cf->is_field=1;
2795      R->cf->is_domain=1;
2796      R->cf->has_simple_Inverse=1;
2797      #endif
2798    }
2799  }
2800  else if (L->m[0].Typ()==LIST_CMD) // something complicated...
2801  {
2802    lists LL=(lists)L->m[0].Data();
2803
2804#ifdef HAVE_RINGS
2805    if (LL->m[0].Typ() == STRING_CMD) // 1st comes a string?
2806    {
2807      rComposeRing(LL, R); // Ring!?
2808    }
2809    else
2810#endif
2811    if (LL->nr < 3)
2812      rComposeC(LL,R); // R, long_R, long_C
2813    else
2814    {
2815      if (LL->m[0].Typ()==INT_CMD)
2816      {
2817        int ch = (int)(long)LL->m[0].Data();
2818        while ((ch!=fftable[is_gf_char]) && (fftable[is_gf_char])) is_gf_char++;
2819        if (fftable[is_gf_char]==0) is_gf_char=-1;
2820
2821        if(is_gf_char!= -1)
2822        {
2823          GFInfo param;
2824
2825          param.GFChar = ch;
2826          param.GFDegree = 1;
2827          param.GFPar_name = (const char*)(((lists)(LL->m[1].Data()))->m[0].Data());
2828
2829          // nfInitChar should be able to handle the case when ch is in fftables!
2830          R->cf = nInitChar(n_GF, (void*)&param);
2831        }
2832      }
2833
2834      if( R->cf == NULL )
2835      {
2836        ring extRing = rCompose((lists)L->m[0].Data(),FALSE,0x7fff);
2837
2838        if (extRing==NULL)
2839        {
2840          WerrorS("could not create the specified coefficient field");
2841          goto rCompose_err;
2842        }
2843
2844        if( extRing->qideal != NULL ) // Algebraic extension
2845        {
2846          AlgExtInfo extParam;
2847
2848          extParam.r = extRing;
2849
2850          R->cf = nInitChar(n_algExt, (void*)&extParam);
2851        }
2852        else // Transcendental extension
2853        {
2854          TransExtInfo extParam;
2855          extParam.r = extRing;
2856          assume( extRing->qideal == NULL );
2857
2858          R->cf = nInitChar(n_transExt, &extParam);
2859        }
2860      }
2861    }
2862  }
2863  else
2864  {
2865    WerrorS("coefficient field must be described by `int` or `list`");
2866    goto rCompose_err;
2867  }
2868
2869  if( R->cf == NULL )
2870  {
2871    WerrorS("could not create coefficient field described by the input!");
2872    goto rCompose_err;
2873  }
2874
2875  // ------------------------- VARS ---------------------------
2876  if (rComposeVar(L,R)) goto rCompose_err;
2877  // ------------------------ ORDER ------------------------------
2878  if (rComposeOrder(L,check_comp,R)) goto rCompose_err;
2879
2880  // ------------------------ ??????? --------------------
2881
2882  if (!isLetterplace) rRenameVars(R);
2883  #ifdef HAVE_SHIFTBBA
2884  else
2885  {
2886    R->isLPring=isLetterplace;
2887    R->ShortOut=FALSE;
2888    R->CanShortOut=FALSE;
2889  }
2890  #endif
2891  if (bitmask!=0x7fff) R->bitmask=bitmask*2;
2892  rComplete(R);
2893
2894  // ------------------------ Q-IDEAL ------------------------
2895
2896  if (L->m[3].Typ()==IDEAL_CMD)
2897  {
2898    ideal q=(ideal)L->m[3].Data();
2899    if (q->m[0]!=NULL)
2900    {
2901      if (R->cf != currRing->cf) //->cf->ch!=currRing->cf->ch)
2902      {
2903      #if 0
2904            WerrorS("coefficient fields must be equal if q-ideal !=0");
2905            goto rCompose_err;
2906      #else
2907        ring orig_ring=currRing;
2908        rChangeCurrRing(R);
2909        int *perm=NULL;
2910        int *par_perm=NULL;
2911        int par_perm_size=0;
2912        nMapFunc nMap;
2913
2914        if ((nMap=nSetMap(orig_ring->cf))==NULL)
2915        {
2916          if (rEqual(orig_ring,currRing))
2917          {
2918            nMap=n_SetMap(currRing->cf, currRing->cf);
2919          }
2920          else
2921          // Allow imap/fetch to be make an exception only for:
2922          if ( (rField_is_Q_a(orig_ring) &&  // Q(a..) -> Q(a..) || Q || Zp || Zp(a)
2923            (rField_is_Q(currRing) || rField_is_Q_a(currRing) ||
2924             rField_is_Zp(currRing) || rField_is_Zp_a(currRing) ||
2925             rField_is_Zn(currRing)))
2926           ||
2927           (rField_is_Zp_a(orig_ring) &&  // Zp(a..) -> Zp(a..) || Zp
2928            (rField_is_Zp(currRing, rInternalChar(orig_ring)) ||
2929             rField_is_Zp_a(currRing, rInternalChar(orig_ring)))) )
2930          {
2931            par_perm_size=rPar(orig_ring);
2932
2933//            if ((orig_ring->minpoly != NULL) || (orig_ring->qideal != NULL))
2934//              naSetChar(rInternalChar(orig_ring),orig_ring);
2935//            else ntSetChar(rInternalChar(orig_ring),orig_ring);
2936
2937            nSetChar(currRing->cf);
2938          }
2939          else
2940          {
2941            WerrorS("coefficient fields must be equal if q-ideal !=0");
2942            goto rCompose_err;
2943          }
2944        }
2945        perm=(int *)omAlloc0((orig_ring->N+1)*sizeof(int));
2946        if (par_perm_size!=0)
2947          par_perm=(int *)omAlloc0(par_perm_size*sizeof(int));
2948        int i;
2949        #if 0
2950        // use imap:
2951        maFindPerm(orig_ring->names,orig_ring->N,orig_ring->parameter,orig_ring->P,
2952          currRing->names,currRing->N,currRing->parameter, currRing->P,
2953          perm,par_perm, currRing->ch);
2954        #else
2955        // use fetch
2956        if ((rPar(orig_ring)>0) && (rPar(currRing)==0))
2957        {
2958          for(i=si_min(rPar(orig_ring),rVar(currRing))-1;i>=0;i--) par_perm[i]=i+1;
2959        }
2960        else if (par_perm_size!=0)
2961          for(i=si_min(rPar(orig_ring),rPar(currRing))-1;i>=0;i--) par_perm[i]=-(i+1);
2962        for(i=si_min(orig_ring->N,rVar(currRing));i>0;i--) perm[i]=i;
2963        #endif
2964        ideal dest_id=idInit(IDELEMS(q),1);
2965        for(i=IDELEMS(q)-1; i>=0; i--)
2966        {
2967          dest_id->m[i]=p_PermPoly(q->m[i],perm,orig_ring, currRing,nMap,
2968                                  par_perm,par_perm_size);
2969          //  PrintS("map:");pWrite(dest_id->m[i]);PrintLn();
2970          pTest(dest_id->m[i]);
2971        }
2972        R->qideal=dest_id;
2973        if (perm!=NULL)
2974          omFreeSize((ADDRESS)perm,(orig_ring->N+1)*sizeof(int));
2975        if (par_perm!=NULL)
2976          omFreeSize((ADDRESS)par_perm,par_perm_size*sizeof(int));
2977        rChangeCurrRing(orig_ring);
2978      #endif
2979      }
2980      else
2981        R->qideal=idrCopyR(q,currRing,R);
2982    }
2983  }
2984  else
2985  {
2986    WerrorS("q-ideal must be given as `ideal`");
2987    goto rCompose_err;
2988  }
2989
2990
2991  // ---------------------------------------------------------------
2992  #ifdef HAVE_PLURAL
2993  if (L->nr==5)
2994  {
2995    if (nc_CallPlural((matrix)L->m[4].Data(),
2996                      (matrix)L->m[5].Data(),
2997                      NULL,NULL,
2998                      R,
2999                      true, // !!!
3000                      true, false,
3001                      currRing, FALSE)) goto rCompose_err;
3002    // takes care about non-comm. quotient! i.e. calls "nc_SetupQuotient" due to last true
3003  }
3004  #endif
3005  return R;
3006
3007rCompose_err:
3008  if (R->N>0)
3009  {
3010    int i;
3011    if (R->names!=NULL)
3012    {
3013      i=R->N-1;
3014      while (i>=0) { omfree(R->names[i]); i--; }
3015      omFree(R->names);
3016    }
3017  }
3018  omfree(R->order);
3019  omfree(R->block0);
3020  omfree(R->block1);
3021  omfree(R->wvhdl);
3022  omFree(R);
3023  return NULL;
3024}
3025
3026// from matpol.cc
3027
3028/*2
3029* compute the jacobi matrix of an ideal
3030*/
3031BOOLEAN mpJacobi(leftv res,leftv a)
3032{
3033  int     i,j;
3034  matrix result;
3035  ideal id=(ideal)a->Data();
3036
3037  result =mpNew(IDELEMS(id),rVar(currRing));
3038  for (i=1; i<=IDELEMS(id); i++)
3039  {
3040    for (j=1; j<=rVar(currRing); j++)
3041    {
3042      MATELEM(result,i,j) = pDiff(id->m[i-1],j);
3043    }
3044  }
3045  res->data=(char *)result;
3046  return FALSE;
3047}
3048
3049/*2
3050* returns the Koszul-matrix of degree d of a vectorspace with dimension n
3051* uses the first n entrees of id, if id <> NULL
3052*/
3053BOOLEAN mpKoszul(leftv res,leftv c/*ip*/, leftv b/*in*/, leftv id)
3054{
3055  int n=(int)(long)b->Data();
3056  int d=(int)(long)c->Data();
3057  int     k,l,sign,row,col;
3058  matrix  result;
3059  ideal temp;
3060  BOOLEAN bo;
3061  poly    p;
3062
3063  if ((d>n) || (d<1) || (n<1))
3064  {
3065    res->data=(char *)mpNew(1,1);
3066    return FALSE;
3067  }
3068  int *choise = (int*)omAlloc(d*sizeof(int));
3069  if (id==NULL)
3070    temp=idMaxIdeal(1);
3071  else
3072    temp=(ideal)id->Data();
3073
3074  k = binom(n,d);
3075  l = k*d;
3076  l /= n-d+1;
3077  result =mpNew(l,k);
3078  col = 1;
3079  idInitChoise(d,1,n,&bo,choise);
3080  while (!bo)
3081  {
3082    sign = 1;
3083    for (l=1;l<=d;l++)
3084    {
3085      if (choise[l-1]<=IDELEMS(temp))
3086      {
3087        p = pCopy(temp->m[choise[l-1]-1]);
3088        if (sign == -1) p = pNeg(p);
3089        sign *= -1;
3090        row = idGetNumberOfChoise(l-1,d,1,n,choise);
3091        MATELEM(result,row,col) = p;
3092      }
3093    }
3094    col++;
3095    idGetNextChoise(d,n,&bo,choise);
3096  }
3097  omFreeSize(choise,d*sizeof(int));
3098  if (id==NULL) idDelete(&temp);
3099
3100  res->data=(char *)result;
3101  return FALSE;
3102}
3103
3104// from syz1.cc
3105/*2
3106* read out the Betti numbers from resolution
3107* (interpreter interface)
3108*/
3109BOOLEAN syBetti2(leftv res, leftv u, leftv w)
3110{
3111  syStrategy syzstr=(syStrategy)u->Data();
3112
3113  BOOLEAN minim=(int)(long)w->Data();
3114  int row_shift=0;
3115  int add_row_shift=0;
3116  intvec *weights=NULL;
3117  intvec *ww=(intvec *)atGet(u,"isHomog",INTVEC_CMD);
3118  if (ww!=NULL)
3119  {
3120     weights=ivCopy(ww);
3121     add_row_shift = ww->min_in();
3122     (*weights) -= add_row_shift;
3123  }
3124
3125  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift,weights);
3126  //row_shift += add_row_shift;
3127  //Print("row_shift=%d, add_row_shift=%d\n",row_shift,add_row_shift);
3128  atSet(res,omStrDup("rowShift"),(void*)(long)add_row_shift,INT_CMD);
3129
3130  return FALSE;
3131}
3132BOOLEAN syBetti1(leftv res, leftv u)
3133{
3134  sleftv tmp;
3135  memset(&tmp,0,sizeof(tmp));
3136  tmp.rtyp=INT_CMD;
3137  tmp.data=(void *)1;
3138  return syBetti2(res,u,&tmp);
3139}
3140
3141/*3
3142* converts a resolution into a list of modules
3143*/
3144lists syConvRes(syStrategy syzstr,BOOLEAN toDel,int add_row_shift)
3145{
3146  resolvente fullres = syzstr->fullres;
3147  resolvente minres = syzstr->minres;
3148
3149  const int length = syzstr->length;
3150
3151  if ((fullres==NULL) && (minres==NULL))
3152  {
3153    if (syzstr->hilb_coeffs==NULL)
3154    { // La Scala
3155      fullres = syReorder(syzstr->res, length, syzstr);
3156    }
3157    else
3158    { // HRES
3159      minres = syReorder(syzstr->orderedRes, length, syzstr);
3160      syKillEmptyEntres(minres, length);
3161    }
3162  }
3163
3164  resolvente tr;
3165  int typ0=IDEAL_CMD;
3166
3167  if (minres!=NULL)
3168    tr = minres;
3169  else
3170    tr = fullres;
3171
3172  resolvente trueres=NULL;
3173  intvec ** w=NULL;
3174
3175  if (length>0)
3176  {
3177    trueres = (resolvente)omAlloc0((length)*sizeof(ideal));
3178    for (int i=length-1;i>=0;i--)
3179    {
3180      if (tr[i]!=NULL)
3181      {
3182        trueres[i] = idCopy(tr[i]);
3183      }
3184    }
3185    if ( id_RankFreeModule(trueres[0], currRing) > 0)
3186      typ0 = MODUL_CMD;
3187    if (syzstr->weights!=NULL)
3188    {
3189      w = (intvec**)omAlloc0(length*sizeof(intvec*));
3190      for (int i=length-1;i>=0;i--)
3191      {
3192        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
3193      }
3194    }
3195  }
3196
3197  lists li = liMakeResolv(trueres, length, syzstr->list_length,typ0,
3198                          w, add_row_shift);
3199
3200  if (toDel)
3201    syKillComputation(syzstr);
3202  else
3203  {
3204    if( fullres != NULL && syzstr->fullres == NULL )
3205      syzstr->fullres = fullres;
3206
3207    if( minres != NULL && syzstr->minres == NULL )
3208      syzstr->minres = minres;
3209  }
3210  return li;
3211}
3212
3213/*3
3214* converts a list of modules into a resolution
3215*/
3216syStrategy syConvList(lists li)
3217{
3218  int typ0;
3219  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
3220
3221  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
3222  if (fr != NULL)
3223  {
3224
3225    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
3226    for (int i=result->length-1;i>=0;i--)
3227    {
3228      if (fr[i]!=NULL)
3229        result->fullres[i] = idCopy(fr[i]);
3230    }
3231    result->list_length=result->length;
3232    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
3233  }
3234  else
3235  {
3236    omFreeSize(result, sizeof(ssyStrategy));
3237    result = NULL;
3238  }
3239  return result;
3240}
3241
3242/*3
3243* converts a list of modules into a minimal resolution
3244*/
3245syStrategy syForceMin(lists li)
3246{
3247  int typ0;
3248  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
3249
3250  resolvente fr = liFindRes(li,&(result->length),&typ0);
3251  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
3252  for (int i=result->length-1;i>=0;i--)
3253  {
3254    if (fr[i]!=NULL)
3255      result->minres[i] = idCopy(fr[i]);
3256  }
3257  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
3258  return result;
3259}
3260// from weight.cc
3261BOOLEAN kWeight(leftv res,leftv id)
3262{
3263  ideal F=(ideal)id->Data();
3264  intvec * iv = new intvec(rVar(currRing));
3265  polyset s;
3266  int  sl, n, i;
3267  int  *x;
3268
3269  res->data=(char *)iv;
3270  s = F->m;
3271  sl = IDELEMS(F) - 1;
3272  n = rVar(currRing);
3273  double wNsqr = (double)2.0 / (double)n;
3274  wFunctional = wFunctionalBuch;
3275  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
3276  wCall(s, sl, x, wNsqr, currRing);
3277  for (i = n; i!=0; i--)
3278    (*iv)[i-1] = x[i + n + 1];
3279  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
3280  return FALSE;
3281}
3282
3283BOOLEAN kQHWeight(leftv res,leftv v)
3284{
3285  res->data=(char *)id_QHomWeight((ideal)v->Data(), currRing);
3286  if (res->data==NULL)
3287    res->data=(char *)new intvec(rVar(currRing));
3288  return FALSE;
3289}
3290/*==============================================================*/
3291// from clapsing.cc
3292#if 0
3293BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
3294{
3295  BOOLEAN b=singclap_factorize((poly)(u->CopyD()), &v, 0);
3296  res->data=(void *)b;
3297}
3298#endif
3299
3300BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
3301{
3302  res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
3303                  (poly)w->CopyD(), currRing);
3304  return errorreported;
3305}
3306
3307BOOLEAN jjCHARSERIES(leftv res, leftv u)
3308{
3309  res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
3310  return (res->data==NULL);
3311}
3312
3313// from semic.cc
3314#ifdef HAVE_SPECTRUM
3315
3316// ----------------------------------------------------------------------------
3317//  Initialize a  spectrum  deep from a  singular  lists
3318// ----------------------------------------------------------------------------
3319
3320void copy_deep( spectrum& spec, lists l )
3321{
3322    spec.mu = (int)(long)(l->m[0].Data( ));
3323    spec.pg = (int)(long)(l->m[1].Data( ));
3324    spec.n  = (int)(long)(l->m[2].Data( ));
3325
3326    spec.copy_new( spec.n );
3327
3328    intvec  *num = (intvec*)l->m[3].Data( );
3329    intvec  *den = (intvec*)l->m[4].Data( );
3330    intvec  *mul = (intvec*)l->m[5].Data( );
3331
3332    for( int i=0; i<spec.n; i++ )
3333    {
3334        spec.s[i] = (Rational)((*num)[i])/(Rational)((*den)[i]);
3335        spec.w[i] = (*mul)[i];
3336    }
3337}
3338
3339// ----------------------------------------------------------------------------
3340//  singular lists  constructor for  spectrum
3341// ----------------------------------------------------------------------------
3342
3343spectrum /*former spectrum::spectrum ( lists l )*/
3344spectrumFromList( lists l )
3345{
3346    spectrum result;
3347    copy_deep( result, l );
3348    return result;
3349}
3350
3351// ----------------------------------------------------------------------------
3352//  generate a Singular  lists  from a spectrum
3353// ----------------------------------------------------------------------------
3354
3355/* former spectrum::thelist ( void )*/
3356lists   getList( spectrum& spec )
3357{
3358    lists   L  = (lists)omAllocBin( slists_bin);
3359
3360    L->Init( 6 );
3361
3362    intvec            *num  = new intvec( spec.n );
3363    intvec            *den  = new intvec( spec.n );
3364    intvec            *mult = new intvec( spec.n );
3365
3366    for( int i=0; i<spec.n; i++ )
3367    {
3368        (*num) [i] = spec.s[i].get_num_si( );
3369        (*den) [i] = spec.s[i].get_den_si( );
3370        (*mult)[i] = spec.w[i];
3371    }
3372
3373    L->m[0].rtyp = INT_CMD;    //  milnor number
3374    L->m[1].rtyp = INT_CMD;    //  geometrical genus
3375    L->m[2].rtyp = INT_CMD;    //  # of spectrum numbers
3376    L->m[3].rtyp = INTVEC_CMD; //  numerators
3377    L->m[4].rtyp = INTVEC_CMD; //  denomiantors
3378    L->m[5].rtyp = INTVEC_CMD; //  multiplicities
3379
3380    L->m[0].data = (void*)(long)spec.mu;
3381    L->m[1].data = (void*)(long)spec.pg;
3382    L->m[2].data = (void*)(long)spec.n;
3383    L->m[3].data = (void*)num;
3384    L->m[4].data = (void*)den;
3385    L->m[5].data = (void*)mult;
3386
3387    return  L;
3388}
3389// from spectrum.cc
3390// ----------------------------------------------------------------------------
3391//  print out an error message for a spectrum list
3392// ----------------------------------------------------------------------------
3393
3394typedef enum
3395{
3396    semicOK,
3397    semicMulNegative,
3398
3399    semicListTooShort,
3400    semicListTooLong,
3401
3402    semicListFirstElementWrongType,
3403    semicListSecondElementWrongType,
3404    semicListThirdElementWrongType,
3405    semicListFourthElementWrongType,
3406    semicListFifthElementWrongType,
3407    semicListSixthElementWrongType,
3408
3409    semicListNNegative,
3410    semicListWrongNumberOfNumerators,
3411    semicListWrongNumberOfDenominators,
3412    semicListWrongNumberOfMultiplicities,
3413
3414    semicListMuNegative,
3415    semicListPgNegative,
3416    semicListNumNegative,
3417    semicListDenNegative,
3418    semicListMulNegative,
3419
3420    semicListNotSymmetric,
3421    semicListNotMonotonous,
3422
3423    semicListMilnorWrong,
3424    semicListPGWrong
3425
3426} semicState;
3427
3428void    list_error( semicState state )
3429{
3430    switch( state )
3431    {
3432        case semicListTooShort:
3433            WerrorS( "the list is too short" );
3434            break;
3435        case semicListTooLong:
3436            WerrorS( "the list is too long" );
3437            break;
3438
3439        case semicListFirstElementWrongType:
3440            WerrorS( "first element of the list should be int" );
3441            break;
3442        case semicListSecondElementWrongType:
3443            WerrorS( "second element of the list should be int" );
3444            break;
3445        case semicListThirdElementWrongType:
3446            WerrorS( "third element of the list should be int" );
3447            break;
3448        case semicListFourthElementWrongType:
3449            WerrorS( "fourth element of the list should be intvec" );
3450            break;
3451        case semicListFifthElementWrongType:
3452            WerrorS( "fifth element of the list should be intvec" );
3453            break;
3454        case semicListSixthElementWrongType:
3455            WerrorS( "sixth element of the list should be intvec" );
3456            break;
3457
3458        case semicListNNegative:
3459            WerrorS( "first element of the list should be positive" );
3460            break;
3461        case semicListWrongNumberOfNumerators:
3462            WerrorS( "wrong number of numerators" );
3463            break;
3464        case semicListWrongNumberOfDenominators:
3465            WerrorS( "wrong number of denominators" );
3466            break;
3467        case semicListWrongNumberOfMultiplicities:
3468            WerrorS( "wrong number of multiplicities" );
3469            break;
3470
3471        case semicListMuNegative:
3472            WerrorS( "the Milnor number should be positive" );
3473            break;
3474        case semicListPgNegative:
3475            WerrorS( "the geometrical genus should be nonnegative" );
3476            break;
3477        case semicListNumNegative:
3478            WerrorS( "all numerators should be positive" );
3479            break;
3480        case semicListDenNegative:
3481            WerrorS( "all denominators should be positive" );
3482            break;
3483        case semicListMulNegative:
3484            WerrorS( "all multiplicities should be positive" );
3485            break;
3486
3487        case semicListNotSymmetric:
3488            WerrorS( "it is not symmetric" );
3489            break;
3490        case semicListNotMonotonous:
3491            WerrorS( "it is not monotonous" );
3492            break;
3493
3494        case semicListMilnorWrong:
3495            WerrorS( "the Milnor number is wrong" );
3496            break;
3497        case semicListPGWrong:
3498            WerrorS( "the geometrical genus is wrong" );
3499            break;
3500
3501        default:
3502            WerrorS( "unspecific error" );
3503            break;
3504    }
3505}
3506// ----------------------------------------------------------------------------
3507//  this is the main spectrum computation function
3508// ----------------------------------------------------------------------------
3509
3510enum    spectrumState
3511{
3512    spectrumOK,
3513    spectrumZero,
3514    spectrumBadPoly,
3515    spectrumNoSingularity,
3516    spectrumNotIsolated,
3517    spectrumDegenerate,
3518    spectrumWrongRing,
3519    spectrumNoHC,
3520    spectrumUnspecErr
3521};
3522
3523// from splist.cc
3524// ----------------------------------------------------------------------------
3525//  Compute the spectrum of a  spectrumPolyList
3526// ----------------------------------------------------------------------------
3527
3528/* former spectrumPolyList::spectrum ( lists*, int) */
3529spectrumState   spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3530{
3531  spectrumPolyNode  **node = &speclist.root;
3532  spectrumPolyNode  *search;
3533
3534  poly              f,tmp;
3535  int               found,cmp;
3536
3537  Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3538                 ( fast==2 ? 2 : 1 ) );
3539
3540  Rational weight_prev( 0,1 );
3541
3542  int     mu = 0;          // the milnor number
3543  int     pg = 0;          // the geometrical genus
3544  int     n  = 0;          // number of different spectral numbers
3545  int     z  = 0;          // number of spectral number equal to smax
3546
3547  while( (*node)!=(spectrumPolyNode*)NULL &&
3548         ( fast==0 || (*node)->weight<=smax ) )
3549  {
3550        // ---------------------------------------
3551        //  determine the first normal form which
3552        //  contains the monomial  node->mon
3553        // ---------------------------------------
3554
3555    found  = FALSE;
3556    search = *node;
3557
3558    while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3559    {
3560      if( search->nf!=(poly)NULL )
3561      {
3562        f = search->nf;
3563
3564        do
3565        {
3566                    // --------------------------------
3567                    //  look for  (*node)->mon  in   f
3568                    // --------------------------------
3569
3570          cmp = pCmp( (*node)->mon,f );
3571
3572          if( cmp<0 )
3573          {
3574            f = pNext( f );
3575          }
3576          else if( cmp==0 )
3577          {
3578                        // -----------------------------
3579                        //  we have found a normal form
3580                        // -----------------------------
3581
3582            found = TRUE;
3583
3584                        //  normalize coefficient
3585
3586            number inv = nInvers( pGetCoeff( f ) );
3587            search->nf=__p_Mult_nn( search->nf,inv,currRing );
3588            nDelete( &inv );
3589
3590                        //  exchange  normal forms
3591
3592            tmp         = (*node)->nf;
3593            (*node)->nf = search->nf;
3594            search->nf  = tmp;
3595          }
3596        }
3597        while( cmp<0 && f!=(poly)NULL );
3598      }
3599      search = search->next;
3600    }
3601
3602    if( found==FALSE )
3603    {
3604            // ------------------------------------------------
3605            //  the weight of  node->mon  is a spectrum number
3606            // ------------------------------------------------
3607
3608      mu++;
3609
3610      if( (*node)->weight<=(Rational)1 )              pg++;
3611      if( (*node)->weight==smax )           z++;
3612      if( (*node)->weight>weight_prev )     n++;
3613
3614      weight_prev = (*node)->weight;
3615      node = &((*node)->next);
3616    }
3617    else
3618    {
3619            // -----------------------------------------------
3620            //  determine all other normal form which contain
3621            //  the monomial  node->mon
3622            //  replace for  node->mon  its normal form
3623            // -----------------------------------------------
3624
3625      while( search!=(spectrumPolyNode*)NULL )
3626      {
3627        if( search->nf!=(poly)NULL )
3628        {
3629          f = search->nf;
3630
3631          do
3632          {
3633                        // --------------------------------
3634                        //  look for  (*node)->mon  in   f
3635                        // --------------------------------
3636
3637            cmp = pCmp( (*node)->mon,f );
3638
3639            if( cmp<0 )
3640            {
3641              f = pNext( f );
3642            }
3643            else if( cmp==0 )
3644            {
3645              search->nf = pSub( search->nf,
3646                                 __pp_Mult_nn( (*node)->nf,pGetCoeff( f ),currRing ) );
3647              pNorm( search->nf );
3648            }
3649          }
3650          while( cmp<0 && f!=(poly)NULL );
3651        }
3652        search = search->next;
3653      }
3654      speclist.delete_node( node );
3655    }
3656
3657  }
3658
3659    // --------------------------------------------------------
3660    //  fast computation exploits the symmetry of the spectrum
3661    // --------------------------------------------------------
3662
3663  if( fast==2 )
3664  {
3665    mu = 2*mu - z;
3666    n  = ( z > 0 ? 2*n - 1 : 2*n );
3667  }
3668
3669    // --------------------------------------------------------
3670    //  compute the spectrum numbers with their multiplicities
3671    // --------------------------------------------------------
3672
3673  intvec            *nom  = new intvec( n );
3674  intvec            *den  = new intvec( n );
3675  intvec            *mult = new intvec( n );
3676
3677  int count         = 0;
3678  int multiplicity  = 1;
3679
3680  for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3681              ( fast==0 || search->weight<=smax );
3682       search=search->next )
3683  {
3684    if( search->next==(spectrumPolyNode*)NULL ||
3685        search->weight<search->next->weight )
3686    {
3687      (*nom) [count] = search->weight.get_num_si( );
3688      (*den) [count] = search->weight.get_den_si( );
3689      (*mult)[count] = multiplicity;
3690
3691      multiplicity=1;
3692      count++;
3693    }
3694    else
3695    {
3696      multiplicity++;
3697    }
3698  }
3699
3700    // --------------------------------------------------------
3701    //  fast computation exploits the symmetry of the spectrum
3702    // --------------------------------------------------------
3703
3704  if( fast==2 )
3705  {
3706    int n1,n2;
3707    for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3708    {
3709      (*nom) [n2] = rVar(currRing)*(*den)[n1]-(*nom)[n1];
3710      (*den) [n2] = (*den)[n1];
3711      (*mult)[n2] = (*mult)[n1];
3712    }
3713  }
3714
3715    // -----------------------------------
3716    //  test if the spectrum is symmetric
3717    // -----------------------------------
3718
3719  if( fast==0 || fast==1 )
3720  {
3721    int symmetric=TRUE;
3722
3723    for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3724    {
3725      if( (*mult)[n1]!=(*mult)[n2] ||
3726          (*den) [n1]!= (*den)[n2] ||
3727          (*nom)[n1]+(*nom)[n2]!=rVar(currRing)*(*den) [n1] )
3728      {
3729        symmetric = FALSE;
3730      }
3731    }
3732
3733    if( symmetric==FALSE )
3734    {
3735            // ---------------------------------------------
3736            //  the spectrum is not symmetric => degenerate
3737            //  principal part
3738            // ---------------------------------------------
3739
3740      *L = (lists)omAllocBin( slists_bin);
3741      (*L)->Init( 1 );
3742      (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3743      (*L)->m[0].data = (void*)(long)mu;
3744
3745      return spectrumDegenerate;
3746    }
3747  }
3748
3749  *L = (lists)omAllocBin( slists_bin);
3750
3751  (*L)->Init( 6 );
3752
3753  (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3754  (*L)->m[1].rtyp = INT_CMD;    //  geometrical genus
3755  (*L)->m[2].rtyp = INT_CMD;    //  number of spectrum values
3756  (*L)->m[3].rtyp = INTVEC_CMD; //  nominators
3757  (*L)->m[4].rtyp = INTVEC_CMD; //  denomiantors
3758  (*L)->m[5].rtyp = INTVEC_CMD; //  multiplicities
3759
3760  (*L)->m[0].data = (void*)(long)mu;
3761  (*L)->m[1].data = (void*)(long)pg;
3762  (*L)->m[2].data = (void*)(long)n;
3763  (*L)->m[3].data = (void*)nom;
3764  (*L)->m[4].data = (void*)den;
3765  (*L)->m[5].data = (void*)mult;
3766
3767  return  spectrumOK;
3768}
3769
3770spectrumState   spectrumCompute( poly h,lists *L,int fast )
3771{
3772  int i;
3773
3774  #ifdef SPECTRUM_DEBUG
3775  #ifdef SPECTRUM_PRINT
3776  #ifdef SPECTRUM_IOSTREAM
3777    cout << "spectrumCompute\n";
3778    if( fast==0 ) cout << "    no optimization" << endl;
3779    if( fast==1 ) cout << "    weight optimization" << endl;
3780    if( fast==2 ) cout << "    symmetry optimization" << endl;
3781  #else
3782    fputs( "spectrumCompute\n",stdout );
3783    if( fast==0 ) fputs( "    no optimization\n", stdout );
3784    if( fast==1 ) fputs( "    weight optimization\n", stdout );
3785    if( fast==2 ) fputs( "    symmetry optimization\n", stdout );
3786  #endif
3787  #endif
3788  #endif
3789
3790  // ----------------------
3791  //  check if  h  is zero
3792  // ----------------------
3793
3794  if( h==(poly)NULL )
3795  {
3796    return  spectrumZero;
3797  }
3798
3799  // ----------------------------------
3800  //  check if  h  has a constant term
3801  // ----------------------------------
3802
3803  if( hasConstTerm( h, currRing ) )
3804  {
3805    return  spectrumBadPoly;
3806  }
3807
3808  // --------------------------------
3809  //  check if  h  has a linear term
3810  // --------------------------------
3811
3812  if( hasLinearTerm( h, currRing ) )
3813  {
3814    *L = (lists)omAllocBin( slists_bin);
3815    (*L)->Init( 1 );
3816    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3817    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3818
3819    return  spectrumNoSingularity;
3820  }
3821
3822  // ----------------------------------
3823  //  compute the jacobi ideal of  (h)
3824  // ----------------------------------
3825
3826  ideal J = NULL;
3827  J = idInit( rVar(currRing),1 );
3828
3829  #ifdef SPECTRUM_DEBUG
3830  #ifdef SPECTRUM_PRINT
3831  #ifdef SPECTRUM_IOSTREAM
3832    cout << "\n   computing the Jacobi ideal...\n";
3833  #else
3834    fputs( "\n   computing the Jacobi ideal...\n",stdout );
3835  #endif
3836  #endif
3837  #endif
3838
3839  for( i=0; i<rVar(currRing); i++ )
3840  {
3841    J->m[i] = pDiff( h,i+1); //j );
3842
3843    #ifdef SPECTRUM_DEBUG
3844    #ifdef SPECTRUM_PRINT
3845    #ifdef SPECTRUM_IOSTREAM
3846      cout << "        ";
3847    #else
3848      fputs("        ", stdout );
3849    #endif
3850      pWrite( J->m[i] );
3851    #endif
3852    #endif
3853  }
3854
3855  // --------------------------------------------
3856  //  compute a standard basis  stdJ  of  jac(h)
3857  // --------------------------------------------
3858
3859  #ifdef SPECTRUM_DEBUG
3860  #ifdef SPECTRUM_PRINT
3861  #ifdef SPECTRUM_IOSTREAM
3862    cout << endl;
3863    cout << "    computing a standard basis..." << endl;
3864  #else
3865    fputs( "\n", stdout );
3866    fputs( "    computing a standard basis...\n", stdout );
3867  #endif
3868  #endif
3869  #endif
3870
3871  ideal stdJ = kStd(J,currRing->qideal,isNotHomog,NULL);
3872  idSkipZeroes( stdJ );
3873
3874  #ifdef SPECTRUM_DEBUG
3875  #ifdef SPECTRUM_PRINT
3876    for( i=0; i<IDELEMS(stdJ); i++ )
3877    {
3878      #ifdef SPECTRUM_IOSTREAM
3879        cout << "        ";
3880      #else
3881        fputs( "        ",stdout );
3882      #endif
3883
3884      pWrite( stdJ->m[i] );
3885    }
3886  #endif
3887  #endif
3888
3889  idDelete( &J );
3890
3891  // ------------------------------------------
3892  //  check if the  h  has a singularity
3893  // ------------------------------------------
3894
3895  if( hasOne( stdJ, currRing ) )
3896  {
3897    // -------------------------------
3898    //  h is smooth in the origin
3899    //  return only the Milnor number
3900    // -------------------------------
3901
3902    *L = (lists)omAllocBin( slists_bin);
3903    (*L)->Init( 1 );
3904    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3905    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3906
3907    return  spectrumNoSingularity;
3908  }
3909
3910  // ------------------------------------------
3911  //  check if the singularity  h  is isolated
3912  // ------------------------------------------
3913
3914  for( i=rVar(currRing); i>0; i-- )
3915  {
3916    if( hasAxis( stdJ,i, currRing )==FALSE )
3917    {
3918      return  spectrumNotIsolated;
3919    }
3920  }
3921
3922  // ------------------------------------------
3923  //  compute the highest corner  hc  of  stdJ
3924  // ------------------------------------------
3925
3926  #ifdef SPECTRUM_DEBUG
3927  #ifdef SPECTRUM_PRINT
3928  #ifdef SPECTRUM_IOSTREAM
3929    cout << "\n    computing the highest corner...\n";
3930  #else
3931    fputs( "\n    computing the highest corner...\n", stdout );
3932  #endif
3933  #endif
3934  #endif
3935
3936  poly hc = (poly)NULL;
3937
3938  scComputeHC( stdJ,currRing->qideal, 0,hc );
3939
3940  if( hc!=(poly)NULL )
3941  {
3942    pGetCoeff(hc) = nInit(1);
3943
3944    for( i=rVar(currRing); i>0; i-- )
3945    {
3946      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3947    }
3948    pSetm( hc );
3949  }
3950  else
3951  {
3952    return  spectrumNoHC;
3953  }
3954
3955  #ifdef SPECTRUM_DEBUG
3956  #ifdef SPECTRUM_PRINT
3957  #ifdef SPECTRUM_IOSTREAM
3958    cout << "       ";
3959  #else
3960    fputs( "       ", stdout );
3961  #endif
3962    pWrite( hc );
3963  #endif
3964  #endif
3965
3966  // ----------------------------------------
3967  //  compute the Newton polygon  nph  of  h
3968  // ----------------------------------------
3969
3970  #ifdef SPECTRUM_DEBUG
3971  #ifdef SPECTRUM_PRINT
3972  #ifdef SPECTRUM_IOSTREAM
3973    cout << "\n    computing the newton polygon...\n";
3974  #else
3975    fputs( "\n    computing the newton polygon...\n", stdout );
3976  #endif
3977  #endif
3978  #endif
3979
3980  newtonPolygon nph( h, currRing );
3981
3982  #ifdef SPECTRUM_DEBUG
3983  #ifdef SPECTRUM_PRINT
3984    cout << nph;
3985  #endif
3986  #endif
3987
3988  // -----------------------------------------------
3989  //  compute the weight corner  wc  of  (stdj,nph)
3990  // -----------------------------------------------
3991
3992  #ifdef SPECTRUM_DEBUG
3993  #ifdef SPECTRUM_PRINT
3994  #ifdef SPECTRUM_IOSTREAM
3995    cout << "\n    computing the weight corner...\n";
3996  #else
3997    fputs( "\n    computing the weight corner...\n", stdout );
3998  #endif
3999  #endif
4000  #endif
4001
4002  poly    wc = ( fast==0 ? pCopy( hc ) :
4003               ( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
4004              /* fast==2 */computeWC( nph,
4005                      ((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
4006
4007  #ifdef SPECTRUM_DEBUG
4008  #ifdef SPECTRUM_PRINT
4009  #ifdef SPECTRUM_IOSTREAM
4010    cout << "        ";
4011  #else
4012    fputs( "        ", stdout );
4013  #endif
4014    pWrite( wc );
4015  #endif
4016  #endif
4017
4018  // -------------
4019  //  compute  NF
4020  // -------------
4021
4022  #ifdef SPECTRUM_DEBUG
4023  #ifdef SPECTRUM_PRINT
4024  #ifdef SPECTRUM_IOSTREAM
4025    cout << "\n    computing NF...\n" << endl;
4026  #else
4027    fputs( "\n    computing NF...\n", stdout );
4028  #endif
4029  #endif
4030  #endif
4031
4032  spectrumPolyList NF( &nph );
4033
4034  computeNF( stdJ,hc,wc,&NF, currRing );
4035
4036  #ifdef SPECTRUM_DEBUG
4037  #ifdef SPECTRUM_PRINT
4038    cout << NF;
4039  #ifdef SPECTRUM_IOSTREAM
4040    cout << endl;
4041  #else
4042    fputs( "\n", stdout );
4043  #endif
4044  #endif
4045  #endif
4046
4047  // ----------------------------
4048  //  compute the spectrum of  h
4049  // ----------------------------
4050//  spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
4051
4052  return spectrumStateFromList(NF, L, fast );
4053}
4054
4055// ----------------------------------------------------------------------------
4056//  this procedure is called from the interpreter
4057// ----------------------------------------------------------------------------
4058//  first  = polynomial
4059//  result = list of spectrum numbers
4060// ----------------------------------------------------------------------------
4061
4062void spectrumPrintError(spectrumState state)
4063{
4064  switch( state )
4065  {
4066    case spectrumZero:
4067      WerrorS( "polynomial is zero" );
4068      break;
4069    case spectrumBadPoly:
4070      WerrorS( "polynomial has constant term" );
4071      break;
4072    case spectrumNoSingularity:
4073      WerrorS( "not a singularity" );
4074      break;
4075    case spectrumNotIsolated:
4076      WerrorS( "the singularity is not isolated" );
4077      break;
4078    case spectrumNoHC:
4079      WerrorS( "highest corner cannot be computed" );
4080      break;
4081    case spectrumDegenerate:
4082      WerrorS( "principal part is degenerate" );
4083      break;
4084    case spectrumOK:
4085      break;
4086
4087    default:
4088      WerrorS( "unknown error occurred" );
4089      break;
4090  }
4091}
4092
4093BOOLEAN spectrumProc( leftv result,leftv first )
4094{
4095  spectrumState state = spectrumOK;
4096
4097  // -------------------
4098  //  check consistency
4099  // -------------------
4100
4101  //  check for a local ring
4102
4103  if( !ringIsLocal(currRing ) )
4104  {
4105    WerrorS( "only works for local orderings" );
4106    state = spectrumWrongRing;
4107  }
4108
4109  //  no quotient rings are allowed
4110
4111  else if( currRing->qideal != NULL )
4112  {
4113    WerrorS( "does not work in quotient rings" );
4114    state = spectrumWrongRing;
4115  }
4116  else
4117  {
4118    lists   L    = (lists)NULL;
4119    int     flag = 1; // weight corner optimization is safe
4120
4121    state = spectrumCompute( (poly)first->Data( ),&L,flag );
4122
4123    if( state==spectrumOK )
4124    {
4125      result->rtyp = LIST_CMD;
4126      result->data = (char*)L;
4127    }
4128    else
4129    {
4130      spectrumPrintError(state);
4131    }
4132  }
4133
4134  return  (state!=spectrumOK);
4135}
4136
4137// ----------------------------------------------------------------------------
4138//  this procedure is called from the interpreter
4139// ----------------------------------------------------------------------------
4140//  first  = polynomial
4141//  result = list of spectrum numbers
4142// ----------------------------------------------------------------------------
4143
4144BOOLEAN spectrumfProc( leftv result,leftv first )
4145{
4146  spectrumState state = spectrumOK;
4147
4148  // -------------------
4149  //  check consistency
4150  // -------------------
4151
4152  //  check for a local polynomial ring
4153
4154  if( currRing->OrdSgn != -1 )
4155  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
4156  // or should we use:
4157  //if( !ringIsLocal( ) )
4158  {
4159    WerrorS( "only works for local orderings" );
4160    state = spectrumWrongRing;
4161  }
4162  else if( currRing->qideal != NULL )
4163  {
4164    WerrorS( "does not work in quotient rings" );
4165    state = spectrumWrongRing;
4166  }
4167  else
4168  {
4169    lists   L    = (lists)NULL;
4170    int     flag = 2; // symmetric optimization
4171
4172    state = spectrumCompute( (poly)first->Data( ),&L,flag );
4173
4174    if( state==spectrumOK )
4175    {
4176      result->rtyp = LIST_CMD;
4177      result->data = (char*)L;
4178    }
4179    else
4180    {
4181      spectrumPrintError(state);
4182    }
4183  }
4184
4185  return  (state!=spectrumOK);
4186}
4187
4188// ----------------------------------------------------------------------------
4189//  check if a list is a spectrum
4190//  check for:
4191//      list has 6 elements
4192//      1st element is int (mu=Milnor number)
4193//      2nd element is int (pg=geometrical genus)
4194//      3rd element is int (n =number of different spectrum numbers)
4195//      4th element is intvec (num=numerators)
4196//      5th element is intvec (den=denomiantors)
4197//      6th element is intvec (mul=multiplicities)
4198//      exactly n numerators
4199//      exactly n denominators
4200//      exactly n multiplicities
4201//      mu>0
4202//      pg>=0
4203//      n>0
4204//      num>0
4205//      den>0
4206//      mul>0
4207//      symmetriy with respect to numberofvariables/2
4208//      monotony
4209//      mu = sum of all multiplicities
4210//      pg = sum of all multiplicities where num/den<=1
4211// ----------------------------------------------------------------------------
4212
4213semicState  list_is_spectrum( lists l )
4214{
4215    // -------------------
4216    //  check list length
4217    // -------------------
4218
4219    if( l->nr < 5 )
4220    {
4221        return  semicListTooShort;
4222    }
4223    else if( l->nr > 5 )
4224    {
4225        return  semicListTooLong;
4226    }
4227
4228    // -------------
4229    //  check types
4230    // -------------
4231
4232    if( l->m[0].rtyp != INT_CMD )
4233    {
4234        return  semicListFirstElementWrongType;
4235    }
4236    else if( l->m[1].rtyp != INT_CMD )
4237    {
4238        return  semicListSecondElementWrongType;
4239    }
4240    else if( l->m[2].rtyp != INT_CMD )
4241    {
4242        return  semicListThirdElementWrongType;
4243    }
4244    else if( l->m[3].rtyp != INTVEC_CMD )
4245    {
4246        return  semicListFourthElementWrongType;
4247    }
4248    else if( l->m[4].rtyp != INTVEC_CMD )
4249    {
4250        return  semicListFifthElementWrongType;
4251    }
4252    else if( l->m[5].rtyp != INTVEC_CMD )
4253    {
4254        return  semicListSixthElementWrongType;
4255    }
4256
4257    // -------------------------
4258    //  check number of entries
4259    // -------------------------
4260
4261    int     mu = (int)(long)(l->m[0].Data( ));
4262    int     pg = (int)(long)(l->m[1].Data( ));
4263    int     n  = (int)(long)(l->m[2].Data( ));
4264
4265    if( n <= 0 )
4266    {
4267        return  semicListNNegative;
4268    }
4269
4270    intvec  *num = (intvec*)l->m[3].Data( );
4271    intvec  *den = (intvec*)l->m[4].Data( );
4272    intvec  *mul = (intvec*)l->m[5].Data( );
4273
4274    if( n != num->length( ) )
4275    {
4276        return  semicListWrongNumberOfNumerators;
4277    }
4278    else if( n != den->length( ) )
4279    {
4280        return  semicListWrongNumberOfDenominators;
4281    }
4282    else if( n != mul->length( ) )
4283    {
4284        return  semicListWrongNumberOfMultiplicities;
4285    }
4286
4287    // --------
4288    //  values
4289    // --------
4290
4291    if( mu <= 0 )
4292    {
4293        return  semicListMuNegative;
4294    }
4295    if( pg < 0 )
4296    {
4297        return  semicListPgNegative;
4298    }
4299
4300    int i;
4301
4302    for( i=0; i<n; i++ )
4303    {
4304        if( (*num)[i] <= 0 )
4305        {
4306            return  semicListNumNegative;
4307        }
4308        if( (*den)[i] <= 0 )
4309        {
4310            return  semicListDenNegative;
4311        }
4312        if( (*mul)[i] <= 0 )
4313        {
4314            return  semicListMulNegative;
4315        }
4316    }
4317
4318    // ----------------
4319    //  check symmetry
4320    // ----------------
4321
4322    int     j;
4323
4324    for( i=0, j=n-1; i<=j; i++,j-- )
4325    {
4326        if( (*num)[i] != rVar(currRing)*((*den)[i]) - (*num)[j] ||
4327            (*den)[i] != (*den)[j] ||
4328            (*mul)[i] != (*mul)[j] )
4329        {
4330            return  semicListNotSymmetric;
4331        }
4332    }
4333
4334    // ----------------
4335    //  check monotony
4336    // ----------------
4337
4338    for( i=0, j=1; i<n/2; i++,j++ )
4339    {
4340        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
4341        {
4342            return  semicListNotMonotonous;
4343        }
4344    }
4345
4346    // ---------------------
4347    //  check Milnor number
4348    // ---------------------
4349
4350    for( mu=0, i=0; i<n; i++ )
4351    {
4352        mu += (*mul)[i];
4353    }
4354
4355    if( mu != (int)(long)(l->m[0].Data( )) )
4356    {
4357        return  semicListMilnorWrong;
4358    }
4359
4360    // -------------------------
4361    //  check geometrical genus
4362    // -------------------------
4363
4364    for( pg=0, i=0; i<n; i++ )
4365    {
4366        if( (*num)[i]<=(*den)[i] )
4367        {
4368            pg += (*mul)[i];
4369        }
4370    }
4371
4372    if( pg != (int)(long)(l->m[1].Data( )) )
4373    {
4374        return  semicListPGWrong;
4375    }
4376
4377    return  semicOK;
4378}
4379
4380// ----------------------------------------------------------------------------
4381//  this procedure is called from the interpreter
4382// ----------------------------------------------------------------------------
4383//  first  = list of spectrum numbers
4384//  second = list of spectrum numbers
4385//  result = sum of the two lists
4386// ----------------------------------------------------------------------------
4387
4388BOOLEAN spaddProc( leftv result,leftv first,leftv second )
4389{
4390    semicState  state;
4391
4392    // -----------------
4393    //  check arguments
4394    // -----------------
4395
4396    lists l1 = (lists)first->Data( );
4397    lists l2 = (lists)second->Data( );
4398
4399    if( (state=list_is_spectrum( l1 )) != semicOK )
4400    {
4401        WerrorS( "first argument is not a spectrum:" );
4402        list_error( state );
4403    }
4404    else if( (state=list_is_spectrum( l2 )) != semicOK )
4405    {
4406        WerrorS( "second argument is not a spectrum:" );
4407        list_error( state );
4408    }
4409    else
4410    {
4411        spectrum s1= spectrumFromList ( l1 );
4412        spectrum s2= spectrumFromList ( l2 );
4413        spectrum sum( s1+s2 );
4414
4415        result->rtyp = LIST_CMD;
4416        result->data = (char*)(getList(sum));
4417    }
4418
4419    return  (state!=semicOK);
4420}
4421
4422// ----------------------------------------------------------------------------
4423//  this procedure is called from the interpreter
4424// ----------------------------------------------------------------------------
4425//  first  = list of spectrum numbers
4426//  second = integer
4427//  result = the multiple of the first list by the second factor
4428// ----------------------------------------------------------------------------
4429
4430BOOLEAN spmulProc( leftv result,leftv first,leftv second )
4431{
4432    semicState  state;
4433
4434    // -----------------
4435    //  check arguments
4436    // -----------------
4437
4438    lists   l = (lists)first->Data( );
4439    int     k = (int)(long)second->Data( );
4440
4441    if( (state=list_is_spectrum( l ))!=semicOK )
4442    {
4443        WerrorS( "first argument is not a spectrum" );
4444        list_error( state );
4445    }
4446    else if( k < 0 )
4447    {
4448        WerrorS( "second argument should be positive" );
4449        state = semicMulNegative;
4450    }
4451    else
4452    {
4453        spectrum s= spectrumFromList( l );
4454        spectrum product( k*s );
4455
4456        result->rtyp = LIST_CMD;
4457        result->data = (char*)getList(product);
4458    }
4459
4460    return  (state!=semicOK);
4461}
4462
4463// ----------------------------------------------------------------------------
4464//  this procedure is called from the interpreter
4465// ----------------------------------------------------------------------------
4466//  first  = list of spectrum numbers
4467//  second = list of spectrum numbers
4468//  result = semicontinuity index
4469// ----------------------------------------------------------------------------
4470
4471BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
4472{
4473  semicState  state;
4474  BOOLEAN qh=(((int)(long)w->Data())==1);
4475
4476  // -----------------
4477  //  check arguments
4478  // -----------------
4479
4480  lists l1 = (lists)u->Data( );
4481  lists l2 = (lists)v->Data( );
4482
4483  if( (state=list_is_spectrum( l1 ))!=semicOK )
4484  {
4485    WerrorS( "first argument is not a spectrum" );
4486    list_error( state );
4487  }
4488  else if( (state=list_is_spectrum( l2 ))!=semicOK )
4489  {
4490    WerrorS( "second argument is not a spectrum" );
4491    list_error( state );
4492  }
4493  else
4494  {
4495    spectrum s1= spectrumFromList( l1 );
4496    spectrum s2= spectrumFromList( l2 );
4497
4498    res->rtyp = INT_CMD;
4499    if (qh)
4500      res->data = (void*)(long)(s1.mult_spectrumh( s2 ));
4501    else
4502      res->data = (void*)(long)(s1.mult_spectrum( s2 ));
4503  }
4504
4505  // -----------------
4506  //  check status
4507  // -----------------
4508
4509  return  (state!=semicOK);
4510}
4511BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
4512{
4513  sleftv tmp;
4514  memset(&tmp,0,sizeof(tmp));
4515  tmp.rtyp=INT_CMD;
4516  /* tmp.data = (void *)0;  -- done by memset */
4517
4518  return  semicProc3(res,u,v,&tmp);
4519}
4520
4521#endif
4522
4523BOOLEAN loNewtonP( leftv res, leftv arg1 )
4524{
4525  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4526  return FALSE;
4527}
4528
4529BOOLEAN loSimplex( leftv res, leftv args )
4530{
4531  if ( !(rField_is_long_R(currRing)) )
4532  {
4533    WerrorS("Ground field not implemented!");
4534    return TRUE;
4535  }
4536
4537  simplex * LP;
4538  matrix m;
4539
4540  leftv v= args;
4541  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4542    return TRUE;
4543  else
4544    m= (matrix)(v->CopyD());
4545
4546  LP = new simplex(MATROWS(m),MATCOLS(m));
4547  LP->mapFromMatrix(m);
4548
4549  v= v->next;
4550  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
4551    return TRUE;
4552  else
4553    LP->m= (int)(long)(v->Data());
4554
4555  v= v->next;
4556  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
4557    return TRUE;
4558  else
4559    LP->n= (int)(long)(v->Data());
4560
4561  v= v->next;
4562  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
4563    return TRUE;
4564  else
4565    LP->m1= (int)(long)(v->Data());
4566
4567  v= v->next;
4568  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
4569    return TRUE;
4570  else
4571    LP->m2= (int)(long)(v->Data());
4572
4573  v= v->next;
4574  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
4575    return TRUE;
4576  else
4577    LP->m3= (int)(long)(v->Data());
4578
4579#ifdef mprDEBUG_PROT
4580  Print("m (constraints) %d\n",LP->m);
4581  Print("n (columns) %d\n",LP->n);
4582  Print("m1 (<=) %d\n",LP->m1);
4583  Print("m2 (>=) %d\n",LP->m2);
4584  Print("m3 (==) %d\n",LP->m3);
4585#endif
4586
4587  LP->compute();
4588
4589  lists lres= (lists)omAlloc( sizeof(slists) );
4590  lres->Init( 6 );
4591
4592  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4593  lres->m[0].data=(void*)LP->mapToMatrix(m);
4594
4595  lres->m[1].rtyp= INT_CMD;   // found a solution?
4596  lres->m[1].data=(void*)(long)LP->icase;
4597
4598  lres->m[2].rtyp= INTVEC_CMD;
4599  lres->m[2].data=(void*)LP->posvToIV();
4600
4601  lres->m[3].rtyp= INTVEC_CMD;
4602  lres->m[3].data=(void*)LP->zrovToIV();
4603
4604  lres->m[4].rtyp= INT_CMD;
4605  lres->m[4].data=(void*)(long)LP->m;
4606
4607  lres->m[5].rtyp= INT_CMD;
4608  lres->m[5].data=(void*)(long)LP->n;
4609
4610  res->data= (void*)lres;
4611
4612  return FALSE;
4613}
4614
4615BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4616{
4617  ideal gls = (ideal)(arg1->Data());
4618  int imtype= (int)(long)arg2->Data();
4619
4620  uResultant::resMatType mtype= determineMType( imtype );
4621
4622  // check input ideal ( = polynomial system )
4623  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4624  {
4625    return TRUE;
4626  }
4627
4628  uResultant *resMat= new uResultant( gls, mtype, false );
4629  if (resMat!=NULL)
4630  {
4631    res->rtyp = MODUL_CMD;
4632    res->data= (void*)resMat->accessResMat()->getMatrix();
4633    if (!errorreported) delete resMat;
4634  }
4635  return errorreported;
4636}
4637
4638BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4639{
4640
4641  poly gls;
4642  gls= (poly)(arg1->Data());
4643  int howclean= (int)(long)arg3->Data();
4644
4645  if ( !(rField_is_R(currRing) ||
4646         rField_is_Q(currRing) ||
4647         rField_is_long_R(currRing) ||
4648         rField_is_long_C(currRing)) )
4649  {
4650    WerrorS("Ground field not implemented!");
4651    return TRUE;
4652  }
4653
4654  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4655                          rField_is_long_C(currRing)) )
4656  {
4657    unsigned long int ii = (unsigned long int)arg2->Data();
4658    setGMPFloatDigits( ii, ii );
4659  }
4660
4661  if ( gls == NULL || pIsConstant( gls ) )
4662  {
4663    WerrorS("Input polynomial is constant!");
4664    return TRUE;
4665  }
4666
4667  int ldummy;
4668  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4669  int i,vpos=0;
4670  poly piter;
4671  lists elist;
4672  lists rlist;
4673
4674  elist= (lists)omAlloc( sizeof(slists) );
4675  elist->Init( 0 );
4676
4677  if ( rVar(currRing) > 1 )
4678  {
4679    piter= gls;
4680    for ( i= 1; i <= rVar(currRing); i++ )
4681      if ( pGetExp( piter, i ) )
4682      {
4683        vpos= i;
4684        break;
4685      }
4686    while ( piter )
4687    {
4688      for ( i= 1; i <= rVar(currRing); i++ )
4689        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4690        {
4691          WerrorS("The input polynomial must be univariate!");
4692          return TRUE;
4693        }
4694      pIter( piter );
4695    }
4696  }
4697
4698  rootContainer * roots= new rootContainer();
4699  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4700  piter= gls;
4701  for ( i= deg; i >= 0; i-- )
4702  {
4703    if ( piter && pTotaldegree(piter) == i )
4704    {
4705      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4706      //nPrint( pcoeffs[i] );PrintS("  ");
4707      pIter( piter );
4708    }
4709    else
4710    {
4711      pcoeffs[i]= nInit(0);
4712    }
4713  }
4714
4715#ifdef mprDEBUG_PROT
4716  for (i=deg; i >= 0; i--)
4717  {
4718    nPrint( pcoeffs[i] );PrintS("  ");
4719  }
4720  PrintLn();
4721#endif
4722
4723  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4724  roots->solver( howclean );
4725
4726  int elem= roots->getAnzRoots();
4727  char *dummy;
4728  int j;
4729
4730  rlist= (lists)omAlloc( sizeof(slists) );
4731  rlist->Init( elem );
4732
4733  if (rField_is_long_C(currRing))
4734  {
4735    for ( j= 0; j < elem; j++ )
4736    {
4737      rlist->m[j].rtyp=NUMBER_CMD;
4738      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4739      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4740    }
4741  }
4742  else
4743  {
4744    for ( j= 0; j < elem; j++ )
4745    {
4746      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4747      rlist->m[j].rtyp=STRING_CMD;
4748      rlist->m[j].data=(void *)dummy;
4749    }
4750  }
4751
4752  elist->Clean();
4753  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4754
4755  // this is (via fillContainer) the same data as in root
4756  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4757  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4758
4759  delete roots;
4760
4761  res->rtyp= LIST_CMD;
4762  res->data= (void*)rlist;
4763
4764  return FALSE;
4765}
4766
4767BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4768{
4769  int i;
4770  ideal p,w;
4771  p= (ideal)arg1->Data();
4772  w= (ideal)arg2->Data();
4773
4774  // w[0] = f(p^0)
4775  // w[1] = f(p^1)
4776  // ...
4777  // p can be a vector of numbers (multivariate polynom)
4778  //   or one number (univariate polynom)
4779  // tdg = deg(f)
4780
4781  int n= IDELEMS( p );
4782  int m= IDELEMS( w );
4783  int tdg= (int)(long)arg3->Data();
4784
4785  res->data= (void*)NULL;
4786
4787  // check the input
4788  if ( tdg < 1 )
4789  {
4790    WerrorS("Last input parameter must be > 0!");
4791    return TRUE;
4792  }
4793  if ( n != rVar(currRing) )
4794  {
4795    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4796    return TRUE;
4797  }
4798  if ( m != (int)pow((double)tdg+1,(double)n) )
4799  {
4800    Werror("Size of second input ideal must be equal to %d!",
4801      (int)pow((double)tdg+1,(double)n));
4802    return TRUE;
4803  }
4804  if ( !(rField_is_Q(currRing) /* ||
4805         rField_is_R() || rField_is_long_R() ||
4806         rField_is_long_C()*/ ) )
4807         {
4808    WerrorS("Ground field not implemented!");
4809    return TRUE;
4810  }
4811
4812  number tmp;
4813  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4814  for ( i= 0; i < n; i++ )
4815  {
4816    pevpoint[i]=nInit(0);
4817    if (  (p->m)[i] )
4818    {
4819      tmp = pGetCoeff( (p->m)[i] );
4820      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4821      {
4822        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4823        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4824        return TRUE;
4825      }
4826    } else tmp= NULL;
4827    if ( !nIsZero(tmp) )
4828    {
4829      if ( !pIsConstant((p->m)[i]))
4830      {
4831        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4832        WerrorS("Elements of first input ideal must be numbers!");
4833        return TRUE;
4834      }
4835      pevpoint[i]= nCopy( tmp );
4836    }
4837  }
4838
4839  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4840  for ( i= 0; i < m; i++ )
4841  {
4842    wresults[i]= nInit(0);
4843    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4844    {
4845      if ( !pIsConstant((w->m)[i]))
4846      {
4847        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4848        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4849        WerrorS("Elements of second input ideal must be numbers!");
4850        return TRUE;
4851      }
4852      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4853    }
4854  }
4855
4856  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4857  number *ncpoly= vm.interpolateDense( wresults );
4858  // do not free ncpoly[]!!
4859  poly rpoly= vm.numvec2poly( ncpoly );
4860
4861  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4862  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4863
4864  res->data= (void*)rpoly;
4865  return FALSE;
4866}
4867
4868BOOLEAN nuUResSolve( leftv res, leftv args )
4869{
4870  leftv v= args;
4871
4872  ideal gls;
4873  int imtype;
4874  int howclean;
4875
4876  // get ideal
4877  if ( v->Typ() != IDEAL_CMD )
4878    return TRUE;
4879  else gls= (ideal)(v->Data());
4880  v= v->next;
4881
4882  // get resultant matrix type to use (0,1)
4883  if ( v->Typ() != INT_CMD )
4884    return TRUE;
4885  else imtype= (int)(long)v->Data();
4886  v= v->next;
4887
4888  if (imtype==0)
4889  {
4890    ideal test_id=idInit(1,1);
4891    int j;
4892    for(j=IDELEMS(gls)-1;j>=0;j--)
4893    {
4894      if (gls->m[j]!=NULL)
4895      {
4896        test_id->m[0]=gls->m[j];
4897        intvec *dummy_w=id_QHomWeight(test_id, currRing);
4898        if (dummy_w!=NULL)
4899        {
4900          WerrorS("Newton polytope not of expected dimension");
4901          delete dummy_w;
4902          return TRUE;
4903        }
4904      }
4905    }
4906  }
4907
4908  // get and set precision in digits ( > 0 )
4909  if ( v->Typ() != INT_CMD )
4910    return TRUE;
4911  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4912                          rField_is_long_C(currRing)) )
4913  {
4914    unsigned long int ii=(unsigned long int)v->Data();
4915    setGMPFloatDigits( ii, ii );
4916  }
4917  v= v->next;
4918
4919  // get interpolation steps (0,1,2)
4920  if ( v->Typ() != INT_CMD )
4921    return TRUE;
4922  else howclean= (int)(long)v->Data();
4923
4924  uResultant::resMatType mtype= determineMType( imtype );
4925  int i,count;
4926  lists listofroots= NULL;
4927  number smv= NULL;
4928  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4929
4930  //emptylist= (lists)omAlloc( sizeof(slists) );
4931  //emptylist->Init( 0 );
4932
4933  //res->rtyp = LIST_CMD;
4934  //res->data= (void *)emptylist;
4935
4936  // check input ideal ( = polynomial system )
4937  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4938  {
4939    return TRUE;
4940  }
4941
4942  uResultant * ures;
4943  rootContainer ** iproots;
4944  rootContainer ** muiproots;
4945  rootArranger * arranger;
4946
4947  // main task 1: setup of resultant matrix
4948  ures= new uResultant( gls, mtype );
4949  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4950  {
4951    WerrorS("Error occurred during matrix setup!");
4952    return TRUE;
4953  }
4954
4955  // if dense resultant, check if minor nonsingular
4956  if ( mtype == uResultant::denseResMat )
4957  {
4958    smv= ures->accessResMat()->getSubDet();
4959#ifdef mprDEBUG_PROT
4960    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4961#endif
4962    if ( nIsZero(smv) )
4963    {
4964      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4965      return TRUE;
4966    }
4967  }
4968
4969  // main task 2: Interpolate specialized resultant polynomials
4970  if ( interpolate_det )
4971    iproots= ures->interpolateDenseSP( false, smv );
4972  else
4973    iproots= ures->specializeInU( false, smv );
4974
4975  // main task 3: Interpolate specialized resultant polynomials
4976  if ( interpolate_det )
4977    muiproots= ures->interpolateDenseSP( true, smv );
4978  else
4979    muiproots= ures->specializeInU( true, smv );
4980
4981#ifdef mprDEBUG_PROT
4982  int c= iproots[0]->getAnzElems();
4983  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4984  c= muiproots[0]->getAnzElems();
4985  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4986#endif
4987
4988  // main task 4: Compute roots of specialized polys and match them up
4989  arranger= new rootArranger( iproots, muiproots, howclean );
4990  arranger->solve_all();
4991
4992  // get list of roots
4993  if ( arranger->success() )
4994  {
4995    arranger->arrange();
4996    listofroots= listOfRoots(arranger, gmp_output_digits );
4997  }
4998  else
4999  {
5000    WerrorS("Solver was unable to find any roots!");
5001    return TRUE;
5002  }
5003
5004  // free everything
5005  count= iproots[0]->getAnzElems();
5006  for (i=0; i < count; i++) delete iproots[i];
5007  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5008  count= muiproots[0]->getAnzElems();
5009  for (i=0; i < count; i++) delete muiproots[i];
5010  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5011
5012  delete ures;
5013  delete arranger;
5014  nDelete( &smv );
5015
5016  res->data= (void *)listofroots;
5017
5018  //emptylist->Clean();
5019  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5020
5021  return FALSE;
5022}
5023
5024// from mpr_numeric.cc
5025lists listOfRoots( rootArranger* self, const unsigned int oprec )
5026{
5027  int i,j;
5028  int count= self->roots[0]->getAnzRoots(); // number of roots
5029  int elem= self->roots[0]->getAnzElems();  // number of koordinates per root
5030
5031  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
5032
5033  if ( self->found_roots )
5034  {
5035    listofroots->Init( count );
5036
5037    for (i=0; i < count; i++)
5038    {
5039      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
5040      onepoint->Init(elem);
5041      for ( j= 0; j < elem; j++ )
5042      {
5043        if ( !rField_is_long_C(currRing) )
5044        {
5045          onepoint->m[j].rtyp=STRING_CMD;
5046          onepoint->m[j].data=(void *)complexToStr((*self->roots[j])[i],oprec, currRing->cf);
5047        }
5048        else
5049        {
5050          onepoint->m[j].rtyp=NUMBER_CMD;
5051          onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
5052        }
5053        onepoint->m[j].next= NULL;
5054        onepoint->m[j].name= NULL;
5055      }
5056      listofroots->m[i].rtyp=LIST_CMD;
5057      listofroots->m[i].data=(void *)onepoint;
5058      listofroots->m[j].next= NULL;
5059      listofroots->m[j].name= NULL;
5060    }
5061
5062  }
5063  else
5064  {
5065    listofroots->Init( 0 );
5066  }
5067
5068  return listofroots;
5069}
5070
5071// from ring.cc
5072void rSetHdl(idhdl h)
5073{
5074  ring rg = NULL;
5075  if (h!=NULL)
5076  {
5077//   Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
5078    rg = IDRING(h);
5079    if (rg==NULL) return; //id <>NULL, ring==NULL
5080    omCheckAddrSize((ADDRESS)h,sizeof(idrec));
5081    if (IDID(h))  // OB: ????
5082      omCheckAddr((ADDRESS)IDID(h));
5083    rTest(rg);
5084  }
5085  else return;
5086
5087  // clean up history
5088  if (currRing!=NULL)
5089  {
5090    if(sLastPrinted.RingDependend())
5091    {
5092      sLastPrinted.CleanUp();
5093      //memset(&sLastPrinted,0,sizeof(sleftv)); // done by Cleanup,Init
5094    }
5095
5096    if (rg!=currRing)/*&&(currRing!=NULL)*/
5097    {
5098      if (rg->cf!=currRing->cf)
5099      {
5100        denominator_list dd=DENOMINATOR_LIST;
5101        if (DENOMINATOR_LIST!=NULL)
5102        {
5103          if (TEST_V_ALLWARN)
5104            Warn("deleting denom_list for ring change to %s",IDID(h));
5105          do
5106          {
5107            n_Delete(&(dd->n),currRing->cf);
5108            dd=dd->next;
5109            omFree(DENOMINATOR_LIST);
5110            DENOMINATOR_LIST=dd;
5111          } while(DENOMINATOR_LIST!=NULL);
5112        }
5113      }
5114    }
5115  }
5116
5117  // test for valid "currRing":
5118  if ((rg!=NULL) && (rg->idroot==NULL))
5119  {
5120    ring old=rg;
5121    rg=rAssure_HasComp(rg);
5122    if (old!=rg)
5123    {
5124      rKill(old);
5125      IDRING(h)=rg;
5126    }
5127  }
5128   /*------------ change the global ring -----------------------*/
5129  rChangeCurrRing(rg);
5130  currRingHdl = h;
5131}
5132
5133static leftv rOptimizeOrdAsSleftv(leftv ord)
5134{
5135  // change some bad orderings/combination into better ones
5136  leftv h=ord;
5137  while(h!=NULL)
5138  {
5139    BOOLEAN change=FALSE;
5140    intvec *iv = (intvec *)(h->data);
5141 // ws(-i) -> wp(i)
5142    if ((*iv)[1]==ringorder_ws)
5143    {
5144      BOOLEAN neg=TRUE;
5145      for(int i=2;i<iv->length();i++)
5146        if((*iv)[i]>=0) { neg=FALSE; break; }
5147      if (neg)
5148      {
5149        (*iv)[1]=ringorder_wp;
5150        for(int i=2;i<iv->length();i++)
5151          (*iv)[i]= - (*iv)[i];
5152        change=TRUE;
5153      }
5154    }
5155 // Ws(-i) -> Wp(i)
5156    if ((*iv)[1]==ringorder_Ws)
5157    {
5158      BOOLEAN neg=TRUE;
5159      for(int i=2;i<iv->length();i++)
5160        if((*iv)[i]>=0) { neg=FALSE; break; }
5161      if (neg)
5162      {
5163        (*iv)[1]=ringorder_Wp;
5164        for(int i=2;i<iv->length();i++)
5165          (*iv)[i]= -(*iv)[i];
5166        change=TRUE;
5167      }
5168    }
5169 // wp(1) -> dp
5170    if ((*iv)[1]==ringorder_wp)
5171    {
5172      BOOLEAN all_one=TRUE;
5173      for(int i=2;i<iv->length();i++)
5174        if((*iv)[i]!=1) { all_one=FALSE; break; }
5175      if (all_one)
5176      {
5177        intvec *iv2=new intvec(3);
5178        (*iv2)[0]=1;
5179        (*iv2)[1]=ringorder_dp;
5180        (*iv2)[2]=iv->length()-2;
5181        delete iv;
5182        iv=iv2;
5183        h->data=iv2;
5184        change=TRUE;
5185      }
5186    }
5187 // Wp(1) -> Dp
5188    if ((*iv)[1]==ringorder_Wp)
5189    {
5190      BOOLEAN all_one=TRUE;
5191      for(int i=2;i<iv->length();i++)
5192        if((*iv)[i]!=1) { all_one=FALSE; break; }
5193      if (all_one)
5194      {
5195        intvec *iv2=new intvec(3);
5196        (*iv2)[0]=1;
5197        (*iv2)[1]=ringorder_Dp;
5198        (*iv2)[2]=iv->length()-2;
5199        delete iv;
5200        iv=iv2;
5201        h->data=iv2;
5202        change=TRUE;
5203      }
5204    }
5205 // dp(1)/Dp(1)/rp(1) -> lp(1)
5206    if (((*iv)[1]==ringorder_dp)
5207    || ((*iv)[1]==ringorder_Dp)
5208    || ((*iv)[1]==ringorder_rp))
5209    {
5210      if (iv->length()==3)
5211      {
5212        if ((*iv)[2]==1)
5213        {
5214          if(h->next!=NULL)
5215          {
5216            intvec *iv2 = (intvec *)(h->next->data);
5217            if ((*iv2)[1]==ringorder_lp)
5218            {
5219              (*iv)[1]=ringorder_lp;
5220              change=TRUE;
5221            }
5222          }
5223        }
5224      }
5225    }
5226 // lp(i),lp(j) -> lp(i+j)
5227    if(((*iv)[1]==ringorder_lp)
5228    && (h->next!=NULL))
5229    {
5230      intvec *iv2 = (intvec *)(h->next->data);
5231      if ((*iv2)[1]==ringorder_lp)
5232      {
5233        leftv hh=h->next;
5234        h->next=hh->next;
5235        hh->next=NULL;
5236        if ((*iv2)[0]==1)
5237          (*iv)[2] += 1; // last block unspecified, at least 1
5238        else
5239          (*iv)[2] += (*iv2)[2];
5240        hh->CleanUp();
5241        omFree(hh);
5242        change=TRUE;
5243      }
5244    }
5245   // -------------------
5246    if (!change) h=h->next;
5247 }
5248 return ord;
5249}
5250
5251
5252BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
5253{
5254  int last = 0, o=0, n = 1, i=0, typ = 1, j;
5255  ord=rOptimizeOrdAsSleftv(ord);
5256  sleftv *sl = ord;
5257
5258  // determine nBlocks
5259  while (sl!=NULL)
5260  {
5261    intvec *iv = (intvec *)(sl->data);
5262    if (((*iv)[1]==ringorder_c)||((*iv)[1]==ringorder_C))
5263      i++;
5264    else if ((*iv)[1]==ringorder_L)
5265    {
5266      R->bitmask=(*iv)[2]*2+1;
5267      n--;
5268    }
5269    else if (((*iv)[1]!=ringorder_a)
5270    && ((*iv)[1]!=ringorder_a64)
5271    && ((*iv)[1]!=ringorder_am))
5272      o++;
5273    n++;
5274    sl=sl->next;
5275  }
5276  // check whether at least one real ordering
5277  if (o==0)
5278  {
5279    WerrorS("invalid combination of orderings");
5280    return TRUE;
5281  }
5282  // if no c/C ordering is given, increment n
5283  if (i==0) n++;
5284  else if (i != 1)
5285  {
5286    // throw error if more than one is given
5287    WerrorS("more than one ordering c/C specified");
5288    return TRUE;
5289  }
5290
5291  // initialize fields of R
5292  R->order=(rRingOrder_t *)omAlloc0(n*sizeof(rRingOrder_t));
5293  R->block0=(int *)omAlloc0(n*sizeof(int));
5294  R->block1=(int *)omAlloc0(n*sizeof(int));
5295  R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
5296
5297  int *weights=(int*)omAlloc0((R->N+1)*sizeof(int));
5298
5299  // init order, so that rBlocks works correctly
5300  for (j=0; j < n-1; j++)
5301    R->order[j] = ringorder_unspec;
5302  // set last _C order, if no c/C order was given
5303  if (i == 0) R->order[n-2] = ringorder_C;
5304
5305  /* init orders */
5306  sl=ord;
5307  n=-1;
5308  while (sl!=NULL)
5309  {
5310    intvec *iv;
5311    iv = (intvec *)(sl->data);
5312    if ((*iv)[1]!=ringorder_L)
5313    {
5314      n++;
5315
5316      /* the format of an ordering:
5317       *  iv[0]: factor
5318       *  iv[1]: ordering
5319       *  iv[2..end]: weights
5320       */
5321      R->order[n] = (rRingOrder_t)((*iv)[1]);
5322      typ=1;
5323      switch ((*iv)[1])
5324      {
5325          case ringorder_ws:
5326          case ringorder_Ws:
5327            typ=-1;
5328          case ringorder_wp:
5329          case ringorder_Wp:
5330            R->wvhdl[n]=(int*)omAlloc((iv->length()-1)*sizeof(int));
5331            R->block0[n] = last+1;
5332            for (i=2; i<iv->length(); i++)
5333            {
5334              R->wvhdl[n][i-2] = (*iv)[i];
5335              last++;
5336              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5337            }
5338            R->block1[n] = si_min(last,R->N);
5339            break;
5340          case ringorder_ls:
5341          case ringorder_ds:
5342          case ringorder_Ds:
5343          case ringorder_rs:
5344            typ=-1;
5345          case ringorder_lp:
5346          case ringorder_dp:
5347          case ringorder_Dp:
5348          case ringorder_rp:
5349            R->block0[n] = last+1;
5350            if (iv->length() == 3) last+=(*iv)[2];
5351            else last += (*iv)[0];
5352            R->block1[n] = si_min(last,R->N);
5353            if (rCheckIV(iv)) return TRUE;
5354            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
5355            {
5356              if (weights[i]==0) weights[i]=typ;
5357            }
5358            break;
5359
5360          case ringorder_s: // no 'rank' params!
5361          {
5362
5363            if(iv->length() > 3)
5364              return TRUE;
5365
5366            if(iv->length() == 3)
5367            {
5368              const int s = (*iv)[2];
5369              R->block0[n] = s;
5370              R->block1[n] = s;
5371            }
5372            break;
5373          }
5374          case ringorder_IS:
5375          {
5376            if(iv->length() != 3) return TRUE;
5377
5378            const int s = (*iv)[2];
5379
5380            if( 1 < s || s < -1 ) return TRUE;
5381
5382            R->block0[n] = s;
5383            R->block1[n] = s;
5384            break;
5385          }
5386          case ringorder_S:
5387          case ringorder_c:
5388          case ringorder_C:
5389          {
5390            if (rCheckIV(iv)) return TRUE;
5391            break;
5392          }
5393          case ringorder_aa:
5394          case ringorder_a:
5395          {
5396            R->block0[n] = last+1;
5397            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5398            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int));
5399            for (i=2; i<iv->length(); i++)
5400            {
5401              R->wvhdl[n][i-2]=(*iv)[i];
5402              last++;
5403              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5404            }
5405            last=R->block0[n]-1;
5406            break;
5407          }
5408          case ringorder_am:
5409          {
5410            R->block0[n] = last+1;
5411            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5412            R->wvhdl[n] = (int*)omAlloc(iv->length()*sizeof(int));
5413            if (R->block1[n]- R->block0[n]+2>=iv->length())
5414               WarnS("missing module weights");
5415            for (i=2; i<=(R->block1[n]-R->block0[n]+2); i++)
5416            {
5417              R->wvhdl[n][i-2]=(*iv)[i];
5418              last++;
5419              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5420            }
5421            R->wvhdl[n][i-2]=iv->length() -3 -(R->block1[n]- R->block0[n]);
5422            for (; i<iv->length(); i++)
5423            {
5424              R->wvhdl[n][i-1]=(*iv)[i];
5425            }
5426            last=R->block0[n]-1;
5427            break;
5428          }
5429          case ringorder_a64:
5430          {
5431            R->block0[n] = last+1;
5432            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5433            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int64));
5434            int64 *w=(int64 *)R->wvhdl[n];
5435            for (i=2; i<iv->length(); i++)
5436            {
5437              w[i-2]=(*iv)[i];
5438              last++;
5439              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5440            }
5441            last=R->block0[n]-1;
5442            break;
5443          }
5444          case ringorder_M:
5445          {
5446            int Mtyp=rTypeOfMatrixOrder(iv);
5447            if (Mtyp==0) return TRUE;
5448            if (Mtyp==-1) typ = -1;
5449
5450            R->wvhdl[n] =( int *)omAlloc((iv->length()-1)*sizeof(int));
5451            for (i=2; i<iv->length();i++)
5452              R->wvhdl[n][i-2]=(*iv)[i];
5453
5454            R->block0[n] = last+1;
5455            last += (int)sqrt((double)(iv->length()-2));
5456            R->block1[n] = si_min(last,R->N);
5457            for(i=R->block1[n];i>=R->block0[n];i--)
5458            {
5459              if (weights[i]==0) weights[i]=typ;
5460            }
5461            break;
5462          }
5463
5464          case ringorder_no:
5465            R->order[n] = ringorder_unspec;
5466            return TRUE;
5467
5468          default:
5469            Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
5470            R->order[n] = ringorder_unspec;
5471            return TRUE;
5472      }
5473    }
5474    if (last>R->N)
5475    {
5476      Werror("mismatch of number of vars (%d) and ordering (>=%d vars)",
5477             R->N,last);
5478      return TRUE;
5479    }
5480    sl=sl->next;
5481  }
5482  // find OrdSgn:
5483  R->OrdSgn = 1;
5484  for(i=1;i<=R->N;i++)
5485  { if (weights[i]<0) { R->OrdSgn=-1;break; }}
5486  omFree(weights);
5487
5488  // check for complete coverage
5489  while ( n >= 0 && (
5490          (R->order[n]==ringorder_c)
5491      ||  (R->order[n]==ringorder_C)
5492      ||  (R->order[n]==ringorder_s)
5493      ||  (R->order[n]==ringorder_S)
5494      ||  (R->order[n]==ringorder_IS)
5495                    )) n--;
5496
5497  assume( n >= 0 );
5498
5499  if (R->block1[n] != R->N)
5500  {
5501    if (((R->order[n]==ringorder_dp) ||
5502         (R->order[n]==ringorder_ds) ||
5503         (R->order[n]==ringorder_Dp) ||
5504         (R->order[n]==ringorder_Ds) ||
5505         (R->order[n]==ringorder_rp) ||
5506         (R->order[n]==ringorder_rs) ||
5507         (R->order[n]==ringorder_lp) ||
5508         (R->order[n]==ringorder_ls))
5509        &&
5510        R->block0[n] <= R->N)
5511    {
5512      R->block1[n] = R->N;
5513    }
5514    else
5515    {
5516      Werror("mismatch of number of vars (%d) and ordering (%d vars)",
5517             R->N,R->block1[n]);
5518      return TRUE;
5519    }
5520  }
5521  return FALSE;
5522}
5523
5524static BOOLEAN rSleftvList2StringArray(leftv sl, char** p)
5525{
5526
5527  while(sl!=NULL)
5528  {
5529    if ((sl->rtyp == IDHDL)||(sl->rtyp==ALIAS_CMD))
5530    {
5531      *p = omStrDup(sl->Name());
5532    }
5533    else if (sl->name!=NULL)
5534    {
5535      *p = (char*)sl->name;
5536      sl->name=NULL;
5537    }
5538    else if (sl->rtyp==POLY_CMD)
5539    {
5540      sleftv s_sl;
5541      iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
5542      if (s_sl.name != NULL)
5543      {
5544        *p = (char*)s_sl.name; s_sl.name=NULL;
5545      }
5546      else
5547        *p = NULL;
5548      sl->next = s_sl.next;
5549      s_sl.next = NULL;
5550      s_sl.CleanUp();
5551      if (*p == NULL) return TRUE;
5552    }
5553    else return TRUE;
5554    p++;
5555    sl=sl->next;
5556  }
5557  return FALSE;
5558}
5559
5560const short MAX_SHORT = 32767; // (1 << (sizeof(short)*8)) - 1;
5561
5562////////////////////
5563//
5564// rInit itself:
5565//
5566// INPUT:  pn: ch & parameter (names), rv: variable (names)
5567//         ord: ordering (all !=NULL)
5568// RETURN: currRingHdl on success
5569//         NULL        on error
5570// NOTE:   * makes new ring to current ring, on success
5571//         * considers input sleftv's as read-only
5572ring rInit(leftv pn, leftv rv, leftv ord)
5573{
5574  int float_len=0;
5575  int float_len2=0;
5576  ring R = NULL;
5577  //BOOLEAN ffChar=FALSE;
5578
5579  /* ch -------------------------------------------------------*/
5580  // get ch of ground field
5581
5582  // allocated ring
5583  R = (ring) omAlloc0Bin(sip_sring_bin);
5584
5585  coeffs cf = NULL;
5586
5587  assume( pn != NULL );
5588  const int P = pn->listLength();
5589
5590  if (pn->Typ()==CRING_CMD)
5591  {
5592    cf=(coeffs)pn->CopyD();
5593    leftv pnn=pn;
5594    if(P>1) /*parameter*/
5595    {
5596      pnn = pnn->next;
5597      const int pars = pnn->listLength();
5598      assume( pars > 0 );
5599      char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5600
5601      if (rSleftvList2StringArray(pnn, names))
5602      {
5603        WerrorS("parameter expected");
5604        goto rInitError;
5605      }
5606
5607      TransExtInfo extParam;
5608
5609      extParam.r = rDefault( cf, pars, names); // Q/Zp [ p_1, ... p_pars ]
5610      for(int i=pars-1; i>=0;i--)
5611      {
5612        omFree(names[i]);
5613      }
5614      omFree(names);
5615
5616      cf = nInitChar(n_transExt, &extParam);
5617    }
5618    assume( cf != NULL );
5619  }
5620  else if (pn->Typ()==INT_CMD)
5621  {
5622    int ch = (int)(long)pn->Data();
5623    leftv pnn=pn;
5624
5625    /* parameter? -------------------------------------------------------*/
5626    pnn = pnn->next;
5627
5628    if (pnn == NULL) // no params!?
5629    {
5630      if (ch!=0)
5631      {
5632        int ch2=IsPrime(ch);
5633        if ((ch<2)||(ch!=ch2))
5634        {
5635          Warn("%d is invalid as characteristic of the ground field. 32003 is used.", ch);
5636          ch=32003;
5637        }
5638        #ifndef TEST_ZN_AS_ZP
5639        cf = nInitChar(n_Zp, (void*)(long)ch);
5640        #else
5641        mpz_t modBase;
5642        mpz_init_set_ui(modBase, (long)ch);
5643        ZnmInfo info;
5644        info.base= modBase;
5645        info.exp= 1;
5646        cf=nInitChar(n_Zn,(void*) &info);
5647        cf->is_field=1;
5648        cf->is_domain=1;
5649        cf->has_simple_Inverse=1;
5650        #endif
5651      }
5652      else
5653        cf = nInitChar(n_Q, (void*)(long)ch);
5654    }
5655    else
5656    {
5657      const int pars = pnn->listLength();
5658
5659      assume( pars > 0 );
5660
5661      // predefined finite field: (p^k, a)
5662      if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
5663      {
5664        GFInfo param;
5665
5666        param.GFChar = ch;
5667        param.GFDegree = 1;
5668        param.GFPar_name = pnn->name;
5669
5670        cf = nInitChar(n_GF, &param);
5671      }
5672      else // (0/p, a, b, ..., z)
5673      {
5674        if ((ch!=0) && (ch!=IsPrime(ch)))
5675        {
5676          WerrorS("too many parameters");
5677          goto rInitError;
5678        }
5679
5680        char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5681
5682        if (rSleftvList2StringArray(pnn, names))
5683        {
5684          WerrorS("parameter expected");
5685          goto rInitError;
5686        }
5687
5688        TransExtInfo extParam;
5689
5690        extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5691        for(int i=pars-1; i>=0;i--)
5692        {
5693          omFree(names[i]);
5694        }
5695        omFree(names);
5696
5697        cf = nInitChar(n_transExt, &extParam);
5698      }
5699    }
5700
5701    //if (cf==NULL) ->Error: Invalid ground field specification
5702  }
5703  else if ((pn->name != NULL)
5704  && ((strcmp(pn->name,"real")==0) || (strcmp(pn->name,"complex")==0)))
5705  {
5706    leftv pnn=pn->next;
5707    BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5708    if ((pnn!=NULL) && (pnn->Typ()==INT_CMD))
5709    {
5710      float_len=(int)(long)pnn->Data();
5711      float_len2=float_len;
5712      pnn=pnn->next;
5713      if ((pnn!=NULL) && (pnn->Typ()==INT_CMD))
5714      {
5715        float_len2=(int)(long)pnn->Data();
5716        pnn=pnn->next;
5717      }
5718    }
5719
5720    if (!complex_flag)
5721      complex_flag= (pnn!=NULL) && (pnn->name!=NULL);
5722    if( !complex_flag && (float_len2 <= (short)SHORT_REAL_LENGTH))
5723       cf=nInitChar(n_R, NULL);
5724    else // longR or longC?
5725    {
5726       LongComplexInfo param;
5727
5728       param.float_len = si_min (float_len, 32767);
5729       param.float_len2 = si_min (float_len2, 32767);
5730
5731       // set the parameter name
5732       if (complex_flag)
5733       {
5734         if (param.float_len < SHORT_REAL_LENGTH)
5735         {
5736           param.float_len= SHORT_REAL_LENGTH;
5737           param.float_len2= SHORT_REAL_LENGTH;
5738         }
5739         if ((pnn == NULL) || (pnn->name == NULL))
5740           param.par_name=(const char*)"i"; //default to i
5741         else
5742           param.par_name = (const char*)pnn->name;
5743       }
5744
5745       cf = nInitChar(complex_flag ? n_long_C: n_long_R, (void*)&param);
5746    }
5747    assume( cf != NULL );
5748  }
5749#ifdef HAVE_RINGS
5750  else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5751  {
5752    // TODO: change to use coeffs_BIGINT!?
5753    mpz_t modBase;
5754    unsigned int modExponent = 1;
5755    mpz_init_set_si(modBase, 0);
5756    if (pn->next!=NULL)
5757    {
5758      leftv pnn=pn;
5759      if (pnn->next->Typ()==INT_CMD)
5760      {
5761        pnn=pnn->next;
5762        mpz_set_ui(modBase, (long) pnn->Data());
5763        if ((pnn->next!=NULL) && (pnn->next->Typ()==INT_CMD))
5764        {
5765          pnn=pnn->next;
5766          modExponent = (long) pnn->Data();
5767        }
5768        while ((pnn->next!=NULL) && (pnn->next->Typ()==INT_CMD))
5769        {
5770          pnn=pnn->next;
5771          mpz_mul_ui(modBase, modBase, (int)(long) pnn->Data());
5772        }
5773      }
5774      else if (pnn->next->Typ()==BIGINT_CMD)
5775      {
5776        number p=(number)pnn->next->CopyD();
5777        nlGMP(p,modBase,coeffs_BIGINT); // TODO? // extern void   nlGMP(number &i, mpz_t n, const coeffs r); // FIXME: n_MPZ( modBase, p, coeffs_BIGINT); ?
5778        n_Delete(&p,coeffs_BIGINT);
5779      }
5780    }
5781    else
5782      cf=nInitChar(n_Z,NULL);
5783
5784    if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_sgn1(modBase) < 0))
5785    {
5786      WerrorS("Wrong ground ring specification (module is 1)");
5787      goto rInitError;
5788    }
5789    if (modExponent < 1)
5790    {
5791      WerrorS("Wrong ground ring specification (exponent smaller than 1");
5792      goto rInitError;
5793    }
5794    // module is 0 ---> integers ringtype = 4;
5795    // we have an exponent
5796    if (modExponent > 1 && cf == NULL)
5797    {
5798      if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
5799      {
5800        /* this branch should be active for modExponent = 2..32 resp. 2..64,
5801           depending on the size of a long on the respective platform */
5802        //ringtype = 1;       // Use Z/2^ch
5803        cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5804      }
5805      else
5806      {
5807        if (mpz_sgn1(modBase)==0)
5808        {
5809          WerrorS("modulus must not be 0 or parameter not allowed");
5810          goto rInitError;
5811        }
5812        //ringtype = 3;
5813        ZnmInfo info;
5814        info.base= modBase;
5815        info.exp= modExponent;
5816        cf=nInitChar(n_Znm,(void*) &info); //exponent is missing
5817      }
5818    }
5819    // just a module m > 1
5820    else if (cf == NULL)
5821    {
5822      if (mpz_sgn1(modBase)==0)
5823      {
5824        WerrorS("modulus must not be 0 or parameter not allowed");
5825        goto rInitError;
5826      }
5827      //ringtype = 2;
5828      ZnmInfo info;
5829      info.base= modBase;
5830      info.exp= modExponent;
5831      cf=nInitChar(n_Zn,(void*) &info);
5832    }
5833    assume( cf != NULL );
5834    mpz_clear(modBase);
5835  }
5836#endif
5837  // ring NEW = OLD, (), (); where OLD is a polynomial ring...
5838  else if ((pn->Typ()==RING_CMD) && (P == 1))
5839  {
5840    TransExtInfo extParam;
5841    extParam.r = (ring)pn->Data();
5842    cf = nInitChar(n_transExt, &extParam);
5843  }
5844  //else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5845  //{
5846  //  AlgExtInfo extParam;
5847  //  extParam.r = (ring)pn->Data();
5848
5849  //  cf = nInitChar(n_algExt, &extParam);   // Q[a]/<minideal>
5850  //}
5851  else
5852  {
5853    WerrorS("Wrong or unknown ground field specification");
5854#if 0
5855// debug stuff for unknown cf descriptions:
5856    sleftv* p = pn;
5857    while (p != NULL)
5858    {
5859      Print( "pn[%p]: type: %d [%s]: %p, name: %s", (void*)p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5860      PrintLn();
5861      p = p->next;
5862    }
5863#endif
5864    goto rInitError;
5865  }
5866
5867  /*every entry in the new ring is initialized to 0*/
5868
5869  /* characteristic -----------------------------------------------*/
5870  /* input: 0 ch=0 : Q     parameter=NULL    ffChar=FALSE   float_len
5871   *         0    1 : Q(a,...)        *names         FALSE
5872   *         0   -1 : R               NULL           FALSE  0
5873   *         0   -1 : R               NULL           FALSE  prec. >6
5874   *         0   -1 : C               *names         FALSE  prec. 0..?
5875   *         p    p : Fp              NULL           FALSE
5876   *         p   -p : Fp(a)           *names         FALSE
5877   *         q    q : GF(q=p^n)       *names         TRUE
5878  */
5879  if (cf==NULL)
5880  {
5881    WerrorS("Invalid ground field specification");
5882    goto rInitError;
5883//    const int ch=32003;
5884//    cf=nInitChar(n_Zp, (void*)(long)ch);
5885  }
5886
5887  assume( R != NULL );
5888
5889  R->cf = cf;
5890
5891  /* names and number of variables-------------------------------------*/
5892  {
5893    int l=rv->listLength();
5894
5895    if (l>MAX_SHORT)
5896    {
5897      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5898       goto rInitError;
5899    }
5900    R->N = l; /*rv->listLength();*/
5901  }
5902  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5903  if (rSleftvList2StringArray(rv, R->names))
5904  {
5905    WerrorS("name of ring variable expected");
5906    goto rInitError;
5907  }
5908
5909  /* check names and parameters for conflicts ------------------------- */
5910  rRenameVars(R); // conflicting variables will be renamed
5911  /* ordering -------------------------------------------------------------*/
5912  if (rSleftvOrdering2Ordering(ord, R))
5913    goto rInitError;
5914
5915  // Complete the initialization
5916  if (rComplete(R,1))
5917    goto rInitError;
5918
5919/*#ifdef HAVE_RINGS
5920// currently, coefficients which are ring elements require a global ordering:
5921  if (rField_is_Ring(R) && (R->OrdSgn==-1))
5922  {
5923    WerrorS("global ordering required for these coefficients");