source: git/Singular/ipshell.cc @ 77bb59

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