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

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