source: git/Singular/ipshell.cc @ a8941bc

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