source: git/Singular/ipshell.cc @ e5ba1d

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