source: git/Singular/ipshell.cc @ d143376

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