source: git/Singular/ipshell.cc @ 9e26458

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