source: git/Singular/ipshell.cc @ 5abb79f

fieker-DuValspielwiese
Last change on this file since 5abb79f was 8fd62c, checked in by Hans Schoenemann <hannes@…>, 6 years ago
add rPlusVar, rMinusVar (ring manipulation: R+"y","y"+R,R-"y")
  • Property mode set to 100644
File size: 161.6 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
76leftv iiCurrArgs=NULL;
77idhdl iiCurrProc=NULL;
78const char *lastreserved=NULL;
79
80static BOOLEAN iiNoKeepRing=TRUE;
81
82/*0 implementation*/
83
84const char * iiTwoOps(int t)
85{
86  if (t<127)
87  {
88    static 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
985BOOLEAN 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");
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  // 1:
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!=2)
2262    || (LL->m[0].rtyp!=INT_CMD)
2263    || (LL->m[1].rtyp!=INT_CMD))
2264  && ((LL->nr!=1)
2265    || (LL->m[0].rtyp!=INT_CMD)))
2266  {
2267    WerrorS("invalid coeff. field description list");
2268    return;
2269  }
2270  int r1=(int)(long)LL->m[0].data;
2271  int r2=(int)(long)LL->m[1].data;
2272  if (L->nr==2) // complex
2273    R->cf = nInitChar(n_long_C, NULL);
2274  else if ((r1<=SHORT_REAL_LENGTH)
2275  && (r2<=SHORT_REAL_LENGTH))
2276    R->cf = nInitChar(n_R, NULL);
2277  else
2278  {
2279    LongComplexInfo* p = (LongComplexInfo *)omAlloc0(sizeof(LongComplexInfo));
2280    p->float_len=r1;
2281    p->float_len2=r2;
2282    R->cf = nInitChar(n_long_R, p);
2283  }
2284
2285  if ((r1<=SHORT_REAL_LENGTH)   // should go into nInitChar
2286  && (r2<=SHORT_REAL_LENGTH))
2287  {
2288    R->cf->float_len=SHORT_REAL_LENGTH/2;
2289    R->cf->float_len2=SHORT_REAL_LENGTH;
2290  }
2291  else
2292  {
2293    R->cf->float_len=si_min(r1,32767);
2294    R->cf->float_len2=si_min(r2,32767);
2295  }
2296  // ----------------------------------------
2297  // 2: list (par)
2298  if (L->nr==2)
2299  {
2300    //R->cf->extRing->N=1;
2301    if (L->m[2].rtyp!=STRING_CMD)
2302    {
2303      WerrorS("invalid coeff. field description, expecting parameter name");
2304      return;
2305    }
2306    //(rParameter(R))=(char**)omAlloc0(rPar(R)*sizeof(char_ptr));
2307    rParameter(R)[0]=omStrDup((char *)L->m[2].data);
2308  }
2309  // ----------------------------------------
2310}
2311
2312#ifdef HAVE_RINGS
2313void rComposeRing(lists L, ring R)
2314/* field is R or C */
2315{
2316  // ----------------------------------------
2317  // 0: string: integer
2318  // no further entries --> Z
2319  mpz_t modBase;
2320  unsigned int modExponent = 1;
2321
2322  if (L->nr == 0)
2323  {
2324    mpz_init_set_ui(modBase,0);
2325    modExponent = 1;
2326  }
2327  // ----------------------------------------
2328  // 1:
2329  else
2330  {
2331    if (L->m[1].rtyp!=LIST_CMD) WerrorS("invalid data, expecting list of numbers");
2332    lists LL=(lists)L->m[1].data;
2333    if ((LL->nr >= 0) && LL->m[0].rtyp == BIGINT_CMD)
2334    {
2335      number tmp= (number) LL->m[0].data; // never use CopyD() on list elements
2336                                    // assume that tmp is integer, not rational
2337      mpz_init(modBase);
2338      n_MPZ (modBase, tmp, coeffs_BIGINT);
2339    }
2340    else if (LL->nr >= 0 && LL->m[0].rtyp == INT_CMD)
2341    {
2342      mpz_init_set_ui(modBase,(unsigned long) LL->m[0].data);
2343    }
2344    else
2345    {
2346      mpz_init_set_ui(modBase,0);
2347    }
2348    if (LL->nr >= 1)
2349    {
2350      modExponent = (unsigned long) LL->m[1].data;
2351    }
2352    else
2353    {
2354      modExponent = 1;
2355    }
2356  }
2357  // ----------------------------------------
2358  if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_sgn1(modBase) < 0))
2359  {
2360    WerrorS("Wrong ground ring specification (module is 1)");
2361    return;
2362  }
2363  if (modExponent < 1)
2364  {
2365    WerrorS("Wrong ground ring specification (exponent smaller than 1)");
2366    return;
2367  }
2368  // module is 0 ---> integers
2369  if (mpz_sgn1(modBase) == 0)
2370  {
2371    R->cf=nInitChar(n_Z,NULL);
2372  }
2373  // we have an exponent
2374  else if (modExponent > 1)
2375  {
2376    //R->cf->ch = R->cf->modExponent;
2377    if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
2378    {
2379      /* this branch should be active for modExponent = 2..32 resp. 2..64,
2380           depending on the size of a long on the respective platform */
2381      R->cf=nInitChar(n_Z2m,(void*)(long)modExponent);       // Use Z/2^ch
2382    }
2383    else
2384    {
2385      //ringtype 3
2386      ZnmInfo info;
2387      info.base= modBase;
2388      info.exp= modExponent;
2389      R->cf=nInitChar(n_Znm,(void*) &info);
2390    }
2391  }
2392  // just a module m > 1
2393  else
2394  {
2395    //ringtype = 2;
2396    //const int ch = mpz_get_ui(modBase);
2397    ZnmInfo info;
2398    info.base= modBase;
2399    info.exp= modExponent;
2400    R->cf=nInitChar(n_Zn,(void*) &info);
2401  }
2402  mpz_clear(modBase);
2403}
2404#endif
2405
2406static void rRenameVars(ring R)
2407{
2408  int i,j;
2409  BOOLEAN ch;
2410  do
2411  {
2412    ch=0;
2413    for(i=0;i<R->N-1;i++)
2414    {
2415      for(j=i+1;j<R->N;j++)
2416      {
2417        if (strcmp(R->names[i],R->names[j])==0)
2418        {
2419          ch=TRUE;
2420          Warn("name conflict var(%d) and var(%d): `%s`, rename to `@%s`",i+1,j+1,R->names[i],R->names[i]);
2421          omFree(R->names[j]);
2422          R->names[j]=(char *)omAlloc(2+strlen(R->names[i]));
2423          sprintf(R->names[j],"@%s",R->names[i]);
2424        }
2425      }
2426    }
2427  }
2428  while (ch);
2429  for(i=0;i<rPar(R); i++)
2430  {
2431    for(j=0;j<R->N;j++)
2432    {
2433      if (strcmp(rParameter(R)[i],R->names[j])==0)
2434      {
2435        Warn("name conflict par(%d) and var(%d): `%s`, renaming the VARIABLE to `@@(%d)`",i+1,j+1,R->names[j],i+1);
2436//        omFree(rParameter(R)[i]);
2437//        rParameter(R)[i]=(char *)omAlloc(10);
2438//        sprintf(rParameter(R)[i],"@@(%d)",i+1);
2439        omFree(R->names[j]);
2440        R->names[j]=(char *)omAlloc(10);
2441        sprintf(R->names[j],"@@(%d)",i+1);
2442      }
2443    }
2444  }
2445}
2446
2447static inline BOOLEAN rComposeVar(const lists  L, ring R)
2448{
2449  assume(R!=NULL);
2450  if (L->m[1].Typ()==LIST_CMD)
2451  {
2452    lists v=(lists)L->m[1].Data();
2453    R->N = v->nr+1;
2454    if (R->N<=0)
2455    {
2456      WerrorS("no ring variables");
2457      return TRUE;
2458    }
2459    R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
2460    int i;
2461    for(i=0;i<R->N;i++)
2462    {
2463      if (v->m[i].Typ()==STRING_CMD)
2464        R->names[i]=omStrDup((char *)v->m[i].Data());
2465      else if (v->m[i].Typ()==POLY_CMD)
2466      {
2467        poly p=(poly)v->m[i].Data();
2468        int nr=pIsPurePower(p);
2469        if (nr>0)
2470          R->names[i]=omStrDup(currRing->names[nr-1]);
2471        else
2472        {
2473          Werror("var name %d must be a string or a ring variable",i+1);
2474          return TRUE;
2475        }
2476      }
2477      else
2478      {
2479        Werror("var name %d must be `string` (not %d)",i+1, v->m[i].Typ());
2480        return TRUE;
2481      }
2482    }
2483  }
2484  else
2485  {
2486    WerrorS("variable must be given as `list`");
2487    return TRUE;
2488  }
2489  return FALSE;
2490}
2491
2492static inline BOOLEAN rComposeOrder(const lists  L, const BOOLEAN check_comp, ring R)
2493{
2494  assume(R!=NULL);
2495  long bitmask=0L;
2496  if (L->m[2].Typ()==LIST_CMD)
2497  {
2498    lists v=(lists)L->m[2].Data();
2499    int n= v->nr+2;
2500    int j_in_R,j_in_L;
2501    // do we have an entry "L",... ?: set bitmask
2502    for (int j=0; j < n-1; j++)
2503    {
2504      if (v->m[j].Typ()==LIST_CMD)
2505      {
2506        lists vv=(lists)v->m[j].Data();
2507        if ((vv->nr==1)
2508        &&(vv->m[0].Typ()==STRING_CMD)
2509        &&(strcmp((char*)vv->m[0].Data(),"L")==0))
2510        {
2511          number nn=(number)vv->m[1].Data();
2512          if (vv->m[1].Typ()==BIGINT_CMD)
2513            bitmask=n_Int(nn,coeffs_BIGINT);
2514          else if (vv->m[1].Typ()==INT_CMD)
2515            bitmask=(long)nn;
2516          else
2517          {
2518            Werror("illegal argument for pseudo ordering L: %d",vv->m[1].Typ());
2519            return TRUE;
2520          }
2521          break;
2522        }
2523      }
2524    }
2525    if (bitmask!=0) n--;
2526
2527    // initialize fields of R
2528    R->order=(rRingOrder_t *)omAlloc0(n*sizeof(rRingOrder_t));
2529    R->block0=(int *)omAlloc0(n*sizeof(int));
2530    R->block1=(int *)omAlloc0(n*sizeof(int));
2531    R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
2532    // init order, so that rBlocks works correctly
2533    for (j_in_R= n-2; j_in_R>=0; j_in_R--)
2534      R->order[j_in_R] = ringorder_unspec;
2535    // orderings
2536    for(j_in_R=0,j_in_L=0;j_in_R<n-1;j_in_R++,j_in_L++)
2537    {
2538    // todo: a(..), M
2539      if (v->m[j_in_L].Typ()!=LIST_CMD)
2540      {
2541        WerrorS("ordering must be list of lists");
2542        return TRUE;
2543      }
2544      lists vv=(lists)v->m[j_in_L].Data();
2545      if ((vv->nr==1)
2546      && (vv->m[0].Typ()==STRING_CMD))
2547      {
2548        if (strcmp((char*)vv->m[0].Data(),"L")==0)
2549        {
2550          j_in_R--;
2551          continue;
2552        }
2553        if ((vv->m[1].Typ()!=INTVEC_CMD) && (vv->m[1].Typ()!=INT_CMD))
2554        {
2555          PrintS(lString(vv));
2556          WerrorS("ordering name must be a (string,intvec)(1)");
2557          return TRUE;
2558        }
2559        R->order[j_in_R]=rOrderName(omStrDup((char*)vv->m[0].Data())); // assume STRING
2560
2561        if (j_in_R==0) R->block0[0]=1;
2562        else
2563        {
2564          int jj=j_in_R-1;
2565          while((jj>=0)
2566          && ((R->order[jj]== ringorder_a)
2567             || (R->order[jj]== ringorder_aa)
2568             || (R->order[jj]== ringorder_am)
2569             || (R->order[jj]== ringorder_c)
2570             || (R->order[jj]== ringorder_C)
2571             || (R->order[jj]== ringorder_s)
2572             || (R->order[jj]== ringorder_S)
2573          ))
2574          {
2575            //Print("jj=%, skip %s\n",rSimpleOrdStr(R->order[jj]));
2576            jj--;
2577          }
2578          if (jj<0) R->block0[j_in_R]=1;
2579          else      R->block0[j_in_R]=R->block1[jj]+1;
2580        }
2581        intvec *iv;
2582        if (vv->m[1].Typ()==INT_CMD)
2583          iv=new intvec((int)(long)vv->m[1].Data(),(int)(long)vv->m[1].Data());
2584        else
2585          iv=ivCopy((intvec*)vv->m[1].Data()); //assume INTVEC
2586        int iv_len=iv->length();
2587        if ((R->order[j_in_R]!=ringorder_s)
2588        &&(R->order[j_in_R]!=ringorder_c)
2589        &&(R->order[j_in_R]!=ringorder_C))
2590        {
2591          R->block1[j_in_R]=si_max(R->block0[j_in_R],R->block0[j_in_R]+iv_len-1);
2592          if (R->block1[j_in_R]>R->N)
2593          {
2594            if (R->block0[j_in_R]>R->N)
2595            {
2596              Werror("not enough variables for ordering %d (%s)",j_in_R,rSimpleOrdStr(R->order[j_in_R]));
2597              return TRUE;
2598            }
2599            R->block1[j_in_R]=R->N;
2600            iv_len=R->block1[j_in_R]-R->block0[j_in_R]+1;
2601          }
2602          //Print("block %d from %d to %d\n",j,R->block0[j], R->block1[j]);
2603        }
2604        int i;
2605        switch (R->order[j_in_R])
2606        {
2607           case ringorder_ws:
2608           case ringorder_Ws:
2609              R->OrdSgn=-1;
2610           case ringorder_aa:
2611           case ringorder_a:
2612           case ringorder_wp:
2613           case ringorder_Wp:
2614             R->wvhdl[j_in_R] =( int *)omAlloc(iv_len*sizeof(int));
2615             for (i=0; i<iv_len;i++)
2616             {
2617               R->wvhdl[j_in_R][i]=(*iv)[i];
2618             }
2619             break;
2620           case ringorder_am:
2621             R->wvhdl[j_in_R] =( int *)omAlloc((iv->length()+1)*sizeof(int));
2622             for (i=0; i<iv_len;i++)
2623             {
2624               R->wvhdl[j_in_R][i]=(*iv)[i];
2625             }
2626             R->wvhdl[j_in_R][i]=iv->length() - iv_len;
2627             //printf("ivlen:%d,iv->len:%d,mod:%d\n",iv_len,iv->length(),R->wvhdl[j][i]);
2628             for (; i<iv->length(); i++)
2629             {
2630                R->wvhdl[j_in_R][i+1]=(*iv)[i];
2631             }
2632             break;
2633           case ringorder_M:
2634             R->wvhdl[j_in_R] =( int *)omAlloc((iv->length())*sizeof(int));
2635             for (i=0; i<iv->length();i++) R->wvhdl[j_in_R][i]=(*iv)[i];
2636             R->block1[j_in_R]=si_max(R->block0[j_in_R],R->block0[j_in_R]+(int)sqrt((double)(iv->length()-1)));
2637             if (R->block1[j_in_R]>R->N)
2638             {
2639               WerrorS("ordering matrix too big");
2640               return TRUE;
2641             }
2642             break;
2643           case ringorder_ls:
2644           case ringorder_ds:
2645           case ringorder_Ds:
2646           case ringorder_rs:
2647             R->OrdSgn=-1;
2648           case ringorder_lp:
2649           case ringorder_dp:
2650           case ringorder_Dp:
2651           case ringorder_rp:
2652             break;
2653           case ringorder_S:
2654             break;
2655           case ringorder_c:
2656           case ringorder_C:
2657             R->block1[j_in_R]=R->block0[j_in_R]=0;
2658             break;
2659
2660           case ringorder_s:
2661             R->block1[j_in_R]=R->block0[j_in_R]=(*iv)[0];
2662             rSetSyzComp(R->block0[j_in_R],R);
2663             break;
2664
2665           case ringorder_IS:
2666           {
2667             R->block1[j_in_R] = R->block0[j_in_R] = 0;
2668             if( iv->length() > 0 )
2669             {
2670               const int s = (*iv)[0];
2671               assume( -2 < s && s < 2 );
2672               R->block1[j_in_R] = R->block0[j_in_R] = s;
2673             }
2674             break;
2675           }
2676           case 0:
2677           case ringorder_unspec:
2678             break;
2679           case ringorder_L: /* cannot happen */
2680           case ringorder_a64: /*not implemented */
2681             WerrorS("ring order not implemented");
2682             return TRUE;
2683        }
2684        delete iv;
2685      }
2686      else
2687      {
2688        PrintS(lString(vv));
2689        WerrorS("ordering name must be a (string,intvec)");
2690        return TRUE;
2691      }
2692    }
2693    // sanity check
2694    j_in_R=n-2;
2695    if ((R->order[j_in_R]==ringorder_c)
2696    || (R->order[j_in_R]==ringorder_C)
2697    || (R->order[j_in_R]==ringorder_unspec)) j_in_R--;
2698    if (R->block1[j_in_R] != R->N)
2699    {
2700      if (((R->order[j_in_R]==ringorder_dp) ||
2701           (R->order[j_in_R]==ringorder_ds) ||
2702           (R->order[j_in_R]==ringorder_Dp) ||
2703           (R->order[j_in_R]==ringorder_Ds) ||
2704           (R->order[j_in_R]==ringorder_rp) ||
2705           (R->order[j_in_R]==ringorder_rs) ||
2706           (R->order[j_in_R]==ringorder_lp) ||
2707           (R->order[j_in_R]==ringorder_ls))
2708          &&
2709            R->block0[j_in_R] <= R->N)
2710      {
2711        R->block1[j_in_R] = R->N;
2712      }
2713      else
2714      {
2715        Werror("ordering incomplete: size (%d) should be %d",R->block1[j_in_R],R->N);
2716        return TRUE;
2717      }
2718    }
2719    if (R->block0[j_in_R]>R->N)
2720    {
2721      Werror("not enough variables (%d) for ordering block %d, scanned so far:",R->N,j_in_R+1);
2722      for(int ii=0;ii<=j_in_R;ii++)
2723        Werror("ord[%d]: %s from v%d to v%d",ii+1,rSimpleOrdStr(R->order[ii]),R->block0[ii],R->block1[ii]);
2724      return TRUE;
2725    }
2726    if (check_comp)
2727    {
2728      BOOLEAN comp_order=FALSE;
2729      int jj;
2730      for(jj=0;jj<n;jj++)
2731      {
2732        if ((R->order[jj]==ringorder_c) ||
2733            (R->order[jj]==ringorder_C)) { comp_order=TRUE; break; }
2734      }
2735      if (!comp_order)
2736      {
2737        R->order=(rRingOrder_t*)omRealloc0Size(R->order,n*sizeof(rRingOrder_t),(n+1)*sizeof(rRingOrder_t));
2738        R->block0=(int*)omRealloc0Size(R->block0,n*sizeof(int),(n+1)*sizeof(int));
2739        R->block1=(int*)omRealloc0Size(R->block1,n*sizeof(int),(n+1)*sizeof(int));
2740        R->wvhdl=(int**)omRealloc0Size(R->wvhdl,n*sizeof(int_ptr),(n+1)*sizeof(int_ptr));
2741        R->order[n-1]=ringorder_C;
2742        R->block0[n-1]=0;
2743        R->block1[n-1]=0;
2744        R->wvhdl[n-1]=NULL;
2745        n++;
2746      }
2747    }
2748  }
2749  else
2750  {
2751    WerrorS("ordering must be given as `list`");
2752    return TRUE;
2753  }
2754  if (bitmask!=0) R->bitmask=bitmask*2;
2755  return FALSE;
2756}
2757
2758ring rCompose(const lists  L, const BOOLEAN check_comp, const long bitmask,const int isLetterplace)
2759{
2760  if ((L->nr!=3)
2761#ifdef HAVE_PLURAL
2762  &&(L->nr!=5)
2763#endif
2764  )
2765    return NULL;
2766  int is_gf_char=0;
2767  // 0: char/ cf - ring
2768  // 1: list (var)
2769  // 2: list (ord)
2770  // 3: qideal
2771  // possibly:
2772  // 4: C
2773  // 5: D
2774
2775  ring R = (ring) omAlloc0Bin(sip_sring_bin);
2776
2777  // ------------------------------------------------------------------
2778  // 0: char:
2779  if (L->m[0].Typ()==CRING_CMD)
2780  {
2781    R->cf=(coeffs)L->m[0].Data();
2782    R->cf->ref++;
2783  }
2784  else if (L->m[0].Typ()==INT_CMD)
2785  {
2786    int ch = (int)(long)L->m[0].Data();
2787    assume( ch >= 0 );
2788
2789    if (ch == 0) // Q?
2790      R->cf = nInitChar(n_Q, NULL);
2791    else
2792    {
2793      int l = IsPrime(ch); // Zp?
2794      if( l != ch )
2795      {
2796        Warn("%d is invalid characteristic of ground field. %d is used.", ch, l);
2797        ch = l;
2798      }
2799      #ifndef TEST_ZN_AS_ZP
2800      R->cf = nInitChar(n_Zp, (void*)(long)ch);
2801      #else
2802      mpz_t modBase;
2803      mpz_init_set_ui(modBase,(long) ch);
2804      ZnmInfo info;
2805      info.base= modBase;
2806      info.exp= 1;
2807      R->cf=nInitChar(n_Zn,(void*) &info); //exponent is missing
2808      R->cf->is_field=1;
2809      R->cf->is_domain=1;
2810      R->cf->has_simple_Inverse=1;
2811      #endif
2812    }
2813  }
2814  else if (L->m[0].Typ()==LIST_CMD) // something complicated...
2815  {
2816    lists LL=(lists)L->m[0].Data();
2817
2818#ifdef HAVE_RINGS
2819    if (LL->m[0].Typ() == STRING_CMD) // 1st comes a string?
2820    {
2821      rComposeRing(LL, R); // Ring!?
2822    }
2823    else
2824#endif
2825    if (LL->nr < 3)
2826      rComposeC(LL,R); // R, long_R, long_C
2827    else
2828    {
2829      if (LL->m[0].Typ()==INT_CMD)
2830      {
2831        int ch = (int)(long)LL->m[0].Data();
2832        while ((ch!=fftable[is_gf_char]) && (fftable[is_gf_char])) is_gf_char++;
2833        if (fftable[is_gf_char]==0) is_gf_char=-1;
2834
2835        if(is_gf_char!= -1)
2836        {
2837          GFInfo param;
2838
2839          param.GFChar = ch;
2840          param.GFDegree = 1;
2841          param.GFPar_name = (const char*)(((lists)(LL->m[1].Data()))->m[0].Data());
2842
2843          // nfInitChar should be able to handle the case when ch is in fftables!
2844          R->cf = nInitChar(n_GF, (void*)&param);
2845        }
2846      }
2847
2848      if( R->cf == NULL )
2849      {
2850        ring extRing = rCompose((lists)L->m[0].Data(),FALSE,0x7fff);
2851
2852        if (extRing==NULL)
2853        {
2854          WerrorS("could not create the specified coefficient field");
2855          goto rCompose_err;
2856        }
2857
2858        if( extRing->qideal != NULL ) // Algebraic extension
2859        {
2860          AlgExtInfo extParam;
2861
2862          extParam.r = extRing;
2863
2864          R->cf = nInitChar(n_algExt, (void*)&extParam);
2865        }
2866        else // Transcendental extension
2867        {
2868          TransExtInfo extParam;
2869          extParam.r = extRing;
2870          assume( extRing->qideal == NULL );
2871
2872          R->cf = nInitChar(n_transExt, &extParam);
2873        }
2874      }
2875    }
2876  }
2877  else
2878  {
2879    WerrorS("coefficient field must be described by `int` or `list`");
2880    goto rCompose_err;
2881  }
2882
2883  if( R->cf == NULL )
2884  {
2885    WerrorS("could not create coefficient field described by the input!");
2886    goto rCompose_err;
2887  }
2888
2889  // ------------------------- VARS ---------------------------
2890  if (rComposeVar(L,R)) goto rCompose_err;
2891  // ------------------------ ORDER ------------------------------
2892  if (rComposeOrder(L,check_comp,R)) goto rCompose_err;
2893
2894  // ------------------------ ??????? --------------------
2895
2896  if (!isLetterplace) rRenameVars(R);
2897  #ifdef HAVE_SHIFTBBA
2898  else
2899  {
2900    R->isLPring=isLetterplace;
2901    R->ShortOut=FALSE;
2902    R->CanShortOut=FALSE;
2903  }
2904  #endif
2905  if (bitmask!=0x7fff) R->bitmask=bitmask*2;
2906  rComplete(R);
2907
2908  // ------------------------ Q-IDEAL ------------------------
2909
2910  if (L->m[3].Typ()==IDEAL_CMD)
2911  {
2912    ideal q=(ideal)L->m[3].Data();
2913    if (q->m[0]!=NULL)
2914    {
2915      if (R->cf != currRing->cf) //->cf->ch!=currRing->cf->ch)
2916      {
2917      #if 0
2918            WerrorS("coefficient fields must be equal if q-ideal !=0");
2919            goto rCompose_err;
2920      #else
2921        ring orig_ring=currRing;
2922        rChangeCurrRing(R);
2923        int *perm=NULL;
2924        int *par_perm=NULL;
2925        int par_perm_size=0;
2926        nMapFunc nMap;
2927
2928        if ((nMap=nSetMap(orig_ring->cf))==NULL)
2929        {
2930          if (rEqual(orig_ring,currRing))
2931          {
2932            nMap=n_SetMap(currRing->cf, currRing->cf);
2933          }
2934          else
2935          // Allow imap/fetch to be make an exception only for:
2936          if ( (rField_is_Q_a(orig_ring) &&  // Q(a..) -> Q(a..) || Q || Zp || Zp(a)
2937            (rField_is_Q(currRing) || rField_is_Q_a(currRing) ||
2938             rField_is_Zp(currRing) || rField_is_Zp_a(currRing) ||
2939             rField_is_Zn(currRing)))
2940           ||
2941           (rField_is_Zp_a(orig_ring) &&  // Zp(a..) -> Zp(a..) || Zp
2942            (rField_is_Zp(currRing, rInternalChar(orig_ring)) ||
2943             rField_is_Zp_a(currRing, rInternalChar(orig_ring)))) )
2944          {
2945            par_perm_size=rPar(orig_ring);
2946
2947//            if ((orig_ring->minpoly != NULL) || (orig_ring->qideal != NULL))
2948//              naSetChar(rInternalChar(orig_ring),orig_ring);
2949//            else ntSetChar(rInternalChar(orig_ring),orig_ring);
2950
2951            nSetChar(currRing->cf);
2952          }
2953          else
2954          {
2955            WerrorS("coefficient fields must be equal if q-ideal !=0");
2956            goto rCompose_err;
2957          }
2958        }
2959        perm=(int *)omAlloc0((orig_ring->N+1)*sizeof(int));
2960        if (par_perm_size!=0)
2961          par_perm=(int *)omAlloc0(par_perm_size*sizeof(int));
2962        int i;
2963        #if 0
2964        // use imap:
2965        maFindPerm(orig_ring->names,orig_ring->N,orig_ring->parameter,orig_ring->P,
2966          currRing->names,currRing->N,currRing->parameter, currRing->P,
2967          perm,par_perm, currRing->ch);
2968        #else
2969        // use fetch
2970        if ((rPar(orig_ring)>0) && (rPar(currRing)==0))
2971        {
2972          for(i=si_min(rPar(orig_ring),rVar(currRing))-1;i>=0;i--) par_perm[i]=i+1;
2973        }
2974        else if (par_perm_size!=0)
2975          for(i=si_min(rPar(orig_ring),rPar(currRing))-1;i>=0;i--) par_perm[i]=-(i+1);
2976        for(i=si_min(orig_ring->N,rVar(currRing));i>0;i--) perm[i]=i;
2977        #endif
2978        ideal dest_id=idInit(IDELEMS(q),1);
2979        for(i=IDELEMS(q)-1; i>=0; i--)
2980        {
2981          dest_id->m[i]=p_PermPoly(q->m[i],perm,orig_ring, currRing,nMap,
2982                                  par_perm,par_perm_size);
2983          //  PrintS("map:");pWrite(dest_id->m[i]);PrintLn();
2984          pTest(dest_id->m[i]);
2985        }
2986        R->qideal=dest_id;
2987        if (perm!=NULL)
2988          omFreeSize((ADDRESS)perm,(orig_ring->N+1)*sizeof(int));
2989        if (par_perm!=NULL)
2990          omFreeSize((ADDRESS)par_perm,par_perm_size*sizeof(int));
2991        rChangeCurrRing(orig_ring);
2992      #endif
2993      }
2994      else
2995        R->qideal=idrCopyR(q,currRing,R);
2996    }
2997  }
2998  else
2999  {
3000    WerrorS("q-ideal must be given as `ideal`");
3001    goto rCompose_err;
3002  }
3003
3004
3005  // ---------------------------------------------------------------
3006  #ifdef HAVE_PLURAL
3007  if (L->nr==5)
3008  {
3009    if (nc_CallPlural((matrix)L->m[4].Data(),
3010                      (matrix)L->m[5].Data(),
3011                      NULL,NULL,
3012                      R,
3013                      true, // !!!
3014                      true, false,
3015                      currRing, FALSE)) goto rCompose_err;
3016    // takes care about non-comm. quotient! i.e. calls "nc_SetupQuotient" due to last true
3017  }
3018  #endif
3019  return R;
3020
3021rCompose_err:
3022  if (R->N>0)
3023  {
3024    int i;
3025    if (R->names!=NULL)
3026    {
3027      i=R->N-1;
3028      while (i>=0) { omfree(R->names[i]); i--; }
3029      omFree(R->names);
3030    }
3031  }
3032  omfree(R->order);
3033  omfree(R->block0);
3034  omfree(R->block1);
3035  omfree(R->wvhdl);
3036  omFree(R);
3037  return NULL;
3038}
3039
3040// from matpol.cc
3041
3042/*2
3043* compute the jacobi matrix of an ideal
3044*/
3045BOOLEAN mpJacobi(leftv res,leftv a)
3046{
3047  int     i,j;
3048  matrix result;
3049  ideal id=(ideal)a->Data();
3050
3051  result =mpNew(IDELEMS(id),rVar(currRing));
3052  for (i=1; i<=IDELEMS(id); i++)
3053  {
3054    for (j=1; j<=rVar(currRing); j++)
3055    {
3056      MATELEM(result,i,j) = pDiff(id->m[i-1],j);
3057    }
3058  }
3059  res->data=(char *)result;
3060  return FALSE;
3061}
3062
3063/*2
3064* returns the Koszul-matrix of degree d of a vectorspace with dimension n
3065* uses the first n entrees of id, if id <> NULL
3066*/
3067BOOLEAN mpKoszul(leftv res,leftv c/*ip*/, leftv b/*in*/, leftv id)
3068{
3069  int n=(int)(long)b->Data();
3070  int d=(int)(long)c->Data();
3071  int     k,l,sign,row,col;
3072  matrix  result;
3073  ideal temp;
3074  BOOLEAN bo;
3075  poly    p;
3076
3077  if ((d>n) || (d<1) || (n<1))
3078  {
3079    res->data=(char *)mpNew(1,1);
3080    return FALSE;
3081  }
3082  int *choise = (int*)omAlloc(d*sizeof(int));
3083  if (id==NULL)
3084    temp=idMaxIdeal(1);
3085  else
3086    temp=(ideal)id->Data();
3087
3088  k = binom(n,d);
3089  l = k*d;
3090  l /= n-d+1;
3091  result =mpNew(l,k);
3092  col = 1;
3093  idInitChoise(d,1,n,&bo,choise);
3094  while (!bo)
3095  {
3096    sign = 1;
3097    for (l=1;l<=d;l++)
3098    {
3099      if (choise[l-1]<=IDELEMS(temp))
3100      {
3101        p = pCopy(temp->m[choise[l-1]-1]);
3102        if (sign == -1) p = pNeg(p);
3103        sign *= -1;
3104        row = idGetNumberOfChoise(l-1,d,1,n,choise);
3105        MATELEM(result,row,col) = p;
3106      }
3107    }
3108    col++;
3109    idGetNextChoise(d,n,&bo,choise);
3110  }
3111  omFreeSize(choise,d*sizeof(int));
3112  if (id==NULL) idDelete(&temp);
3113
3114  res->data=(char *)result;
3115  return FALSE;
3116}
3117
3118// from syz1.cc
3119/*2
3120* read out the Betti numbers from resolution
3121* (interpreter interface)
3122*/
3123BOOLEAN syBetti2(leftv res, leftv u, leftv w)
3124{
3125  syStrategy syzstr=(syStrategy)u->Data();
3126
3127  BOOLEAN minim=(int)(long)w->Data();
3128  int row_shift=0;
3129  int add_row_shift=0;
3130  intvec *weights=NULL;
3131  intvec *ww=(intvec *)atGet(u,"isHomog",INTVEC_CMD);
3132  if (ww!=NULL)
3133  {
3134     weights=ivCopy(ww);
3135     add_row_shift = ww->min_in();
3136     (*weights) -= add_row_shift;
3137  }
3138
3139  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift,weights);
3140  //row_shift += add_row_shift;
3141  //Print("row_shift=%d, add_row_shift=%d\n",row_shift,add_row_shift);
3142  atSet(res,omStrDup("rowShift"),(void*)(long)add_row_shift,INT_CMD);
3143
3144  return FALSE;
3145}
3146BOOLEAN syBetti1(leftv res, leftv u)
3147{
3148  sleftv tmp;
3149  memset(&tmp,0,sizeof(tmp));
3150  tmp.rtyp=INT_CMD;
3151  tmp.data=(void *)1;
3152  return syBetti2(res,u,&tmp);
3153}
3154
3155/*3
3156* converts a resolution into a list of modules
3157*/
3158lists syConvRes(syStrategy syzstr,BOOLEAN toDel,int add_row_shift)
3159{
3160  resolvente fullres = syzstr->fullres;
3161  resolvente minres = syzstr->minres;
3162
3163  const int length = syzstr->length;
3164
3165  if ((fullres==NULL) && (minres==NULL))
3166  {
3167    if (syzstr->hilb_coeffs==NULL)
3168    { // La Scala
3169      fullres = syReorder(syzstr->res, length, syzstr);
3170    }
3171    else
3172    { // HRES
3173      minres = syReorder(syzstr->orderedRes, length, syzstr);
3174      syKillEmptyEntres(minres, length);
3175    }
3176  }
3177
3178  resolvente tr;
3179  int typ0=IDEAL_CMD;
3180
3181  if (minres!=NULL)
3182    tr = minres;
3183  else
3184    tr = fullres;
3185
3186  resolvente trueres=NULL;
3187  intvec ** w=NULL;
3188
3189  if (length>0)
3190  {
3191    trueres = (resolvente)omAlloc0((length)*sizeof(ideal));
3192    for (int i=length-1;i>=0;i--)
3193    {
3194      if (tr[i]!=NULL)
3195      {
3196        trueres[i] = idCopy(tr[i]);
3197      }
3198    }
3199    if ( id_RankFreeModule(trueres[0], currRing) > 0)
3200      typ0 = MODUL_CMD;
3201    if (syzstr->weights!=NULL)
3202    {
3203      w = (intvec**)omAlloc0(length*sizeof(intvec*));
3204      for (int i=length-1;i>=0;i--)
3205      {
3206        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
3207      }
3208    }
3209  }
3210
3211  lists li = liMakeResolv(trueres, length, syzstr->list_length,typ0,
3212                          w, add_row_shift);
3213
3214  if (toDel)
3215    syKillComputation(syzstr);
3216  else
3217  {
3218    if( fullres != NULL && syzstr->fullres == NULL )
3219      syzstr->fullres = fullres;
3220
3221    if( minres != NULL && syzstr->minres == NULL )
3222      syzstr->minres = minres;
3223  }
3224  return li;
3225}
3226
3227/*3
3228* converts a list of modules into a resolution
3229*/
3230syStrategy syConvList(lists li)
3231{
3232  int typ0;
3233  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
3234
3235  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
3236  if (fr != NULL)
3237  {
3238
3239    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
3240    for (int i=result->length-1;i>=0;i--)
3241    {
3242      if (fr[i]!=NULL)
3243        result->fullres[i] = idCopy(fr[i]);
3244    }
3245    result->list_length=result->length;
3246    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
3247  }
3248  else
3249  {
3250    omFreeSize(result, sizeof(ssyStrategy));
3251    result = NULL;
3252  }
3253  return result;
3254}
3255
3256/*3
3257* converts a list of modules into a minimal resolution
3258*/
3259syStrategy syForceMin(lists li)
3260{
3261  int typ0;
3262  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
3263
3264  resolvente fr = liFindRes(li,&(result->length),&typ0);
3265  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
3266  for (int i=result->length-1;i>=0;i--)
3267  {
3268    if (fr[i]!=NULL)
3269      result->minres[i] = idCopy(fr[i]);
3270  }
3271  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
3272  return result;
3273}
3274// from weight.cc
3275BOOLEAN kWeight(leftv res,leftv id)
3276{
3277  ideal F=(ideal)id->Data();
3278  intvec * iv = new intvec(rVar(currRing));
3279  polyset s;
3280  int  sl, n, i;
3281  int  *x;
3282
3283  res->data=(char *)iv;
3284  s = F->m;
3285  sl = IDELEMS(F) - 1;
3286  n = rVar(currRing);
3287  double wNsqr = (double)2.0 / (double)n;
3288  wFunctional = wFunctionalBuch;
3289  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
3290  wCall(s, sl, x, wNsqr, currRing);
3291  for (i = n; i!=0; i--)
3292    (*iv)[i-1] = x[i + n + 1];
3293  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
3294  return FALSE;
3295}
3296
3297BOOLEAN kQHWeight(leftv res,leftv v)
3298{
3299  res->data=(char *)id_QHomWeight((ideal)v->Data(), currRing);
3300  if (res->data==NULL)
3301    res->data=(char *)new intvec(rVar(currRing));
3302  return FALSE;
3303}
3304/*==============================================================*/
3305// from clapsing.cc
3306#if 0
3307BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
3308{
3309  BOOLEAN b=singclap_factorize((poly)(u->CopyD()), &v, 0);
3310  res->data=(void *)b;
3311}
3312#endif
3313
3314BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
3315{
3316  res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
3317                  (poly)w->CopyD(), currRing);
3318  return errorreported;
3319}
3320
3321BOOLEAN jjCHARSERIES(leftv res, leftv u)
3322{
3323  res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
3324  return (res->data==NULL);
3325}
3326
3327// from semic.cc
3328#ifdef HAVE_SPECTRUM
3329
3330// ----------------------------------------------------------------------------
3331//  Initialize a  spectrum  deep from a  singular  lists
3332// ----------------------------------------------------------------------------
3333
3334void copy_deep( spectrum& spec, lists l )
3335{
3336    spec.mu = (int)(long)(l->m[0].Data( ));
3337    spec.pg = (int)(long)(l->m[1].Data( ));
3338    spec.n  = (int)(long)(l->m[2].Data( ));
3339
3340    spec.copy_new( spec.n );
3341
3342    intvec  *num = (intvec*)l->m[3].Data( );
3343    intvec  *den = (intvec*)l->m[4].Data( );
3344    intvec  *mul = (intvec*)l->m[5].Data( );
3345
3346    for( int i=0; i<spec.n; i++ )
3347    {
3348        spec.s[i] = (Rational)((*num)[i])/(Rational)((*den)[i]);
3349        spec.w[i] = (*mul)[i];
3350    }
3351}
3352
3353// ----------------------------------------------------------------------------
3354//  singular lists  constructor for  spectrum
3355// ----------------------------------------------------------------------------
3356
3357spectrum /*former spectrum::spectrum ( lists l )*/
3358spectrumFromList( lists l )
3359{
3360    spectrum result;
3361    copy_deep( result, l );
3362    return result;
3363}
3364
3365// ----------------------------------------------------------------------------
3366//  generate a Singular  lists  from a spectrum
3367// ----------------------------------------------------------------------------
3368
3369/* former spectrum::thelist ( void )*/
3370lists   getList( spectrum& spec )
3371{
3372    lists   L  = (lists)omAllocBin( slists_bin);
3373
3374    L->Init( 6 );
3375
3376    intvec            *num  = new intvec( spec.n );
3377    intvec            *den  = new intvec( spec.n );
3378    intvec            *mult = new intvec( spec.n );
3379
3380    for( int i=0; i<spec.n; i++ )
3381    {
3382        (*num) [i] = spec.s[i].get_num_si( );
3383        (*den) [i] = spec.s[i].get_den_si( );
3384        (*mult)[i] = spec.w[i];
3385    }
3386
3387    L->m[0].rtyp = INT_CMD;    //  milnor number
3388    L->m[1].rtyp = INT_CMD;    //  geometrical genus
3389    L->m[2].rtyp = INT_CMD;    //  # of spectrum numbers
3390    L->m[3].rtyp = INTVEC_CMD; //  numerators
3391    L->m[4].rtyp = INTVEC_CMD; //  denomiantors
3392    L->m[5].rtyp = INTVEC_CMD; //  multiplicities
3393
3394    L->m[0].data = (void*)(long)spec.mu;
3395    L->m[1].data = (void*)(long)spec.pg;
3396    L->m[2].data = (void*)(long)spec.n;
3397    L->m[3].data = (void*)num;
3398    L->m[4].data = (void*)den;
3399    L->m[5].data = (void*)mult;
3400
3401    return  L;
3402}
3403// from spectrum.cc
3404// ----------------------------------------------------------------------------
3405//  print out an error message for a spectrum list
3406// ----------------------------------------------------------------------------
3407
3408typedef enum
3409{
3410    semicOK,
3411    semicMulNegative,
3412
3413    semicListTooShort,
3414    semicListTooLong,
3415
3416    semicListFirstElementWrongType,
3417    semicListSecondElementWrongType,
3418    semicListThirdElementWrongType,
3419    semicListFourthElementWrongType,
3420    semicListFifthElementWrongType,
3421    semicListSixthElementWrongType,
3422
3423    semicListNNegative,
3424    semicListWrongNumberOfNumerators,
3425    semicListWrongNumberOfDenominators,
3426    semicListWrongNumberOfMultiplicities,
3427
3428    semicListMuNegative,
3429    semicListPgNegative,
3430    semicListNumNegative,
3431    semicListDenNegative,
3432    semicListMulNegative,
3433
3434    semicListNotSymmetric,
3435    semicListNotMonotonous,
3436
3437    semicListMilnorWrong,
3438    semicListPGWrong
3439
3440} semicState;
3441
3442void    list_error( semicState state )
3443{
3444    switch( state )
3445    {
3446        case semicListTooShort:
3447            WerrorS( "the list is too short" );
3448            break;
3449        case semicListTooLong:
3450            WerrorS( "the list is too long" );
3451            break;
3452
3453        case semicListFirstElementWrongType:
3454            WerrorS( "first element of the list should be int" );
3455            break;
3456        case semicListSecondElementWrongType:
3457            WerrorS( "second element of the list should be int" );
3458            break;
3459        case semicListThirdElementWrongType:
3460            WerrorS( "third element of the list should be int" );
3461            break;
3462        case semicListFourthElementWrongType:
3463            WerrorS( "fourth element of the list should be intvec" );
3464            break;
3465        case semicListFifthElementWrongType:
3466            WerrorS( "fifth element of the list should be intvec" );
3467            break;
3468        case semicListSixthElementWrongType:
3469            WerrorS( "sixth element of the list should be intvec" );
3470            break;
3471
3472        case semicListNNegative:
3473            WerrorS( "first element of the list should be positive" );
3474            break;
3475        case semicListWrongNumberOfNumerators:
3476            WerrorS( "wrong number of numerators" );
3477            break;
3478        case semicListWrongNumberOfDenominators:
3479            WerrorS( "wrong number of denominators" );
3480            break;
3481        case semicListWrongNumberOfMultiplicities:
3482            WerrorS( "wrong number of multiplicities" );
3483            break;
3484
3485        case semicListMuNegative:
3486            WerrorS( "the Milnor number should be positive" );
3487            break;
3488        case semicListPgNegative:
3489            WerrorS( "the geometrical genus should be nonnegative" );
3490            break;
3491        case semicListNumNegative:
3492            WerrorS( "all numerators should be positive" );
3493            break;
3494        case semicListDenNegative:
3495            WerrorS( "all denominators should be positive" );
3496            break;
3497        case semicListMulNegative:
3498            WerrorS( "all multiplicities should be positive" );
3499            break;
3500
3501        case semicListNotSymmetric:
3502            WerrorS( "it is not symmetric" );
3503            break;
3504        case semicListNotMonotonous:
3505            WerrorS( "it is not monotonous" );
3506            break;
3507
3508        case semicListMilnorWrong:
3509            WerrorS( "the Milnor number is wrong" );
3510            break;
3511        case semicListPGWrong:
3512            WerrorS( "the geometrical genus is wrong" );
3513            break;
3514
3515        default:
3516            WerrorS( "unspecific error" );
3517            break;
3518    }
3519}
3520// ----------------------------------------------------------------------------
3521//  this is the main spectrum computation function
3522// ----------------------------------------------------------------------------
3523
3524enum    spectrumState
3525{
3526    spectrumOK,
3527    spectrumZero,
3528    spectrumBadPoly,
3529    spectrumNoSingularity,
3530    spectrumNotIsolated,
3531    spectrumDegenerate,
3532    spectrumWrongRing,
3533    spectrumNoHC,
3534    spectrumUnspecErr
3535};
3536
3537// from splist.cc
3538// ----------------------------------------------------------------------------
3539//  Compute the spectrum of a  spectrumPolyList
3540// ----------------------------------------------------------------------------
3541
3542/* former spectrumPolyList::spectrum ( lists*, int) */
3543spectrumState   spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3544{
3545  spectrumPolyNode  **node = &speclist.root;
3546  spectrumPolyNode  *search;
3547
3548  poly              f,tmp;
3549  int               found,cmp;
3550
3551  Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3552                 ( fast==2 ? 2 : 1 ) );
3553
3554  Rational weight_prev( 0,1 );
3555
3556  int     mu = 0;          // the milnor number
3557  int     pg = 0;          // the geometrical genus
3558  int     n  = 0;          // number of different spectral numbers
3559  int     z  = 0;          // number of spectral number equal to smax
3560
3561  while( (*node)!=(spectrumPolyNode*)NULL &&
3562         ( fast==0 || (*node)->weight<=smax ) )
3563  {
3564        // ---------------------------------------
3565        //  determine the first normal form which
3566        //  contains the monomial  node->mon
3567        // ---------------------------------------
3568
3569    found  = FALSE;
3570    search = *node;
3571
3572    while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3573    {
3574      if( search->nf!=(poly)NULL )
3575      {
3576        f = search->nf;
3577
3578        do
3579        {
3580                    // --------------------------------
3581                    //  look for  (*node)->mon  in   f
3582                    // --------------------------------
3583
3584          cmp = pCmp( (*node)->mon,f );
3585
3586          if( cmp<0 )
3587          {
3588            f = pNext( f );
3589          }
3590          else if( cmp==0 )
3591          {
3592                        // -----------------------------
3593                        //  we have found a normal form
3594                        // -----------------------------
3595
3596            found = TRUE;
3597
3598                        //  normalize coefficient
3599
3600            number inv = nInvers( pGetCoeff( f ) );
3601            search->nf=__p_Mult_nn( search->nf,inv,currRing );
3602            nDelete( &inv );
3603
3604                        //  exchange  normal forms
3605
3606            tmp         = (*node)->nf;
3607            (*node)->nf = search->nf;
3608            search->nf  = tmp;
3609          }
3610        }
3611        while( cmp<0 && f!=(poly)NULL );
3612      }
3613      search = search->next;
3614    }
3615
3616    if( found==FALSE )
3617    {
3618            // ------------------------------------------------
3619            //  the weight of  node->mon  is a spectrum number
3620            // ------------------------------------------------
3621
3622      mu++;
3623
3624      if( (*node)->weight<=(Rational)1 )              pg++;
3625      if( (*node)->weight==smax )           z++;
3626      if( (*node)->weight>weight_prev )     n++;
3627
3628      weight_prev = (*node)->weight;
3629      node = &((*node)->next);
3630    }
3631    else
3632    {
3633            // -----------------------------------------------
3634            //  determine all other normal form which contain
3635            //  the monomial  node->mon
3636            //  replace for  node->mon  its normal form
3637            // -----------------------------------------------
3638
3639      while( search!=(spectrumPolyNode*)NULL )
3640      {
3641        if( search->nf!=(poly)NULL )
3642        {
3643          f = search->nf;
3644
3645          do
3646          {
3647                        // --------------------------------
3648                        //  look for  (*node)->mon  in   f
3649                        // --------------------------------
3650
3651            cmp = pCmp( (*node)->mon,f );
3652
3653            if( cmp<0 )
3654            {
3655              f = pNext( f );
3656            }
3657            else if( cmp==0 )
3658            {
3659              search->nf = pSub( search->nf,
3660                                 __pp_Mult_nn( (*node)->nf,pGetCoeff( f ),currRing ) );
3661              pNorm( search->nf );
3662            }
3663          }
3664          while( cmp<0 && f!=(poly)NULL );
3665        }
3666        search = search->next;
3667      }
3668      speclist.delete_node( node );
3669    }
3670
3671  }
3672
3673    // --------------------------------------------------------
3674    //  fast computation exploits the symmetry of the spectrum
3675    // --------------------------------------------------------
3676
3677  if( fast==2 )
3678  {
3679    mu = 2*mu - z;
3680    n  = ( z > 0 ? 2*n - 1 : 2*n );
3681  }
3682
3683    // --------------------------------------------------------
3684    //  compute the spectrum numbers with their multiplicities
3685    // --------------------------------------------------------
3686
3687  intvec            *nom  = new intvec( n );
3688  intvec            *den  = new intvec( n );
3689  intvec            *mult = new intvec( n );
3690
3691  int count         = 0;
3692  int multiplicity  = 1;
3693
3694  for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3695              ( fast==0 || search->weight<=smax );
3696       search=search->next )
3697  {
3698    if( search->next==(spectrumPolyNode*)NULL ||
3699        search->weight<search->next->weight )
3700    {
3701      (*nom) [count] = search->weight.get_num_si( );
3702      (*den) [count] = search->weight.get_den_si( );
3703      (*mult)[count] = multiplicity;
3704
3705      multiplicity=1;
3706      count++;
3707    }
3708    else
3709    {
3710      multiplicity++;
3711    }
3712  }
3713
3714    // --------------------------------------------------------
3715    //  fast computation exploits the symmetry of the spectrum
3716    // --------------------------------------------------------
3717
3718  if( fast==2 )
3719  {
3720    int n1,n2;
3721    for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3722    {
3723      (*nom) [n2] = rVar(currRing)*(*den)[n1]-(*nom)[n1];
3724      (*den) [n2] = (*den)[n1];
3725      (*mult)[n2] = (*mult)[n1];
3726    }
3727  }
3728
3729    // -----------------------------------
3730    //  test if the spectrum is symmetric
3731    // -----------------------------------
3732
3733  if( fast==0 || fast==1 )
3734  {
3735    int symmetric=TRUE;
3736
3737    for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3738    {
3739      if( (*mult)[n1]!=(*mult)[n2] ||
3740          (*den) [n1]!= (*den)[n2] ||
3741          (*nom)[n1]+(*nom)[n2]!=rVar(currRing)*(*den) [n1] )
3742      {
3743        symmetric = FALSE;
3744      }
3745    }
3746
3747    if( symmetric==FALSE )
3748    {
3749            // ---------------------------------------------
3750            //  the spectrum is not symmetric => degenerate
3751            //  principal part
3752            // ---------------------------------------------
3753
3754      *L = (lists)omAllocBin( slists_bin);
3755      (*L)->Init( 1 );
3756      (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3757      (*L)->m[0].data = (void*)(long)mu;
3758
3759      return spectrumDegenerate;
3760    }
3761  }
3762
3763  *L = (lists)omAllocBin( slists_bin);
3764
3765  (*L)->Init( 6 );
3766
3767  (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3768  (*L)->m[1].rtyp = INT_CMD;    //  geometrical genus
3769  (*L)->m[2].rtyp = INT_CMD;    //  number of spectrum values
3770  (*L)->m[3].rtyp = INTVEC_CMD; //  nominators
3771  (*L)->m[4].rtyp = INTVEC_CMD; //  denomiantors
3772  (*L)->m[5].rtyp = INTVEC_CMD; //  multiplicities
3773
3774  (*L)->m[0].data = (void*)(long)mu;
3775  (*L)->m[1].data = (void*)(long)pg;
3776  (*L)->m[2].data = (void*)(long)n;
3777  (*L)->m[3].data = (void*)nom;
3778  (*L)->m[4].data = (void*)den;
3779  (*L)->m[5].data = (void*)mult;
3780
3781  return  spectrumOK;
3782}
3783
3784spectrumState   spectrumCompute( poly h,lists *L,int fast )
3785{
3786  int i;
3787
3788  #ifdef SPECTRUM_DEBUG
3789  #ifdef SPECTRUM_PRINT
3790  #ifdef SPECTRUM_IOSTREAM
3791    cout << "spectrumCompute\n";
3792    if( fast==0 ) cout << "    no optimization" << endl;
3793    if( fast==1 ) cout << "    weight optimization" << endl;
3794    if( fast==2 ) cout << "    symmetry optimization" << endl;
3795  #else
3796    fputs( "spectrumCompute\n",stdout );
3797    if( fast==0 ) fputs( "    no optimization\n", stdout );
3798    if( fast==1 ) fputs( "    weight optimization\n", stdout );
3799    if( fast==2 ) fputs( "    symmetry optimization\n", stdout );
3800  #endif
3801  #endif
3802  #endif
3803
3804  // ----------------------
3805  //  check if  h  is zero
3806  // ----------------------
3807
3808  if( h==(poly)NULL )
3809  {
3810    return  spectrumZero;
3811  }
3812
3813  // ----------------------------------
3814  //  check if  h  has a constant term
3815  // ----------------------------------
3816
3817  if( hasConstTerm( h, currRing ) )
3818  {
3819    return  spectrumBadPoly;
3820  }
3821
3822  // --------------------------------
3823  //  check if  h  has a linear term
3824  // --------------------------------
3825
3826  if( hasLinearTerm( h, currRing ) )
3827  {
3828    *L = (lists)omAllocBin( slists_bin);
3829    (*L)->Init( 1 );
3830    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3831    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3832
3833    return  spectrumNoSingularity;
3834  }
3835
3836  // ----------------------------------
3837  //  compute the jacobi ideal of  (h)
3838  // ----------------------------------
3839
3840  ideal J = NULL;
3841  J = idInit( rVar(currRing),1 );
3842
3843  #ifdef SPECTRUM_DEBUG
3844  #ifdef SPECTRUM_PRINT
3845  #ifdef SPECTRUM_IOSTREAM
3846    cout << "\n   computing the Jacobi ideal...\n";
3847  #else
3848    fputs( "\n   computing the Jacobi ideal...\n",stdout );
3849  #endif
3850  #endif
3851  #endif
3852
3853  for( i=0; i<rVar(currRing); i++ )
3854  {
3855    J->m[i] = pDiff( h,i+1); //j );
3856
3857    #ifdef SPECTRUM_DEBUG
3858    #ifdef SPECTRUM_PRINT
3859    #ifdef SPECTRUM_IOSTREAM
3860      cout << "        ";
3861    #else
3862      fputs("        ", stdout );
3863    #endif
3864      pWrite( J->m[i] );
3865    #endif
3866    #endif
3867  }
3868
3869  // --------------------------------------------
3870  //  compute a standard basis  stdJ  of  jac(h)
3871  // --------------------------------------------
3872
3873  #ifdef SPECTRUM_DEBUG
3874  #ifdef SPECTRUM_PRINT
3875  #ifdef SPECTRUM_IOSTREAM
3876    cout << endl;
3877    cout << "    computing a standard basis..." << endl;
3878  #else
3879    fputs( "\n", stdout );
3880    fputs( "    computing a standard basis...\n", stdout );
3881  #endif
3882  #endif
3883  #endif
3884
3885  ideal stdJ = kStd(J,currRing->qideal,isNotHomog,NULL);
3886  idSkipZeroes( stdJ );
3887
3888  #ifdef SPECTRUM_DEBUG
3889  #ifdef SPECTRUM_PRINT
3890    for( i=0; i<IDELEMS(stdJ); i++ )
3891    {
3892      #ifdef SPECTRUM_IOSTREAM
3893        cout << "        ";
3894      #else
3895        fputs( "        ",stdout );
3896      #endif
3897
3898      pWrite( stdJ->m[i] );
3899    }
3900  #endif
3901  #endif
3902
3903  idDelete( &J );
3904
3905  // ------------------------------------------
3906  //  check if the  h  has a singularity
3907  // ------------------------------------------
3908
3909  if( hasOne( stdJ, currRing ) )
3910  {
3911    // -------------------------------
3912    //  h is smooth in the origin
3913    //  return only the Milnor number
3914    // -------------------------------
3915
3916    *L = (lists)omAllocBin( slists_bin);
3917    (*L)->Init( 1 );
3918    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3919    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3920
3921    return  spectrumNoSingularity;
3922  }
3923
3924  // ------------------------------------------
3925  //  check if the singularity  h  is isolated
3926  // ------------------------------------------
3927
3928  for( i=rVar(currRing); i>0; i-- )
3929  {
3930    if( hasAxis( stdJ,i, currRing )==FALSE )
3931    {
3932      return  spectrumNotIsolated;
3933    }
3934  }
3935
3936  // ------------------------------------------
3937  //  compute the highest corner  hc  of  stdJ
3938  // ------------------------------------------
3939
3940  #ifdef SPECTRUM_DEBUG
3941  #ifdef SPECTRUM_PRINT
3942  #ifdef SPECTRUM_IOSTREAM
3943    cout << "\n    computing the highest corner...\n";
3944  #else
3945    fputs( "\n    computing the highest corner...\n", stdout );
3946  #endif
3947  #endif
3948  #endif
3949
3950  poly hc = (poly)NULL;
3951
3952  scComputeHC( stdJ,currRing->qideal, 0,hc );
3953
3954  if( hc!=(poly)NULL )
3955  {
3956    pGetCoeff(hc) = nInit(1);
3957
3958    for( i=rVar(currRing); i>0; i-- )
3959    {
3960      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3961    }
3962    pSetm( hc );
3963  }
3964  else
3965  {
3966    return  spectrumNoHC;
3967  }
3968
3969  #ifdef SPECTRUM_DEBUG
3970  #ifdef SPECTRUM_PRINT
3971  #ifdef SPECTRUM_IOSTREAM
3972    cout << "       ";
3973  #else
3974    fputs( "       ", stdout );
3975  #endif
3976    pWrite( hc );
3977  #endif
3978  #endif
3979
3980  // ----------------------------------------
3981  //  compute the Newton polygon  nph  of  h
3982  // ----------------------------------------
3983
3984  #ifdef SPECTRUM_DEBUG
3985  #ifdef SPECTRUM_PRINT
3986  #ifdef SPECTRUM_IOSTREAM
3987    cout << "\n    computing the newton polygon...\n";
3988  #else
3989    fputs( "\n    computing the newton polygon...\n", stdout );
3990  #endif
3991  #endif
3992  #endif
3993
3994  newtonPolygon nph( h, currRing );
3995
3996  #ifdef SPECTRUM_DEBUG
3997  #ifdef SPECTRUM_PRINT
3998    cout << nph;
3999  #endif
4000  #endif
4001
4002  // -----------------------------------------------
4003  //  compute the weight corner  wc  of  (stdj,nph)
4004  // -----------------------------------------------
4005
4006  #ifdef SPECTRUM_DEBUG
4007  #ifdef SPECTRUM_PRINT
4008  #ifdef SPECTRUM_IOSTREAM
4009    cout << "\n    computing the weight corner...\n";
4010  #else
4011    fputs( "\n    computing the weight corner...\n", stdout );
4012  #endif
4013  #endif
4014  #endif
4015
4016  poly    wc = ( fast==0 ? pCopy( hc ) :
4017               ( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
4018              /* fast==2 */computeWC( nph,
4019                      ((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
4020
4021  #ifdef SPECTRUM_DEBUG
4022  #ifdef SPECTRUM_PRINT
4023  #ifdef SPECTRUM_IOSTREAM
4024    cout << "        ";
4025  #else
4026    fputs( "        ", stdout );
4027  #endif
4028    pWrite( wc );
4029  #endif
4030  #endif
4031
4032  // -------------
4033  //  compute  NF
4034  // -------------
4035
4036  #ifdef SPECTRUM_DEBUG
4037  #ifdef SPECTRUM_PRINT
4038  #ifdef SPECTRUM_IOSTREAM
4039    cout << "\n    computing NF...\n" << endl;
4040  #else
4041    fputs( "\n    computing NF...\n", stdout );
4042  #endif
4043  #endif
4044  #endif
4045
4046  spectrumPolyList NF( &nph );
4047
4048  computeNF( stdJ,hc,wc,&NF, currRing );
4049
4050  #ifdef SPECTRUM_DEBUG
4051  #ifdef SPECTRUM_PRINT
4052    cout << NF;
4053  #ifdef SPECTRUM_IOSTREAM
4054    cout << endl;
4055  #else
4056    fputs( "\n", stdout );
4057  #endif
4058  #endif
4059  #endif
4060
4061  // ----------------------------
4062  //  compute the spectrum of  h
4063  // ----------------------------
4064//  spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
4065
4066  return spectrumStateFromList(NF, L, fast );
4067}
4068
4069// ----------------------------------------------------------------------------
4070//  this procedure is called from the interpreter
4071// ----------------------------------------------------------------------------
4072//  first  = polynomial
4073//  result = list of spectrum numbers
4074// ----------------------------------------------------------------------------
4075
4076void spectrumPrintError(spectrumState state)
4077{
4078  switch( state )
4079  {
4080    case spectrumZero:
4081      WerrorS( "polynomial is zero" );
4082      break;
4083    case spectrumBadPoly:
4084      WerrorS( "polynomial has constant term" );
4085      break;
4086    case spectrumNoSingularity:
4087      WerrorS( "not a singularity" );
4088      break;
4089    case spectrumNotIsolated:
4090      WerrorS( "the singularity is not isolated" );
4091      break;
4092    case spectrumNoHC:
4093      WerrorS( "highest corner cannot be computed" );
4094      break;
4095    case spectrumDegenerate:
4096      WerrorS( "principal part is degenerate" );
4097      break;
4098    case spectrumOK:
4099      break;
4100
4101    default:
4102      WerrorS( "unknown error occurred" );
4103      break;
4104  }
4105}
4106
4107BOOLEAN spectrumProc( leftv result,leftv first )
4108{
4109  spectrumState state = spectrumOK;
4110
4111  // -------------------
4112  //  check consistency
4113  // -------------------
4114
4115  //  check for a local ring
4116
4117  if( !ringIsLocal(currRing ) )
4118  {
4119    WerrorS( "only works for local orderings" );
4120    state = spectrumWrongRing;
4121  }
4122
4123  //  no quotient rings are allowed
4124
4125  else if( currRing->qideal != NULL )
4126  {
4127    WerrorS( "does not work in quotient rings" );
4128    state = spectrumWrongRing;
4129  }
4130  else
4131  {
4132    lists   L    = (lists)NULL;
4133    int     flag = 1; // weight corner optimization is safe
4134
4135    state = spectrumCompute( (poly)first->Data( ),&L,flag );
4136
4137    if( state==spectrumOK )
4138    {
4139      result->rtyp = LIST_CMD;
4140      result->data = (char*)L;
4141    }
4142    else
4143    {
4144      spectrumPrintError(state);
4145    }
4146  }
4147
4148  return  (state!=spectrumOK);
4149}
4150
4151// ----------------------------------------------------------------------------
4152//  this procedure is called from the interpreter
4153// ----------------------------------------------------------------------------
4154//  first  = polynomial
4155//  result = list of spectrum numbers
4156// ----------------------------------------------------------------------------
4157
4158BOOLEAN spectrumfProc( leftv result,leftv first )
4159{
4160  spectrumState state = spectrumOK;
4161
4162  // -------------------
4163  //  check consistency
4164  // -------------------
4165
4166  //  check for a local polynomial ring
4167
4168  if( currRing->OrdSgn != -1 )
4169  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
4170  // or should we use:
4171  //if( !ringIsLocal( ) )
4172  {
4173    WerrorS( "only works for local orderings" );
4174    state = spectrumWrongRing;
4175  }
4176  else if( currRing->qideal != NULL )
4177  {
4178    WerrorS( "does not work in quotient rings" );
4179    state = spectrumWrongRing;
4180  }
4181  else
4182  {
4183    lists   L    = (lists)NULL;
4184    int     flag = 2; // symmetric optimization
4185
4186    state = spectrumCompute( (poly)first->Data( ),&L,flag );
4187
4188    if( state==spectrumOK )
4189    {
4190      result->rtyp = LIST_CMD;
4191      result->data = (char*)L;
4192    }
4193    else
4194    {
4195      spectrumPrintError(state);
4196    }
4197  }
4198
4199  return  (state!=spectrumOK);
4200}
4201
4202// ----------------------------------------------------------------------------
4203//  check if a list is a spectrum
4204//  check for:
4205//      list has 6 elements
4206//      1st element is int (mu=Milnor number)
4207//      2nd element is int (pg=geometrical genus)
4208//      3rd element is int (n =number of different spectrum numbers)
4209//      4th element is intvec (num=numerators)
4210//      5th element is intvec (den=denomiantors)
4211//      6th element is intvec (mul=multiplicities)
4212//      exactly n numerators
4213//      exactly n denominators
4214//      exactly n multiplicities
4215//      mu>0
4216//      pg>=0
4217//      n>0
4218//      num>0
4219//      den>0
4220//      mul>0
4221//      symmetriy with respect to numberofvariables/2
4222//      monotony
4223//      mu = sum of all multiplicities
4224//      pg = sum of all multiplicities where num/den<=1
4225// ----------------------------------------------------------------------------
4226
4227semicState  list_is_spectrum( lists l )
4228{
4229    // -------------------
4230    //  check list length
4231    // -------------------
4232
4233    if( l->nr < 5 )
4234    {
4235        return  semicListTooShort;
4236    }
4237    else if( l->nr > 5 )
4238    {
4239        return  semicListTooLong;
4240    }
4241
4242    // -------------
4243    //  check types
4244    // -------------
4245
4246    if( l->m[0].rtyp != INT_CMD )
4247    {
4248        return  semicListFirstElementWrongType;
4249    }
4250    else if( l->m[1].rtyp != INT_CMD )
4251    {
4252        return  semicListSecondElementWrongType;
4253    }
4254    else if( l->m[2].rtyp != INT_CMD )
4255    {
4256        return  semicListThirdElementWrongType;
4257    }
4258    else if( l->m[3].rtyp != INTVEC_CMD )
4259    {
4260        return  semicListFourthElementWrongType;
4261    }
4262    else if( l->m[4].rtyp != INTVEC_CMD )
4263    {
4264        return  semicListFifthElementWrongType;
4265    }
4266    else if( l->m[5].rtyp != INTVEC_CMD )
4267    {
4268        return  semicListSixthElementWrongType;
4269    }
4270
4271    // -------------------------
4272    //  check number of entries
4273    // -------------------------
4274
4275    int     mu = (int)(long)(l->m[0].Data( ));
4276    int     pg = (int)(long)(l->m[1].Data( ));
4277    int     n  = (int)(long)(l->m[2].Data( ));
4278
4279    if( n <= 0 )
4280    {
4281        return  semicListNNegative;
4282    }
4283
4284    intvec  *num = (intvec*)l->m[3].Data( );
4285    intvec  *den = (intvec*)l->m[4].Data( );
4286    intvec  *mul = (intvec*)l->m[5].Data( );
4287
4288    if( n != num->length( ) )
4289    {
4290        return  semicListWrongNumberOfNumerators;
4291    }
4292    else if( n != den->length( ) )
4293    {
4294        return  semicListWrongNumberOfDenominators;
4295    }
4296    else if( n != mul->length( ) )
4297    {
4298        return  semicListWrongNumberOfMultiplicities;
4299    }
4300
4301    // --------
4302    //  values
4303    // --------
4304
4305    if( mu <= 0 )
4306    {
4307        return  semicListMuNegative;
4308    }
4309    if( pg < 0 )
4310    {
4311        return  semicListPgNegative;
4312    }
4313
4314    int i;
4315
4316    for( i=0; i<n; i++ )
4317    {
4318        if( (*num)[i] <= 0 )
4319        {
4320            return  semicListNumNegative;
4321        }
4322        if( (*den)[i] <= 0 )
4323        {
4324            return  semicListDenNegative;
4325        }
4326        if( (*mul)[i] <= 0 )
4327        {
4328            return  semicListMulNegative;
4329        }
4330    }
4331
4332    // ----------------
4333    //  check symmetry
4334    // ----------------
4335
4336    int     j;
4337
4338    for( i=0, j=n-1; i<=j; i++,j-- )
4339    {
4340        if( (*num)[i] != rVar(currRing)*((*den)[i]) - (*num)[j] ||
4341            (*den)[i] != (*den)[j] ||
4342            (*mul)[i] != (*mul)[j] )
4343        {
4344            return  semicListNotSymmetric;
4345        }
4346    }
4347
4348    // ----------------
4349    //  check monotony
4350    // ----------------
4351
4352    for( i=0, j=1; i<n/2; i++,j++ )
4353    {
4354        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
4355        {
4356            return  semicListNotMonotonous;
4357        }
4358    }
4359
4360    // ---------------------
4361    //  check Milnor number
4362    // ---------------------
4363
4364    for( mu=0, i=0; i<n; i++ )
4365    {
4366        mu += (*mul)[i];
4367    }
4368
4369    if( mu != (int)(long)(l->m[0].Data( )) )
4370    {
4371        return  semicListMilnorWrong;
4372    }
4373
4374    // -------------------------
4375    //  check geometrical genus
4376    // -------------------------
4377
4378    for( pg=0, i=0; i<n; i++ )
4379    {
4380        if( (*num)[i]<=(*den)[i] )
4381        {
4382            pg += (*mul)[i];
4383        }
4384    }
4385
4386    if( pg != (int)(long)(l->m[1].Data( )) )
4387    {
4388        return  semicListPGWrong;
4389    }
4390
4391    return  semicOK;
4392}
4393
4394// ----------------------------------------------------------------------------
4395//  this procedure is called from the interpreter
4396// ----------------------------------------------------------------------------
4397//  first  = list of spectrum numbers
4398//  second = list of spectrum numbers
4399//  result = sum of the two lists
4400// ----------------------------------------------------------------------------
4401
4402BOOLEAN spaddProc( leftv result,leftv first,leftv second )
4403{
4404    semicState  state;
4405
4406    // -----------------
4407    //  check arguments
4408    // -----------------
4409
4410    lists l1 = (lists)first->Data( );
4411    lists l2 = (lists)second->Data( );
4412
4413    if( (state=list_is_spectrum( l1 )) != semicOK )
4414    {
4415        WerrorS( "first argument is not a spectrum:" );
4416        list_error( state );
4417    }
4418    else if( (state=list_is_spectrum( l2 )) != semicOK )
4419    {
4420        WerrorS( "second argument is not a spectrum:" );
4421        list_error( state );
4422    }
4423    else
4424    {
4425        spectrum s1= spectrumFromList ( l1 );
4426        spectrum s2= spectrumFromList ( l2 );
4427        spectrum sum( s1+s2 );
4428
4429        result->rtyp = LIST_CMD;
4430        result->data = (char*)(getList(sum));
4431    }
4432
4433    return  (state!=semicOK);
4434}
4435
4436// ----------------------------------------------------------------------------
4437//  this procedure is called from the interpreter
4438// ----------------------------------------------------------------------------
4439//  first  = list of spectrum numbers
4440//  second = integer
4441//  result = the multiple of the first list by the second factor
4442// ----------------------------------------------------------------------------
4443
4444BOOLEAN spmulProc( leftv result,leftv first,leftv second )
4445{
4446    semicState  state;
4447
4448    // -----------------
4449    //  check arguments
4450    // -----------------
4451
4452    lists   l = (lists)first->Data( );
4453    int     k = (int)(long)second->Data( );
4454
4455    if( (state=list_is_spectrum( l ))!=semicOK )
4456    {
4457        WerrorS( "first argument is not a spectrum" );
4458        list_error( state );
4459    }
4460    else if( k < 0 )
4461    {
4462        WerrorS( "second argument should be positive" );
4463        state = semicMulNegative;
4464    }
4465    else
4466    {
4467        spectrum s= spectrumFromList( l );
4468        spectrum product( k*s );
4469
4470        result->rtyp = LIST_CMD;
4471        result->data = (char*)getList(product);
4472    }
4473
4474    return  (state!=semicOK);
4475}
4476
4477// ----------------------------------------------------------------------------
4478//  this procedure is called from the interpreter
4479// ----------------------------------------------------------------------------
4480//  first  = list of spectrum numbers
4481//  second = list of spectrum numbers
4482//  result = semicontinuity index
4483// ----------------------------------------------------------------------------
4484
4485BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
4486{
4487  semicState  state;
4488  BOOLEAN qh=(((int)(long)w->Data())==1);
4489
4490  // -----------------
4491  //  check arguments
4492  // -----------------
4493
4494  lists l1 = (lists)u->Data( );
4495  lists l2 = (lists)v->Data( );
4496
4497  if( (state=list_is_spectrum( l1 ))!=semicOK )
4498  {
4499    WerrorS( "first argument is not a spectrum" );
4500    list_error( state );
4501  }
4502  else if( (state=list_is_spectrum( l2 ))!=semicOK )
4503  {
4504    WerrorS( "second argument is not a spectrum" );
4505    list_error( state );
4506  }
4507  else
4508  {
4509    spectrum s1= spectrumFromList( l1 );
4510    spectrum s2= spectrumFromList( l2 );
4511
4512    res->rtyp = INT_CMD;
4513    if (qh)
4514      res->data = (void*)(long)(s1.mult_spectrumh( s2 ));
4515    else
4516      res->data = (void*)(long)(s1.mult_spectrum( s2 ));
4517  }
4518
4519  // -----------------
4520  //  check status
4521  // -----------------
4522
4523  return  (state!=semicOK);
4524}
4525BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
4526{
4527  sleftv tmp;
4528  memset(&tmp,0,sizeof(tmp));
4529  tmp.rtyp=INT_CMD;
4530  /* tmp.data = (void *)0;  -- done by memset */
4531
4532  return  semicProc3(res,u,v,&tmp);
4533}
4534
4535#endif
4536
4537BOOLEAN loNewtonP( leftv res, leftv arg1 )
4538{
4539  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4540  return FALSE;
4541}
4542
4543BOOLEAN loSimplex( leftv res, leftv args )
4544{
4545  if ( !(rField_is_long_R(currRing)) )
4546  {
4547    WerrorS("Ground field not implemented!");
4548    return TRUE;
4549  }
4550
4551  simplex * LP;
4552  matrix m;
4553
4554  leftv v= args;
4555  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4556    return TRUE;
4557  else
4558    m= (matrix)(v->CopyD());
4559
4560  LP = new simplex(MATROWS(m),MATCOLS(m));
4561  LP->mapFromMatrix(m);
4562
4563  v= v->next;
4564  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
4565    return TRUE;
4566  else
4567    LP->m= (int)(long)(v->Data());
4568
4569  v= v->next;
4570  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
4571    return TRUE;
4572  else
4573    LP->n= (int)(long)(v->Data());
4574
4575  v= v->next;
4576  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
4577    return TRUE;
4578  else
4579    LP->m1= (int)(long)(v->Data());
4580
4581  v= v->next;
4582  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
4583    return TRUE;
4584  else
4585    LP->m2= (int)(long)(v->Data());
4586
4587  v= v->next;
4588  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
4589    return TRUE;
4590  else
4591    LP->m3= (int)(long)(v->Data());
4592
4593#ifdef mprDEBUG_PROT
4594  Print("m (constraints) %d\n",LP->m);
4595  Print("n (columns) %d\n",LP->n);
4596  Print("m1 (<=) %d\n",LP->m1);
4597  Print("m2 (>=) %d\n",LP->m2);
4598  Print("m3 (==) %d\n",LP->m3);
4599#endif
4600
4601  LP->compute();
4602
4603  lists lres= (lists)omAlloc( sizeof(slists) );
4604  lres->Init( 6 );
4605
4606  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4607  lres->m[0].data=(void*)LP->mapToMatrix(m);
4608
4609  lres->m[1].rtyp= INT_CMD;   // found a solution?
4610  lres->m[1].data=(void*)(long)LP->icase;
4611
4612  lres->m[2].rtyp= INTVEC_CMD;
4613  lres->m[2].data=(void*)LP->posvToIV();
4614
4615  lres->m[3].rtyp= INTVEC_CMD;
4616  lres->m[3].data=(void*)LP->zrovToIV();
4617
4618  lres->m[4].rtyp= INT_CMD;
4619  lres->m[4].data=(void*)(long)LP->m;
4620
4621  lres->m[5].rtyp= INT_CMD;
4622  lres->m[5].data=(void*)(long)LP->n;
4623
4624  res->data= (void*)lres;
4625
4626  return FALSE;
4627}
4628
4629BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4630{
4631  ideal gls = (ideal)(arg1->Data());
4632  int imtype= (int)(long)arg2->Data();
4633
4634  uResultant::resMatType mtype= determineMType( imtype );
4635
4636  // check input ideal ( = polynomial system )
4637  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4638  {
4639    return TRUE;
4640  }
4641
4642  uResultant *resMat= new uResultant( gls, mtype, false );
4643  if (resMat!=NULL)
4644  {
4645    res->rtyp = MODUL_CMD;
4646    res->data= (void*)resMat->accessResMat()->getMatrix();
4647    if (!errorreported) delete resMat;
4648  }
4649  return errorreported;
4650}
4651
4652BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4653{
4654
4655  poly gls;
4656  gls= (poly)(arg1->Data());
4657  int howclean= (int)(long)arg3->Data();
4658
4659  if ( !(rField_is_R(currRing) ||
4660         rField_is_Q(currRing) ||
4661         rField_is_long_R(currRing) ||
4662         rField_is_long_C(currRing)) )
4663  {
4664    WerrorS("Ground field not implemented!");
4665    return TRUE;
4666  }
4667
4668  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4669                          rField_is_long_C(currRing)) )
4670  {
4671    unsigned long int ii = (unsigned long int)arg2->Data();
4672    setGMPFloatDigits( ii, ii );
4673  }
4674
4675  if ( gls == NULL || pIsConstant( gls ) )
4676  {
4677    WerrorS("Input polynomial is constant!");
4678    return TRUE;
4679  }
4680
4681  int ldummy;
4682  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4683  int i,vpos=0;
4684  poly piter;
4685  lists elist;
4686  lists rlist;
4687
4688  elist= (lists)omAlloc( sizeof(slists) );
4689  elist->Init( 0 );
4690
4691  if ( rVar(currRing) > 1 )
4692  {
4693    piter= gls;
4694    for ( i= 1; i <= rVar(currRing); i++ )
4695      if ( pGetExp( piter, i ) )
4696      {
4697        vpos= i;
4698        break;
4699      }
4700    while ( piter )
4701    {
4702      for ( i= 1; i <= rVar(currRing); i++ )
4703        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4704        {
4705          WerrorS("The input polynomial must be univariate!");
4706          return TRUE;
4707        }
4708      pIter( piter );
4709    }
4710  }
4711
4712  rootContainer * roots= new rootContainer();
4713  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4714  piter= gls;
4715  for ( i= deg; i >= 0; i-- )
4716  {
4717    if ( piter && pTotaldegree(piter) == i )
4718    {
4719      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4720      //nPrint( pcoeffs[i] );PrintS("  ");
4721      pIter( piter );
4722    }
4723    else
4724    {
4725      pcoeffs[i]= nInit(0);
4726    }
4727  }
4728
4729#ifdef mprDEBUG_PROT
4730  for (i=deg; i >= 0; i--)
4731  {
4732    nPrint( pcoeffs[i] );PrintS("  ");
4733  }
4734  PrintLn();
4735#endif
4736
4737  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4738  roots->solver( howclean );
4739
4740  int elem= roots->getAnzRoots();
4741  char *dummy;
4742  int j;
4743
4744  rlist= (lists)omAlloc( sizeof(slists) );
4745  rlist->Init( elem );
4746
4747  if (rField_is_long_C(currRing))
4748  {
4749    for ( j= 0; j < elem; j++ )
4750    {
4751      rlist->m[j].rtyp=NUMBER_CMD;
4752      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4753      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4754    }
4755  }
4756  else
4757  {
4758    for ( j= 0; j < elem; j++ )
4759    {
4760      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4761      rlist->m[j].rtyp=STRING_CMD;
4762      rlist->m[j].data=(void *)dummy;
4763    }
4764  }
4765
4766  elist->Clean();
4767  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4768
4769  // this is (via fillContainer) the same data as in root
4770  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4771  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4772
4773  delete roots;
4774
4775  res->rtyp= LIST_CMD;
4776  res->data= (void*)rlist;
4777
4778  return FALSE;
4779}
4780
4781BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4782{
4783  int i;
4784  ideal p,w;
4785  p= (ideal)arg1->Data();
4786  w= (ideal)arg2->Data();
4787
4788  // w[0] = f(p^0)
4789  // w[1] = f(p^1)
4790  // ...
4791  // p can be a vector of numbers (multivariate polynom)
4792  //   or one number (univariate polynom)
4793  // tdg = deg(f)
4794
4795  int n= IDELEMS( p );
4796  int m= IDELEMS( w );
4797  int tdg= (int)(long)arg3->Data();
4798
4799  res->data= (void*)NULL;
4800
4801  // check the input
4802  if ( tdg < 1 )
4803  {
4804    WerrorS("Last input parameter must be > 0!");
4805    return TRUE;
4806  }
4807  if ( n != rVar(currRing) )
4808  {
4809    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4810    return TRUE;
4811  }
4812  if ( m != (int)pow((double)tdg+1,(double)n) )
4813  {
4814    Werror("Size of second input ideal must be equal to %d!",
4815      (int)pow((double)tdg+1,(double)n));
4816    return TRUE;
4817  }
4818  if ( !(rField_is_Q(currRing) /* ||
4819         rField_is_R() || rField_is_long_R() ||
4820         rField_is_long_C()*/ ) )
4821         {
4822    WerrorS("Ground field not implemented!");
4823    return TRUE;
4824  }
4825
4826  number tmp;
4827  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4828  for ( i= 0; i < n; i++ )
4829  {
4830    pevpoint[i]=nInit(0);
4831    if (  (p->m)[i] )
4832    {
4833      tmp = pGetCoeff( (p->m)[i] );
4834      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4835      {
4836        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4837        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4838        return TRUE;
4839      }
4840    } else tmp= NULL;
4841    if ( !nIsZero(tmp) )
4842    {
4843      if ( !pIsConstant((p->m)[i]))
4844      {
4845        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4846        WerrorS("Elements of first input ideal must be numbers!");
4847        return TRUE;
4848      }
4849      pevpoint[i]= nCopy( tmp );
4850    }
4851  }
4852
4853  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4854  for ( i= 0; i < m; i++ )
4855  {
4856    wresults[i]= nInit(0);
4857    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4858    {
4859      if ( !pIsConstant((w->m)[i]))
4860      {
4861        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4862        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4863        WerrorS("Elements of second input ideal must be numbers!");
4864        return TRUE;
4865      }
4866      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4867    }
4868  }
4869
4870  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4871  number *ncpoly= vm.interpolateDense( wresults );
4872  // do not free ncpoly[]!!
4873  poly rpoly= vm.numvec2poly( ncpoly );
4874
4875  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4876  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4877
4878  res->data= (void*)rpoly;
4879  return FALSE;
4880}
4881
4882BOOLEAN nuUResSolve( leftv res, leftv args )
4883{
4884  leftv v= args;
4885
4886  ideal gls;
4887  int imtype;
4888  int howclean;
4889
4890  // get ideal
4891  if ( v->Typ() != IDEAL_CMD )
4892    return TRUE;
4893  else gls= (ideal)(v->Data());
4894  v= v->next;
4895
4896  // get resultant matrix type to use (0,1)
4897  if ( v->Typ() != INT_CMD )
4898    return TRUE;
4899  else imtype= (int)(long)v->Data();
4900  v= v->next;
4901
4902  if (imtype==0)
4903  {
4904    ideal test_id=idInit(1,1);
4905    int j;
4906    for(j=IDELEMS(gls)-1;j>=0;j--)
4907    {
4908      if (gls->m[j]!=NULL)
4909      {
4910        test_id->m[0]=gls->m[j];
4911        intvec *dummy_w=id_QHomWeight(test_id, currRing);
4912        if (dummy_w!=NULL)
4913        {
4914          WerrorS("Newton polytope not of expected dimension");
4915          delete dummy_w;
4916          return TRUE;
4917        }
4918      }
4919    }
4920  }
4921
4922  // get and set precision in digits ( > 0 )
4923  if ( v->Typ() != INT_CMD )
4924    return TRUE;
4925  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4926                          rField_is_long_C(currRing)) )
4927  {
4928    unsigned long int ii=(unsigned long int)v->Data();
4929    setGMPFloatDigits( ii, ii );
4930  }
4931  v= v->next;
4932
4933  // get interpolation steps (0,1,2)
4934  if ( v->Typ() != INT_CMD )
4935    return TRUE;
4936  else howclean= (int)(long)v->Data();
4937
4938  uResultant::resMatType mtype= determineMType( imtype );
4939  int i,count;
4940  lists listofroots= NULL;
4941  number smv= NULL;
4942  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4943
4944  //emptylist= (lists)omAlloc( sizeof(slists) );
4945  //emptylist->Init( 0 );
4946
4947  //res->rtyp = LIST_CMD;
4948  //res->data= (void *)emptylist;
4949
4950  // check input ideal ( = polynomial system )
4951  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4952  {
4953    return TRUE;
4954  }
4955
4956  uResultant * ures;
4957  rootContainer ** iproots;
4958  rootContainer ** muiproots;
4959  rootArranger * arranger;
4960
4961  // main task 1: setup of resultant matrix
4962  ures= new uResultant( gls, mtype );
4963  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4964  {
4965    WerrorS("Error occurred during matrix setup!");
4966    return TRUE;
4967  }
4968
4969  // if dense resultant, check if minor nonsingular
4970  if ( mtype == uResultant::denseResMat )
4971  {
4972    smv= ures->accessResMat()->getSubDet();
4973#ifdef mprDEBUG_PROT
4974    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4975#endif
4976    if ( nIsZero(smv) )
4977    {
4978      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4979      return TRUE;
4980    }
4981  }
4982
4983  // main task 2: Interpolate specialized resultant polynomials
4984  if ( interpolate_det )
4985    iproots= ures->interpolateDenseSP( false, smv );
4986  else
4987    iproots= ures->specializeInU( false, smv );
4988
4989  // main task 3: Interpolate specialized resultant polynomials
4990  if ( interpolate_det )
4991    muiproots= ures->interpolateDenseSP( true, smv );
4992  else
4993    muiproots= ures->specializeInU( true, smv );
4994
4995#ifdef mprDEBUG_PROT
4996  int c= iproots[0]->getAnzElems();
4997  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4998  c= muiproots[0]->getAnzElems();
4999  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5000#endif
5001
5002  // main task 4: Compute roots of specialized polys and match them up
5003  arranger= new rootArranger( iproots, muiproots, howclean );
5004  arranger->solve_all();
5005
5006  // get list of roots
5007  if ( arranger->success() )
5008  {
5009    arranger->arrange();
5010    listofroots= listOfRoots(arranger, gmp_output_digits );
5011  }
5012  else
5013  {
5014    WerrorS("Solver was unable to find any roots!");
5015    return TRUE;
5016  }
5017
5018  // free everything
5019  count= iproots[0]->getAnzElems();
5020  for (i=0; i < count; i++) delete iproots[i];
5021  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5022  count= muiproots[0]->getAnzElems();
5023  for (i=0; i < count; i++) delete muiproots[i];
5024  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5025
5026  delete ures;
5027  delete arranger;
5028  nDelete( &smv );
5029
5030  res->data= (void *)listofroots;
5031
5032  //emptylist->Clean();
5033  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5034
5035  return FALSE;
5036}
5037
5038// from mpr_numeric.cc
5039lists listOfRoots( rootArranger* self, const unsigned int oprec )
5040{
5041  int i,j;
5042  int count= self->roots[0]->getAnzRoots(); // number of roots
5043  int elem= self->roots[0]->getAnzElems();  // number of koordinates per root
5044
5045  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
5046
5047  if ( self->found_roots )
5048  {
5049    listofroots->Init( count );
5050
5051    for (i=0; i < count; i++)
5052    {
5053      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
5054      onepoint->Init(elem);
5055      for ( j= 0; j < elem; j++ )
5056      {
5057        if ( !rField_is_long_C(currRing) )
5058        {
5059          onepoint->m[j].rtyp=STRING_CMD;
5060          onepoint->m[j].data=(void *)complexToStr((*self->roots[j])[i],oprec, currRing->cf);
5061        }
5062        else
5063        {
5064          onepoint->m[j].rtyp=NUMBER_CMD;
5065          onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
5066        }
5067        onepoint->m[j].next= NULL;
5068        onepoint->m[j].name= NULL;
5069      }
5070      listofroots->m[i].rtyp=LIST_CMD;
5071      listofroots->m[i].data=(void *)onepoint;
5072      listofroots->m[j].next= NULL;
5073      listofroots->m[j].name= NULL;
5074    }
5075
5076  }
5077  else
5078  {
5079    listofroots->Init( 0 );
5080  }
5081
5082  return listofroots;
5083}
5084
5085// from ring.cc
5086void rSetHdl(idhdl h)
5087{
5088  ring rg = NULL;
5089  if (h!=NULL)
5090  {
5091//   Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
5092    rg = IDRING(h);
5093    if (rg==NULL) return; //id <>NULL, ring==NULL
5094    omCheckAddrSize((ADDRESS)h,sizeof(idrec));
5095    if (IDID(h))  // OB: ????
5096      omCheckAddr((ADDRESS)IDID(h));
5097    rTest(rg);
5098  }
5099  else return;
5100
5101  // clean up history
5102  if (currRing!=NULL)
5103  {
5104    if(sLastPrinted.RingDependend())
5105    {
5106      sLastPrinted.CleanUp();
5107      //memset(&sLastPrinted,0,sizeof(sleftv)); // done by Cleanup,Init
5108    }
5109
5110    if (rg!=currRing)/*&&(currRing!=NULL)*/
5111    {
5112      if (rg->cf!=currRing->cf)
5113      {
5114        denominator_list dd=DENOMINATOR_LIST;
5115        if (DENOMINATOR_LIST!=NULL)
5116        {
5117          if (TEST_V_ALLWARN)
5118            Warn("deleting denom_list for ring change to %s",IDID(h));
5119          do
5120          {
5121            n_Delete(&(dd->n),currRing->cf);
5122            dd=dd->next;
5123            omFree(DENOMINATOR_LIST);
5124            DENOMINATOR_LIST=dd;
5125          } while(DENOMINATOR_LIST!=NULL);
5126        }
5127      }
5128    }
5129  }
5130
5131  // test for valid "currRing":
5132  if ((rg!=NULL) && (rg->idroot==NULL))
5133  {
5134    ring old=rg;
5135    rg=rAssure_HasComp(rg);
5136    if (old!=rg)
5137    {
5138      rKill(old);
5139      IDRING(h)=rg;
5140    }
5141  }
5142   /*------------ change the global ring -----------------------*/
5143  rChangeCurrRing(rg);
5144  currRingHdl = h;
5145}
5146
5147static leftv rOptimizeOrdAsSleftv(leftv ord)
5148{
5149  // change some bad orderings/combination into better ones
5150  leftv h=ord;
5151  while(h!=NULL)
5152  {
5153    BOOLEAN change=FALSE;
5154    intvec *iv = (intvec *)(h->data);
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 // Ws(-i) -> Wp(i)
5170    if ((*iv)[1]==ringorder_Ws)
5171    {
5172      BOOLEAN neg=TRUE;
5173      for(int i=2;i<iv->length();i++)
5174        if((*iv)[i]>=0) { neg=FALSE; break; }
5175      if (neg)
5176      {
5177        (*iv)[1]=ringorder_Wp;
5178        for(int i=2;i<iv->length();i++)
5179          (*iv)[i]= -(*iv)[i];
5180        change=TRUE;
5181      }
5182    }
5183 // wp(1) -> dp
5184    if ((*iv)[1]==ringorder_wp)
5185    {
5186      BOOLEAN all_one=TRUE;
5187      for(int i=2;i<iv->length();i++)
5188        if((*iv)[i]!=1) { all_one=FALSE; break; }
5189      if (all_one)
5190      {
5191        intvec *iv2=new intvec(3);
5192        (*iv2)[0]=1;
5193        (*iv2)[1]=ringorder_dp;
5194        (*iv2)[2]=iv->length()-2;
5195        delete iv;
5196        iv=iv2;
5197        h->data=iv2;
5198        change=TRUE;
5199      }
5200    }
5201 // Wp(1) -> Dp
5202    if ((*iv)[1]==ringorder_Wp)
5203    {
5204      BOOLEAN all_one=TRUE;
5205      for(int i=2;i<iv->length();i++)
5206        if((*iv)[i]!=1) { all_one=FALSE; break; }
5207      if (all_one)
5208      {
5209        intvec *iv2=new intvec(3);
5210        (*iv2)[0]=1;
5211        (*iv2)[1]=ringorder_Dp;
5212        (*iv2)[2]=iv->length()-2;
5213        delete iv;
5214        iv=iv2;
5215        h->data=iv2;
5216        change=TRUE;
5217      }
5218    }
5219 // dp(1)/Dp(1)/rp(1) -> lp(1)
5220    if (((*iv)[1]==ringorder_dp)
5221    || ((*iv)[1]==ringorder_Dp)
5222    || ((*iv)[1]==ringorder_rp))
5223    {
5224      if (iv->length()==3)
5225      {
5226        if ((*iv)[2]==1)
5227        {
5228          (*iv)[1]=ringorder_lp;
5229          change=TRUE;
5230        }
5231      }
5232    }
5233 // lp(i),lp(j) -> lp(i+j)
5234    if(((*iv)[1]==ringorder_lp)
5235    && (h->next!=NULL))
5236    {
5237      intvec *iv2 = (intvec *)(h->next->data);
5238      if ((*iv2)[1]==ringorder_lp)
5239      {
5240        leftv hh=h->next;
5241        h->next=hh->next;
5242        hh->next=NULL;
5243        if ((*iv2)[0]==1)
5244          (*iv)[2] += 1; // last block unspecified, at least 1
5245        else
5246          (*iv)[2] += (*iv2)[2];
5247        hh->CleanUp();
5248        omFree(hh);
5249        change=TRUE;
5250      }
5251    }
5252   // -------------------
5253    if (!change) h=h->next;
5254 }
5255 return ord;
5256}
5257
5258
5259BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
5260{
5261  int last = 0, o=0, n = 1, i=0, typ = 1, j;
5262  ord=rOptimizeOrdAsSleftv(ord);
5263  sleftv *sl = ord;
5264
5265  // determine nBlocks
5266  while (sl!=NULL)
5267  {
5268    intvec *iv = (intvec *)(sl->data);
5269    if (((*iv)[1]==ringorder_c)||((*iv)[1]==ringorder_C))
5270      i++;
5271    else if ((*iv)[1]==ringorder_L)
5272    {
5273      R->bitmask=(*iv)[2]*2+1;
5274      n--;
5275    }
5276    else if (((*iv)[1]!=ringorder_a)
5277    && ((*iv)[1]!=ringorder_a64)
5278    && ((*iv)[1]!=ringorder_am))
5279      o++;
5280    n++;
5281    sl=sl->next;
5282  }
5283  // check whether at least one real ordering
5284  if (o==0)
5285  {
5286    WerrorS("invalid combination of orderings");
5287    return TRUE;
5288  }
5289  // if no c/C ordering is given, increment n
5290  if (i==0) n++;
5291  else if (i != 1)
5292  {
5293    // throw error if more than one is given
5294    WerrorS("more than one ordering c/C specified");
5295    return TRUE;
5296  }
5297
5298  // initialize fields of R
5299  R->order=(rRingOrder_t *)omAlloc0(n*sizeof(rRingOrder_t));
5300  R->block0=(int *)omAlloc0(n*sizeof(int));
5301  R->block1=(int *)omAlloc0(n*sizeof(int));
5302  R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
5303
5304  int *weights=(int*)omAlloc0((R->N+1)*sizeof(int));
5305
5306  // init order, so that rBlocks works correctly
5307  for (j=0; j < n-1; j++)
5308    R->order[j] = ringorder_unspec;
5309  // set last _C order, if no c/C order was given
5310  if (i == 0) R->order[n-2] = ringorder_C;
5311
5312  /* init orders */
5313  sl=ord;
5314  n=-1;
5315  while (sl!=NULL)
5316  {
5317    intvec *iv;
5318    iv = (intvec *)(sl->data);
5319    if ((*iv)[1]!=ringorder_L)
5320    {
5321      n++;
5322
5323      /* the format of an ordering:
5324       *  iv[0]: factor
5325       *  iv[1]: ordering
5326       *  iv[2..end]: weights
5327       */
5328      R->order[n] = (rRingOrder_t)((*iv)[1]);
5329      typ=1;
5330      switch ((*iv)[1])
5331      {
5332          case ringorder_ws:
5333          case ringorder_Ws:
5334            typ=-1;
5335          case ringorder_wp:
5336          case ringorder_Wp:
5337            R->wvhdl[n]=(int*)omAlloc((iv->length()-1)*sizeof(int));
5338            R->block0[n] = last+1;
5339            for (i=2; i<iv->length(); i++)
5340            {
5341              R->wvhdl[n][i-2] = (*iv)[i];
5342              last++;
5343              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5344            }
5345            R->block1[n] = si_min(last,R->N);
5346            break;
5347          case ringorder_ls:
5348          case ringorder_ds:
5349          case ringorder_Ds:
5350          case ringorder_rs:
5351            typ=-1;
5352          case ringorder_lp:
5353          case ringorder_dp:
5354          case ringorder_Dp:
5355          case ringorder_rp:
5356            R->block0[n] = last+1;
5357            if (iv->length() == 3) last+=(*iv)[2];
5358            else last += (*iv)[0];
5359            R->block1[n] = si_min(last,R->N);
5360            if (rCheckIV(iv)) return TRUE;
5361            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
5362            {
5363              if (weights[i]==0) weights[i]=typ;
5364            }
5365            break;
5366
5367          case ringorder_s: // no 'rank' params!
5368          {
5369
5370            if(iv->length() > 3)
5371              return TRUE;
5372
5373            if(iv->length() == 3)
5374            {
5375              const int s = (*iv)[2];
5376              R->block0[n] = s;
5377              R->block1[n] = s;
5378            }
5379            break;
5380          }
5381          case ringorder_IS:
5382          {
5383            if(iv->length() != 3) return TRUE;
5384
5385            const int s = (*iv)[2];
5386
5387            if( 1 < s || s < -1 ) return TRUE;
5388
5389            R->block0[n] = s;
5390            R->block1[n] = s;
5391            break;
5392          }
5393          case ringorder_S:
5394          case ringorder_c:
5395          case ringorder_C:
5396          {
5397            if (rCheckIV(iv)) return TRUE;
5398            break;
5399          }
5400          case ringorder_aa:
5401          case ringorder_a:
5402          {
5403            R->block0[n] = last+1;
5404            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5405            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int));
5406            for (i=2; i<iv->length(); i++)
5407            {
5408              R->wvhdl[n][i-2]=(*iv)[i];
5409              last++;
5410              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5411            }
5412            last=R->block0[n]-1;
5413            break;
5414          }
5415          case ringorder_am:
5416          {
5417            R->block0[n] = last+1;
5418            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5419            R->wvhdl[n] = (int*)omAlloc(iv->length()*sizeof(int));
5420            if (R->block1[n]- R->block0[n]+2>=iv->length())
5421               WarnS("missing module weights");
5422            for (i=2; i<=(R->block1[n]-R->block0[n]+2); i++)
5423            {
5424              R->wvhdl[n][i-2]=(*iv)[i];
5425              last++;
5426              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5427            }
5428            R->wvhdl[n][i-2]=iv->length() -3 -(R->block1[n]- R->block0[n]);
5429            for (; i<iv->length(); i++)
5430            {
5431              R->wvhdl[n][i-1]=(*iv)[i];
5432            }
5433            last=R->block0[n]-1;
5434            break;
5435          }
5436          case ringorder_a64:
5437          {
5438            R->block0[n] = last+1;
5439            R->block1[n] = si_min(last+iv->length()-2 , R->N);
5440            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int64));
5441            int64 *w=(int64 *)R->wvhdl[n];
5442            for (i=2; i<iv->length(); i++)
5443            {
5444              w[i-2]=(*iv)[i];
5445              last++;
5446              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
5447            }
5448            last=R->block0[n]-1;
5449            break;
5450          }
5451          case ringorder_M:
5452          {
5453            int Mtyp=rTypeOfMatrixOrder(iv);
5454            if (Mtyp==0) return TRUE;
5455            if (Mtyp==-1) typ = -1;
5456
5457            R->wvhdl[n] =( int *)omAlloc((iv->length()-1)*sizeof(int));
5458            for (i=2; i<iv->length();i++)
5459              R->wvhdl[n][i-2]=(*iv)[i];
5460
5461            R->block0[n] = last+1;
5462            last += (int)sqrt((double)(iv->length()-2));
5463            R->block1[n] = si_min(last,R->N);
5464            for(i=R->block1[n];i>=R->block0[n];i--)
5465            {
5466              if (weights[i]==0) weights[i]=typ;
5467            }
5468            break;
5469          }
5470
5471          case ringorder_no:
5472            R->order[n] = ringorder_unspec;
5473            return TRUE;
5474
5475          default:
5476            Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
5477            R->order[n] = ringorder_unspec;
5478            return TRUE;
5479      }
5480    }
5481    if (last>R->N)
5482    {
5483      Werror("mismatch of number of vars (%d) and ordering (>=%d vars)",
5484             R->N,last);
5485      return TRUE;
5486    }
5487    sl=sl->next;
5488  }
5489  // find OrdSgn:
5490  R->OrdSgn = 1;
5491  for(i=1;i<=R->N;i++)
5492  { if (weights[i]<0) { R->OrdSgn=-1;break; }}
5493  omFree(weights);
5494
5495  // check for complete coverage
5496  while ( n >= 0 && (
5497          (R->order[n]==ringorder_c)
5498      ||  (R->order[n]==ringorder_C)
5499      ||  (R->order[n]==ringorder_s)
5500      ||  (R->order[n]==ringorder_S)
5501      ||  (R->order[n]==ringorder_IS)
5502                    )) n--;
5503
5504  assume( n >= 0 );
5505
5506  if (R->block1[n] != R->N)
5507  {
5508    if (((R->order[n]==ringorder_dp) ||
5509         (R->order[n]==ringorder_ds) ||
5510         (R->order[n]==ringorder_Dp) ||
5511         (R->order[n]==ringorder_Ds) ||
5512         (R->order[n]==ringorder_rp) ||
5513         (R->order[n]==ringorder_rs) ||
5514         (R->order[n]==ringorder_lp) ||
5515         (R->order[n]==ringorder_ls))
5516        &&
5517        R->block0[n] <= R->N)
5518    {
5519      R->block1[n] = R->N;
5520    }
5521    else
5522    {
5523      Werror("mismatch of number of vars (%d) and ordering (%d vars)",
5524             R->N,R->block1[n]);
5525      return TRUE;
5526    }
5527  }
5528  return FALSE;
5529}
5530
5531static BOOLEAN rSleftvList2StringArray(leftv sl, char** p)
5532{
5533
5534  while(sl!=NULL)
5535  {
5536    if ((sl->rtyp == IDHDL)||(sl->rtyp==ALIAS_CMD))
5537    {
5538      *p = omStrDup(sl->Name());
5539    }
5540    else if (sl->name!=NULL)
5541    {
5542      *p = (char*)sl->name;
5543      sl->name=NULL;
5544    }
5545    else if (sl->rtyp==POLY_CMD)
5546    {
5547      sleftv s_sl;
5548      iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
5549      if (s_sl.name != NULL)
5550      {
5551        *p = (char*)s_sl.name; s_sl.name=NULL;
5552      }
5553      else
5554        *p = NULL;
5555      sl->next = s_sl.next;
5556      s_sl.next = NULL;
5557      s_sl.CleanUp();
5558      if (*p == NULL) return TRUE;
5559    }
5560    else return TRUE;
5561    p++;
5562    sl=sl->next;
5563  }
5564  return FALSE;
5565}
5566
5567const short MAX_SHORT = 32767; // (1 << (sizeof(short)*8)) - 1;
5568
5569////////////////////
5570//
5571// rInit itself:
5572//
5573// INPUT:  pn: ch & parameter (names), rv: variable (names)
5574//         ord: ordering (all !=NULL)
5575// RETURN: currRingHdl on success
5576//         NULL        on error
5577// NOTE:   * makes new ring to current ring, on success
5578//         * considers input sleftv's as read-only
5579ring rInit(leftv pn, leftv rv, leftv ord)
5580{
5581  int float_len=0;
5582  int float_len2=0;
5583  ring R = NULL;
5584  //BOOLEAN ffChar=FALSE;
5585
5586  /* ch -------------------------------------------------------*/
5587  // get ch of ground field
5588
5589  // allocated ring
5590  R = (ring) omAlloc0Bin(sip_sring_bin);
5591
5592  coeffs cf = NULL;
5593
5594  assume( pn != NULL );
5595  const int P = pn->listLength();
5596
5597  if (pn->Typ()==CRING_CMD)
5598  {
5599    cf=(coeffs)pn->CopyD();
5600    leftv pnn=pn;
5601    if(P>1) /*parameter*/
5602    {
5603      pnn = pnn->next;
5604      const int pars = pnn->listLength();
5605      assume( pars > 0 );
5606      char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5607
5608      if (rSleftvList2StringArray(pnn, names))
5609      {
5610        WerrorS("parameter expected");
5611        goto rInitError;
5612      }
5613
5614      TransExtInfo extParam;
5615
5616      extParam.r = rDefault( cf, pars, names); // Q/Zp [ p_1, ... p_pars ]
5617      for(int i=pars-1; i>=0;i--)
5618      {
5619        omFree(names[i]);
5620      }
5621      omFree(names);
5622
5623      cf = nInitChar(n_transExt, &extParam);
5624    }
5625    assume( cf != NULL );
5626  }
5627  else if (pn->Typ()==INT_CMD)
5628  {
5629    int ch = (int)(long)pn->Data();
5630    leftv pnn=pn;
5631
5632    /* parameter? -------------------------------------------------------*/
5633    pnn = pnn->next;
5634
5635    if (pnn == NULL) // no params!?
5636    {
5637      if (ch!=0)
5638      {
5639        int ch2=IsPrime(ch);
5640        if ((ch<2)||(ch!=ch2))
5641        {
5642          Warn("%d is invalid as characteristic of the ground field. 32003 is used.", ch);
5643          ch=32003;
5644        }
5645        #ifndef TEST_ZN_AS_ZP
5646        cf = nInitChar(n_Zp, (void*)(long)ch);
5647        #else
5648        mpz_t modBase;
5649        mpz_init_set_ui(modBase, (long)ch);
5650        ZnmInfo info;
5651        info.base= modBase;
5652        info.exp= 1;
5653        cf=nInitChar(n_Zn,(void*) &info);
5654        cf->is_field=1;
5655        cf->is_domain=1;
5656        cf->has_simple_Inverse=1;
5657        #endif
5658      }
5659      else
5660        cf = nInitChar(n_Q, (void*)(long)ch);
5661    }
5662    else
5663    {
5664      const int pars = pnn->listLength();
5665
5666      assume( pars > 0 );
5667
5668      // predefined finite field: (p^k, a)
5669      if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
5670      {
5671        GFInfo param;
5672
5673        param.GFChar = ch;
5674        param.GFDegree = 1;
5675        param.GFPar_name = pnn->name;
5676
5677        cf = nInitChar(n_GF, &param);
5678      }
5679      else // (0/p, a, b, ..., z)
5680      {
5681        if ((ch!=0) && (ch!=IsPrime(ch)))
5682        {
5683          WerrorS("too many parameters");
5684          goto rInitError;
5685        }
5686
5687        char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5688
5689        if (rSleftvList2StringArray(pnn, names))
5690        {
5691          WerrorS("parameter expected");
5692          goto rInitError;
5693        }
5694
5695        TransExtInfo extParam;
5696
5697        extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5698        for(int i=pars-1; i>=0;i--)
5699        {
5700          omFree(names[i]);
5701        }
5702        omFree(names);
5703
5704        cf = nInitChar(n_transExt, &extParam);
5705      }
5706    }
5707
5708    //if (cf==NULL) ->Error: Invalid ground field specification
5709  }
5710  else if ((pn->name != NULL)
5711  && ((strcmp(pn->name,"real")==0) || (strcmp(pn->name,"complex")==0)))
5712  {
5713    leftv pnn=pn->next;
5714    BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5715    if ((pnn!=NULL) && (pnn->Typ()==INT_CMD))
5716    {
5717      float_len=(int)(long)pnn->Data();
5718      float_len2=float_len;
5719      pnn=pnn->next;
5720      if ((pnn!=NULL) && (pnn->Typ()==INT_CMD))
5721      {
5722        float_len2=(int)(long)pnn->Data();
5723        pnn=pnn->next;
5724      }
5725    }
5726
5727    if (!complex_flag)
5728      complex_flag= (pnn!=NULL) && (pnn->name!=NULL);
5729    if( !complex_flag && (float_len2 <= (short)SHORT_REAL_LENGTH))
5730       cf=nInitChar(n_R, NULL);
5731    else // longR or longC?
5732    {
5733       LongComplexInfo param;
5734
5735       param.float_len = si_min (float_len, 32767);
5736       param.float_len2 = si_min (float_len2, 32767);
5737
5738       // set the parameter name
5739       if (complex_flag)
5740       {
5741         if (param.float_len < SHORT_REAL_LENGTH)
5742         {
5743           param.float_len= SHORT_REAL_LENGTH;
5744           param.float_len2= SHORT_REAL_LENGTH;
5745         }
5746         if ((pnn == NULL) || (pnn->name == NULL))
5747           param.par_name=(const char*)"i"; //default to i
5748         else
5749           param.par_name = (const char*)pnn->name;
5750       }
5751
5752       cf = nInitChar(complex_flag ? n_long_C: n_long_R, (void*)&param);
5753    }
5754    assume( cf != NULL );
5755  }
5756#ifdef HAVE_RINGS
5757  else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5758  {
5759    // TODO: change to use coeffs_BIGINT!?
5760    mpz_t modBase;
5761    unsigned int modExponent = 1;
5762    mpz_init_set_si(modBase, 0);
5763    if (pn->next!=NULL)
5764    {
5765      leftv pnn=pn;
5766      if (pnn->next->Typ()==INT_CMD)
5767      {
5768        pnn=pnn->next;
5769        mpz_set_ui(modBase, (long) pnn->Data());
5770        if ((pnn->next!=NULL) && (pnn->next->Typ()==INT_CMD))
5771        {
5772          pnn=pnn->next;
5773          modExponent = (long) pnn->Data();
5774        }
5775        while ((pnn->next!=NULL) && (pnn->next->Typ()==INT_CMD))
5776        {
5777          pnn=pnn->next;
5778          mpz_mul_ui(modBase, modBase, (int)(long) pnn->Data());
5779        }
5780      }
5781      else if (pnn->next->Typ()==BIGINT_CMD)
5782      {
5783        number p=(number)pnn->next->CopyD();
5784        nlGMP(p,modBase,coeffs_BIGINT); // TODO? // extern void   nlGMP(number &i, mpz_t n, const coeffs r); // FIXME: n_MPZ( modBase, p, coeffs_BIGINT); ?
5785        n_Delete(&p,coeffs_BIGINT);
5786      }
5787    }
5788    else
5789      cf=nInitChar(n_Z,NULL);
5790
5791    if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_sgn1(modBase) < 0))
5792    {
5793      WerrorS("Wrong ground ring specification (module is 1)");
5794      goto rInitError;
5795    }
5796    if (modExponent < 1)
5797    {
5798      WerrorS("Wrong ground ring specification (exponent smaller than 1");
5799      goto rInitError;
5800    }
5801    // module is 0 ---> integers ringtype = 4;
5802    // we have an exponent
5803    if (modExponent > 1 && cf == NULL)
5804    {
5805      if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(unsigned long)))
5806      {
5807        /* this branch should be active for modExponent = 2..32 resp. 2..64,
5808           depending on the size of a long on the respective platform */
5809        //ringtype = 1;       // Use Z/2^ch
5810        cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5811      }
5812      else
5813      {
5814        if (mpz_sgn1(modBase)==0)
5815        {
5816          WerrorS("modulus must not be 0 or parameter not allowed");
5817          goto rInitError;
5818        }
5819        //ringtype = 3;
5820        ZnmInfo info;
5821        info.base= modBase;
5822        info.exp= modExponent;
5823        cf=nInitChar(n_Znm,(void*) &info); //exponent is missing
5824      }
5825    }
5826    // just a module m > 1
5827    else if (cf == NULL)
5828    {
5829      if (mpz_sgn1(modBase)==0)
5830      {
5831        WerrorS("modulus must not be 0 or parameter not allowed");
5832        goto rInitError;
5833      }
5834      //ringtype = 2;
5835      ZnmInfo info;
5836      info.base= modBase;
5837      info.exp= modExponent;
5838      cf=nInitChar(n_Zn,(void*) &info);
5839    }
5840    assume( cf != NULL );
5841    mpz_clear(modBase);
5842  }
5843#endif
5844  // ring NEW = OLD, (), (); where OLD is a polynomial ring...
5845  else if ((pn->Typ()==RING_CMD) && (P == 1))
5846  {
5847    TransExtInfo extParam;
5848    extParam.r = (ring)pn->Data();
5849    cf = nInitChar(n_transExt, &extParam);
5850  }
5851  //else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5852  //{
5853  //  AlgExtInfo extParam;
5854  //  extParam.r = (ring)pn->Data();
5855
5856  //  cf = nInitChar(n_algExt, &extParam);   // Q[a]/<minideal>
5857  //}
5858  else
5859  {
5860    WerrorS("Wrong or unknown ground field specification");
5861#if 0
5862// debug stuff for unknown cf descriptions:
5863    sleftv* p = pn;
5864    while (p != NULL)
5865    {
5866      Print( "pn[%p]: type: %d [%s]: %p, name: %s", (void*)p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5867      PrintLn();
5868      p = p->next;
5869    }
5870#endif
5871    goto rInitError;
5872  }
5873
5874  /*every entry in the new ring is initialized to 0*/
5875
5876  /* characteristic -----------------------------------------------*/
5877  /* input: 0 ch=0 : Q     parameter=NULL    ffChar=FALSE   float_len
5878   *         0    1 : Q(a,...)        *names         FALSE
5879   *         0   -1 : R               NULL           FALSE  0
5880   *         0   -1 : R               NULL           FALSE  prec. >6
5881   *         0   -1 : C               *names         FALSE  prec. 0..?
5882   *         p    p : Fp              NULL           FALSE
5883   *         p   -p : Fp(a)           *names         FALSE
5884   *         q    q : GF(q=p^n)       *names         TRUE
5885  */
5886  if (cf==NULL)
5887  {
5888    WerrorS("Invalid ground field specification");
5889    goto rInitError;
5890//    const int ch=32003;
5891//    cf=nInitChar(n_Zp, (void*)(long)ch);
5892  }
5893
5894  assume( R != NULL );
5895
5896  R->cf = cf;
5897
5898  /* names and number of variables-------------------------------------*/
5899  {
5900    int l=rv->listLength();
5901
5902    if (l>MAX_SHORT)
5903    {
5904      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5905       goto rInitError;
5906    }
5907    R->N = l; /*rv->listLength();*/
5908  }
5909  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5910  if (rSleftvList2StringArray(rv, R->names))
5911  {
5912    WerrorS("name of ring variable expected");
5913    goto rInitError;
5914  }
5915
5916  /* check names and parameters for conflicts ------------------------- */
5917  rRenameVars(R); // conflicting variables will be renamed
5918  /* ordering -------------------------------------------------------------*/
5919  if (rSleftvOrdering2Ordering(ord, R))
5920    goto rInitError;
5921
5922  // Complete the initialization
5923  if (rComplete(R,1))
5924    goto rInitError;
5925
5926/*#ifdef HAVE_RINGS
5927// currently, coefficients which are ring elements require a global ordering:
5928  if (rField_is_Ring(R) && (R->OrdSgn==-1))
5929  {
5930    WerrorS("global ordering required for these coefficients");
5931    goto rInitError;
5932  }
5933#endif*/
5934
5935  rTest(R);
5936
5937  // try to enter the ring into the name list
5938  // need to clean up sleftv here, before this ring can be set to
5939  // new currRing or currRing can be killed beacuse new ring has
5940  // same name
5941  pn->CleanUp();
5942  rv->CleanUp();
5943  ord->CleanUp();
5944  //if ((tmp = enterid(s, myynest, RING_CMD, &IDROOT))==NULL)
5945  //  goto rInitError;
5946
5947  //memcpy(IDRING(tmp),R,sizeof(*R));
5948  // set current ring
5949  //omFreeBin(R,  ip_sring_bin);
5950  //return tmp;
5951  return R;
5952
5953  // error case:
5954  rInitError:
5955  if  ((R != NULL)&&(R->cf!=NULL)) rDelete(R);
5956  pn->CleanUp();
5957  rv->CleanUp();
5958  ord->CleanUp();
5959  return NULL;
5960}
5961
5962ring rSubring(ring org_ring, sleftv* rv)
5963{
5964  ring R = rCopy0(org_ring);
5965  int *perm=(int *)omAlloc0((org_ring->N+1)*sizeof(int));
5966  int n = rBlocks(org_ring), i=0, j;
5967
5968  /* names and number of variables-------------------------------------*/
5969  {
5970    int l=rv->listLength();
5971    if (l>MAX_SHORT)
5972    {
5973      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5974       goto rInitError;
5975    }
5976    R->N = l; /*rv->listLength();*/
5977  }
5978  omFree(R->names);
5979  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5980  if (rSleftvList2StringArray(rv, R->names))
5981  {
5982    WerrorS("name of ring variable expected");
5983    goto rInitError;
5984  }
5985
5986  /* check names for subring in org_ring ------------------------- */
5987  {
5988    i=0;
5989
5990    for(j=0;j<R->N;j++)
5991    {
5992      for(;i<org_ring->N;i++)
5993      {
5994        if (strcmp(org_ring->names[i],R->names[j])==0)
5995        {
5996          perm[i+1]=j+1;
5997          break;
5998        }
5999      }
6000      if (i>org_ring->N)
6001      {
6002        Werror("variable %d (%s) not in basering",j+1,R->names[j]);
6003        break;
6004      }
6005    }
6006  }
6007  //Print("perm=");
6008  //for(i=1;i<org_ring->N;i++) Print("v%d -> v%d\n",i,perm[i]);
6009  /* ordering -------------------------------------------------------------*/
6010
6011  for(i=0;i<n;i++)
6012  {
6013    int min_var=-1;
6014    int max_var=-1;
6015    for(j=R->block0[i];j<=R->block1[i];j++)
6016    {
6017      if (perm[j]>0)
6018      {
6019        if (min_var==-1) min_var=perm[j];
6020        max_var=perm[j];
6021      }
6022    }
6023    if (min_var!=-1)
6024    {
6025      //Print("block %d: old %d..%d, now:%d..%d\n",
6026      //      i,R->block0[i],R->block1[i],min_var,max_var);
6027      R->block0[i]=min_var;
6028      R->block1[i]=max_var;
6029      if (R->wvhdl[i]!=NULL)
6030      {
6031        omFree(R->wvhdl[i]);
6032        R->wvhdl[i]=(int*)omAlloc0((max_var-min_var+1)*sizeof(int));
6033        for(j=org_ring->block0[i];j<=org_ring->block1[i];j++)
6034        {
6035          if (perm[j]>0)
6036          {
6037            R->wvhdl[i][perm[j]-R->block0[i]]=
6038                org_ring->wvhdl[i][j-org_ring->block0[i]];
6039            //Print("w%d=%d (orig_w%d)\n",perm[j],R->wvhdl[i][perm[j]-R->block0[i]],j);
6040          }
6041        }
6042      }
6043    }
6044    else
6045    {
6046      if(R->block0[i]>0)
6047      {
6048        //Print("skip block %d\n",i);
6049        R->order[i]=ringorder_unspec;
6050        if (R->wvhdl[i] !=NULL) omFree(R->wvhdl[i]);
6051        R->wvhdl[i]=NULL;
6052      }
6053      //else Print("keep block %d\n",i);
6054    }
6055  }
6056  i=n-1;
6057  while(i>0)
6058  {
6059    // removed unneded blocks
6060    if(R->order[i-1]==ringorder_unspec)
6061    {
6062      for(j=i;j<=n;j++)
6063      {
6064        R->order[j-1]=R->order[j];
6065        R->block0[j-1]=R->block0[j];
6066        R->block1[j-1]=R->block1[j];
6067        if (R->wvhdl[j-1] !=NULL) omFree(R->wvhdl[j-1]);
6068        R->wvhdl[j-1]=R->wvhdl[j];
6069      }
6070      R->order[n]=ringorder_unspec;
6071      n--;
6072    }
6073    i--;
6074  }
6075  n=rBlocks(org_ring)-1;
6076  while (R->order[n]==0)  n--;
6077  while (R->order[n]==ringorder_unspec)  n--;
6078  if ((R->order[n]==ringorder_c) ||  (R->order[n]==ringorder_C)) n--;
6079  if (R->block1[n] != R->N)
6080  {
6081    if (((R->order[n]==ringorder_dp) ||
6082         (R->order[n]==ringorder_ds) ||
6083         (R->order[n]==ringorder_Dp) ||
6084         (R->order[n]==ringorder_Ds) ||
6085         (R->order[n]==ringorder_rp) ||
6086         (R->order[n]==ringorder_rs) ||
6087         (R->order[n]==ringorder_lp) ||
6088         (R->order[n]==ringorder_ls))
6089        &&
6090        R->block0[n] <= R->N)
6091    {
6092      R->block1[n] = R->N;
6093    }
6094    else
6095    {
6096      Werror("mismatch of number of vars (%d) and ordering (%d vars) in block %d",
6097             R->N,R->block1[n],n);
6098      return NULL;
6099    }
6100  }
6101  omFree(perm);
6102  // find OrdSgn:
6103  R->OrdSgn = org_ring->OrdSgn; // IMPROVE!
6104  //for(i=1;i<=R->N;i++)
6105  //{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
6106  //omFree(weights);
6107  // Complete the initialization
6108  if (rComplete(R,1))
6109    goto rInitError;
6110
6111  rTest(R);
6112
6113  if (rv != NULL) rv->CleanUp();
6114
6115  return R;
6116
6117  // error case:
6118  rInitError:
6119  if  (R != NULL) rDelete(R);
6120  if (rv != NULL) rv->CleanUp();
6121  return NULL;
6122}
6123
6124void rKill(ring r)
6125{
6126  if ((r->ref<=0)&&(r->order!=NULL))
6127  {
6128#ifdef RDEBUG
6129    if (traceit &TRACE_SHOW_RINGS) Print("kill ring %lx\n",(long)r);
6130#endif
6131    int j;
6132    for (j=0;j<myynest;j++)
6133    {
6134      if (iiLocalRing[j]==r)
6135      {
6136        if (j==0) WarnS("killing the basering for level 0");
6137        iiLocalRing[j]=NULL;
6138      }
6139    }
6140// any variables depending on r ?
6141    while (r->idroot!=NULL)
6142    {
6143      r->idroot->lev=myynest; // avoid warning about kill global objects
6144      killhdl2(r->idroot,&(r->idroot),r);
6145    }
6146    if (r==currRing)
6147    {
6148      // all dependend stuff is done, clean global vars:
6149      if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
6150      if (sLastPrinted.RingDependend())
6151      {
6152        sLastPrinted.CleanUp();
6153      }
6154      //if ((myynest>0) && (iiRETURNEXPR.RingDependend()))
6155      //{
6156      //  WerrorS("return value depends on local ring variable (export missing ?)");
6157      //  iiRETURNEXPR.CleanUp();
6158      //}
6159      currRing=NULL;
6160      currRingHdl=NULL;
6161    }
6162
6163    /* nKillChar(r); will be called from inside of rDelete */
6164    rDelete(r);
6165    return;
6166  }
6167  r->ref--;
6168}
6169
6170void rKill(idhdl h)
6171{
6172  ring r = IDRING(h);
6173  int ref=0;
6174  if (r!=NULL)
6175  {
6176    // avoid, that sLastPrinted is the last reference to the base ring:
6177    // clean up before killing the last "named" refrence:
6178    if ((sLastPrinted.rtyp==RING_CMD)
6179    && (sLastPrinted.data==(void*)r))
6180    {
6181      sLastPrinted.CleanUp(r);
6182    }
6183    ref=r->ref;
6184    if ((ref<=0)&&(r==currRing))
6185    {
6186      // cleanup DENOMINATOR_LIST
6187      if (DENOMINATOR_LIST!=NULL)
6188      {
6189        denominator_list dd=DENOMINATOR_LIST;
6190        if (TEST_V_ALLWARN)
6191          Warn("deleting denom_list for ring change from %s",IDID(h));
6192        do
6193        {
6194          n_Delete(&(dd->n),currRing->cf);
6195          dd=dd->next;
6196          omFree(DENOMINATOR_LIST);
6197          DENOMINATOR_LIST=dd;
6198        } while(DENOMINATOR_LIST!=NULL);
6199      }
6200    }
6201    rKill(r);
6202  }
6203  if (h==currRingHdl)
6204  {
6205    if (ref<=0) { currRing=NULL; currRingHdl=NULL;}
6206    else
6207    {
6208      currRingHdl=rFindHdl(r,currRingHdl);
6209    }
6210  }
6211}
6212
6213idhdl rSimpleFindHdl(ring r, idhdl root, idhdl n)
6214{
6215  idhdl h=root;
6216  while (h!=NULL)
6217  {
6218    if ((IDTYP(h)==RING_CMD)
6219    && (h!=n)
6220    && (IDRING(h)==r)
6221    )
6222    {
6223      return h;
6224    }
6225    h=IDNEXT(h);
6226  }
6227  return NULL;
6228}
6229
6230extern BOOLEAN jjPROC(leftv res, leftv u, leftv v);
6231ideal kGroebner(ideal F, ideal Q)
6232{
6233  //test|=Sy_bit(OPT_PROT);
6234  idhdl save_ringhdl=currRingHdl;
6235  ideal resid;
6236  idhdl new_ring=NULL;
6237  if ((currRingHdl==NULL) || (IDRING(currRingHdl)!=currRing))
6238  {
6239    currRingHdl=enterid(" GROEBNERring",0,RING_CMD,&IDROOT,FALSE);
6240    new_ring=currRingHdl;
6241    IDRING(currRingHdl)=currRing;
6242  }
6243  sleftv v; memset(&v,0,sizeof(v)); v.rtyp=IDEAL_CMD; v.data=(char *) F;
6244  idhdl h=ggetid("groebner");
6245  sleftv u; memset(&u,0,sizeof(u)); u.rtyp=IDHDL; u.data=(char *) h;
6246            u.name=IDID(h);
6247
6248  sleftv res; memset(&res,0,sizeof(res));
6249  if(jjPROC(&res,&u,&v))
6250  {
6251    resid=kStd(F,Q,testHomog,NULL);
6252  }
6253  else
6254  {
6255    //printf("typ:%d\n",res.rtyp);
6256    resid=(ideal)(res.data);
6257  }
6258  // cleanup GROEBNERring, save_ringhdl, u,v,(res )
6259  if (new_ring!=NULL)
6260  {
6261    idhdl h=IDROOT;
6262    if (h==new_ring) IDROOT=h->next;
6263    else
6264    {
6265      while ((h!=NULL) &&(h->next!=new_ring)) h=h->next;
6266      if (h!=NULL) h->next=h->next->next;
6267    }
6268    if (h!=NULL) omFreeSize(h,sizeof(*h));
6269  }
6270  currRingHdl=save_ringhdl;
6271  u.CleanUp();
6272  v.CleanUp();
6273  return resid;
6274}
6275
6276static void jjINT_S_TO_ID(int n,int *e, leftv res)
6277{
6278  if (n==0) n=1;
6279  ideal l=idInit(n,1);
6280  int i;
6281  poly p;
6282  for(i=rVar(currRing);i>0;i--)
6283  {
6284    if (e[i]>0)
6285    {
6286      n--;
6287      p=pOne();
6288      pSetExp(p,i,1);
6289      pSetm(p);
6290      l->m[n]=p;
6291      if (n==0) break;
6292    }
6293  }
6294  res->data=(char*)l;
6295  setFlag(res,FLAG_STD);
6296  omFreeSize((ADDRESS)e,(rVar(currRing)+1)*sizeof(int));
6297}
6298BOOLEAN jjVARIABLES_P(leftv res, leftv u)
6299{
6300  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
6301  int n=pGetVariables((poly)u->Data(),e);
6302  jjINT_S_TO_ID(n,e,res);
6303  return FALSE;
6304}
6305
6306BOOLEAN jjVARIABLES_ID(leftv res, leftv u)
6307{
6308  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
6309  ideal I=(ideal)u->Data();
6310  int i;
6311  int n=0;
6312  for(i=I->nrows*I->ncols-1;i>=0;i--)
6313  {
6314    int n0=pGetVariables(I->m[i],e);
6315    if (n0>n) n=n0;
6316  }
6317  jjINT_S_TO_ID(n,e,res);
6318  return FALSE;
6319}
6320
6321void paPrint(const char *n,package p)
6322{
6323  Print(" %s (",n);
6324  switch (p->language)
6325  {
6326    case LANG_SINGULAR: PrintS("S"); break;
6327    case LANG_C:        PrintS("C"); break;
6328    case LANG_TOP:      PrintS("T"); break;
6329    case LANG_MAX:      PrintS("M"); break;
6330    case LANG_NONE:     PrintS("N"); break;
6331    default:            PrintS("U");
6332  }
6333  if(p->libname!=NULL)
6334  Print(",%s", p->libname);
6335  PrintS(")");
6336}
6337
6338BOOLEAN iiApplyINTVEC(leftv res, leftv a, int op, leftv proc)
6339{
6340  intvec *aa=(intvec*)a->Data();
6341  sleftv tmp_out;
6342  sleftv tmp_in;
6343  leftv curr=res;
6344  BOOLEAN bo=FALSE;
6345  for(int i=0;i<aa->length(); i++)
6346  {
6347    memset(&tmp_in,0,sizeof(tmp_in));
6348    tmp_in.rtyp=INT_CMD;
6349    tmp_in.data=(void*)(long)(*aa)[i];
6350    if (proc==NULL)
6351      bo=iiExprArith1(&tmp_out,&tmp_in,op);
6352    else
6353      bo=jjPROC(&tmp_out,proc,&tmp_in);
6354    if (bo)
6355    {
6356      res->CleanUp(currRing);
6357      Werror("apply fails at index %d",i+1);
6358      return TRUE;
6359    }
6360    if (i==0) { memcpy(res,&tmp_out,sizeof(tmp_out)); }
6361    else
6362    {
6363      curr->next=(leftv)omAllocBin(sleftv_bin);
6364      curr=curr->next;
6365      memcpy(curr,&tmp_out,sizeof(tmp_out));
6366    }
6367  }
6368  return FALSE;
6369}
6370BOOLEAN iiApplyBIGINTMAT(leftv, leftv, int, leftv)
6371{
6372  WerrorS("not implemented");
6373  return TRUE;
6374}
6375BOOLEAN iiApplyIDEAL(leftv, leftv, int, leftv)
6376{
6377  WerrorS("not implemented");
6378  return TRUE;
6379}
6380BOOLEAN iiApplyLIST(leftv res, leftv a, int op, leftv proc)
6381{
6382  lists aa=(lists)a->Data();
6383  sleftv tmp_out;
6384  sleftv tmp_in;
6385  leftv curr=res;
6386  BOOLEAN bo=FALSE;
6387  for(int i=0;i<=aa->nr; i++)
6388  {
6389    memset(&tmp_in,0,sizeof(tmp_in));
6390    tmp_in.Copy(&(aa->m[i]));
6391    if (proc==NULL)
6392      bo=iiExprArith1(&tmp_out,&tmp_in,op);
6393    else
6394      bo=jjPROC(&tmp_out,proc,&tmp_in);
6395    tmp_in.CleanUp();
6396    if (bo)
6397    {
6398      res->CleanUp(currRing);
6399      Werror("apply fails at index %d",i+1);
6400      return TRUE;
6401    }
6402    if (i==0) { memcpy(res,&tmp_out,sizeof(tmp_out)); }
6403    else
6404    {
6405      curr->next=(leftv)omAllocBin(sleftv_bin);
6406      curr=curr->next;
6407      memcpy(curr,&tmp_out,sizeof(tmp_out));
6408    }
6409  }
6410  return FALSE;
6411}
6412BOOLEAN iiApply(leftv res, leftv a, int op, leftv proc)
6413{
6414  memset(res,0,sizeof(sleftv));
6415  res->rtyp=a->Typ();
6416  switch (res->rtyp /*a->Typ()*/)
6417  {
6418    case INTVEC_CMD:
6419    case INTMAT_CMD:
6420        return iiApplyINTVEC(res,a,op,proc);
6421    case BIGINTMAT_CMD:
6422        return iiApplyBIGINTMAT(res,a,op,proc);
6423    case IDEAL_CMD:
6424    case MODUL_CMD:
6425    case MATRIX_CMD:
6426        return iiApplyIDEAL(res,a,op,proc);
6427    case LIST_CMD:
6428        return iiApplyLIST(res,a,op,proc);
6429  }
6430  WerrorS("first argument to `apply` must allow an index");
6431  return TRUE;
6432}
6433
6434BOOLEAN iiTestAssume(leftv a, leftv b)
6435{
6436  // assume a: level
6437  if ((a->Typ()==INT_CMD)&&((long)a->Data()>=0))
6438  {
6439    if ((TEST_V_ALLWARN) && (myynest==0)) WarnS("ASSUME at top level is of no use: see documentation");
6440    char       assume_yylinebuf[80];
6441    strncpy(assume_yylinebuf,my_yylinebuf,79);
6442    int lev=(long)a->Data();
6443    int startlev=0;
6444    idhdl h=ggetid("assumeLevel");
6445    if ((h!=NULL)&&(IDTYP(h)==INT_CMD)) startlev=(long)IDINT(h);
6446    if(lev <=startlev)
6447    {
6448      BOOLEAN bo=b->Eval();
6449      if (bo) { WerrorS("syntax error in ASSUME");return TRUE;}
6450      if (b->Typ()!=INT_CMD) { WerrorS("ASUMME(<level>,<int expr>)");return TRUE; }
6451      if (b->Data()==NULL) { Werror("ASSUME failed:%s",assume_yylinebuf);return TRUE;}
6452    }
6453  }
6454  b->CleanUp();
6455  a->CleanUp();
6456  return FALSE;
6457}
6458
6459#include "libparse.h"
6460
6461BOOLEAN iiARROW(leftv r, char* a, char *s)
6462{
6463  char *ss=(char*)omAlloc(strlen(a)+strlen(s)+30); /* max. 27 currently */
6464  // find end of s:
6465  int end_s=strlen(s);
6466  while ((end_s>0) && ((s[end_s]<=' ')||(s[end_s]==';'))) end_s--;
6467  s[end_s+1]='\0';
6468  char *name=(char *)omAlloc(strlen(a)+strlen(s)+30);
6469  sprintf(name,"%s->%s",a,s);
6470  // find start of last expression
6471  int start_s=end_s-1;
6472  while ((start_s>=0) && (s[start_s]!=';')) start_s--;
6473  if (start_s<0) // ';' not found
6474  {
6475    sprintf(ss,"parameter def %s;return(%s);\n",a,s);
6476  }
6477  else // s[start_s] is ';'
6478  {
6479    s[start_s]='\0';
6480    sprintf(ss,"parameter def %s;%s;return(%s);\n",a,s,s+start_s+1);
6481  }
6482  memset(r,0,sizeof(*r));
6483  // now produce procinfo for PROC_CMD:
6484  r->data = (void *)omAlloc0Bin(procinfo_bin);
6485  ((procinfo *)(r->data))->language=LANG_NONE;
6486  iiInitSingularProcinfo((procinfo *)r->data,"",name,0,0);
6487  ((procinfo *)r->data)->data.s.body=ss;
6488  omFree(name);
6489  r->rtyp=PROC_CMD;
6490  //r->rtyp=STRING_CMD;
6491  //r->data=ss;
6492  return FALSE;
6493}
6494
6495BOOLEAN iiAssignCR(leftv r, leftv arg)
6496{
6497  char* ring_name=omStrDup((char*)r->Name());
6498  int t=arg->Typ();
6499  if (t==RING_CMD)
6500  {
6501    sleftv tmp;
6502    memset(&tmp,0,sizeof(tmp));
6503    tmp.rtyp=IDHDL;
6504    tmp.data=(char*)rDefault(ring_name);
6505    if (tmp.data!=NULL)
6506    {
6507      BOOLEAN b=iiAssign(&tmp,arg);
6508      if (b) return TRUE;
6509      rSetHdl(ggetid(ring_name));
6510      omFree(ring_name);
6511      return FALSE;
6512    }
6513    else
6514      return TRUE;
6515  }
6516  else if (t==CRING_CMD)
6517  {
6518    sleftv tmp;
6519    sleftv n;
6520    memset(&n,0,sizeof(n));
6521    n.name=ring_name;
6522    if (iiDeclCommand(&tmp,&n,myynest,CRING_CMD,&IDROOT)) return TRUE;
6523    if (iiAssign(&tmp,arg)) return TRUE;
6524    //Print("create %s\n",r->Name());
6525    //Print("from %s(%d)\n",Tok2Cmdname(arg->Typ()),arg->Typ());
6526    return FALSE;
6527  }
6528  //Print("create %s\n",r->Name());
6529  //Print("from %s(%d)\n",Tok2Cmdname(arg->Typ()),arg->Typ());
6530  return TRUE;// not handled -> error for now
6531}
6532
6533static void iiReportTypes(int nr,int t,const short *T)
6534{
6535  char buf[250];
6536  buf[0]='\0';
6537  if (nr==0)
6538    sprintf(buf,"wrong length of parameters(%d), expected ",t);
6539  else
6540    sprintf(buf,"par. %d is of type `%s`, expected ",nr,Tok2Cmdname(t));
6541  for(int i=1;i<=T[0];i++)
6542  {
6543    strcat(buf,"`");
6544    strcat(buf,Tok2Cmdname(T[i]));
6545    strcat(buf,"`");
6546    if (i<T[0]) strcat(buf,",");
6547  }
6548  WerrorS(buf);
6549}
6550
6551BOOLEAN iiCheckTypes(leftv args, const short *type_list, int report)
6552{
6553  int l=0;
6554  if (args==NULL)
6555  {
6556    if (type_list[0]==0) return TRUE;
6557  }
6558  else l=args->listLength();
6559  if (l!=(int)type_list[0])
6560  {
6561    if (report) iiReportTypes(0,l,type_list);
6562    return FALSE;
6563  }
6564  for(int i=1;i<=l;i++,args=args->next)
6565  {
6566    short t=type_list[i];
6567    if (t!=ANY_TYPE)
6568    {
6569      if (((t==IDHDL)&&(args->rtyp!=IDHDL))
6570      || (t!=args->Typ()))
6571      {
6572        if (report) iiReportTypes(i,args->Typ(),type_list);
6573        return FALSE;
6574      }
6575    }
6576  }
6577  return TRUE;
6578}
Note: See TracBrowser for help on using the repository browser.