source: git/Singular/ipshell.cc @ 6179c1

spielwiese
Last change on this file since 6179c1 was 6179c1, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix rInit: ch must be prime (Tst/Old/ch.tst)
  • Property mode set to 100644
File size: 135.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT:
7*/
8#include <kernel/mod2.h>
9#include <misc/auxiliary.h>
10
11
12#include <misc/options.h>
13#include <misc/mylimits.h>
14
15#ifdef HAVE_FACTORY
16#define SI_DONT_HAVE_GLOBAL_VARS
17#include <factory/factory.h>
18#endif
19
20#include <Singular/maps_ip.h>
21#include <Singular/tok.h>
22#include <misc/options.h>
23#include <Singular/ipid.h>
24#include <misc/intvec.h>
25#include <omalloc/omalloc.h>
26#include <kernel/febase.h>
27#include <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(strlen(IDPROC(h)->libname)>0)
172                     Print(" from %s",IDPROC(h)->libname);
173                   if(IDPROC(h)->is_static)
174                     PrintS(" (static)");
175                   break;
176    case STRING_CMD:
177                   {
178                     char *s;
179                     l=strlen(IDSTRING(h));
180                     memset(buffer,0,22);
181                     strncpy(buffer,IDSTRING(h),si_min(l,20));
182                     if ((s=strchr(buffer,'\n'))!=NULL)
183                     {
184                       *s='\0';
185                     }
186                     PrintS(" ");
187                     PrintS(buffer);
188                     if((s!=NULL) ||(l>20))
189                     {
190                       Print("..., %d char(s)",l);
191                     }
192                     break;
193                   }
194    case LIST_CMD: Print(", size: %d",IDLIST(h)->nr+1);
195                   break;
196    case QRING_CMD:
197    case RING_CMD:
198                   if ((IDRING(h)==currRing) && (currRingHdl!=h))
199                     PrintS("(*)"); /* this is an alias to currRing */
200#ifdef RDEBUG
201                   if (traceit &TRACE_SHOW_RINGS)
202                     Print(" <%lx>",(long)(IDRING(h)));
203#endif
204                   break;
205    /*default:     break;*/
206  }
207  PrintLn();
208}
209
210void type_cmd(leftv v)
211{
212  BOOLEAN oldShortOut = FALSE;
213
214  if (currRing != NULL)
215  {
216    oldShortOut = currRing->ShortOut;
217    currRing->ShortOut = 1;
218  }
219  int t=v->Typ();
220  Print("// %s %s ",v->Name(),Tok2Cmdname(t));
221  switch (t)
222  {
223    case MAP_CMD:Print(" from %s\n",((map)(v->Data()))->preimage); break;
224    case INTMAT_CMD: Print(" %d x %d\n",((intvec*)(v->Data()))->rows(),
225                                      ((intvec*)(v->Data()))->cols()); break;
226    case MATRIX_CMD:Print(" %u x %u\n" ,
227       MATROWS((matrix)(v->Data())),
228       MATCOLS((matrix)(v->Data())));break;
229    case MODUL_CMD: Print(", rk %d\n", (int)(((ideal)(v->Data()))->rank));break;
230    case LIST_CMD: Print(", size %d\n",((lists)(v->Data()))->nr+1); break;
231
232    case PROC_CMD:
233    case RING_CMD:
234    case IDEAL_CMD:
235    case QRING_CMD: PrintLn(); break;
236
237    //case INT_CMD:
238    //case STRING_CMD:
239    //case INTVEC_CMD:
240    //case POLY_CMD:
241    //case VECTOR_CMD:
242    //case PACKAGE_CMD:
243
244    default:
245      break;
246  }
247  v->Print();
248  if (currRing != NULL)
249    currRing->ShortOut = oldShortOut;
250}
251
252static void killlocals0(int v, idhdl * localhdl, const ring r)
253{
254  idhdl h = *localhdl;
255  while (h!=NULL)
256  {
257    int vv;
258    //Print("consider %s, lev: %d:",IDID(h),IDLEV(h));
259    if ((vv=IDLEV(h))>0)
260    {
261      if (vv < v)
262      {
263        if (iiNoKeepRing)
264        {
265          //PrintS(" break\n");
266          return;
267        }
268        h = IDNEXT(h);
269        //PrintLn();
270      }
271      else //if (vv >= v)
272      {
273        idhdl nexth = IDNEXT(h);
274        killhdl2(h,localhdl,r);
275        h = nexth;
276        //PrintS("kill\n");
277      }
278    }
279    else
280    {
281      h = IDNEXT(h);
282      //PrintLn();
283    }
284  }
285}
286void killlocals_list(lists l,int v)
287{
288  int i;
289  for(i=l->nr; i>=0; i--)
290  {
291    if (l->m[i].rtyp == LIST_CMD)
292      killlocals_list((lists)l->m[i].data,v);
293    else if ((l->m[i].rtyp == RING_CMD) || (l->m[i].rtyp == QRING_CMD))
294      killlocals0(v,&(((ring)(l->m[i].data))->idroot),currRing);
295  }
296}
297void killlocals_rec(idhdl *root,int v, ring r)
298{
299  idhdl h=*root;
300  while (h!=NULL)
301  {
302    if (IDLEV(h)>=v)
303    {
304//      Print("kill %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
305      idhdl n=IDNEXT(h);
306      killhdl2(h,root,r);
307      h=n;
308    }
309    else if (IDTYP(h)==PACKAGE_CMD)
310    {
311 //     Print("into pack %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
312      if (IDPACKAGE(h)!=basePack)
313        killlocals_rec(&(IDRING(h)->idroot),v,r);
314      h=IDNEXT(h);
315    }
316    else if ((IDTYP(h)==RING_CMD)
317    ||(IDTYP(h)==QRING_CMD))
318    {
319      if ((IDRING(h)!=NULL) && (IDRING(h)->idroot!=NULL))
320      // we have to test IDRING(h)!=NULL: qring Q=groebner(...): killlocals
321      {
322  //    Print("into ring %s, lev %d for lev %d\n",IDID(h),IDLEV(h),v);
323        killlocals_rec(&(IDRING(h)->idroot),v,IDRING(h));
324      }
325      h=IDNEXT(h);
326    }
327    else
328    {
329//      Print("skip %s lev %d for lev %d\n",IDID(h),IDLEV(h),v);
330      h=IDNEXT(h);
331    }
332  }
333}
334BOOLEAN killlocals_list(int v, lists L)
335{
336  if (L==NULL) return FALSE;
337  BOOLEAN changed=FALSE;
338  int n=L->nr;
339  for(;n>=0;n--)
340  {
341    leftv h=&(L->m[n]);
342    void *d=h->data;
343    if (((h->rtyp==RING_CMD) || (h->rtyp==QRING_CMD))
344    && (((ring)d)->idroot!=NULL))
345    {
346      if (d!=currRing) {changed=TRUE;rChangeCurrRing((ring)d);}
347      killlocals0(v,&(((ring)h->data)->idroot),(ring)h->data);
348    }
349    else if (h->rtyp==LIST_CMD)
350      changed|=killlocals_list(v,(lists)d);
351  }
352  return changed;
353}
354void killlocals(int v)
355{
356  BOOLEAN changed=FALSE;
357  idhdl sh=currRingHdl;
358  ring cr=currRing;
359  if (sh!=NULL) changed=((IDLEV(sh)<v) || (IDRING(sh)->ref>0));
360  //if (changed) Print("currRing=%s(%x), lev=%d,ref=%d\n",IDID(sh),IDRING(sh),IDLEV(sh),IDRING(sh)->ref);
361
362  killlocals_rec(&(basePack->idroot),v,currRing);
363
364  if (iiRETURNEXPR_len > myynest)
365  {
366    int t=iiRETURNEXPR[myynest].Typ();
367    if ((/*iiRETURNEXPR[myynest].Typ()*/ t==RING_CMD)
368    || (/*iiRETURNEXPR[myynest].Typ()*/ t==QRING_CMD))
369    {
370      leftv h=&iiRETURNEXPR[myynest];
371      if (((ring)h->data)->idroot!=NULL)
372        killlocals0(v,&(((ring)h->data)->idroot),(ring)h->data);
373    }
374    else if (/*iiRETURNEXPR[myynest].Typ()*/ t==LIST_CMD)
375    {
376      leftv h=&iiRETURNEXPR[myynest];
377      changed |=killlocals_list(v,(lists)h->data);
378    }
379  }
380  if (changed)
381  {
382    currRingHdl=rFindHdl(cr,NULL,NULL);
383    if (currRingHdl==NULL)
384      currRing=NULL;
385    else
386      rChangeCurrRing(cr);
387  }
388
389  if (myynest<=1) iiNoKeepRing=TRUE;
390  //Print("end killlocals  >= %d\n",v);
391  //listall();
392}
393
394void list_cmd(int typ, const char* what, const char *prefix,BOOLEAN iterate, BOOLEAN fullname)
395{
396  idhdl h,start;
397  BOOLEAN all = typ<0;
398  BOOLEAN really_all=FALSE;
399  BOOLEAN do_packages=FALSE;
400
401  if ( typ == -1 ) do_packages=TRUE;
402  if ( typ==0 )
403  {
404    if (strcmp(what,"all")==0)
405    {
406      really_all=TRUE;
407      h=basePack->idroot;
408    }
409    else
410    {
411      h = ggetid(what);
412      if (h!=NULL)
413      {
414        if (iterate) list1(prefix,h,TRUE,fullname);
415        if (IDTYP(h)==ALIAS_CMD) PrintS("A");
416        if ((IDTYP(h)==RING_CMD)
417            || (IDTYP(h)==QRING_CMD)
418            //|| (IDTYP(h)==PACKE_CMD)
419        )
420        {
421          h=IDRING(h)->idroot;
422        }
423        else if((IDTYP(h)==PACKAGE_CMD) || (IDTYP(h)==POINTER_CMD))
424        {
425          //Print("list_cmd:package or pointer\n");
426          all=TRUE;typ=PROC_CMD;fullname=TRUE;really_all=TRUE;
427          h=IDPACKAGE(h)->idroot;
428        }
429        else
430          return;
431      }
432      else
433      {
434        Werror("%s is undefined",what);
435        return;
436      }
437    }
438    all=TRUE;
439  }
440  else if (RingDependend(typ))
441  {
442    h = currRing->idroot;
443  }
444  else
445    h = IDROOT;
446  start=h;
447  while (h!=NULL)
448  {
449    if ((all && (IDTYP(h)!=PROC_CMD) &&(IDTYP(h)!=PACKAGE_CMD))
450    || (typ == IDTYP(h))
451    || ((IDTYP(h)==QRING_CMD) && (typ==RING_CMD)))
452    {
453      list1(prefix,h,start==currRingHdl, fullname);
454      if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
455        && (really_all || (all && (h==currRingHdl)))
456        && ((IDLEV(h)==0)||(IDLEV(h)==myynest)))
457      {
458        list_cmd(0,IDID(h),"//      ",FALSE);
459      }
460      if (IDTYP(h)==PACKAGE_CMD && really_all)
461      {
462        package save_p=currPack;
463        currPack=IDPACKAGE(h);
464        list_cmd(0,IDID(h),"//      ",FALSE);
465        currPack=save_p;
466      }
467    }
468    h = IDNEXT(h);
469  }
470}
471
472void test_cmd(int i)
473{
474  int ii;
475
476  if (i<0)
477  {
478    ii= -i;
479    if (ii < 32)
480    {
481      test &= ~Sy_bit(ii);
482    }
483    else if (ii < 64)
484    {
485      verbose &= ~Sy_bit(ii-32);
486    }
487    else
488      WerrorS("out of bounds\n");
489  }
490  else if (i<32)
491  {
492    ii=i;
493    if (Sy_bit(ii) & kOptions)
494    {
495      Warn("Gerhard, use the option command");
496      test |= Sy_bit(ii);
497    }
498    else if (Sy_bit(ii) & validOpts)
499      test |= Sy_bit(ii);
500  }
501  else if (i<64)
502  {
503    ii=i-32;
504    verbose |= Sy_bit(ii);
505  }
506  else
507    WerrorS("out of bounds\n");
508}
509
510int exprlist_length(leftv v)
511{
512  int rc = 0;
513  while (v!=NULL)
514  {
515    switch (v->Typ())
516    {
517      case INT_CMD:
518      case POLY_CMD:
519      case VECTOR_CMD:
520      case NUMBER_CMD:
521        rc++;
522        break;
523      case INTVEC_CMD:
524      case INTMAT_CMD:
525        rc += ((intvec *)(v->Data()))->length();
526        break;
527      case MATRIX_CMD:
528      case IDEAL_CMD:
529      case MODUL_CMD:
530        {
531          matrix mm = (matrix)(v->Data());
532          rc += mm->rows() * mm->cols();
533        }
534        break;
535      case LIST_CMD:
536        rc+=((lists)v->Data())->nr+1;
537        break;
538      default:
539        rc++;
540    }
541    v = v->next;
542  }
543  return rc;
544}
545
546int iiIsPrime0(unsigned p)  /* brute force !!!! */
547{
548  unsigned i,j=0 /*only to avoid compiler warnings*/;
549#ifdef HAVE_FACTORY
550  if (p<=32749) // max. small prime in factory
551  {
552    int a=0;
553    int e=cf_getNumSmallPrimes()-1;
554    i=e/2;
555    do
556    {
557      j=cf_getSmallPrime(i);
558      if (p==j) return p;
559      if (p<j) e=i-1;
560      else     a=i+1;
561      i=a+(e-a)/2;
562    } while ( a<= e);
563    if (p>j) return j;
564    else     return cf_getSmallPrime(i-1);
565  }
566#endif
567#ifdef HAVE_FACTORY
568  unsigned end_i=cf_getNumSmallPrimes()-1;
569#else
570  unsigned end_i=p/2;
571#endif
572  unsigned end_p=(unsigned)sqrt((double)p);
573restart:
574  for (i=0; i<end_i; i++)
575  {
576#ifdef HAVE_FACTORY
577    j=cf_getSmallPrime(i);
578#else
579    if (i==0) j=2;
580    else j=2*i-1;
581#endif
582    if ((p%j) == 0)
583    {
584    #ifdef HAVE_FACTORY
585      if (p<=32751) return iiIsPrime0(p-2);
586    #endif
587      p-=2;
588      goto restart;
589    }
590    if (j > end_p) return p;
591  }
592#ifdef HAVE_FACTORY
593  if (i>=end_i)
594  {
595    while(j<=end_p)
596    {
597      j+=2;
598      if ((p%j) == 0)
599      {
600        if (p<=32751) return iiIsPrime0(p-2);
601        p-=2;
602        goto restart;
603      }
604    }
605  }
606#endif
607  return p;
608}
609int IsPrime(int p)  /* brute force !!!! */
610{
611  int i,j;
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 res,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 (currRingHdl==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 w)
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->float_len,SHORT_REAL_LENGTH/2);
1688    LL->m[1].rtyp=INT_CMD;
1689    LL->m[1].data=(void *)si_max(R->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)[0]);
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 1 /* TODO */
1759  if (rField_is_numeric(r))
1760  {
1761    rDecomposeC(&(L->m[0]),r);
1762  }
1763#ifdef HAVE_RINGS
1764  else if (rField_is_Ring(r))
1765  {
1766    rDecomposeRing(&(L->m[0]),r);
1767  }
1768#endif
1769  else if (rIsExtension(r))
1770  {
1771    if (nCoeff_is_algExt(r->cf))
1772      rDecomposeCF(&(L->m[0]),r->cf->extRing,r);
1773    else
1774    {
1775      lists Lc=(lists)omAlloc0Bin(slists_bin);
1776      Lc->Init(4);
1777      // char:
1778      Lc->m[0].rtyp=INT_CMD;
1779      Lc->m[0].data=(void*)r->cf->ch;
1780      // var:
1781      lists Lv=(lists)omAlloc0Bin(slists_bin);
1782      Lv->Init(1);
1783      Lv->m[0].rtyp=STRING_CMD;
1784      Lv->m[0].data=(void *)omStrDup(rParameter(r)[0]);
1785      Lc->m[1].rtyp=LIST_CMD;
1786      Lc->m[1].data=(void*)Lv;
1787      // ord:
1788      lists Lo=(lists)omAlloc0Bin(slists_bin);
1789      Lo->Init(1);
1790      lists Loo=(lists)omAlloc0Bin(slists_bin);
1791      Loo->Init(2);
1792      Loo->m[0].rtyp=STRING_CMD;
1793      Loo->m[0].data=(void *)omStrDup(rSimpleOrdStr(ringorder_lp));
1794      intvec *iv=new intvec(1); (*iv)[0]=1;
1795      Loo->m[1].rtyp=INTVEC_CMD;
1796      Loo->m[1].data=(void *)iv;
1797      Lo->m[0].rtyp=LIST_CMD;
1798      Lo->m[0].data=(void*)Loo;
1799
1800      Lc->m[2].rtyp=LIST_CMD;
1801      Lc->m[2].data=(void*)Lo;
1802      // q-ideal:
1803      Lc->m[3].rtyp=IDEAL_CMD;
1804      Lc->m[3].data=(void *)idInit(1,1);
1805      // ----------------------
1806      L->m[0].rtyp=LIST_CMD;
1807      L->m[0].data=(void*)Lc;
1808    }
1809    if (L->m[0].rtyp==0)
1810    {
1811      //omFreeBin(slists_bin,(void *)L);
1812      return NULL;
1813    }
1814  }
1815  else
1816  #endif
1817  {
1818    L->m[0].rtyp=INT_CMD;
1819    L->m[0].data=(void *)r->cf->ch;
1820  }
1821  // ----------------------------------------
1822  // 1: list (var)
1823  lists LL=(lists)omAlloc0Bin(slists_bin);
1824  LL->Init(r->N);
1825  int i;
1826  for(i=0; i<r->N; i++)
1827  {
1828    LL->m[i].rtyp=STRING_CMD;
1829    LL->m[i].data=(void *)omStrDup(r->names[i]);
1830  }
1831  L->m[1].rtyp=LIST_CMD;
1832  L->m[1].data=(void *)LL;
1833  // ----------------------------------------
1834  // 2: list (ord)
1835  LL=(lists)omAlloc0Bin(slists_bin);
1836  i=rBlocks(r)-1;
1837  LL->Init(i);
1838  i--;
1839  lists LLL;
1840  for(; i>=0; i--)
1841  {
1842    intvec *iv;
1843    int j;
1844    LL->m[i].rtyp=LIST_CMD;
1845    LLL=(lists)omAlloc0Bin(slists_bin);
1846    LLL->Init(2);
1847    LLL->m[0].rtyp=STRING_CMD;
1848    LLL->m[0].data=(void *)omStrDup(rSimpleOrdStr(r->order[i]));
1849
1850    if(r->order[i] == ringorder_IS) //  || r->order[i] == ringorder_s || r->order[i] == ringorder_S)
1851    {
1852      assume( r->block0[i] == r->block1[i] );
1853      const int s = r->block0[i];
1854      assume( -2 < s && s < 2);
1855
1856      iv=new intvec(1);
1857      (*iv)[0] = s;
1858    } else
1859    if (r->block1[i]-r->block0[i] >=0 )
1860    {
1861      j=r->block1[i]-r->block0[i];
1862      if (r->order[i]==ringorder_M)  j=(j+1)*(j+1)-1;
1863      iv=new intvec(j+1);
1864      if ((r->wvhdl!=NULL) && (r->wvhdl[i]!=NULL))
1865      {
1866        for(;j>=0; j--) (*iv)[j]=r->wvhdl[i][j];
1867      }
1868      else switch (r->order[i])
1869      {
1870        case ringorder_dp:
1871        case ringorder_Dp:
1872        case ringorder_ds:
1873        case ringorder_Ds:
1874        case ringorder_lp:
1875          for(;j>=0; j--) (*iv)[j]=1;
1876          break;
1877        default: /* do nothing */;
1878      }
1879    }
1880    else
1881    {
1882      iv=new intvec(1);
1883    }
1884    LLL->m[1].rtyp=INTVEC_CMD;
1885    LLL->m[1].data=(void *)iv;
1886    LL->m[i].data=(void *)LLL;
1887  }
1888  L->m[2].rtyp=LIST_CMD;
1889  L->m[2].data=(void *)LL;
1890  // ----------------------------------------
1891  // 3: qideal
1892  L->m[3].rtyp=IDEAL_CMD;
1893  if (r->qideal==NULL)
1894    L->m[3].data=(void *)idInit(1,1);
1895  else
1896    L->m[3].data=(void *)idCopy(r->qideal);
1897  // ----------------------------------------
1898  #ifdef HAVE_PLURAL // NC! in rDecompose
1899  if (rIsPluralRing(r))
1900  {
1901    L->m[4].rtyp=MATRIX_CMD;
1902    L->m[4].data=(void *)mp_Copy(r->GetNC()->C, r, r);
1903    L->m[5].rtyp=MATRIX_CMD;
1904    L->m[5].data=(void *)mp_Copy(r->GetNC()->D, r, r);
1905  }
1906  #endif
1907  return L;
1908}
1909
1910void rComposeC(lists L, ring R)
1911/* field is R or C */
1912{
1913  // ----------------------------------------
1914  // 0: char/ cf - ring
1915  if ((L->m[0].rtyp!=INT_CMD) || (L->m[0].data!=(char *)0))
1916  {
1917    Werror("invald coeff. field description, expecting 0");
1918    return;
1919  }
1920  R->cf->ch=-1;
1921  // ----------------------------------------
1922  // 1:
1923  if (L->m[1].rtyp!=LIST_CMD)
1924    Werror("invald coeff. field description, expecting precision list");
1925  lists LL=(lists)L->m[1].data;
1926  int r1=(int)(long)LL->m[0].data;
1927  int r2=(int)(long)LL->m[1].data;
1928  if ((r1<=SHORT_REAL_LENGTH)
1929  && (r2=SHORT_REAL_LENGTH))
1930  {
1931    R->float_len=SHORT_REAL_LENGTH/2;
1932    R->float_len2=SHORT_REAL_LENGTH;
1933  }
1934  else
1935  {
1936    R->float_len=si_min(r1,32767);
1937    R->float_len2=si_min(r2,32767);
1938  }
1939  // ----------------------------------------
1940  // 2: list (par)
1941  if (L->nr==2)
1942  {
1943    R->cf->extRing->N=1;
1944    if (L->m[2].rtyp!=STRING_CMD)
1945    {
1946      Werror("invald coeff. field description, expecting parameter name");
1947      return;
1948    }
1949    R->cf->extRing->names=(char**)omAlloc0(rPar(R)*sizeof(char_ptr));
1950    R->cf->extRing->names[0]=omStrDup((char *)L->m[2].data);
1951  }
1952  // ----------------------------------------
1953}
1954
1955#ifdef HAVE_RINGS
1956void rComposeRing(lists L, ring R)
1957/* field is R or C */
1958{
1959  // ----------------------------------------
1960  // 0: string: integer
1961  // no further entries --> Z
1962  R->cf->modBase = (int_number) omAlloc(sizeof(mpz_t));
1963  if (L->nr == 0)
1964  {
1965    mpz_init_set_ui(R->cf->modBase,0);
1966    R->cf->modExponent = 1;
1967  }
1968  // ----------------------------------------
1969  // 1:
1970  else
1971  {
1972    if (L->m[1].rtyp!=LIST_CMD) Werror("invald data, expecting list of numbers");
1973    lists LL=(lists)L->m[1].data;
1974    mpz_init(R->cf->modBase);
1975    if ((LL->nr >= 0) && LL->m[0].rtyp == BIGINT_CMD)
1976    {
1977      number modBase = (number) LL->m[0].data;
1978      nlGMP(modBase, (number) R->cf->modBase, R->cf);
1979      LL->m[0].data = (void *)modBase;
1980    }
1981    else if ((LL->nr >= 0) && LL->m[0].rtyp == INT_CMD)
1982    {
1983      mpz_set_ui(R->cf->modBase,(unsigned long) LL->m[0].data);
1984    }
1985    else
1986    {
1987      mpz_set_ui(R->cf->modBase,0);
1988    }
1989    if (LL->nr >= 1)
1990    {
1991      R->cf->modExponent = (unsigned long) LL->m[1].data;
1992    }
1993    else
1994    {
1995      R->cf->modExponent = 1;
1996    }
1997  }
1998  // ----------------------------------------
1999  if ((mpz_cmp_ui(R->cf->modBase, 1) == 0) && (mpz_cmp_ui(R->cf->modBase, 0) < 0))
2000  {
2001    Werror("Wrong ground ring specification (module is 1)");
2002    return;
2003  }
2004  if (R->cf->modExponent < 1)
2005  {
2006    Werror("Wrong ground ring specification (exponent smaller than 1");
2007    return;
2008  }
2009  // module is 0 ---> integers
2010  if (mpz_cmp_ui(R->cf->modBase, 0) == 0)
2011  {
2012    R->cf->ch = 0;
2013    R->cf->ringtype = 4;
2014  }
2015  // we have an exponent
2016  else if (R->cf->modExponent > 1)
2017  {
2018    R->cf->ch = R->cf->modExponent;
2019    if ((mpz_cmp_ui(R->cf->modBase, 2) == 0) && (R->cf->modExponent <= 8*sizeof(NATNUMBER)))
2020    {
2021      /* this branch should be active for modExponent = 2..32 resp. 2..64,
2022           depending on the size of a long on the respective platform */
2023      R->cf->ringtype = 1;       // Use Z/2^ch
2024    }
2025    else
2026    {
2027      R->cf->ringtype = 3;
2028    }
2029  }
2030  // just a module m > 1
2031  else
2032  {
2033    R->cf->ringtype = 2;
2034    R->cf->ch = mpz_get_ui(R->cf->modBase);
2035  }
2036}
2037#endif
2038
2039static void rRenameVars(ring R)
2040{
2041  int i,j;
2042  for(i=0;i<R->N-1;i++)
2043  {
2044    for(j=i+1;j<R->N;j++)
2045    {
2046      if (strcmp(R->names[i],R->names[j])==0)
2047      {
2048        Warn("name conflict var(%d) and var(%d): `%s`, rename to `@(%d)`",i+1,j+1,R->names[i],j+1);
2049        omFree(R->names[j]);
2050        R->names[j]=(char *)omAlloc(10);
2051        sprintf(R->names[j],"@(%d)",j+1);
2052      }
2053    }
2054  }
2055  for(i=0;i<rPar(R); i++)
2056  {
2057    for(j=0;j<R->N;j++)
2058    {
2059      if (strcmp(rParameter(R)[i],R->names[j])==0)
2060      {
2061        Warn("name conflict par(%d) and var(%d): `%s`, rename to `@@(%d)`",i+1,j+1,R->names[j],i+1);
2062        omFree(rParameter(R)[i]);
2063        rParameter(R)[i]=(char *)omAlloc(10);
2064        sprintf(rParameter(R)[i],"@@(%d)",i+1);
2065      }
2066    }
2067  }
2068}
2069
2070ring rCompose(const lists  L)
2071{
2072  if ((L->nr!=3)
2073#ifdef HAVE_PLURAL
2074  &&(L->nr!=5)
2075#endif
2076  )
2077    return NULL;
2078  // 0: char/ cf - ring
2079  // 1: list (var)
2080  // 2: list (ord)
2081  // 3: qideal
2082  // possibly:
2083  // 4: C
2084  // 5: D
2085 
2086  ring R = (ring) omAlloc0Bin(sip_sring_bin); // why
2087   
2088   
2089  assume( R->cf == NULL );
2090 
2091  // ------------------------------------------------------------------
2092  // 0: char:
2093  if (L->m[0].Typ()==INT_CMD)
2094  {
2095    int ch = (int)(long)L->m[0].Data();
2096    assume( ch >= 0 );
2097
2098    if (ch == 0) // Q?
2099      R->cf = nInitChar(n_Q, NULL);
2100    else
2101    {
2102      int l = IsPrime(ch); // Zp?
2103      if( l != ch )
2104      {
2105        Warn("%d is invalid characteristic of ground field. %d is used.", ch, l);
2106        ch = l;
2107      }     
2108      R->cf = nInitChar(n_Zp, (void*)ch);       
2109    }
2110  }
2111  else if (L->m[0].Typ()==LIST_CMD) // something complicated...
2112  {
2113    lists LL=(lists)L->m[0].Data();
2114   
2115#ifdef HAVE_RINGS
2116    if (LL->m[0].Typ() == STRING_CMD) // 1st comes a string?
2117    {
2118      rComposeRing(LL, R); // Ring!?
2119    }
2120    else
2121#endif
2122    if (LL->nr < 3)
2123      rComposeC(LL,R); // R, long_R, long_C
2124    else
2125    {     
2126      if (LL->m[0].Typ()==INT_CMD)
2127      {       
2128        int ch = (int)(long)LL->m[0].Data();
2129
2130        if(( ch != 0 ) && (ch!=IsPrime(ch))) // TODO: GF-Test ch!
2131        {
2132          GFInfo param;
2133       
2134          param.GFChar = ch; 
2135          param.GFDegree = 1;
2136          param.GFPar_name = (const char*)(((lists)(LL->m[1].Data()))->m[0].Data());
2137
2138          // nfInitChar should be able to handle the case when ch is in fftables!
2139          R->cf = nInitChar(n_GF, (void*)&param);
2140        }
2141      }
2142
2143      if( R->cf == NULL )
2144      {
2145        ring extRing = rCompose((lists)L->m[0].Data());
2146       
2147        if (extRing==NULL)
2148        {
2149          WerrorS("could not create the specified coefficient field");
2150          goto rCompose_err;
2151        }
2152
2153        if( extRing->qideal != NULL ) // Algebraic extension
2154        {
2155          AlgExtInfo extParam;
2156         
2157          extParam.i = extRing->qideal;
2158          extParam.r = extRing;
2159
2160          extRing->qideal = NULL; // ???
2161
2162          R->cf = nInitChar(n_algExt, (void*)&extParam); 
2163        } else // Transcendental extension
2164        {
2165          TransExtInfo extParam;
2166          extParam.r = extRing;
2167
2168          R->cf = nInitChar(n_transExt, &extParam);
2169        }
2170       
2171      }     
2172    }
2173  }
2174  else
2175  {
2176    WerrorS("coefficient field must be described by `int` or `list`");
2177    goto rCompose_err;
2178  }
2179
2180  if( R->cf == NULL )
2181  {
2182    WerrorS("could not create coefficient field described by the input!");
2183    goto rCompose_err;
2184  }
2185
2186  // ------------------------- VARS ---------------------------
2187  if (L->m[1].Typ()==LIST_CMD)
2188  {
2189    lists v=(lists)L->m[1].Data();
2190    R->N = v->nr+1;
2191    R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
2192    int i;
2193    for(i=0;i<R->N;i++)
2194    {
2195      if (v->m[i].Typ()==STRING_CMD)
2196        R->names[i]=omStrDup((char *)v->m[i].Data());
2197      else if (v->m[i].Typ()==POLY_CMD)
2198      {
2199        poly p=(poly)v->m[i].Data();
2200        int nr=pIsPurePower(p);
2201        if (nr>0)
2202          R->names[i]=omStrDup(currRing->names[nr-1]);
2203        else
2204        {
2205          Werror("var name %d must be a string or a ring variable",i+1);
2206          goto rCompose_err;
2207        }
2208      }
2209      else
2210      {
2211        Werror("var name %d must be `string`",i+1);
2212        goto rCompose_err;
2213      }
2214    }
2215  }
2216  else
2217  {
2218    WerrorS("variable must be given as `list`");
2219    goto rCompose_err;
2220  }
2221  // ------------------------ ORDER ------------------------------
2222  if (L->m[2].Typ()==LIST_CMD)
2223  {
2224    lists v=(lists)L->m[2].Data();
2225    int n= v->nr+2;
2226    int j;
2227    // initialize fields of R
2228    R->order=(int *)omAlloc0(n*sizeof(int));
2229    R->block0=(int *)omAlloc0(n*sizeof(int));
2230    R->block1=(int *)omAlloc0(n*sizeof(int));
2231    R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
2232    // init order, so that rBlocks works correctly
2233    for (j=0; j < n-1; j++)
2234      R->order[j] = (int) ringorder_unspec;
2235    // orderings
2236    R->OrdSgn=1;
2237    for(j=0;j<n-1;j++)
2238    {
2239    // todo: a(..), M
2240      if (v->m[j].Typ()!=LIST_CMD)
2241      {
2242        WerrorS("ordering must be list of lists");
2243        goto rCompose_err;
2244      }
2245      lists vv=(lists)v->m[j].Data();
2246      if ((vv->nr!=1)
2247      || (vv->m[0].Typ()!=STRING_CMD)
2248      || ((vv->m[1].Typ()!=INTVEC_CMD) && (vv->m[1].Typ()!=INT_CMD)))
2249      {
2250        WerrorS("ordering name must be a (string,intvec)");
2251        goto rCompose_err;
2252      }
2253      R->order[j]=rOrderName(omStrDup((char*)vv->m[0].Data())); // assume STRING
2254
2255      if (j==0) R->block0[0]=1;
2256      else
2257      {
2258         int jj=j-1;
2259         while((jj>=0)
2260         && ((R->order[jj]== ringorder_a)
2261            || (R->order[jj]== ringorder_aa)
2262            || (R->order[jj]== ringorder_c)
2263            || (R->order[jj]== ringorder_C)
2264            || (R->order[jj]== ringorder_s)
2265            || (R->order[jj]== ringorder_S)
2266         ))
2267         {
2268           //Print("jj=%, skip %s\n",rSimpleOrdStr(R->order[jj]));
2269           jj--;
2270         }
2271         if (jj<0) R->block0[j]=1;
2272         else       R->block0[j]=R->block1[jj]+1;
2273      }
2274      intvec *iv;
2275      if (vv->m[1].Typ()==INT_CMD)
2276        iv=new intvec((int)(long)vv->m[1].Data(),(int)(long)vv->m[1].Data());
2277      else
2278        iv=ivCopy((intvec*)vv->m[1].Data()); //assume INTVEC
2279      int iv_len=iv->length();
2280      R->block1[j]=si_max(R->block0[j],R->block0[j]+iv_len-1);
2281      if (R->block1[j]>R->N)
2282      {
2283        R->block1[j]=R->N;
2284        iv_len=R->block1[j]-R->block0[j]+1;
2285      }
2286      //Print("block %d from %d to %d\n",j,R->block0[j], R->block1[j]);
2287      int i;
2288      switch (R->order[j])
2289      {
2290         case ringorder_ws:
2291         case ringorder_Ws:
2292            R->OrdSgn=-1;
2293         case ringorder_aa:
2294         case ringorder_a:
2295         case ringorder_wp:
2296         case ringorder_Wp:
2297           R->wvhdl[j] =( int *)omAlloc(iv_len*sizeof(int));
2298           for (i=0; i<iv_len;i++)
2299           {
2300             R->wvhdl[j][i]=(*iv)[i];
2301           }
2302           break;
2303         case ringorder_M:
2304           R->wvhdl[j] =( int *)omAlloc((iv->length())*sizeof(int));
2305           for (i=0; i<iv->length();i++) R->wvhdl[j][i]=(*iv)[i];
2306           R->block1[j]=si_max(R->block0[j],R->block0[j]+(int)sqrt((double)(iv->length()-1)));
2307           if (R->block1[j]>R->N)
2308           {
2309             WerrorS("ordering matrix too big");
2310             goto rCompose_err;
2311           }
2312           break;
2313         case ringorder_ls:
2314         case ringorder_ds:
2315         case ringorder_Ds:
2316         case ringorder_rs:
2317           R->OrdSgn=-1;
2318         case ringorder_lp:
2319         case ringorder_dp:
2320         case ringorder_Dp:
2321         case ringorder_rp:
2322           break;
2323         case ringorder_S:
2324           break;
2325         case ringorder_c:
2326         case ringorder_C:
2327           R->block1[j]=R->block0[j]=0;
2328           break;
2329
2330         case ringorder_s:
2331           break;
2332
2333         case ringorder_IS:
2334         {
2335           R->block1[j] = R->block0[j] = 0;
2336           if( iv->length() > 0 )
2337           {
2338             const int s = (*iv)[0];
2339             assume( -2 < s && s < 2 );
2340             R->block1[j] = R->block0[j] = s;
2341           }
2342           break;
2343         }
2344         case 0:
2345         case ringorder_unspec:
2346           break;
2347      }
2348      delete iv;
2349    }
2350    // sanity check
2351    j=n-2;
2352    if ((R->order[j]==ringorder_c)
2353    || (R->order[j]==ringorder_C)
2354    || (R->order[j]==ringorder_unspec)) j--;
2355    if (R->block1[j] != R->N)
2356    {
2357      if (((R->order[j]==ringorder_dp) ||
2358           (R->order[j]==ringorder_ds) ||
2359           (R->order[j]==ringorder_Dp) ||
2360           (R->order[j]==ringorder_Ds) ||
2361           (R->order[j]==ringorder_rp) ||
2362           (R->order[j]==ringorder_rs) ||
2363           (R->order[j]==ringorder_lp) ||
2364           (R->order[j]==ringorder_ls))
2365          &&
2366            R->block0[j] <= R->N)
2367      {
2368        R->block1[j] = R->N;
2369      }
2370      else
2371      {
2372        Werror("ordering incomplete: size (%d) should be %d",R->block1[j],R->N);
2373        goto rCompose_err;
2374      }
2375    }
2376  }
2377  else
2378  {
2379    WerrorS("ordering must be given as `list`");
2380    goto rCompose_err;
2381  }
2382
2383  // ------------------------ ??????? --------------------
2384
2385  rRenameVars(R);
2386  rComplete(R);
2387
2388#ifdef HAVE_RINGS
2389  // This was a BUG IN SINGULAR: There is no HABE_RINGS!!!
2390
2391// currently, coefficients which are ring elements require a global ordering:
2392  if (rField_is_Ring(R) && (R->OrdSgn==-1))
2393  {
2394    WerrorS("global ordering required for these coefficients");
2395    goto rCompose_err;
2396  }
2397#endif
2398
2399
2400  // ------------------------ Q-IDEAL ------------------------
2401
2402  if (L->m[3].Typ()==IDEAL_CMD)
2403  {
2404    ideal q=(ideal)L->m[3].Data();
2405    if (q->m[0]!=NULL)
2406    {
2407      if (R->cf->ch!=currRing->cf->ch)
2408      {
2409      #if 0
2410            WerrorS("coefficient fields must be equal if q-ideal !=0");
2411            goto rCompose_err;
2412      #else
2413        ring orig_ring=currRing;
2414        rChangeCurrRing(R);
2415        int *perm=NULL;
2416        int *par_perm=NULL;
2417        int par_perm_size=0;
2418        nMapFunc nMap;
2419        BOOLEAN bo;
2420
2421        if ((nMap=nSetMap(orig_ring->cf))==NULL)
2422        {
2423          if (rEqual(orig_ring,currRing))
2424          {
2425            nMap=n_SetMap(currRing->cf, currRing->cf);
2426          }
2427          else
2428          // Allow imap/fetch to be make an exception only for:
2429          if ( (rField_is_Q_a(orig_ring) &&  // Q(a..) -> Q(a..) || Q || Zp || Zp(a)
2430            (rField_is_Q(currRing) || rField_is_Q_a(currRing) ||
2431             (rField_is_Zp(currRing) || rField_is_Zp_a(currRing))))
2432           ||
2433           (rField_is_Zp_a(orig_ring) &&  // Zp(a..) -> Zp(a..) || Zp
2434            (rField_is_Zp(currRing, rInternalChar(orig_ring)) ||
2435             rField_is_Zp_a(currRing, rInternalChar(orig_ring)))) )
2436          {
2437            par_perm_size=rPar(orig_ring);
2438            BITSET save_test=test;
2439           
2440//            if ((orig_ring->minpoly != NULL) || (orig_ring->minideal != NULL))
2441//              naSetChar(rInternalChar(orig_ring),orig_ring);
2442//            else ntSetChar(rInternalChar(orig_ring),orig_ring);
2443           
2444            nSetChar(currRing->cf);
2445            test=save_test;
2446          }
2447          else
2448          {
2449            WerrorS("coefficient fields must be equal if q-ideal !=0");
2450            goto rCompose_err;
2451          }
2452        }
2453        perm=(int *)omAlloc0((orig_ring->N+1)*sizeof(int));
2454        if (par_perm_size!=0)
2455          par_perm=(int *)omAlloc0(par_perm_size*sizeof(int));
2456        int i;
2457        #if 0
2458        // use imap:
2459        maFindPerm(orig_ring->names,orig_ring->N,orig_ring->parameter,orig_ring->P,
2460          currRing->names,currRing->N,currRing->parameter, currRing->P,
2461          perm,par_perm, currRing->ch);
2462        #else
2463        // use fetch
2464        if ((rPar(orig_ring)>0) && (rPar(currRing)==0))
2465        {
2466          for(i=si_min(rPar(orig_ring),rVar(currRing))-1;i>=0;i--) par_perm[i]=i+1;
2467        }
2468        else if (par_perm_size!=0)
2469          for(i=si_min(rPar(orig_ring),rPar(currRing))-1;i>=0;i--) par_perm[i]=-(i+1);
2470        for(i=si_min(orig_ring->N,rVar(currRing));i>0;i--) perm[i]=i;
2471        #endif
2472        ideal dest_id=idInit(IDELEMS(q),1);
2473        for(i=IDELEMS(q)-1; i>=0; i--)
2474        {
2475          dest_id->m[i]=p_PermPoly(q->m[i],perm,orig_ring, currRing,nMap,
2476                                  par_perm,par_perm_size);
2477          //  PrintS("map:");pWrite(dest_id->m[i]);PrintLn();
2478          pTest(dest_id->m[i]);
2479        }
2480        R->qideal=dest_id;
2481        if (perm!=NULL)
2482          omFreeSize((ADDRESS)perm,(orig_ring->N+1)*sizeof(int));
2483        if (par_perm!=NULL)
2484          omFreeSize((ADDRESS)par_perm,par_perm_size*sizeof(int));
2485        rChangeCurrRing(orig_ring);
2486      #endif
2487      }
2488      else
2489        R->qideal=idrCopyR(q,currRing,R);
2490    }
2491  }
2492  else
2493  {
2494    WerrorS("q-ideal must be given as `ideal`");
2495    goto rCompose_err;
2496  }
2497
2498
2499  // ---------------------------------------------------------------
2500  #ifdef HAVE_PLURAL
2501  if (L->nr==5)
2502  {
2503    if (nc_CallPlural((matrix)L->m[4].Data(),
2504                      (matrix)L->m[5].Data(),
2505                      NULL,NULL,
2506                      R,
2507                      true, // !!!
2508                      true, false,
2509                      currRing, FALSE)) goto rCompose_err;
2510    // takes care about non-comm. quotient! i.e. calls "nc_SetupQuotient" due to last true
2511  }
2512  #endif
2513  return R;
2514
2515rCompose_err:
2516  if (R->N>0)
2517  {
2518    int i;
2519    if (R->names!=NULL)
2520    {
2521      i=R->N-1;
2522      while (i>=0) { if (R->names[i]!=NULL) omFree(R->names[i]); i--; }
2523      omFree(R->names);
2524    }
2525  }
2526  if (R->order!=NULL) omFree(R->order);
2527  if (R->block0!=NULL) omFree(R->block0);
2528  if (R->block1!=NULL) omFree(R->block1);
2529  if (R->wvhdl!=NULL) omFree(R->wvhdl);
2530  omFree(R);
2531  return NULL;
2532}
2533
2534// from matpol.cc
2535
2536/*2
2537* compute the jacobi matrix of an ideal
2538*/
2539BOOLEAN mpJacobi(leftv res,leftv a)
2540{
2541  int     i,j;
2542  matrix result;
2543  ideal id=(ideal)a->Data();
2544
2545  result =mpNew(IDELEMS(id),rVar(currRing));
2546  for (i=1; i<=IDELEMS(id); i++)
2547  {
2548    for (j=1; j<=rVar(currRing); j++)
2549    {
2550      MATELEM(result,i,j) = pDiff(id->m[i-1],j);
2551    }
2552  }
2553  res->data=(char *)result;
2554  return FALSE;
2555}
2556
2557/*2
2558* returns the Koszul-matrix of degree d of a vectorspace with dimension n
2559* uses the first n entrees of id, if id <> NULL
2560*/
2561BOOLEAN mpKoszul(leftv res,leftv c/*ip*/, leftv b/*in*/, leftv id)
2562{
2563  int n=(int)(long)b->Data();
2564  int d=(int)(long)c->Data();
2565  int     k,l,sign,row,col;
2566  matrix  result;
2567  ideal temp;
2568  BOOLEAN bo;
2569  poly    p;
2570
2571  if ((d>n) || (d<1) || (n<1))
2572  {
2573    res->data=(char *)mpNew(1,1);
2574    return FALSE;
2575  }
2576  int *choise = (int*)omAlloc(d*sizeof(int));
2577  if (id==NULL)
2578    temp=idMaxIdeal(1);
2579  else
2580    temp=(ideal)id->Data();
2581
2582  k = binom(n,d);
2583  l = k*d;
2584  l /= n-d+1;
2585  result =mpNew(l,k);
2586  col = 1;
2587  idInitChoise(d,1,n,&bo,choise);
2588  while (!bo)
2589  {
2590    sign = 1;
2591    for (l=1;l<=d;l++)
2592    {
2593      if (choise[l-1]<=IDELEMS(temp))
2594      {
2595        p = pCopy(temp->m[choise[l-1]-1]);
2596        if (sign == -1) p = pNeg(p);
2597        sign *= -1;
2598        row = idGetNumberOfChoise(l-1,d,1,n,choise);
2599        MATELEM(result,row,col) = p;
2600      }
2601    }
2602    col++;
2603    idGetNextChoise(d,n,&bo,choise);
2604  }
2605  if (id==NULL) idDelete(&temp);
2606
2607  res->data=(char *)result;
2608  return FALSE;
2609}
2610
2611// from syz1.cc
2612/*2
2613* read out the Betti numbers from resolution
2614* (interpreter interface)
2615*/
2616BOOLEAN syBetti2(leftv res, leftv u, leftv w)
2617{
2618  syStrategy syzstr=(syStrategy)u->Data();
2619
2620  BOOLEAN minim=(int)(long)w->Data();
2621  int row_shift=0;
2622  int add_row_shift=0;
2623  intvec *weights=NULL;
2624  intvec *ww=(intvec *)atGet(u,"isHomog",INTVEC_CMD);
2625  if (ww!=NULL)
2626  {
2627     weights=ivCopy(ww);
2628     add_row_shift = ww->min_in();
2629     (*weights) -= add_row_shift;
2630  }
2631
2632  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift,weights);
2633  //row_shift += add_row_shift;
2634  //Print("row_shift=%d, add_row_shift=%d\n",row_shift,add_row_shift);
2635  atSet(res,omStrDup("rowShift"),(void*)add_row_shift,INT_CMD);
2636
2637  return FALSE;
2638}
2639BOOLEAN syBetti1(leftv res, leftv u)
2640{
2641  sleftv tmp;
2642  memset(&tmp,0,sizeof(tmp));
2643  tmp.rtyp=INT_CMD;
2644  tmp.data=(void *)1;
2645  return syBetti2(res,u,&tmp);
2646}
2647
2648/*3
2649* converts a resolution into a list of modules
2650*/
2651lists syConvRes(syStrategy syzstr,BOOLEAN toDel,int add_row_shift)
2652{
2653  resolvente fullres = syzstr->fullres;
2654  resolvente minres = syzstr->minres;
2655
2656  const int length = syzstr->length;
2657
2658  if ((fullres==NULL) && (minres==NULL))
2659  {
2660    if (syzstr->hilb_coeffs==NULL)
2661    { // La Scala
2662      fullres = syReorder(syzstr->res, length, syzstr);
2663    }
2664    else
2665    { // HRES
2666      minres = syReorder(syzstr->orderedRes, length, syzstr);
2667      syKillEmptyEntres(minres, length);
2668    }
2669  }
2670
2671  resolvente tr;
2672  int typ0=IDEAL_CMD;
2673
2674  if (minres!=NULL)
2675    tr = minres;
2676  else
2677    tr = fullres;
2678
2679  resolvente trueres=NULL; intvec ** w=NULL;
2680
2681  if (length>0)
2682  {
2683    trueres = (resolvente)omAlloc0((length)*sizeof(ideal));
2684    for (int i=(length)-1;i>=0;i--)
2685    {
2686      if (tr[i]!=NULL)
2687      {
2688        trueres[i] = idCopy(tr[i]);
2689      }
2690    }
2691    if (/* id_RankFreeModule(trueres[0], currRing) */ trueres[0] > 0)
2692      typ0 = MODUL_CMD;
2693    if (syzstr->weights!=NULL)
2694    {
2695      w = (intvec**)omAlloc0(length*sizeof(intvec*));
2696      for (int i=length-1;i>=0;i--)
2697      {
2698        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
2699      }
2700    }
2701  }
2702
2703  lists li = liMakeResolv(trueres, length, syzstr->list_length,typ0,
2704                          w, add_row_shift);
2705
2706  if (w != NULL) omFreeSize(w, length*sizeof(intvec*));
2707
2708  if (toDel)
2709    syKillComputation(syzstr);
2710  else
2711  {
2712    if( fullres != NULL && syzstr->fullres == NULL )
2713      syzstr->fullres = fullres;
2714
2715    if( minres != NULL && syzstr->minres == NULL )
2716      syzstr->minres = minres;
2717  }
2718
2719  return li;
2720
2721
2722}
2723
2724/*3
2725* converts a list of modules into a resolution
2726*/
2727syStrategy syConvList(lists li,BOOLEAN toDel)
2728{
2729  int typ0;
2730  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2731
2732  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
2733  if (fr != NULL)
2734  {
2735
2736    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2737    for (int i=result->length-1;i>=0;i--)
2738    {
2739      if (fr[i]!=NULL)
2740        result->fullres[i] = idCopy(fr[i]);
2741    }
2742    result->list_length=result->length;
2743    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2744  }
2745  else
2746  {
2747    omFreeSize(result, sizeof(ssyStrategy));
2748    result = NULL;
2749  }
2750  if (toDel) li->Clean();
2751  return result;
2752}
2753
2754/*3
2755* converts a list of modules into a minimal resolution
2756*/
2757syStrategy syForceMin(lists li)
2758{
2759  int typ0;
2760  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2761
2762  resolvente fr = liFindRes(li,&(result->length),&typ0);
2763  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2764  for (int i=result->length-1;i>=0;i--)
2765  {
2766    if (fr[i]!=NULL)
2767      result->minres[i] = idCopy(fr[i]);
2768  }
2769  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2770  return result;
2771}
2772// from weight.cc
2773BOOLEAN kWeight(leftv res,leftv id)
2774{
2775  ideal F=(ideal)id->Data();
2776  intvec * iv = new intvec(rVar(currRing));
2777  polyset s;
2778  int  sl, n, i;
2779  int  *x;
2780
2781  res->data=(char *)iv;
2782  s = F->m;
2783  sl = IDELEMS(F) - 1;
2784  n = rVar(currRing);
2785  double wNsqr = (double)2.0 / (double)n;
2786  wFunctional = wFunctionalBuch;
2787  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
2788  wCall(s, sl, x, wNsqr, currRing);
2789  for (i = n; i!=0; i--)
2790    (*iv)[i-1] = x[i + n + 1];
2791  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
2792  return FALSE;
2793}
2794
2795BOOLEAN kQHWeight(leftv res,leftv v)
2796{
2797  res->data=(char *)id_QHomWeight((ideal)v->Data(), currRing);
2798  if (res->data==NULL)
2799    res->data=(char *)new intvec(rVar(currRing));
2800  return FALSE;
2801}
2802/*==============================================================*/
2803// from clapsing.cc
2804#if 0
2805BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
2806{
2807  BOOLEAN b=singclap_factorize((poly)(u->CopyD()), &v, 0);
2808  res->data=(void *)b;
2809}
2810#endif
2811
2812#ifdef HAVE_FACTORY
2813BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
2814{
2815  res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
2816                  (poly)w->CopyD(), currRing);
2817  return errorreported;
2818}
2819BOOLEAN jjCHARSERIES(leftv res, leftv u)
2820{
2821#if 1
2822  Werror("Sorry: not yet re-factored: see libpolys/polys/clapsing.cc");
2823  return FALSE;
2824#else   
2825  res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
2826  return (res->data==NULL);
2827#endif
2828}
2829#endif
2830
2831// from semic.cc
2832#ifdef HAVE_SPECTRUM
2833
2834// ----------------------------------------------------------------------------
2835//  Initialize a  spectrum  deep from a  singular  lists
2836// ----------------------------------------------------------------------------
2837
2838void copy_deep( spectrum& spec, lists l )
2839{
2840    spec.mu = (int)(long)(l->m[0].Data( ));
2841    spec.pg = (int)(long)(l->m[1].Data( ));
2842    spec.n  = (int)(long)(l->m[2].Data( ));
2843
2844    spec.copy_new( spec.n );
2845
2846    intvec  *num = (intvec*)l->m[3].Data( );
2847    intvec  *den = (intvec*)l->m[4].Data( );
2848    intvec  *mul = (intvec*)l->m[5].Data( );
2849
2850    for( int i=0; i<spec.n; i++ )
2851    {
2852        spec.s[i] = (Rational)((*num)[i])/(Rational)((*den)[i]);
2853        spec.w[i] = (*mul)[i];
2854    }
2855}
2856
2857// ----------------------------------------------------------------------------
2858//  singular lists  constructor for  spectrum
2859// ----------------------------------------------------------------------------
2860
2861spectrum /*former spectrum::spectrum ( lists l )*/
2862spectrumFromList( lists l )
2863{
2864    spectrum result;
2865    copy_deep( result, l );
2866    return result;
2867}
2868
2869// ----------------------------------------------------------------------------
2870//  generate a Singular  lists  from a spectrum
2871// ----------------------------------------------------------------------------
2872
2873/* former spectrum::thelist ( void )*/
2874lists   getList( spectrum& spec )
2875{
2876    lists   L  = (lists)omAllocBin( slists_bin);
2877
2878    L->Init( 6 );
2879
2880    intvec            *num  = new intvec( spec.n );
2881    intvec            *den  = new intvec( spec.n );
2882    intvec            *mult = new intvec( spec.n );
2883
2884    for( int i=0; i<spec.n; i++ )
2885    {
2886        (*num) [i] = spec.s[i].get_num_si( );
2887        (*den) [i] = spec.s[i].get_den_si( );
2888        (*mult)[i] = spec.w[i];
2889    }
2890
2891    L->m[0].rtyp = INT_CMD;    //  milnor number
2892    L->m[1].rtyp = INT_CMD;    //  geometrical genus
2893    L->m[2].rtyp = INT_CMD;    //  # of spectrum numbers
2894    L->m[3].rtyp = INTVEC_CMD; //  numerators
2895    L->m[4].rtyp = INTVEC_CMD; //  denomiantors
2896    L->m[5].rtyp = INTVEC_CMD; //  multiplicities
2897
2898    L->m[0].data = (void*)spec.mu;
2899    L->m[1].data = (void*)spec.pg;
2900    L->m[2].data = (void*)spec.n;
2901    L->m[3].data = (void*)num;
2902    L->m[4].data = (void*)den;
2903    L->m[5].data = (void*)mult;
2904
2905    return  L;
2906}
2907// from spectrum.cc
2908// ----------------------------------------------------------------------------
2909//  print out an error message for a spectrum list
2910// ----------------------------------------------------------------------------
2911
2912typedef enum
2913{
2914    semicOK,
2915    semicMulNegative,
2916
2917    semicListTooShort,
2918    semicListTooLong,
2919
2920    semicListFirstElementWrongType,
2921    semicListSecondElementWrongType,
2922    semicListThirdElementWrongType,
2923    semicListFourthElementWrongType,
2924    semicListFifthElementWrongType,
2925    semicListSixthElementWrongType,
2926
2927    semicListNNegative,
2928    semicListWrongNumberOfNumerators,
2929    semicListWrongNumberOfDenominators,
2930    semicListWrongNumberOfMultiplicities,
2931
2932    semicListMuNegative,
2933    semicListPgNegative,
2934    semicListNumNegative,
2935    semicListDenNegative,
2936    semicListMulNegative,
2937
2938    semicListNotSymmetric,
2939    semicListNotMonotonous,
2940
2941    semicListMilnorWrong,
2942    semicListPGWrong
2943
2944} semicState;
2945
2946void    list_error( semicState state )
2947{
2948    switch( state )
2949    {
2950        case semicListTooShort:
2951            WerrorS( "the list is too short" );
2952            break;
2953        case semicListTooLong:
2954            WerrorS( "the list is too long" );
2955            break;
2956
2957        case semicListFirstElementWrongType:
2958            WerrorS( "first element of the list should be int" );
2959            break;
2960        case semicListSecondElementWrongType:
2961            WerrorS( "second element of the list should be int" );
2962            break;
2963        case semicListThirdElementWrongType:
2964            WerrorS( "third element of the list should be int" );
2965            break;
2966        case semicListFourthElementWrongType:
2967            WerrorS( "fourth element of the list should be intvec" );
2968            break;
2969        case semicListFifthElementWrongType:
2970            WerrorS( "fifth element of the list should be intvec" );
2971            break;
2972        case semicListSixthElementWrongType:
2973            WerrorS( "sixth element of the list should be intvec" );
2974            break;
2975
2976        case semicListNNegative:
2977            WerrorS( "first element of the list should be positive" );
2978            break;
2979        case semicListWrongNumberOfNumerators:
2980            WerrorS( "wrong number of numerators" );
2981            break;
2982        case semicListWrongNumberOfDenominators:
2983            WerrorS( "wrong number of denominators" );
2984            break;
2985        case semicListWrongNumberOfMultiplicities:
2986            WerrorS( "wrong number of multiplicities" );
2987            break;
2988
2989        case semicListMuNegative:
2990            WerrorS( "the Milnor number should be positive" );
2991            break;
2992        case semicListPgNegative:
2993            WerrorS( "the geometrical genus should be nonnegative" );
2994            break;
2995        case semicListNumNegative:
2996            WerrorS( "all numerators should be positive" );
2997            break;
2998        case semicListDenNegative:
2999            WerrorS( "all denominators should be positive" );
3000            break;
3001        case semicListMulNegative:
3002            WerrorS( "all multiplicities should be positive" );
3003            break;
3004
3005        case semicListNotSymmetric:
3006            WerrorS( "it is not symmetric" );
3007            break;
3008        case semicListNotMonotonous:
3009            WerrorS( "it is not monotonous" );
3010            break;
3011
3012        case semicListMilnorWrong:
3013            WerrorS( "the Milnor number is wrong" );
3014            break;
3015        case semicListPGWrong:
3016            WerrorS( "the geometrical genus is wrong" );
3017            break;
3018
3019        default:
3020            WerrorS( "unspecific error" );
3021            break;
3022    }
3023}
3024// ----------------------------------------------------------------------------
3025//  this is the main spectrum computation function
3026// ----------------------------------------------------------------------------
3027
3028enum    spectrumState
3029{
3030    spectrumOK,
3031    spectrumZero,
3032    spectrumBadPoly,
3033    spectrumNoSingularity,
3034    spectrumNotIsolated,
3035    spectrumDegenerate,
3036    spectrumWrongRing,
3037    spectrumNoHC,
3038    spectrumUnspecErr
3039};
3040
3041// from splist.cc
3042// ----------------------------------------------------------------------------
3043//  Compute the spectrum of a  spectrumPolyList
3044// ----------------------------------------------------------------------------
3045
3046/* former spectrumPolyList::spectrum ( lists*, int) */
3047spectrumState   spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3048{
3049  spectrumPolyNode  **node = &speclist.root;
3050  spectrumPolyNode  *search;
3051
3052  poly              f,tmp;
3053  int               found,cmp;
3054
3055  Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3056                 ( fast==2 ? 2 : 1 ) );
3057
3058  Rational weight_prev( 0,1 );
3059
3060  int     mu = 0;          // the milnor number
3061  int     pg = 0;          // the geometrical genus
3062  int     n  = 0;          // number of different spectral numbers
3063  int     z  = 0;          // number of spectral number equal to smax
3064
3065  int     k = 0;
3066
3067  while( (*node)!=(spectrumPolyNode*)NULL &&
3068         ( fast==0 || (*node)->weight<=smax ) )
3069  {
3070        // ---------------------------------------
3071        //  determine the first normal form which
3072        //  contains the monomial  node->mon
3073        // ---------------------------------------
3074
3075    found  = FALSE;
3076    search = *node;
3077
3078    while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3079    {
3080      if( search->nf!=(poly)NULL )
3081      {
3082        f = search->nf;
3083
3084        do
3085        {
3086                    // --------------------------------
3087                    //  look for  (*node)->mon  in   f
3088                    // --------------------------------
3089
3090          cmp = pCmp( (*node)->mon,f );
3091
3092          if( cmp<0 )
3093          {
3094            f = pNext( f );
3095          }
3096          else if( cmp==0 )
3097          {
3098                        // -----------------------------
3099                        //  we have found a normal form
3100                        // -----------------------------
3101
3102            found = TRUE;
3103
3104                        //  normalize coefficient
3105
3106            number inv = nInvers( pGetCoeff( f ) );
3107            pMult_nn( search->nf,inv );
3108            nDelete( &inv );
3109
3110                        //  exchange  normal forms
3111
3112            tmp         = (*node)->nf;
3113            (*node)->nf = search->nf;
3114            search->nf  = tmp;
3115          }
3116        }
3117        while( cmp<0 && f!=(poly)NULL );
3118      }
3119      search = search->next;
3120    }
3121
3122    if( found==FALSE )
3123    {
3124            // ------------------------------------------------
3125            //  the weight of  node->mon  is a spectrum number
3126            // ------------------------------------------------
3127
3128      mu++;
3129
3130      if( (*node)->weight<=(Rational)1 )              pg++;
3131      if( (*node)->weight==smax )           z++;
3132      if( (*node)->weight>weight_prev )     n++;
3133
3134      weight_prev = (*node)->weight;
3135      node = &((*node)->next);
3136    }
3137    else
3138    {
3139            // -----------------------------------------------
3140            //  determine all other normal form which contain
3141            //  the monomial  node->mon
3142            //  replace for  node->mon  its normal form
3143            // -----------------------------------------------
3144
3145      while( search!=(spectrumPolyNode*)NULL )
3146      {
3147        if( search->nf!=(poly)NULL )
3148        {
3149          f = search->nf;
3150
3151          do
3152          {
3153                        // --------------------------------
3154                        //  look for  (*node)->mon  in   f
3155                        // --------------------------------
3156
3157            cmp = pCmp( (*node)->mon,f );
3158
3159            if( cmp<0 )
3160            {
3161              f = pNext( f );
3162            }
3163            else if( cmp==0 )
3164            {
3165              search->nf = pSub( search->nf,
3166                                 ppMult_nn( (*node)->nf,pGetCoeff( f ) ) );
3167              pNorm( search->nf );
3168            }
3169          }
3170          while( cmp<0 && f!=(poly)NULL );
3171        }
3172        search = search->next;
3173      }
3174      speclist.delete_node( node );
3175    }
3176
3177  }
3178
3179    // --------------------------------------------------------
3180    //  fast computation exploits the symmetry of the spectrum
3181    // --------------------------------------------------------
3182
3183  if( fast==2 )
3184  {
3185    mu = 2*mu - z;
3186    n  = ( z > 0 ? 2*n - 1 : 2*n );
3187  }
3188
3189    // --------------------------------------------------------
3190    //  compute the spectrum numbers with their multiplicities
3191    // --------------------------------------------------------
3192
3193  intvec            *nom  = new intvec( n );
3194  intvec            *den  = new intvec( n );
3195  intvec            *mult = new intvec( n );
3196
3197  int count         = 0;
3198  int multiplicity  = 1;
3199
3200  for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3201              ( fast==0 || search->weight<=smax );
3202       search=search->next )
3203  {
3204    if( search->next==(spectrumPolyNode*)NULL ||
3205        search->weight<search->next->weight )
3206    {
3207      (*nom) [count] = search->weight.get_num_si( );
3208      (*den) [count] = search->weight.get_den_si( );
3209      (*mult)[count] = multiplicity;
3210
3211      multiplicity=1;
3212      count++;
3213    }
3214    else
3215    {
3216      multiplicity++;
3217    }
3218  }
3219
3220    // --------------------------------------------------------
3221    //  fast computation exploits the symmetry of the spectrum
3222    // --------------------------------------------------------
3223
3224  if( fast==2 )
3225  {
3226    int n1,n2;
3227    for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3228    {
3229      (*nom) [n2] = rVar(currRing)*(*den)[n1]-(*nom)[n1];
3230      (*den) [n2] = (*den)[n1];
3231      (*mult)[n2] = (*mult)[n1];
3232    }
3233  }
3234
3235    // -----------------------------------
3236    //  test if the spectrum is symmetric
3237    // -----------------------------------
3238
3239  if( fast==0 || fast==1 )
3240  {
3241    int symmetric=TRUE;
3242
3243    for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3244    {
3245      if( (*mult)[n1]!=(*mult)[n2] ||
3246          (*den) [n1]!= (*den)[n2] ||
3247          (*nom)[n1]+(*nom)[n2]!=rVar(currRing)*(*den) [n1] )
3248      {
3249        symmetric = FALSE;
3250      }
3251    }
3252
3253    if( symmetric==FALSE )
3254    {
3255            // ---------------------------------------------
3256            //  the spectrum is not symmetric => degenerate
3257            //  principal part
3258            // ---------------------------------------------
3259
3260      *L = (lists)omAllocBin( slists_bin);
3261      (*L)->Init( 1 );
3262      (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3263      (*L)->m[0].data = (void*)mu;
3264
3265      return spectrumDegenerate;
3266    }
3267  }
3268
3269  *L = (lists)omAllocBin( slists_bin);
3270
3271  (*L)->Init( 6 );
3272
3273  (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3274  (*L)->m[1].rtyp = INT_CMD;    //  geometrical genus
3275  (*L)->m[2].rtyp = INT_CMD;    //  number of spectrum values
3276  (*L)->m[3].rtyp = INTVEC_CMD; //  nominators
3277  (*L)->m[4].rtyp = INTVEC_CMD; //  denomiantors
3278  (*L)->m[5].rtyp = INTVEC_CMD; //  multiplicities
3279
3280  (*L)->m[0].data = (void*)mu;
3281  (*L)->m[1].data = (void*)pg;
3282  (*L)->m[2].data = (void*)n;
3283  (*L)->m[3].data = (void*)nom;
3284  (*L)->m[4].data = (void*)den;
3285  (*L)->m[5].data = (void*)mult;
3286
3287  return  spectrumOK;
3288}
3289
3290spectrumState   spectrumCompute( poly h,lists *L,int fast )
3291{
3292  int i,j;
3293
3294  #ifdef SPECTRUM_DEBUG
3295  #ifdef SPECTRUM_PRINT
3296  #ifdef SPECTRUM_IOSTREAM
3297    cout << "spectrumCompute\n";
3298    if( fast==0 ) cout << "    no optimization" << endl;
3299    if( fast==1 ) cout << "    weight optimization" << endl;
3300    if( fast==2 ) cout << "    symmetry optimization" << endl;
3301  #else
3302    fprintf( stdout,"spectrumCompute\n" );
3303    if( fast==0 ) fprintf( stdout,"    no optimization\n" );
3304    if( fast==1 ) fprintf( stdout,"    weight optimization\n" );
3305    if( fast==2 ) fprintf( stdout,"    symmetry optimization\n" );
3306  #endif
3307  #endif
3308  #endif
3309
3310  // ----------------------
3311  //  check if  h  is zero
3312  // ----------------------
3313
3314  if( h==(poly)NULL )
3315  {
3316    return  spectrumZero;
3317  }
3318
3319  // ----------------------------------
3320  //  check if  h  has a constant term
3321  // ----------------------------------
3322
3323  if( hasConstTerm( h, currRing ) )
3324  {
3325    return  spectrumBadPoly;
3326  }
3327
3328  // --------------------------------
3329  //  check if  h  has a linear term
3330  // --------------------------------
3331
3332  if( hasLinearTerm( h, currRing ) )
3333  {
3334    *L = (lists)omAllocBin( slists_bin);
3335    (*L)->Init( 1 );
3336    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3337    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3338
3339    return  spectrumNoSingularity;
3340  }
3341
3342  // ----------------------------------
3343  //  compute the jacobi ideal of  (h)
3344  // ----------------------------------
3345
3346  ideal J = NULL;
3347  J = idInit( rVar(currRing),1 );
3348
3349  #ifdef SPECTRUM_DEBUG
3350  #ifdef SPECTRUM_PRINT
3351  #ifdef SPECTRUM_IOSTREAM
3352    cout << "\n   computing the Jacobi ideal...\n";
3353  #else
3354    fprintf( stdout,"\n   computing the Jacobi ideal...\n" );
3355  #endif
3356  #endif
3357  #endif
3358
3359  for( i=0; i<rVar(currRing); i++ )
3360  {
3361    J->m[i] = pDiff( h,i+1); //j );
3362
3363    #ifdef SPECTRUM_DEBUG
3364    #ifdef SPECTRUM_PRINT
3365    #ifdef SPECTRUM_IOSTREAM
3366      cout << "        ";
3367    #else
3368      fprintf( stdout,"        " );
3369    #endif
3370      pWrite( J->m[i] );
3371    #endif
3372    #endif
3373  }
3374
3375  // --------------------------------------------
3376  //  compute a standard basis  stdJ  of  jac(h)
3377  // --------------------------------------------
3378
3379  #ifdef SPECTRUM_DEBUG
3380  #ifdef SPECTRUM_PRINT
3381  #ifdef SPECTRUM_IOSTREAM
3382    cout << endl;
3383    cout << "    computing a standard basis..." << endl;
3384  #else
3385    fprintf( stdout,"\n" );
3386    fprintf( stdout,"    computing a standard basis...\n" );
3387  #endif
3388  #endif
3389  #endif
3390
3391  ideal stdJ = kStd(J,currQuotient,isNotHomog,NULL);
3392  idSkipZeroes( stdJ );
3393
3394  #ifdef SPECTRUM_DEBUG
3395  #ifdef SPECTRUM_PRINT
3396    for( i=0; i<IDELEMS(stdJ); i++ )
3397    {
3398      #ifdef SPECTRUM_IOSTREAM
3399        cout << "        ";
3400      #else
3401        fprintf( stdout,"        " );
3402      #endif
3403
3404      pWrite( stdJ->m[i] );
3405    }
3406  #endif
3407  #endif
3408
3409  idDelete( &J );
3410
3411  // ------------------------------------------
3412  //  check if the  h  has a singularity
3413  // ------------------------------------------
3414
3415  if( hasOne( stdJ, currRing ) )
3416  {
3417    // -------------------------------
3418    //  h is smooth in the origin
3419    //  return only the Milnor number
3420    // -------------------------------
3421
3422    *L = (lists)omAllocBin( slists_bin);
3423    (*L)->Init( 1 );
3424    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3425    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3426
3427    return  spectrumNoSingularity;
3428  }
3429
3430  // ------------------------------------------
3431  //  check if the singularity  h  is isolated
3432  // ------------------------------------------
3433
3434  for( i=rVar(currRing); i>0; i-- )
3435  {
3436    if( hasAxis( stdJ,i, currRing )==FALSE )
3437    {
3438      return  spectrumNotIsolated;
3439    }
3440  }
3441
3442  // ------------------------------------------
3443  //  compute the highest corner  hc  of  stdJ
3444  // ------------------------------------------
3445
3446  #ifdef SPECTRUM_DEBUG
3447  #ifdef SPECTRUM_PRINT
3448  #ifdef SPECTRUM_IOSTREAM
3449    cout << "\n    computing the highest corner...\n";
3450  #else
3451    fprintf( stdout,"\n    computing the highest corner...\n" );
3452  #endif
3453  #endif
3454  #endif
3455
3456  poly hc = (poly)NULL;
3457
3458  scComputeHC( stdJ,currQuotient, 0,hc );
3459
3460  if( hc!=(poly)NULL )
3461  {
3462    pGetCoeff(hc) = nInit(1);
3463
3464    for( i=rVar(currRing); i>0; i-- )
3465    {
3466      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3467    }
3468    pSetm( hc );
3469  }
3470  else
3471  {
3472    return  spectrumNoHC;
3473  }
3474
3475  #ifdef SPECTRUM_DEBUG
3476  #ifdef SPECTRUM_PRINT
3477  #ifdef SPECTRUM_IOSTREAM
3478    cout << "       ";
3479  #else
3480    fprintf( stdout,"       " );
3481  #endif
3482    pWrite( hc );
3483  #endif
3484  #endif
3485
3486  // ----------------------------------------
3487  //  compute the Newton polygon  nph  of  h
3488  // ----------------------------------------
3489
3490  #ifdef SPECTRUM_DEBUG
3491  #ifdef SPECTRUM_PRINT
3492  #ifdef SPECTRUM_IOSTREAM
3493    cout << "\n    computing the newton polygon...\n";
3494  #else
3495    fprintf( stdout,"\n    computing the newton polygon...\n" );
3496  #endif
3497  #endif
3498  #endif
3499
3500  newtonPolygon nph( h, currRing );
3501
3502  #ifdef SPECTRUM_DEBUG
3503  #ifdef SPECTRUM_PRINT
3504    cout << nph;
3505  #endif
3506  #endif
3507
3508  // -----------------------------------------------
3509  //  compute the weight corner  wc  of  (stdj,nph)
3510  // -----------------------------------------------
3511
3512  #ifdef SPECTRUM_DEBUG
3513  #ifdef SPECTRUM_PRINT
3514  #ifdef SPECTRUM_IOSTREAM
3515    cout << "\n    computing the weight corner...\n";
3516  #else
3517    fprintf( stdout,"\n    computing the weight corner...\n" );
3518  #endif
3519  #endif
3520  #endif
3521
3522  poly    wc = ( fast==0 ? pCopy( hc ) :
3523               ( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
3524              /* fast==2 */computeWC( nph,
3525                      ((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
3526
3527  #ifdef SPECTRUM_DEBUG
3528  #ifdef SPECTRUM_PRINT
3529  #ifdef SPECTRUM_IOSTREAM
3530    cout << "        ";
3531  #else
3532    fprintf( stdout,"        " );
3533  #endif
3534    pWrite( wc );
3535  #endif
3536  #endif
3537
3538  // -------------
3539  //  compute  NF
3540  // -------------
3541
3542  #ifdef SPECTRUM_DEBUG
3543  #ifdef SPECTRUM_PRINT
3544  #ifdef SPECTRUM_IOSTREAM
3545    cout << "\n    computing NF...\n" << endl;
3546  #else
3547    fprintf( stdout,"\n    computing NF...\n" );
3548  #endif
3549  #endif
3550  #endif
3551
3552  spectrumPolyList NF( &nph );
3553
3554  computeNF( stdJ,hc,wc,&NF, currRing );
3555
3556  #ifdef SPECTRUM_DEBUG
3557  #ifdef SPECTRUM_PRINT
3558    cout << NF;
3559  #ifdef SPECTRUM_IOSTREAM
3560    cout << endl;
3561  #else
3562    fprintf( stdout,"\n" );
3563  #endif
3564  #endif
3565  #endif
3566
3567  // ----------------------------
3568  //  compute the spectrum of  h
3569  // ----------------------------
3570//  spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
3571
3572  return spectrumStateFromList(NF, L, fast );
3573}
3574
3575// ----------------------------------------------------------------------------
3576//  this procedure is called from the interpreter
3577// ----------------------------------------------------------------------------
3578//  first  = polynomial
3579//  result = list of spectrum numbers
3580// ----------------------------------------------------------------------------
3581
3582void spectrumPrintError(spectrumState state)
3583{
3584  switch( state )
3585  {
3586    case spectrumZero:
3587      WerrorS( "polynomial is zero" );
3588      break;
3589    case spectrumBadPoly:
3590      WerrorS( "polynomial has constant term" );
3591      break;
3592    case spectrumNoSingularity:
3593      WerrorS( "not a singularity" );
3594      break;
3595    case spectrumNotIsolated:
3596      WerrorS( "the singularity is not isolated" );
3597      break;
3598    case spectrumNoHC:
3599      WerrorS( "highest corner cannot be computed" );
3600      break;
3601    case spectrumDegenerate:
3602      WerrorS( "principal part is degenerate" );
3603      break;
3604    case spectrumOK:
3605      break;
3606
3607    default:
3608      WerrorS( "unknown error occurred" );
3609      break;
3610  }
3611}
3612
3613BOOLEAN spectrumProc( leftv result,leftv first )
3614{
3615  spectrumState state = spectrumOK;
3616
3617  // -------------------
3618  //  check consistency
3619  // -------------------
3620
3621  //  check for a local ring
3622
3623  if( !ringIsLocal(currRing ) )
3624  {
3625    WerrorS( "only works for local orderings" );
3626    state = spectrumWrongRing;
3627  }
3628
3629  //  no quotient rings are allowed
3630
3631  else if( currRing->qideal != NULL )
3632  {
3633    WerrorS( "does not work in quotient rings" );
3634    state = spectrumWrongRing;
3635  }
3636  else
3637  {
3638    lists   L    = (lists)NULL;
3639    int     flag = 1; // weight corner optimization is safe
3640
3641    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3642
3643    if( state==spectrumOK )
3644    {
3645      result->rtyp = LIST_CMD;
3646      result->data = (char*)L;
3647    }
3648    else
3649    {
3650      spectrumPrintError(state);
3651    }
3652  }
3653
3654  return  (state!=spectrumOK);
3655}
3656
3657// ----------------------------------------------------------------------------
3658//  this procedure is called from the interpreter
3659// ----------------------------------------------------------------------------
3660//  first  = polynomial
3661//  result = list of spectrum numbers
3662// ----------------------------------------------------------------------------
3663
3664BOOLEAN spectrumfProc( leftv result,leftv first )
3665{
3666  spectrumState state = spectrumOK;
3667
3668  // -------------------
3669  //  check consistency
3670  // -------------------
3671
3672  //  check for a local polynomial ring
3673
3674  if( currRing->OrdSgn != -1 )
3675  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
3676  // or should we use:
3677  //if( !ringIsLocal( ) )
3678  {
3679    WerrorS( "only works for local orderings" );
3680    state = spectrumWrongRing;
3681  }
3682  else if( currRing->qideal != NULL )
3683  {
3684    WerrorS( "does not work in quotient rings" );
3685    state = spectrumWrongRing;
3686  }
3687  else
3688  {
3689    lists   L    = (lists)NULL;
3690    int     flag = 2; // symmetric optimization
3691
3692    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3693
3694    if( state==spectrumOK )
3695    {
3696      result->rtyp = LIST_CMD;
3697      result->data = (char*)L;
3698    }
3699    else
3700    {
3701      spectrumPrintError(state);
3702    }
3703  }
3704
3705  return  (state!=spectrumOK);
3706}
3707
3708// ----------------------------------------------------------------------------
3709//  check if a list is a spectrum
3710//  check for:
3711//      list has 6 elements
3712//      1st element is int (mu=Milnor number)
3713//      2nd element is int (pg=geometrical genus)
3714//      3rd element is int (n =number of different spectrum numbers)
3715//      4th element is intvec (num=numerators)
3716//      5th element is intvec (den=denomiantors)
3717//      6th element is intvec (mul=multiplicities)
3718//      exactly n numerators
3719//      exactly n denominators
3720//      exactly n multiplicities
3721//      mu>0
3722//      pg>=0
3723//      n>0
3724//      num>0
3725//      den>0
3726//      mul>0
3727//      symmetriy with respect to numberofvariables/2
3728//      monotony
3729//      mu = sum of all multiplicities
3730//      pg = sum of all multiplicities where num/den<=1
3731// ----------------------------------------------------------------------------
3732
3733semicState  list_is_spectrum( lists l )
3734{
3735    // -------------------
3736    //  check list length
3737    // -------------------
3738
3739    if( l->nr < 5 )
3740    {
3741        return  semicListTooShort;
3742    }
3743    else if( l->nr > 5 )
3744    {
3745        return  semicListTooLong;
3746    }
3747
3748    // -------------
3749    //  check types
3750    // -------------
3751
3752    if( l->m[0].rtyp != INT_CMD )
3753    {
3754        return  semicListFirstElementWrongType;
3755    }
3756    else if( l->m[1].rtyp != INT_CMD )
3757    {
3758        return  semicListSecondElementWrongType;
3759    }
3760    else if( l->m[2].rtyp != INT_CMD )
3761    {
3762        return  semicListThirdElementWrongType;
3763    }
3764    else if( l->m[3].rtyp != INTVEC_CMD )
3765    {
3766        return  semicListFourthElementWrongType;
3767    }
3768    else if( l->m[4].rtyp != INTVEC_CMD )
3769    {
3770        return  semicListFifthElementWrongType;
3771    }
3772    else if( l->m[5].rtyp != INTVEC_CMD )
3773    {
3774        return  semicListSixthElementWrongType;
3775    }
3776
3777    // -------------------------
3778    //  check number of entries
3779    // -------------------------
3780
3781    int     mu = (int)(long)(l->m[0].Data( ));
3782    int     pg = (int)(long)(l->m[1].Data( ));
3783    int     n  = (int)(long)(l->m[2].Data( ));
3784
3785    if( n <= 0 )
3786    {
3787        return  semicListNNegative;
3788    }
3789
3790    intvec  *num = (intvec*)l->m[3].Data( );
3791    intvec  *den = (intvec*)l->m[4].Data( );
3792    intvec  *mul = (intvec*)l->m[5].Data( );
3793
3794    if( n != num->length( ) )
3795    {
3796        return  semicListWrongNumberOfNumerators;
3797    }
3798    else if( n != den->length( ) )
3799    {
3800        return  semicListWrongNumberOfDenominators;
3801    }
3802    else if( n != mul->length( ) )
3803    {
3804        return  semicListWrongNumberOfMultiplicities;
3805    }
3806
3807    // --------
3808    //  values
3809    // --------
3810
3811    if( mu <= 0 )
3812    {
3813        return  semicListMuNegative;
3814    }
3815    if( pg < 0 )
3816    {
3817        return  semicListPgNegative;
3818    }
3819
3820    int i;
3821
3822    for( i=0; i<n; i++ )
3823    {
3824        if( (*num)[i] <= 0 )
3825        {
3826            return  semicListNumNegative;
3827        }
3828        if( (*den)[i] <= 0 )
3829        {
3830            return  semicListDenNegative;
3831        }
3832        if( (*mul)[i] <= 0 )
3833        {
3834            return  semicListMulNegative;
3835        }
3836    }
3837
3838    // ----------------
3839    //  check symmetry
3840    // ----------------
3841
3842    int     j;
3843
3844    for( i=0, j=n-1; i<=j; i++,j-- )
3845    {
3846        if( (*num)[i] != rVar(currRing)*((*den)[i]) - (*num)[j] ||
3847            (*den)[i] != (*den)[j] ||
3848            (*mul)[i] != (*mul)[j] )
3849        {
3850            return  semicListNotSymmetric;
3851        }
3852    }
3853
3854    // ----------------
3855    //  check monotony
3856    // ----------------
3857
3858    for( i=0, j=1; i<n/2; i++,j++ )
3859    {
3860        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
3861        {
3862            return  semicListNotMonotonous;
3863        }
3864    }
3865
3866    // ---------------------
3867    //  check Milnor number
3868    // ---------------------
3869
3870    for( mu=0, i=0; i<n; i++ )
3871    {
3872        mu += (*mul)[i];
3873    }
3874
3875    if( mu != (int)(long)(l->m[0].Data( )) )
3876    {
3877        return  semicListMilnorWrong;
3878    }
3879
3880    // -------------------------
3881    //  check geometrical genus
3882    // -------------------------
3883
3884    for( pg=0, i=0; i<n; i++ )
3885    {
3886        if( (*num)[i]<=(*den)[i] )
3887        {
3888            pg += (*mul)[i];
3889        }
3890    }
3891
3892    if( pg != (int)(long)(l->m[1].Data( )) )
3893    {
3894        return  semicListPGWrong;
3895    }
3896
3897    return  semicOK;
3898}
3899
3900// ----------------------------------------------------------------------------
3901//  this procedure is called from the interpreter
3902// ----------------------------------------------------------------------------
3903//  first  = list of spectrum numbers
3904//  second = list of spectrum numbers
3905//  result = sum of the two lists
3906// ----------------------------------------------------------------------------
3907
3908BOOLEAN spaddProc( leftv result,leftv first,leftv second )
3909{
3910    semicState  state;
3911
3912    // -----------------
3913    //  check arguments
3914    // -----------------
3915
3916    lists l1 = (lists)first->Data( );
3917    lists l2 = (lists)second->Data( );
3918
3919    if( (state=list_is_spectrum( l1 )) != semicOK )
3920    {
3921        WerrorS( "first argument is not a spectrum:" );
3922        list_error( state );
3923    }
3924    else if( (state=list_is_spectrum( l2 )) != semicOK )
3925    {
3926        WerrorS( "second argument is not a spectrum:" );
3927        list_error( state );
3928    }
3929    else
3930    {
3931        spectrum s1= spectrumFromList ( l1 );
3932        spectrum s2= spectrumFromList ( l2 );
3933        spectrum sum( s1+s2 );
3934
3935        result->rtyp = LIST_CMD;
3936        result->data = (char*)(getList(sum));
3937    }
3938
3939    return  (state!=semicOK);
3940}
3941
3942// ----------------------------------------------------------------------------
3943//  this procedure is called from the interpreter
3944// ----------------------------------------------------------------------------
3945//  first  = list of spectrum numbers
3946//  second = integer
3947//  result = the multiple of the first list by the second factor
3948// ----------------------------------------------------------------------------
3949
3950BOOLEAN spmulProc( leftv result,leftv first,leftv second )
3951{
3952    semicState  state;
3953
3954    // -----------------
3955    //  check arguments
3956    // -----------------
3957
3958    lists   l = (lists)first->Data( );
3959    int     k = (int)(long)second->Data( );
3960
3961    if( (state=list_is_spectrum( l ))!=semicOK )
3962    {
3963        WerrorS( "first argument is not a spectrum" );
3964        list_error( state );
3965    }
3966    else if( k < 0 )
3967    {
3968        WerrorS( "second argument should be positive" );
3969        state = semicMulNegative;
3970    }
3971    else
3972    {
3973        spectrum s= spectrumFromList( l );
3974        spectrum product( k*s );
3975
3976        result->rtyp = LIST_CMD;
3977        result->data = (char*)getList(product);
3978    }
3979
3980    return  (state!=semicOK);
3981}
3982
3983// ----------------------------------------------------------------------------
3984//  this procedure is called from the interpreter
3985// ----------------------------------------------------------------------------
3986//  first  = list of spectrum numbers
3987//  second = list of spectrum numbers
3988//  result = semicontinuity index
3989// ----------------------------------------------------------------------------
3990
3991BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
3992{
3993  semicState  state;
3994  BOOLEAN qh=(((int)(long)w->Data())==1);
3995
3996  // -----------------
3997  //  check arguments
3998  // -----------------
3999
4000  lists l1 = (lists)u->Data( );
4001  lists l2 = (lists)v->Data( );
4002
4003  if( (state=list_is_spectrum( l1 ))!=semicOK )
4004  {
4005    WerrorS( "first argument is not a spectrum" );
4006    list_error( state );
4007  }
4008  else if( (state=list_is_spectrum( l2 ))!=semicOK )
4009  {
4010    WerrorS( "second argument is not a spectrum" );
4011    list_error( state );
4012  }
4013  else
4014  {
4015    spectrum s1= spectrumFromList( l1 );
4016    spectrum s2= spectrumFromList( l2 );
4017
4018    res->rtyp = INT_CMD;
4019    if (qh)
4020      res->data = (void*)(s1.mult_spectrumh( s2 ));
4021    else
4022      res->data = (void*)(s1.mult_spectrum( s2 ));
4023  }
4024
4025  // -----------------
4026  //  check status
4027  // -----------------
4028
4029  return  (state!=semicOK);
4030}
4031BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
4032{
4033  sleftv tmp;
4034  memset(&tmp,0,sizeof(tmp));
4035  tmp.rtyp=INT_CMD;
4036  /* tmp.data = (void *)0;  -- done by memset */
4037
4038  return  semicProc3(res,u,v,&tmp);
4039}
4040
4041#endif
4042
4043//from mpr_inout.cc
4044extern void nPrint(number n);
4045
4046BOOLEAN loNewtonP( leftv res, leftv arg1 )
4047{
4048  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4049  return FALSE;
4050}
4051
4052BOOLEAN loSimplex( leftv res, leftv args )
4053{
4054  if ( !(rField_is_long_R(currRing)) )
4055  {
4056    WerrorS("Ground field not implemented!");
4057    return TRUE;
4058  }
4059
4060  simplex * LP;
4061  matrix m;
4062
4063  leftv v= args;
4064  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4065    return TRUE;
4066  else
4067    m= (matrix)(v->CopyD());
4068
4069  LP = new simplex(MATROWS(m),MATCOLS(m));
4070  LP->mapFromMatrix(m);
4071
4072  v= v->next;
4073  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
4074    return TRUE;
4075  else
4076    LP->m= (int)(long)(v->Data());
4077
4078  v= v->next;
4079  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
4080    return TRUE;
4081  else
4082    LP->n= (int)(long)(v->Data());
4083
4084  v= v->next;
4085  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
4086    return TRUE;
4087  else
4088    LP->m1= (int)(long)(v->Data());
4089
4090  v= v->next;
4091  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
4092    return TRUE;
4093  else
4094    LP->m2= (int)(long)(v->Data());
4095
4096  v= v->next;
4097  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
4098    return TRUE;
4099  else
4100    LP->m3= (int)(long)(v->Data());
4101
4102#ifdef mprDEBUG_PROT
4103  Print("m (constraints) %d\n",LP->m);
4104  Print("n (columns) %d\n",LP->n);
4105  Print("m1 (<=) %d\n",LP->m1);
4106  Print("m2 (>=) %d\n",LP->m2);
4107  Print("m3 (==) %d\n",LP->m3);
4108#endif
4109
4110  LP->compute();
4111
4112  lists lres= (lists)omAlloc( sizeof(slists) );
4113  lres->Init( 6 );
4114
4115  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4116  lres->m[0].data=(void*)LP->mapToMatrix(m);
4117
4118  lres->m[1].rtyp= INT_CMD;   // found a solution?
4119  lres->m[1].data=(void*)LP->icase;
4120
4121  lres->m[2].rtyp= INTVEC_CMD;
4122  lres->m[2].data=(void*)LP->posvToIV();
4123
4124  lres->m[3].rtyp= INTVEC_CMD;
4125  lres->m[3].data=(void*)LP->zrovToIV();
4126
4127  lres->m[4].rtyp= INT_CMD;
4128  lres->m[4].data=(void*)LP->m;
4129
4130  lres->m[5].rtyp= INT_CMD;
4131  lres->m[5].data=(void*)LP->n;
4132
4133  res->data= (void*)lres;
4134
4135  return FALSE;
4136}
4137
4138BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4139{
4140  ideal gls = (ideal)(arg1->Data());
4141  int imtype= (int)(long)arg2->Data();
4142
4143  uResultant::resMatType mtype= determineMType( imtype );
4144
4145  // check input ideal ( = polynomial system )
4146  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4147  {
4148    return TRUE;
4149  }
4150
4151  uResultant *resMat= new uResultant( gls, mtype, false );
4152  if (resMat!=NULL)
4153  {
4154    res->rtyp = MODUL_CMD;
4155    res->data= (void*)resMat->accessResMat()->getMatrix();
4156    if (!errorreported) delete resMat;
4157  }
4158  return errorreported;
4159}
4160
4161BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4162{
4163
4164  poly gls;
4165  gls= (poly)(arg1->Data());
4166  int howclean= (int)(long)arg3->Data();
4167
4168  if ( !(rField_is_R(currRing) ||
4169         rField_is_Q(currRing) ||
4170         rField_is_long_R(currRing) ||
4171         rField_is_long_C(currRing)) )
4172  {
4173    WerrorS("Ground field not implemented!");
4174    return TRUE;
4175  }
4176
4177  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4178                          rField_is_long_C(currRing)) )
4179  {
4180    unsigned long int ii = (unsigned long int)arg2->Data();
4181    setGMPFloatDigits( ii, ii );
4182  }
4183
4184  if ( gls == NULL || pIsConstant( gls ) )
4185  {
4186    WerrorS("Input polynomial is constant!");
4187    return TRUE;
4188  }
4189
4190  int ldummy;
4191  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4192  //  int deg= pDeg( gls );
4193  int len= pLength( gls );
4194  int i,vpos=0;
4195  poly piter;
4196  lists elist;
4197  lists rlist;
4198
4199  elist= (lists)omAlloc( sizeof(slists) );
4200  elist->Init( 0 );
4201
4202  if ( rVar(currRing) > 1 )
4203  {
4204    piter= gls;
4205    for ( i= 1; i <= rVar(currRing); i++ )
4206      if ( pGetExp( piter, i ) )
4207      {
4208        vpos= i;
4209        break;
4210      }
4211    while ( piter )
4212    {
4213      for ( i= 1; i <= rVar(currRing); i++ )
4214        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4215        {
4216          WerrorS("The input polynomial must be univariate!");
4217          return TRUE;
4218        }
4219      pIter( piter );
4220    }
4221  }
4222
4223  rootContainer * roots= new rootContainer();
4224  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4225  piter= gls;
4226  for ( i= deg; i >= 0; i-- )
4227  {
4228    //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4229    if ( piter && pTotaldegree(piter) == i )
4230    {
4231      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4232      //nPrint( pcoeffs[i] );PrintS("  ");
4233      pIter( piter );
4234    }
4235    else
4236    {
4237      pcoeffs[i]= nInit(0);
4238    }
4239  }
4240
4241#ifdef mprDEBUG_PROT
4242  for (i=deg; i >= 0; i--)
4243  {
4244    nPrint( pcoeffs[i] );PrintS("  ");
4245  }
4246  PrintLn();
4247#endif
4248
4249  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4250  roots->solver( howclean );
4251
4252  int elem= roots->getAnzRoots();
4253  char *out;
4254  char *dummy;
4255  int j;
4256
4257  rlist= (lists)omAlloc( sizeof(slists) );
4258  rlist->Init( elem );
4259
4260  if (rField_is_long_C(currRing))
4261  {
4262    for ( j= 0; j < elem; j++ )
4263    {
4264      rlist->m[j].rtyp=NUMBER_CMD;
4265      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4266      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4267    }
4268  }
4269  else
4270  {
4271    for ( j= 0; j < elem; j++ )
4272    {
4273      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4274      rlist->m[j].rtyp=STRING_CMD;
4275      rlist->m[j].data=(void *)dummy;
4276    }
4277  }
4278
4279  elist->Clean();
4280  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4281
4282  // this is (via fillContainer) the same data as in root
4283  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4284  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4285
4286  delete roots;
4287
4288  res->rtyp= LIST_CMD;
4289  res->data= (void*)rlist;
4290
4291  return FALSE;
4292}
4293
4294BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4295{
4296  int i;
4297  ideal p,w;
4298  p= (ideal)arg1->Data();
4299  w= (ideal)arg2->Data();
4300
4301  // w[0] = f(p^0)
4302  // w[1] = f(p^1)
4303  // ...
4304  // p can be a vector of numbers (multivariate polynom)
4305  //   or one number (univariate polynom)
4306  // tdg = deg(f)
4307
4308  int n= IDELEMS( p );
4309  int m= IDELEMS( w );
4310  int tdg= (int)(long)arg3->Data();
4311
4312  res->data= (void*)NULL;
4313
4314  // check the input
4315  if ( tdg < 1 )
4316  {
4317    WerrorS("Last input parameter must be > 0!");
4318    return TRUE;
4319  }
4320  if ( n != rVar(currRing) )
4321  {
4322    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4323    return TRUE;
4324  }
4325  if ( m != (int)pow((double)tdg+1,(double)n) )
4326  {
4327    Werror("Size of second input ideal must be equal to %d!",
4328      (int)pow((double)tdg+1,(double)n));
4329    return TRUE;
4330  }
4331  if ( !(rField_is_Q(currRing) /* ||
4332         rField_is_R() || rField_is_long_R() ||
4333         rField_is_long_C()*/ ) )
4334         {
4335    WerrorS("Ground field not implemented!");
4336    return TRUE;
4337  }
4338
4339  number tmp;
4340  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4341  for ( i= 0; i < n; i++ )
4342  {
4343    pevpoint[i]=nInit(0);
4344    if (  (p->m)[i] )
4345    {
4346      tmp = pGetCoeff( (p->m)[i] );
4347      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4348      {
4349        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4350        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4351        return TRUE;
4352      }
4353    } else tmp= NULL;
4354    if ( !nIsZero(tmp) )
4355    {
4356      if ( !pIsConstant((p->m)[i]))
4357      {
4358        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4359        WerrorS("Elements of first input ideal must be numbers!");
4360        return TRUE;
4361      }
4362      pevpoint[i]= nCopy( tmp );
4363    }
4364  }
4365
4366  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4367  for ( i= 0; i < m; i++ )
4368  {
4369    wresults[i]= nInit(0);
4370    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4371    {
4372      if ( !pIsConstant((w->m)[i]))
4373      {
4374        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4375        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4376        WerrorS("Elements of second input ideal must be numbers!");
4377        return TRUE;
4378      }
4379      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4380    }
4381  }
4382
4383  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4384  number *ncpoly= vm.interpolateDense( wresults );
4385  // do not free ncpoly[]!!
4386  poly rpoly= vm.numvec2poly( ncpoly );
4387
4388  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4389  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4390
4391  res->data= (void*)rpoly;
4392  return FALSE;
4393}
4394
4395BOOLEAN nuUResSolve( leftv res, leftv args )
4396{
4397  leftv v= args;
4398
4399  ideal gls;
4400  int imtype;
4401  int howclean;
4402
4403  // get ideal
4404  if ( v->Typ() != IDEAL_CMD )
4405    return TRUE;
4406  else gls= (ideal)(v->Data());
4407  v= v->next;
4408
4409  // get resultant matrix type to use (0,1)
4410  if ( v->Typ() != INT_CMD )
4411    return TRUE;
4412  else imtype= (int)(long)v->Data();
4413  v= v->next;
4414
4415  if (imtype==0)
4416  {
4417    ideal test_id=idInit(1,1);
4418    int j;
4419    for(j=IDELEMS(gls)-1;j>=0;j--)
4420    {
4421      if (gls->m[j]!=NULL)
4422      {
4423        test_id->m[0]=gls->m[j];
4424        intvec *dummy_w=id_QHomWeight(test_id, currRing);
4425        if (dummy_w!=NULL)
4426        {
4427          WerrorS("Newton polytope not of expected dimension");
4428          delete dummy_w;
4429          return TRUE;
4430        }
4431      }
4432    }
4433  }
4434
4435  // get and set precision in digits ( > 0 )
4436  if ( v->Typ() != INT_CMD )
4437    return TRUE;
4438  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4439                          rField_is_long_C(currRing)) )
4440  {
4441    unsigned long int ii=(unsigned long int)v->Data();
4442    setGMPFloatDigits( ii, ii );
4443  }
4444  v= v->next;
4445
4446  // get interpolation steps (0,1,2)
4447  if ( v->Typ() != INT_CMD )
4448    return TRUE;
4449  else howclean= (int)(long)v->Data();
4450
4451  uResultant::resMatType mtype= determineMType( imtype );
4452  int i,c,count;
4453  lists listofroots= NULL;
4454  lists emptylist;
4455  number smv= NULL;
4456  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4457
4458  //emptylist= (lists)omAlloc( sizeof(slists) );
4459  //emptylist->Init( 0 );
4460
4461  //res->rtyp = LIST_CMD;
4462  //res->data= (void *)emptylist;
4463
4464  // check input ideal ( = polynomial system )
4465  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4466  {
4467    return TRUE;
4468  }
4469
4470  uResultant * ures;
4471  rootContainer ** iproots;
4472  rootContainer ** muiproots;
4473  rootArranger * arranger;
4474
4475  // main task 1: setup of resultant matrix
4476  ures= new uResultant( gls, mtype );
4477  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4478  {
4479    WerrorS("Error occurred during matrix setup!");
4480    return TRUE;
4481  }
4482
4483  // if dense resultant, check if minor nonsingular
4484  if ( mtype == uResultant::denseResMat )
4485  {
4486    smv= ures->accessResMat()->getSubDet();
4487#ifdef mprDEBUG_PROT
4488    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4489#endif
4490    if ( nIsZero(smv) )
4491    {
4492      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4493      return TRUE;
4494    }
4495  }
4496
4497  // main task 2: Interpolate specialized resultant polynomials
4498  if ( interpolate_det )
4499    iproots= ures->interpolateDenseSP( false, smv );
4500  else
4501    iproots= ures->specializeInU( false, smv );
4502
4503  // main task 3: Interpolate specialized resultant polynomials
4504  if ( interpolate_det )
4505    muiproots= ures->interpolateDenseSP( true, smv );
4506  else
4507    muiproots= ures->specializeInU( true, smv );
4508
4509#ifdef mprDEBUG_PROT
4510  c= iproots[0]->getAnzElems();
4511  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4512  c= muiproots[0]->getAnzElems();
4513  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4514#endif
4515
4516  // main task 4: Compute roots of specialized polys and match them up
4517  arranger= new rootArranger( iproots, muiproots, howclean );
4518  arranger->solve_all();
4519
4520  // get list of roots
4521  if ( arranger->success() )
4522  {
4523    arranger->arrange();
4524    listofroots= listOfRoots(arranger, gmp_output_digits );
4525  }
4526  else
4527  {
4528    WerrorS("Solver was unable to find any roots!");
4529    return TRUE;
4530  }
4531
4532  // free everything
4533  count= iproots[0]->getAnzElems();
4534  for (i=0; i < count; i++) delete iproots[i];
4535  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4536  count= muiproots[0]->getAnzElems();
4537  for (i=0; i < count; i++) delete muiproots[i];
4538  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4539
4540  delete ures;
4541  delete arranger;
4542  nDelete( &smv );
4543
4544  res->data= (void *)listofroots;
4545
4546  //emptylist->Clean();
4547  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4548
4549  return FALSE;
4550}
4551
4552// from mpr_numeric.cc
4553lists listOfRoots( rootArranger* self, const unsigned int oprec )
4554{
4555  int i,j,tr;
4556  int count= self->roots[0]->getAnzRoots(); // number of roots
4557  int elem= self->roots[0]->getAnzElems();  // number of koordinates per root
4558
4559  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
4560
4561  if ( self->found_roots )
4562  {
4563    listofroots->Init( count );
4564
4565    for (i=0; i < count; i++)
4566    {
4567      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
4568      onepoint->Init(elem);
4569      for ( j= 0; j < elem; j++ )
4570      {
4571        if ( !rField_is_long_C(currRing) )
4572        {
4573          onepoint->m[j].rtyp=STRING_CMD;
4574          onepoint->m[j].data=(void *)complexToStr((*self->roots[j])[i],oprec, currRing->cf);
4575        }
4576        else
4577        {
4578          onepoint->m[j].rtyp=NUMBER_CMD;
4579          onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
4580        }
4581        onepoint->m[j].next= NULL;
4582        onepoint->m[j].name= NULL;
4583      }
4584      listofroots->m[i].rtyp=LIST_CMD;
4585      listofroots->m[i].data=(void *)onepoint;
4586      listofroots->m[j].next= NULL;
4587      listofroots->m[j].name= NULL;
4588    }
4589
4590  }
4591  else
4592  {
4593    listofroots->Init( 0 );
4594  }
4595
4596  return listofroots;
4597}
4598
4599// from ring.cc
4600void rSetHdl(idhdl h)
4601{
4602  int i;
4603  ring rg = NULL;
4604  if (h!=NULL)
4605  {
4606//   Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
4607    rg = IDRING(h);
4608    if (rg==NULL) return; //id <>NULL, ring==NULL
4609    omCheckAddrSize((ADDRESS)h,sizeof(idrec));
4610    if (IDID(h))  // OB: ????
4611      omCheckAddr((ADDRESS)IDID(h));
4612    rTest(rg);
4613  }
4614
4615  // clean up history
4616  if (sLastPrinted.RingDependend())
4617  {
4618    sLastPrinted.CleanUp();
4619    memset(&sLastPrinted,0,sizeof(sleftv));
4620  }
4621
4622  // test for valid "currRing":
4623  if ((rg!=NULL) && (rg->idroot==NULL))
4624  {
4625    ring old=rg;
4626    rg=rAssure_HasComp(rg);
4627    if (old!=rg)
4628    {
4629      rKill(old);
4630      IDRING(h)=rg;
4631    }
4632  }
4633   /*------------ change the global ring -----------------------*/
4634  rChangeCurrRing(rg);
4635  currRingHdl = h;
4636}
4637
4638BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
4639{
4640  int last = 0, o=0, n = 1, i=0, typ = 1, j;
4641  sleftv *sl = ord;
4642
4643  // determine nBlocks
4644  while (sl!=NULL)
4645  {
4646    intvec *iv = (intvec *)(sl->data);
4647    if (((*iv)[1]==ringorder_c)||((*iv)[1]==ringorder_C))
4648      i++;
4649    else if ((*iv)[1]==ringorder_L)
4650    {
4651      R->bitmask=(*iv)[2];
4652      n--;
4653    }
4654    else if (((*iv)[1]!=ringorder_a)
4655    && ((*iv)[1]!=ringorder_a64))
4656      o++;
4657    n++;
4658    sl=sl->next;
4659  }
4660  // check whether at least one real ordering
4661  if (o==0)
4662  {
4663    WerrorS("invalid combination of orderings");
4664    return TRUE;
4665  }
4666  // if no c/C ordering is given, increment n
4667  if (i==0) n++;
4668  else if (i != 1)
4669  {
4670    // throw error if more than one is given
4671    WerrorS("more than one ordering c/C specified");
4672    return TRUE;
4673  }
4674
4675  // initialize fields of R
4676  R->order=(int *)omAlloc0(n*sizeof(int));
4677  R->block0=(int *)omAlloc0(n*sizeof(int));
4678  R->block1=(int *)omAlloc0(n*sizeof(int));
4679  R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
4680
4681  int *weights=(int*)omAlloc0((R->N+1)*sizeof(int));
4682
4683  // init order, so that rBlocks works correctly
4684  for (j=0; j < n-1; j++)
4685    R->order[j] = (int) ringorder_unspec;
4686  // set last _C order, if no c/C order was given
4687  if (i == 0) R->order[n-2] = ringorder_C;
4688
4689  /* init orders */
4690  sl=ord;
4691  n=-1;
4692  while (sl!=NULL)
4693  {
4694    intvec *iv;
4695    iv = (intvec *)(sl->data);
4696    if ((*iv)[1]!=ringorder_L)
4697    {
4698      n++;
4699
4700      /* the format of an ordering:
4701       *  iv[0]: factor
4702       *  iv[1]: ordering
4703       *  iv[2..end]: weights
4704       */
4705      R->order[n] = (*iv)[1];
4706      typ=1;
4707      switch ((*iv)[1])
4708      {
4709          case ringorder_ws:
4710          case ringorder_Ws:
4711            typ=-1;
4712          case ringorder_wp:
4713          case ringorder_Wp:
4714            R->wvhdl[n]=(int*)omAlloc((iv->length()-1)*sizeof(int));
4715            R->block0[n] = last+1;
4716            for (i=2; i<iv->length(); i++)
4717            {
4718              R->wvhdl[n][i-2] = (*iv)[i];
4719              last++;
4720              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4721            }
4722            R->block1[n] = last;
4723            break;
4724          case ringorder_ls:
4725          case ringorder_ds:
4726          case ringorder_Ds:
4727          case ringorder_rs:
4728            typ=-1;
4729          case ringorder_lp:
4730          case ringorder_dp:
4731          case ringorder_Dp:
4732          case ringorder_rp:
4733            R->block0[n] = last+1;
4734            if (iv->length() == 3) last+=(*iv)[2];
4735            else last += (*iv)[0];
4736            R->block1[n] = last;
4737            //if ((R->block0[n]>R->block1[n])
4738            //|| (R->block1[n]>rVar(R)))
4739            //{
4740            //  R->block1[n]=rVar(R);
4741            //  //WerrorS("ordering larger than number of variables");
4742            //  break;
4743            //}
4744            if (rCheckIV(iv)) return TRUE;
4745            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4746            {
4747              if (weights[i]==0) weights[i]=typ;
4748            }
4749            break;
4750
4751          case ringorder_s: // no 'rank' params!
4752          {
4753
4754            if(iv->length() > 3)
4755              return TRUE;
4756
4757            if(iv->length() == 3)
4758            {
4759              const int s = (*iv)[2];
4760              R->block0[n] = s;
4761              R->block1[n] = s;
4762            }
4763            break;
4764          }
4765          case ringorder_IS:
4766          {
4767            if(iv->length() != 3) return TRUE;
4768
4769            const int s = (*iv)[2];
4770
4771            if( 1 < s || s < -1 ) return TRUE;
4772
4773            R->block0[n] = s;
4774            R->block1[n] = s;
4775            break;
4776          }
4777          case ringorder_S:
4778          case ringorder_c:
4779          case ringorder_C:
4780          {
4781            if (rCheckIV(iv)) return TRUE;
4782            break;
4783          }
4784          case ringorder_aa:
4785          case ringorder_a:
4786          {
4787            R->block0[n] = last+1;
4788            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4789            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int));
4790            for (i=2; i<iv->length(); i++)
4791            {
4792              R->wvhdl[n][i-2]=(*iv)[i];
4793              last++;
4794              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4795            }
4796            last=R->block0[n]-1;
4797            break;
4798          }
4799          case ringorder_a64:
4800          {
4801            R->block0[n] = last+1;
4802            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4803            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int64));
4804            int64 *w=(int64 *)R->wvhdl[n];
4805            for (i=2; i<iv->length(); i++)
4806            {
4807              w[i-2]=(*iv)[i];
4808              last++;
4809              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4810            }
4811            last=R->block0[n]-1;
4812            break;
4813          }
4814          case ringorder_M:
4815          {
4816            int Mtyp=rTypeOfMatrixOrder(iv);
4817            if (Mtyp==0) return TRUE;
4818            if (Mtyp==-1) typ = -1;
4819
4820            R->wvhdl[n] =( int *)omAlloc((iv->length()-1)*sizeof(int));
4821            for (i=2; i<iv->length();i++)
4822              R->wvhdl[n][i-2]=(*iv)[i];
4823
4824            R->block0[n] = last+1;
4825            last += (int)sqrt((double)(iv->length()-2));
4826            R->block1[n] = last;
4827            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4828            {
4829              if (weights[i]==0) weights[i]=typ;
4830            }
4831            break;
4832          }
4833
4834          case ringorder_no:
4835            R->order[n] = ringorder_unspec;
4836            return TRUE;
4837
4838          default:
4839            Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
4840            R->order[n] = ringorder_unspec;
4841            return TRUE;
4842      }
4843    }
4844    sl=sl->next;
4845  }
4846
4847  // check for complete coverage
4848  while ( n >= 0 && (
4849          (R->order[n]==ringorder_c)
4850      ||  (R->order[n]==ringorder_C)
4851      ||  (R->order[n]==ringorder_s)
4852      ||  (R->order[n]==ringorder_S)
4853      ||  (R->order[n]==ringorder_IS)
4854                    )) n--;
4855
4856  assume( n >= 0 );
4857
4858  if (R->block1[n] != R->N)
4859  {
4860    if (((R->order[n]==ringorder_dp) ||
4861         (R->order[n]==ringorder_ds) ||
4862         (R->order[n]==ringorder_Dp) ||
4863         (R->order[n]==ringorder_Ds) ||
4864         (R->order[n]==ringorder_rp) ||
4865         (R->order[n]==ringorder_rs) ||
4866         (R->order[n]==ringorder_lp) ||
4867         (R->order[n]==ringorder_ls))
4868        &&
4869        R->block0[n] <= R->N)
4870    {
4871      R->block1[n] = R->N;
4872    }
4873    else
4874    {
4875      Werror("mismatch of number of vars (%d) and ordering (%d vars)",
4876             R->N,R->block1[n]);
4877      return TRUE;
4878    }
4879  }
4880  // find OrdSgn:
4881  R->OrdSgn = 1;
4882  for(i=1;i<=R->N;i++)
4883  { if (weights[i]<0) { R->OrdSgn=-1;break; }}
4884  omFree(weights);
4885  return FALSE;
4886}
4887
4888BOOLEAN rSleftvList2StringArray(sleftv* sl, char** p)
4889{
4890
4891  while(sl!=NULL)
4892  {
4893    if (sl->Name() == sNoName)
4894    {
4895      if (sl->Typ()==POLY_CMD)
4896      {
4897        sleftv s_sl;
4898        iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
4899        if (s_sl.Name() != sNoName)
4900          *p = omStrDup(s_sl.Name());
4901        else
4902          *p = NULL;
4903        sl->next = s_sl.next;
4904        s_sl.next = NULL;
4905        s_sl.CleanUp();
4906        if (*p == NULL) return TRUE;
4907      }
4908      else
4909        return TRUE;
4910    }
4911    else
4912      *p = omStrDup(sl->Name());
4913    p++;
4914    sl=sl->next;
4915  }
4916  return FALSE;
4917}
4918
4919const short MAX_SHORT = SHRT_MAX; // (1 << (sizeof(short)*8)) - 1;
4920
4921////////////////////
4922//
4923// rInit itself:
4924//
4925// INPUT:  s: name, pn: ch & parameter (names), rv: variable (names)
4926//         ord: ordering
4927// RETURN: currRingHdl on success
4928//         NULL        on error
4929// NOTE:   * makes new ring to current ring, on success
4930//         * considers input sleftv's as read-only
4931//idhdl rInit(char *s, sleftv* pn, sleftv* rv, sleftv* ord)
4932ring rInit(sleftv* pn, sleftv* rv, sleftv* ord)
4933{
4934#ifdef HAVE_RINGS
4935  unsigned int ringtype = 0;
4936  int_number modBase = NULL;
4937  unsigned int modExponent = 1;
4938#endif
4939  int float_len=0;
4940  int float_len2=0;
4941  ring R = NULL;
4942  idhdl tmp = NULL;
4943  BOOLEAN ffChar=FALSE;
4944  int typ = 1;
4945
4946  /* ch -------------------------------------------------------*/
4947  // get ch of ground field
4948  int numberOfAllocatedBlocks;
4949 
4950  // allocated ring
4951  R = (ring) omAlloc0Bin(sip_sring_bin);
4952
4953  coeffs cf = NULL;
4954
4955  assume( pn != NULL );
4956  const int P = pn->listLength();
4957 
4958  if (pn->Typ()==INT_CMD)
4959  {
4960    int ch = (int)(long)pn->Data();
4961   
4962    /* parameter? -------------------------------------------------------*/
4963    pn = pn->next;
4964   
4965    if (pn == NULL) // no params!?
4966    {
4967      if (ch!=0)
4968      {
4969        ch=IsPrime(ch);
4970        cf = nInitChar(n_Zp, (void*)(long)ch);
4971      }
4972      else
4973        cf = nInitChar(n_Q, (void*)(long)ch);
4974    }
4975    else
4976    {
4977      const int pars = pn->listLength();
4978
4979      assume( pars > 0 );
4980
4981      // predefined finite field: (p^k, a)
4982      if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
4983      {
4984        GFInfo param;
4985         
4986        param.GFChar = ch; 
4987        param.GFDegree = 1;     
4988        param.GFPar_name = pn->name;
4989       
4990        cf = nInitChar(n_GF, &param);       
4991      }
4992      else // (0/p, a, b, ..., z)
4993      {
4994        assume( (ch == 0) || (ch==IsPrime(ch)) );
4995
4996//         if ((pars > 1) && (ffChar))
4997//         {
4998//           WerrorS("too many parameters");
4999//           goto rInitError;
5000//         }
5001       
5002        char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5003
5004        if (rSleftvList2StringArray(pn, names))
5005        {
5006          WerrorS("parameter expected");
5007          goto rInitError;
5008        }
5009
5010        TransExtInfo extParam;
5011       
5012        extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5013       
5014        cf = nInitChar(n_transExt, &extParam);
5015      }
5016    }
5017   
5018//    if (cf==NULL) goto rInitError;
5019    assume( cf != NULL );
5020  }
5021  else if ((pn->name != NULL)
5022  && ((strcmp(pn->name,"real")==0) || (strcmp(pn->name,"complex")==0)))
5023  {
5024    BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5025    const int ch=0;
5026    if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5027    {
5028      WarnS("not implemented: size for real/complex");
5029      float_len=(int)(long)pn->next->Data();
5030      float_len2=float_len;
5031      pn=pn->next;
5032      if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5033      {
5034        float_len2=(int)(long)pn->next->Data();
5035        WarnS("not implemented: size for real/complex");
5036        pn=pn->next;
5037      }
5038    }
5039    if ((pn->next==NULL) && complex_flag)
5040    {
5041      pn->next=(leftv)omAlloc0Bin(sleftv_bin);
5042      pn->next->name=omStrDup("i");
5043    }
5044    else
5045      WarnS("not implemented: name for i (complex)");
5046    cf=nInitChar(complex_flag ? n_long_C: n_long_R,NULL);
5047
5048    assume( cf != NULL );
5049  }
5050#ifdef HAVE_RINGS
5051  else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5052  {
5053    modBase = (int_number) omAlloc(sizeof(mpz_t));
5054    mpz_init_set_si(modBase, 0);
5055    if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5056    {
5057      mpz_set_ui(modBase, (int)(long) pn->next->Data());
5058      pn=pn->next;
5059      if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5060      {
5061        modExponent = (long) pn->next->Data();
5062        pn=pn->next;
5063      }
5064      while ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5065      {
5066        mpz_mul_ui(modBase, modBase, (int)(long) pn->next->Data());
5067        pn=pn->next;
5068      }
5069    }
5070    else
5071      cf=nInitChar(n_Z,NULL);
5072   
5073    if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_cmp_ui(modBase, 0) < 0))
5074    {
5075      Werror("Wrong ground ring specification (module is 1)");
5076      goto rInitError;
5077    }
5078    if (modExponent < 1)
5079    {
5080      Werror("Wrong ground ring specification (exponent smaller than 1");
5081      goto rInitError;
5082    }
5083    // module is 0 ---> integers ringtype = 4;
5084    // we have an exponent
5085    if (modExponent > 1)
5086    {
5087      const int ch = modExponent;
5088      if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(NATNUMBER)))
5089      {
5090        /* this branch should be active for modExponent = 2..32 resp. 2..64,
5091           depending on the size of a long on the respective platform */
5092        ringtype = 1;       // Use Z/2^ch
5093        cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5094      }
5095      else
5096      {
5097        ringtype = 3;
5098        cf=nInitChar(n_Zn,(void*)(long)modBase);
5099      }
5100    }
5101    // just a module m > 1
5102    else
5103    {
5104      ringtype = 2;
5105      const int ch = mpz_get_ui(modBase);
5106      cf=nInitChar(n_Zn,(void*)(long)ch);
5107    }
5108    assume( cf != NULL );
5109  }
5110#endif
5111  // ring NEW = OLD, (), (); where OLD is a polynomial ring...
5112  else if ((pn->Typ()==RING_CMD) && (P == 1))
5113  {
5114    TransExtInfo extParam;
5115    extParam.r = (ring)pn->Data();
5116    cf = nInitChar(n_transExt, &extParam);
5117  }
5118  else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5119  {
5120    AlgExtInfo extParam;
5121    extParam.r = (ring)pn->Data();
5122    extParam.i = (extParam.r->qideal);
5123   
5124    cf = nInitChar(n_algExt, &extParam);   // Q[a]/<minideal>
5125  }
5126  else
5127  {
5128    sleftv* p = pn;
5129    Werror("Wrong or unknown ground field specification");
5130#ifndef NDEBUG
5131    while (p != NULL)
5132    {
5133      Print( "pn[%p]: type: %d [%s]: %p, name: %s", p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5134      PrintLn();
5135      p = p->next;
5136    }
5137#endif   
5138    goto rInitError;
5139  }
5140//  pn=pn->next;
5141
5142  int l, last;
5143  sleftv * sl;
5144  /*every entry in the new ring is initialized to 0*/
5145
5146  /* characteristic -----------------------------------------------*/
5147  /* input: 0 ch=0 : Q     parameter=NULL    ffChar=FALSE   float_len
5148   *         0    1 : Q(a,...)        *names         FALSE
5149   *         0   -1 : R               NULL           FALSE  0
5150   *         0   -1 : R               NULL           FALSE  prec. >6
5151   *         0   -1 : C               *names         FALSE  prec. 0..?
5152   *         p    p : Fp              NULL           FALSE
5153   *         p   -p : Fp(a)           *names         FALSE
5154   *         q    q : GF(q=p^n)       *names         TRUE
5155  */
5156  l = 0;
5157
5158  if (cf==NULL)
5159  {
5160    const int ch=32003;
5161    Warn("Invalid ground field specification: using the default field: Z_{%d}", ch);
5162    cf=nInitChar(n_Zp, (void*)(long)ch);
5163  }
5164
5165  assume( R != NULL );
5166 
5167  R->cf = cf;
5168 
5169#ifdef HAVE_RINGS
5170  // the following should have beed set already into cf, right?!
5171//  R->cf->ringtype = ringtype;
5172//  R->cf->modBase = modBase;
5173//  R->cf->modExponent = modExponent;
5174#endif
5175
5176  /* names and number of variables-------------------------------------*/
5177  {
5178    int l=rv->listLength();
5179
5180    if (l>MAX_SHORT)
5181    {
5182      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5183       goto rInitError;
5184    }
5185    R->N = l; /*rv->listLength();*/
5186  }
5187  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5188  if (rSleftvList2StringArray(rv, R->names))
5189  {
5190    WerrorS("name of ring variable expected");
5191    goto rInitError;
5192  }
5193
5194  /* check names and parameters for conflicts ------------------------- */
5195  rRenameVars(R); // conflicting variables will be renamed
5196  /* ordering -------------------------------------------------------------*/
5197  if (rSleftvOrdering2Ordering(ord, R))
5198    goto rInitError;
5199
5200  // Complete the initialization
5201  if (rComplete(R,1))
5202    goto rInitError;
5203
5204#ifdef HABE_RINGS
5205// currently, coefficients which are ring elements require a global ordering:
5206  if (rField_is_Ring(R) && (R->OrdSgn==-1))
5207  {
5208    WerrorS("global ordering required for these coefficients");
5209    goto rInitError;
5210  }
5211#endif
5212
5213  rTest(R);
5214
5215  // try to enter the ring into the name list
5216  // need to clean up sleftv here, before this ring can be set to
5217  // new currRing or currRing can be killed beacuse new ring has
5218  // same name
5219  if (pn != NULL) pn->CleanUp();
5220  if (rv != NULL) rv->CleanUp();
5221  if (ord != NULL) ord->CleanUp();
5222  //if ((tmp = enterid(s, myynest, RING_CMD, &IDROOT))==NULL)
5223  //  goto rInitError;
5224
5225  //memcpy(IDRING(tmp),R,sizeof(*R));
5226  // set current ring
5227  //omFreeBin(R,  ip_sring_bin);
5228  //return tmp;
5229  return R;
5230
5231  // error case:
5232  rInitError:
5233  if  ((R != NULL)&&(R->cf!=NULL)) rDelete(R);
5234  if (pn != NULL) pn->CleanUp();
5235  if (rv != NULL) rv->CleanUp();
5236  if (ord != NULL) ord->CleanUp();
5237  return NULL;
5238}
5239
5240ring rSubring(ring org_ring, sleftv* rv)
5241{
5242  ring R = rCopy0(org_ring);
5243  int *perm=(int *)omAlloc0((org_ring->N+1)*sizeof(int));
5244  int last = 0, o=0, n = rBlocks(org_ring), i=0, typ = 1, j;
5245
5246  /* names and number of variables-------------------------------------*/
5247  {
5248    int l=rv->listLength();
5249    if (l>MAX_SHORT)
5250    {
5251      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5252       goto rInitError;
5253    }
5254    R->N = l; /*rv->listLength();*/
5255  }
5256  omFree(R->names);
5257  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5258  if (rSleftvList2StringArray(rv, R->names))
5259  {
5260    WerrorS("name of ring variable expected");
5261    goto rInitError;
5262  }
5263
5264  /* check names for subring in org_ring ------------------------- */
5265  {
5266    i=0;
5267
5268    for(j=0;j<R->N;j++)
5269    {
5270      for(;i<org_ring->N;i++)
5271      {
5272        if (strcmp(org_ring->names[i],R->names[j])==0)
5273        {
5274          perm[i+1]=j+1;
5275          break;
5276        }
5277      }
5278      if (i>org_ring->N)
5279      {
5280        Werror("variable %d (%s) not in basering",j+1,R->names[j]);
5281        break;
5282      }
5283    }
5284  }
5285  //Print("perm=");
5286  //for(i=1;i<org_ring->N;i++) Print("v%d -> v%d\n",i,perm[i]);
5287  /* ordering -------------------------------------------------------------*/
5288
5289  for(i=0;i<n;i++)
5290  {
5291    int min_var=-1;
5292    int max_var=-1;
5293    for(j=R->block0[i];j<=R->block1[i];j++)
5294    {
5295      if (perm[j]>0)
5296      {
5297        if (min_var==-1) min_var=perm[j];
5298        max_var=perm[j];
5299      }
5300    }
5301    if (min_var!=-1)
5302    {
5303      //Print("block %d: old %d..%d, now:%d..%d\n",
5304      //      i,R->block0[i],R->block1[i],min_var,max_var);
5305      R->block0[i]=min_var;
5306      R->block1[i]=max_var;
5307      if (R->wvhdl[i]!=NULL)
5308      {
5309        omFree(R->wvhdl[i]);
5310        R->wvhdl[i]=(int*)omAlloc0((max_var-min_var+1)*sizeof(int));
5311        for(j=org_ring->block0[i];j<=org_ring->block1[i];j++)
5312        {
5313          if (perm[j]>0)
5314          {
5315            R->wvhdl[i][perm[j]-R->block0[i]]=
5316                org_ring->wvhdl[i][j-org_ring->block0[i]];
5317            //Print("w%d=%d (orig_w%d)\n",perm[j],R->wvhdl[i][perm[j]-R->block0[i]],j);
5318          }
5319        }
5320      }
5321    }
5322    else
5323    {
5324      if(R->block0[i]>0)
5325      {
5326        //Print("skip block %d\n",i);
5327        R->order[i]=ringorder_unspec;
5328        if (R->wvhdl[i] !=NULL) omFree(R->wvhdl[i]);
5329        R->wvhdl[i]=NULL;
5330      }
5331      //else Print("keep block %d\n",i);
5332    }
5333  }
5334  i=n-1;
5335  while(i>0)
5336  {
5337    // removed unneded blocks
5338    if(R->order[i-1]==ringorder_unspec)
5339    {
5340      for(j=i;j<=n;j++)
5341      {
5342        R->order[j-1]=R->order[j];
5343        R->block0[j-1]=R->block0[j];
5344        R->block1[j-1]=R->block1[j];
5345        if (R->wvhdl[j-1] !=NULL) omFree(R->wvhdl[j-1]);
5346        R->wvhdl[j-1]=R->wvhdl[j];
5347      }
5348      R->order[n]=ringorder_unspec;
5349      n--;
5350    }
5351    i--;
5352  }
5353  n=rBlocks(org_ring)-1;
5354  while (R->order[n]==0)  n--;
5355  while (R->order[n]==ringorder_unspec)  n--;
5356  if ((R->order[n]==ringorder_c) ||  (R->order[n]==ringorder_C)) n--;
5357  if (R->block1[n] != R->N)
5358  {
5359    if (((R->order[n]==ringorder_dp) ||
5360         (R->order[n]==ringorder_ds) ||
5361         (R->order[n]==ringorder_Dp) ||
5362         (R->order[n]==ringorder_Ds) ||
5363         (R->order[n]==ringorder_rp) ||
5364         (R->order[n]==ringorder_rs) ||
5365         (R->order[n]==ringorder_lp) ||
5366         (R->order[n]==ringorder_ls))
5367        &&
5368        R->block0[n] <= R->N)
5369    {
5370      R->block1[n] = R->N;
5371    }
5372    else
5373    {
5374      Werror("mismatch of number of vars (%d) and ordering (%d vars) in block %d",
5375             R->N,R->block1[n],n);
5376      return NULL;
5377    }
5378  }
5379  omFree(perm);
5380  // find OrdSgn:
5381  R->OrdSgn = org_ring->OrdSgn; // IMPROVE!
5382  //for(i=1;i<=R->N;i++)
5383  //{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
5384  //omFree(weights);
5385  // Complete the initialization
5386  if (rComplete(R,1))
5387    goto rInitError;
5388
5389  rTest(R);
5390
5391  if (rv != NULL) rv->CleanUp();
5392
5393  return R;
5394
5395  // error case:
5396  rInitError:
5397  if  (R != NULL) rDelete(R);
5398  if (rv != NULL) rv->CleanUp();
5399  return NULL;
5400}
5401
5402void rKill(ring r)
5403{
5404  if ((r->ref<=0)&&(r->order!=NULL))
5405  {
5406#ifdef RDEBUG
5407    if (traceit &TRACE_SHOW_RINGS) Print("kill ring %lx\n",(long)r);
5408#endif
5409    if (r->qideal!=NULL)
5410    {
5411      id_Delete(&r->qideal, r);
5412      r->qideal = NULL;
5413    }
5414    int i=1;
5415    int j;
5416    int *pi=r->order;
5417#ifdef USE_IILOCALRING
5418    for (j=0;j<iiRETURNEXPR_len;j++)
5419    {
5420      if (iiLocalRing[j]==r)
5421      {
5422        if (j<myynest) Warn("killing the basering for level %d",j);
5423        iiLocalRing[j]=NULL;
5424      }
5425    }
5426#else /* USE_IILOCALRING */
5427//#endif /* USE_IILOCALRING */
5428    {
5429      proclevel * nshdl = procstack;
5430      int lev=myynest-1;
5431
5432      for(; nshdl != NULL; nshdl = nshdl->next)
5433      {
5434        if (nshdl->cRing==r)
5435        {
5436          Warn("killing the basering for level %d",lev);
5437          nshdl->cRing=NULL;
5438          nshdl->cRingHdl=NULL;
5439        }
5440      }
5441    }
5442#endif /* USE_IILOCALRING */
5443// any variables depending on r ?
5444    while (r->idroot!=NULL)
5445    {
5446      killhdl2(r->idroot,&(r->idroot),r);
5447    }
5448    if (r==currRing)
5449    {
5450      // all dependend stuff is done, clean global vars:
5451      if (r->qideal!=NULL)
5452      {
5453        currQuotient=NULL;
5454      }
5455      if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
5456      if (sLastPrinted.RingDependend())
5457      {
5458        sLastPrinted.CleanUp();
5459      }
5460      if ((myynest>0) && (iiRETURNEXPR[myynest].RingDependend()))
5461      {
5462        WerrorS("return value depends on local ring variable (export missing ?)");
5463        iiRETURNEXPR[myynest].CleanUp();
5464      }
5465      currRing=NULL;
5466      currRingHdl=NULL;
5467    }
5468
5469    /* nKillChar(r); will be called from inside of rDelete */
5470    rDelete(r);
5471    return;
5472  }
5473  r->ref--;
5474}
5475
5476void rKill(idhdl h)
5477{
5478  ring r = IDRING(h);
5479  int ref=0;
5480  if (r!=NULL)
5481  {
5482    ref=r->ref;
5483    rKill(r);
5484  }
5485  if (h==currRingHdl)
5486  {
5487    if (ref<=0) { currRing=NULL; currRingHdl=NULL;}
5488    else
5489    {
5490      currRingHdl=rFindHdl(r,currRingHdl,NULL);
5491    }
5492  }
5493}
5494
5495idhdl rSimpleFindHdl(ring r, idhdl root, idhdl n)
5496{
5497  //idhdl next_best=NULL;
5498  idhdl h=root;
5499  while (h!=NULL)
5500  {
5501    if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
5502    && (h!=n)
5503    && (IDRING(h)==r)
5504    )
5505    {
5506   //   if (IDLEV(h)==myynest)
5507   //     return h;
5508   //   if ((IDLEV(h)==0) || (next_best==NULL))
5509   //     next_best=h;
5510   //   else if (IDLEV(next_best)<IDLEV(h))
5511   //     next_best=h;
5512      return h;
5513    }
5514    h=IDNEXT(h);
5515  }
5516  //return next_best;
5517  return NULL;
5518}
5519
5520extern BOOLEAN jjPROC(leftv res, leftv u, leftv v);
5521ideal kGroebner(ideal F, ideal Q)
5522{
5523  //test|=Sy_bit(OPT_PROT);
5524  idhdl save_ringhdl=currRingHdl;
5525  ideal resid;
5526  idhdl new_ring=NULL;
5527  if ((currRingHdl==NULL) || (IDRING(currRingHdl)!=currRing))
5528  {
5529    currRingHdl=enterid(omStrDup(" GROEBNERring"),0,RING_CMD,&IDROOT,FALSE);
5530    new_ring=currRingHdl;
5531    IDRING(currRingHdl)=currRing;
5532  }
5533  sleftv v; memset(&v,0,sizeof(v)); v.rtyp=IDEAL_CMD; v.data=(char *) F;
5534  idhdl h=ggetid("groebner");
5535  sleftv u; memset(&u,0,sizeof(u)); u.rtyp=IDHDL; u.data=(char *) h;
5536            u.name=IDID(h);
5537
5538  sleftv res; memset(&res,0,sizeof(res));
5539  if(jjPROC(&res,&u,&v))
5540  {
5541    resid=kStd(F,Q,testHomog,NULL);
5542  }
5543  else
5544  {
5545    //printf("typ:%d\n",res.rtyp);
5546    resid=(ideal)(res.data);
5547  }
5548  // cleanup GROEBNERring, save_ringhdl, u,v,(res )
5549  if (new_ring!=NULL)
5550  {
5551    idhdl h=IDROOT;
5552    if (h==new_ring) IDROOT=h->next;
5553    else
5554    {
5555      while ((h!=NULL) &&(h->next!=new_ring)) h=h->next;
5556      if (h!=NULL) h->next=h->next->next;
5557    }
5558    if (h!=NULL) omFreeSize(h,sizeof(*h));
5559  }
5560  currRingHdl=save_ringhdl;
5561  u.CleanUp();
5562  v.CleanUp();
5563  return resid;
5564}
5565
5566static void jjINT_S_TO_ID(int n,int *e, leftv res)
5567{
5568  if (n==0) n=1;
5569  ideal l=idInit(n,1);
5570  int i;
5571  poly p;
5572  for(i=rVar(currRing);i>0;i--)
5573  {
5574    if (e[i]>0)
5575    {
5576      n--;
5577      p=pOne();
5578      pSetExp(p,i,1);
5579      pSetm(p);
5580      l->m[n]=p;
5581      if (n==0) break;
5582    }
5583  }
5584  res->data=(char*)l;
5585  setFlag(res,FLAG_STD);
5586  omFreeSize((ADDRESS)e,(rVar(currRing)+1)*sizeof(int));
5587}
5588BOOLEAN jjVARIABLES_P(leftv res, leftv u)
5589{
5590  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5591  int n=pGetVariables((poly)u->Data(),e);
5592  jjINT_S_TO_ID(n,e,res);
5593  return FALSE;
5594}
5595
5596BOOLEAN jjVARIABLES_ID(leftv res, leftv u)
5597{
5598  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5599  ideal I=(ideal)u->Data();
5600  int i;
5601  int n=0;
5602  for(i=I->nrows*I->ncols-1;i>=0;i--)
5603  {
5604    n=pGetVariables(I->m[i],e);
5605  }
5606  jjINT_S_TO_ID(n,e,res);
5607  return FALSE;
5608}
Note: See TracBrowser for help on using the repository browser.