source: git/Singular/ipshell.cc @ ae4fd2a

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