source: git/Singular/ipshell.cc @ 69d7df

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