source: git/Singular/ipshell.cc @ 70b3ae

spielwiese
Last change on this file since 70b3ae was 70b3ae, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: ring constructions with wrong coeffs (tr.517)
  • Property mode set to 100644
File size: 141.6 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 (check_comp)
2444    {
2445      BOOLEAN comp_order=FALSE;
2446      int jj;
2447      for(jj=0;jj<n;jj++)
2448      {
2449        if ((R->order[jj]==ringorder_c) ||
2450            (R->order[jj]==ringorder_C)) { comp_order=TRUE; break; }
2451      }
2452      if (!comp_order)
2453      {
2454        R->order=(int*)omRealloc0Size(R->order,n*sizeof(int),(n+1)*sizeof(int));
2455        R->block0=(int*)omRealloc0Size(R->block0,n*sizeof(int),(n+1)*sizeof(int));
2456        R->block1=(int*)omRealloc0Size(R->block1,n*sizeof(int),(n+1)*sizeof(int));
2457        R->wvhdl=(int**)omRealloc0Size(R->wvhdl,n*sizeof(int_ptr),(n+1)*sizeof(int_ptr));
2458        R->order[n-1]=ringorder_C;
2459        R->block0[n-1]=0;
2460        R->block1[n-1]=0;
2461        R->wvhdl[n-1]=NULL;
2462        n++;
2463      }
2464    }
2465  }
2466  else
2467  {
2468    WerrorS("ordering must be given as `list`");
2469    goto rCompose_err;
2470  }
2471
2472  // ------------------------ ??????? --------------------
2473
2474  rRenameVars(R);
2475  rComplete(R);
2476
2477/*#ifdef HAVE_RINGS
2478// currently, coefficients which are ring elements require a global ordering:
2479  if (rField_is_Ring(R) && (R->OrdSgn==-1))
2480  {
2481    WerrorS("global ordering required for these coefficients");
2482    goto rCompose_err;
2483  }
2484#endif*/
2485
2486
2487  // ------------------------ Q-IDEAL ------------------------
2488
2489  if (L->m[3].Typ()==IDEAL_CMD)
2490  {
2491    ideal q=(ideal)L->m[3].Data();
2492    if (q->m[0]!=NULL)
2493    {
2494      if (R->cf != currRing->cf) //->cf->ch!=currRing->cf->ch)
2495      {
2496      #if 0
2497            WerrorS("coefficient fields must be equal if q-ideal !=0");
2498            goto rCompose_err;
2499      #else
2500        ring orig_ring=currRing;
2501        rChangeCurrRing(R);
2502        int *perm=NULL;
2503        int *par_perm=NULL;
2504        int par_perm_size=0;
2505        nMapFunc nMap;
2506
2507        if ((nMap=nSetMap(orig_ring->cf))==NULL)
2508        {
2509          if (rEqual(orig_ring,currRing))
2510          {
2511            nMap=n_SetMap(currRing->cf, currRing->cf);
2512          }
2513          else
2514          // Allow imap/fetch to be make an exception only for:
2515          if ( (rField_is_Q_a(orig_ring) &&  // Q(a..) -> Q(a..) || Q || Zp || Zp(a)
2516            (rField_is_Q(currRing) || rField_is_Q_a(currRing) ||
2517             (rField_is_Zp(currRing) || rField_is_Zp_a(currRing))))
2518           ||
2519           (rField_is_Zp_a(orig_ring) &&  // Zp(a..) -> Zp(a..) || Zp
2520            (rField_is_Zp(currRing, rInternalChar(orig_ring)) ||
2521             rField_is_Zp_a(currRing, rInternalChar(orig_ring)))) )
2522          {
2523            par_perm_size=rPar(orig_ring);
2524
2525//            if ((orig_ring->minpoly != NULL) || (orig_ring->qideal != NULL))
2526//              naSetChar(rInternalChar(orig_ring),orig_ring);
2527//            else ntSetChar(rInternalChar(orig_ring),orig_ring);
2528
2529            nSetChar(currRing->cf);
2530          }
2531          else
2532          {
2533            WerrorS("coefficient fields must be equal if q-ideal !=0");
2534            goto rCompose_err;
2535          }
2536        }
2537        perm=(int *)omAlloc0((orig_ring->N+1)*sizeof(int));
2538        if (par_perm_size!=0)
2539          par_perm=(int *)omAlloc0(par_perm_size*sizeof(int));
2540        int i;
2541        #if 0
2542        // use imap:
2543        maFindPerm(orig_ring->names,orig_ring->N,orig_ring->parameter,orig_ring->P,
2544          currRing->names,currRing->N,currRing->parameter, currRing->P,
2545          perm,par_perm, currRing->ch);
2546        #else
2547        // use fetch
2548        if ((rPar(orig_ring)>0) && (rPar(currRing)==0))
2549        {
2550          for(i=si_min(rPar(orig_ring),rVar(currRing))-1;i>=0;i--) par_perm[i]=i+1;
2551        }
2552        else if (par_perm_size!=0)
2553          for(i=si_min(rPar(orig_ring),rPar(currRing))-1;i>=0;i--) par_perm[i]=-(i+1);
2554        for(i=si_min(orig_ring->N,rVar(currRing));i>0;i--) perm[i]=i;
2555        #endif
2556        ideal dest_id=idInit(IDELEMS(q),1);
2557        for(i=IDELEMS(q)-1; i>=0; i--)
2558        {
2559          dest_id->m[i]=p_PermPoly(q->m[i],perm,orig_ring, currRing,nMap,
2560                                  par_perm,par_perm_size);
2561          //  PrintS("map:");pWrite(dest_id->m[i]);PrintLn();
2562          pTest(dest_id->m[i]);
2563        }
2564        R->qideal=dest_id;
2565        if (perm!=NULL)
2566          omFreeSize((ADDRESS)perm,(orig_ring->N+1)*sizeof(int));
2567        if (par_perm!=NULL)
2568          omFreeSize((ADDRESS)par_perm,par_perm_size*sizeof(int));
2569        rChangeCurrRing(orig_ring);
2570      #endif
2571      }
2572      else
2573        R->qideal=idrCopyR(q,currRing,R);
2574    }
2575  }
2576  else
2577  {
2578    WerrorS("q-ideal must be given as `ideal`");
2579    goto rCompose_err;
2580  }
2581
2582
2583  // ---------------------------------------------------------------
2584  #ifdef HAVE_PLURAL
2585  if (L->nr==5)
2586  {
2587    if (nc_CallPlural((matrix)L->m[4].Data(),
2588                      (matrix)L->m[5].Data(),
2589                      NULL,NULL,
2590                      R,
2591                      true, // !!!
2592                      true, false,
2593                      currRing, FALSE)) goto rCompose_err;
2594    // takes care about non-comm. quotient! i.e. calls "nc_SetupQuotient" due to last true
2595  }
2596  #endif
2597  return R;
2598
2599rCompose_err:
2600  if (R->N>0)
2601  {
2602    int i;
2603    if (R->names!=NULL)
2604    {
2605      i=R->N-1;
2606      while (i>=0) { if (R->names[i]!=NULL) omFree(R->names[i]); i--; }
2607      omFree(R->names);
2608    }
2609  }
2610  if (R->order!=NULL) omFree(R->order);
2611  if (R->block0!=NULL) omFree(R->block0);
2612  if (R->block1!=NULL) omFree(R->block1);
2613  if (R->wvhdl!=NULL) omFree(R->wvhdl);
2614  omFree(R);
2615  return NULL;
2616}
2617
2618// from matpol.cc
2619
2620/*2
2621* compute the jacobi matrix of an ideal
2622*/
2623BOOLEAN mpJacobi(leftv res,leftv a)
2624{
2625  int     i,j;
2626  matrix result;
2627  ideal id=(ideal)a->Data();
2628
2629  result =mpNew(IDELEMS(id),rVar(currRing));
2630  for (i=1; i<=IDELEMS(id); i++)
2631  {
2632    for (j=1; j<=rVar(currRing); j++)
2633    {
2634      MATELEM(result,i,j) = pDiff(id->m[i-1],j);
2635    }
2636  }
2637  res->data=(char *)result;
2638  return FALSE;
2639}
2640
2641/*2
2642* returns the Koszul-matrix of degree d of a vectorspace with dimension n
2643* uses the first n entrees of id, if id <> NULL
2644*/
2645BOOLEAN mpKoszul(leftv res,leftv c/*ip*/, leftv b/*in*/, leftv id)
2646{
2647  int n=(int)(long)b->Data();
2648  int d=(int)(long)c->Data();
2649  int     k,l,sign,row,col;
2650  matrix  result;
2651  ideal temp;
2652  BOOLEAN bo;
2653  poly    p;
2654
2655  if ((d>n) || (d<1) || (n<1))
2656  {
2657    res->data=(char *)mpNew(1,1);
2658    return FALSE;
2659  }
2660  int *choise = (int*)omAlloc(d*sizeof(int));
2661  if (id==NULL)
2662    temp=idMaxIdeal(1);
2663  else
2664    temp=(ideal)id->Data();
2665
2666  k = binom(n,d);
2667  l = k*d;
2668  l /= n-d+1;
2669  result =mpNew(l,k);
2670  col = 1;
2671  idInitChoise(d,1,n,&bo,choise);
2672  while (!bo)
2673  {
2674    sign = 1;
2675    for (l=1;l<=d;l++)
2676    {
2677      if (choise[l-1]<=IDELEMS(temp))
2678      {
2679        p = pCopy(temp->m[choise[l-1]-1]);
2680        if (sign == -1) p = pNeg(p);
2681        sign *= -1;
2682        row = idGetNumberOfChoise(l-1,d,1,n,choise);
2683        MATELEM(result,row,col) = p;
2684      }
2685    }
2686    col++;
2687    idGetNextChoise(d,n,&bo,choise);
2688  }
2689  if (id==NULL) idDelete(&temp);
2690
2691  res->data=(char *)result;
2692  return FALSE;
2693}
2694
2695// from syz1.cc
2696/*2
2697* read out the Betti numbers from resolution
2698* (interpreter interface)
2699*/
2700BOOLEAN syBetti2(leftv res, leftv u, leftv w)
2701{
2702  syStrategy syzstr=(syStrategy)u->Data();
2703
2704  BOOLEAN minim=(int)(long)w->Data();
2705  int row_shift=0;
2706  int add_row_shift=0;
2707  intvec *weights=NULL;
2708  intvec *ww=(intvec *)atGet(u,"isHomog",INTVEC_CMD);
2709  if (ww!=NULL)
2710  {
2711     weights=ivCopy(ww);
2712     add_row_shift = ww->min_in();
2713     (*weights) -= add_row_shift;
2714  }
2715
2716  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift,weights);
2717  //row_shift += add_row_shift;
2718  //Print("row_shift=%d, add_row_shift=%d\n",row_shift,add_row_shift);
2719  atSet(res,omStrDup("rowShift"),(void*)(long)add_row_shift,INT_CMD);
2720
2721  return FALSE;
2722}
2723BOOLEAN syBetti1(leftv res, leftv u)
2724{
2725  sleftv tmp;
2726  memset(&tmp,0,sizeof(tmp));
2727  tmp.rtyp=INT_CMD;
2728  tmp.data=(void *)1;
2729  return syBetti2(res,u,&tmp);
2730}
2731
2732/*3
2733* converts a resolution into a list of modules
2734*/
2735lists syConvRes(syStrategy syzstr,BOOLEAN toDel,int add_row_shift)
2736{
2737  resolvente fullres = syzstr->fullres;
2738  resolvente minres = syzstr->minres;
2739
2740  const int length = syzstr->length;
2741
2742  if ((fullres==NULL) && (minres==NULL))
2743  {
2744    if (syzstr->hilb_coeffs==NULL)
2745    { // La Scala
2746      fullres = syReorder(syzstr->res, length, syzstr);
2747    }
2748    else
2749    { // HRES
2750      minres = syReorder(syzstr->orderedRes, length, syzstr);
2751      syKillEmptyEntres(minres, length);
2752    }
2753  }
2754
2755  resolvente tr;
2756  int typ0=IDEAL_CMD;
2757
2758  if (minres!=NULL)
2759    tr = minres;
2760  else
2761    tr = fullres;
2762
2763  resolvente trueres=NULL; intvec ** w=NULL;
2764
2765  if (length>0)
2766  {
2767    trueres = (resolvente)omAlloc0((length)*sizeof(ideal));
2768    for (int i=(length)-1;i>=0;i--)
2769    {
2770      if (tr[i]!=NULL)
2771      {
2772        trueres[i] = idCopy(tr[i]);
2773      }
2774    }
2775    if ( id_RankFreeModule(trueres[0], currRing) > 0)
2776      typ0 = MODUL_CMD;
2777    if (syzstr->weights!=NULL)
2778    {
2779      w = (intvec**)omAlloc0(length*sizeof(intvec*));
2780      for (int i=length-1;i>=0;i--)
2781      {
2782        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
2783      }
2784    }
2785  }
2786
2787  lists li = liMakeResolv(trueres, length, syzstr->list_length,typ0,
2788                          w, add_row_shift);
2789
2790  if (w != NULL) omFreeSize(w, length*sizeof(intvec*));
2791
2792  if (toDel)
2793    syKillComputation(syzstr);
2794  else
2795  {
2796    if( fullres != NULL && syzstr->fullres == NULL )
2797      syzstr->fullres = fullres;
2798
2799    if( minres != NULL && syzstr->minres == NULL )
2800      syzstr->minres = minres;
2801  }
2802
2803  return li;
2804
2805
2806}
2807
2808/*3
2809* converts a list of modules into a resolution
2810*/
2811syStrategy syConvList(lists li,BOOLEAN toDel)
2812{
2813  int typ0;
2814  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2815
2816  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
2817  if (fr != NULL)
2818  {
2819
2820    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2821    for (int i=result->length-1;i>=0;i--)
2822    {
2823      if (fr[i]!=NULL)
2824        result->fullres[i] = idCopy(fr[i]);
2825    }
2826    result->list_length=result->length;
2827    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2828  }
2829  else
2830  {
2831    omFreeSize(result, sizeof(ssyStrategy));
2832    result = NULL;
2833  }
2834  if (toDel) li->Clean();
2835  return result;
2836}
2837
2838/*3
2839* converts a list of modules into a minimal resolution
2840*/
2841syStrategy syForceMin(lists li)
2842{
2843  int typ0;
2844  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2845
2846  resolvente fr = liFindRes(li,&(result->length),&typ0);
2847  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2848  for (int i=result->length-1;i>=0;i--)
2849  {
2850    if (fr[i]!=NULL)
2851      result->minres[i] = idCopy(fr[i]);
2852  }
2853  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2854  return result;
2855}
2856// from weight.cc
2857BOOLEAN kWeight(leftv res,leftv id)
2858{
2859  ideal F=(ideal)id->Data();
2860  intvec * iv = new intvec(rVar(currRing));
2861  polyset s;
2862  int  sl, n, i;
2863  int  *x;
2864
2865  res->data=(char *)iv;
2866  s = F->m;
2867  sl = IDELEMS(F) - 1;
2868  n = rVar(currRing);
2869  double wNsqr = (double)2.0 / (double)n;
2870  wFunctional = wFunctionalBuch;
2871  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
2872  wCall(s, sl, x, wNsqr, currRing);
2873  for (i = n; i!=0; i--)
2874    (*iv)[i-1] = x[i + n + 1];
2875  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
2876  return FALSE;
2877}
2878
2879BOOLEAN kQHWeight(leftv res,leftv v)
2880{
2881  res->data=(char *)id_QHomWeight((ideal)v->Data(), currRing);
2882  if (res->data==NULL)
2883    res->data=(char *)new intvec(rVar(currRing));
2884  return FALSE;
2885}
2886/*==============================================================*/
2887// from clapsing.cc
2888#if 0
2889BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
2890{
2891  BOOLEAN b=singclap_factorize((poly)(u->CopyD()), &v, 0);
2892  res->data=(void *)b;
2893}
2894#endif
2895
2896BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
2897{
2898#ifdef HAVE_FACTORY
2899  res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
2900                  (poly)w->CopyD(), currRing);
2901  return errorreported;
2902#else
2903  Werror("Sorry: not yet re-factored: see libpolys/polys/clapsing.cc");
2904  return FALSE;
2905#endif
2906}
2907
2908BOOLEAN jjCHARSERIES(leftv res, leftv u)
2909{
2910#if defined(HAVE_FACTORY) && defined(HAVE_LIBFAC)
2911  res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
2912  return (res->data==NULL);
2913#else
2914  Werror("Sorry: not yet re-factored: see libpolys/polys/clapsing.cc");
2915  return FALSE;
2916#endif
2917}
2918
2919// from semic.cc
2920#ifdef HAVE_SPECTRUM
2921
2922// ----------------------------------------------------------------------------
2923//  Initialize a  spectrum  deep from a  singular  lists
2924// ----------------------------------------------------------------------------
2925
2926void copy_deep( spectrum& spec, lists l )
2927{
2928    spec.mu = (int)(long)(l->m[0].Data( ));
2929    spec.pg = (int)(long)(l->m[1].Data( ));
2930    spec.n  = (int)(long)(l->m[2].Data( ));
2931
2932    spec.copy_new( spec.n );
2933
2934    intvec  *num = (intvec*)l->m[3].Data( );
2935    intvec  *den = (intvec*)l->m[4].Data( );
2936    intvec  *mul = (intvec*)l->m[5].Data( );
2937
2938    for( int i=0; i<spec.n; i++ )
2939    {
2940        spec.s[i] = (Rational)((*num)[i])/(Rational)((*den)[i]);
2941        spec.w[i] = (*mul)[i];
2942    }
2943}
2944
2945// ----------------------------------------------------------------------------
2946//  singular lists  constructor for  spectrum
2947// ----------------------------------------------------------------------------
2948
2949spectrum /*former spectrum::spectrum ( lists l )*/
2950spectrumFromList( lists l )
2951{
2952    spectrum result;
2953    copy_deep( result, l );
2954    return result;
2955}
2956
2957// ----------------------------------------------------------------------------
2958//  generate a Singular  lists  from a spectrum
2959// ----------------------------------------------------------------------------
2960
2961/* former spectrum::thelist ( void )*/
2962lists   getList( spectrum& spec )
2963{
2964    lists   L  = (lists)omAllocBin( slists_bin);
2965
2966    L->Init( 6 );
2967
2968    intvec            *num  = new intvec( spec.n );
2969    intvec            *den  = new intvec( spec.n );
2970    intvec            *mult = new intvec( spec.n );
2971
2972    for( int i=0; i<spec.n; i++ )
2973    {
2974        (*num) [i] = spec.s[i].get_num_si( );
2975        (*den) [i] = spec.s[i].get_den_si( );
2976        (*mult)[i] = spec.w[i];
2977    }
2978
2979    L->m[0].rtyp = INT_CMD;    //  milnor number
2980    L->m[1].rtyp = INT_CMD;    //  geometrical genus
2981    L->m[2].rtyp = INT_CMD;    //  # of spectrum numbers
2982    L->m[3].rtyp = INTVEC_CMD; //  numerators
2983    L->m[4].rtyp = INTVEC_CMD; //  denomiantors
2984    L->m[5].rtyp = INTVEC_CMD; //  multiplicities
2985
2986    L->m[0].data = (void*)(long)spec.mu;
2987    L->m[1].data = (void*)(long)spec.pg;
2988    L->m[2].data = (void*)(long)spec.n;
2989    L->m[3].data = (void*)num;
2990    L->m[4].data = (void*)den;
2991    L->m[5].data = (void*)mult;
2992
2993    return  L;
2994}
2995// from spectrum.cc
2996// ----------------------------------------------------------------------------
2997//  print out an error message for a spectrum list
2998// ----------------------------------------------------------------------------
2999
3000typedef enum
3001{
3002    semicOK,
3003    semicMulNegative,
3004
3005    semicListTooShort,
3006    semicListTooLong,
3007
3008    semicListFirstElementWrongType,
3009    semicListSecondElementWrongType,
3010    semicListThirdElementWrongType,
3011    semicListFourthElementWrongType,
3012    semicListFifthElementWrongType,
3013    semicListSixthElementWrongType,
3014
3015    semicListNNegative,
3016    semicListWrongNumberOfNumerators,
3017    semicListWrongNumberOfDenominators,
3018    semicListWrongNumberOfMultiplicities,
3019
3020    semicListMuNegative,
3021    semicListPgNegative,
3022    semicListNumNegative,
3023    semicListDenNegative,
3024    semicListMulNegative,
3025
3026    semicListNotSymmetric,
3027    semicListNotMonotonous,
3028
3029    semicListMilnorWrong,
3030    semicListPGWrong
3031
3032} semicState;
3033
3034void    list_error( semicState state )
3035{
3036    switch( state )
3037    {
3038        case semicListTooShort:
3039            WerrorS( "the list is too short" );
3040            break;
3041        case semicListTooLong:
3042            WerrorS( "the list is too long" );
3043            break;
3044
3045        case semicListFirstElementWrongType:
3046            WerrorS( "first element of the list should be int" );
3047            break;
3048        case semicListSecondElementWrongType:
3049            WerrorS( "second element of the list should be int" );
3050            break;
3051        case semicListThirdElementWrongType:
3052            WerrorS( "third element of the list should be int" );
3053            break;
3054        case semicListFourthElementWrongType:
3055            WerrorS( "fourth element of the list should be intvec" );
3056            break;
3057        case semicListFifthElementWrongType:
3058            WerrorS( "fifth element of the list should be intvec" );
3059            break;
3060        case semicListSixthElementWrongType:
3061            WerrorS( "sixth element of the list should be intvec" );
3062            break;
3063
3064        case semicListNNegative:
3065            WerrorS( "first element of the list should be positive" );
3066            break;
3067        case semicListWrongNumberOfNumerators:
3068            WerrorS( "wrong number of numerators" );
3069            break;
3070        case semicListWrongNumberOfDenominators:
3071            WerrorS( "wrong number of denominators" );
3072            break;
3073        case semicListWrongNumberOfMultiplicities:
3074            WerrorS( "wrong number of multiplicities" );
3075            break;
3076
3077        case semicListMuNegative:
3078            WerrorS( "the Milnor number should be positive" );
3079            break;
3080        case semicListPgNegative:
3081            WerrorS( "the geometrical genus should be nonnegative" );
3082            break;
3083        case semicListNumNegative:
3084            WerrorS( "all numerators should be positive" );
3085            break;
3086        case semicListDenNegative:
3087            WerrorS( "all denominators should be positive" );
3088            break;
3089        case semicListMulNegative:
3090            WerrorS( "all multiplicities should be positive" );
3091            break;
3092
3093        case semicListNotSymmetric:
3094            WerrorS( "it is not symmetric" );
3095            break;
3096        case semicListNotMonotonous:
3097            WerrorS( "it is not monotonous" );
3098            break;
3099
3100        case semicListMilnorWrong:
3101            WerrorS( "the Milnor number is wrong" );
3102            break;
3103        case semicListPGWrong:
3104            WerrorS( "the geometrical genus is wrong" );
3105            break;
3106
3107        default:
3108            WerrorS( "unspecific error" );
3109            break;
3110    }
3111}
3112// ----------------------------------------------------------------------------
3113//  this is the main spectrum computation function
3114// ----------------------------------------------------------------------------
3115
3116enum    spectrumState
3117{
3118    spectrumOK,
3119    spectrumZero,
3120    spectrumBadPoly,
3121    spectrumNoSingularity,
3122    spectrumNotIsolated,
3123    spectrumDegenerate,
3124    spectrumWrongRing,
3125    spectrumNoHC,
3126    spectrumUnspecErr
3127};
3128
3129// from splist.cc
3130// ----------------------------------------------------------------------------
3131//  Compute the spectrum of a  spectrumPolyList
3132// ----------------------------------------------------------------------------
3133
3134/* former spectrumPolyList::spectrum ( lists*, int) */
3135spectrumState   spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3136{
3137  spectrumPolyNode  **node = &speclist.root;
3138  spectrumPolyNode  *search;
3139
3140  poly              f,tmp;
3141  int               found,cmp;
3142
3143  Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3144                 ( fast==2 ? 2 : 1 ) );
3145
3146  Rational weight_prev( 0,1 );
3147
3148  int     mu = 0;          // the milnor number
3149  int     pg = 0;          // the geometrical genus
3150  int     n  = 0;          // number of different spectral numbers
3151  int     z  = 0;          // number of spectral number equal to smax
3152
3153  while( (*node)!=(spectrumPolyNode*)NULL &&
3154         ( fast==0 || (*node)->weight<=smax ) )
3155  {
3156        // ---------------------------------------
3157        //  determine the first normal form which
3158        //  contains the monomial  node->mon
3159        // ---------------------------------------
3160
3161    found  = FALSE;
3162    search = *node;
3163
3164    while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3165    {
3166      if( search->nf!=(poly)NULL )
3167      {
3168        f = search->nf;
3169
3170        do
3171        {
3172                    // --------------------------------
3173                    //  look for  (*node)->mon  in   f
3174                    // --------------------------------
3175
3176          cmp = pCmp( (*node)->mon,f );
3177
3178          if( cmp<0 )
3179          {
3180            f = pNext( f );
3181          }
3182          else if( cmp==0 )
3183          {
3184                        // -----------------------------
3185                        //  we have found a normal form
3186                        // -----------------------------
3187
3188            found = TRUE;
3189
3190                        //  normalize coefficient
3191
3192            number inv = nInvers( pGetCoeff( f ) );
3193            pMult_nn( search->nf,inv );
3194            nDelete( &inv );
3195
3196                        //  exchange  normal forms
3197
3198            tmp         = (*node)->nf;
3199            (*node)->nf = search->nf;
3200            search->nf  = tmp;
3201          }
3202        }
3203        while( cmp<0 && f!=(poly)NULL );
3204      }
3205      search = search->next;
3206    }
3207
3208    if( found==FALSE )
3209    {
3210            // ------------------------------------------------
3211            //  the weight of  node->mon  is a spectrum number
3212            // ------------------------------------------------
3213
3214      mu++;
3215
3216      if( (*node)->weight<=(Rational)1 )              pg++;
3217      if( (*node)->weight==smax )           z++;
3218      if( (*node)->weight>weight_prev )     n++;
3219
3220      weight_prev = (*node)->weight;
3221      node = &((*node)->next);
3222    }
3223    else
3224    {
3225            // -----------------------------------------------
3226            //  determine all other normal form which contain
3227            //  the monomial  node->mon
3228            //  replace for  node->mon  its normal form
3229            // -----------------------------------------------
3230
3231      while( search!=(spectrumPolyNode*)NULL )
3232      {
3233        if( search->nf!=(poly)NULL )
3234        {
3235          f = search->nf;
3236
3237          do
3238          {
3239                        // --------------------------------
3240                        //  look for  (*node)->mon  in   f
3241                        // --------------------------------
3242
3243            cmp = pCmp( (*node)->mon,f );
3244
3245            if( cmp<0 )
3246            {
3247              f = pNext( f );
3248            }
3249            else if( cmp==0 )
3250            {
3251              search->nf = pSub( search->nf,
3252                                 ppMult_nn( (*node)->nf,pGetCoeff( f ) ) );
3253              pNorm( search->nf );
3254            }
3255          }
3256          while( cmp<0 && f!=(poly)NULL );
3257        }
3258        search = search->next;
3259      }
3260      speclist.delete_node( node );
3261    }
3262
3263  }
3264
3265    // --------------------------------------------------------
3266    //  fast computation exploits the symmetry of the spectrum
3267    // --------------------------------------------------------
3268
3269  if( fast==2 )
3270  {
3271    mu = 2*mu - z;
3272    n  = ( z > 0 ? 2*n - 1 : 2*n );
3273  }
3274
3275    // --------------------------------------------------------
3276    //  compute the spectrum numbers with their multiplicities
3277    // --------------------------------------------------------
3278
3279  intvec            *nom  = new intvec( n );
3280  intvec            *den  = new intvec( n );
3281  intvec            *mult = new intvec( n );
3282
3283  int count         = 0;
3284  int multiplicity  = 1;
3285
3286  for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3287              ( fast==0 || search->weight<=smax );
3288       search=search->next )
3289  {
3290    if( search->next==(spectrumPolyNode*)NULL ||
3291        search->weight<search->next->weight )
3292    {
3293      (*nom) [count] = search->weight.get_num_si( );
3294      (*den) [count] = search->weight.get_den_si( );
3295      (*mult)[count] = multiplicity;
3296
3297      multiplicity=1;
3298      count++;
3299    }
3300    else
3301    {
3302      multiplicity++;
3303    }
3304  }
3305
3306    // --------------------------------------------------------
3307    //  fast computation exploits the symmetry of the spectrum
3308    // --------------------------------------------------------
3309
3310  if( fast==2 )
3311  {
3312    int n1,n2;
3313    for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3314    {
3315      (*nom) [n2] = rVar(currRing)*(*den)[n1]-(*nom)[n1];
3316      (*den) [n2] = (*den)[n1];
3317      (*mult)[n2] = (*mult)[n1];
3318    }
3319  }
3320
3321    // -----------------------------------
3322    //  test if the spectrum is symmetric
3323    // -----------------------------------
3324
3325  if( fast==0 || fast==1 )
3326  {
3327    int symmetric=TRUE;
3328
3329    for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3330    {
3331      if( (*mult)[n1]!=(*mult)[n2] ||
3332          (*den) [n1]!= (*den)[n2] ||
3333          (*nom)[n1]+(*nom)[n2]!=rVar(currRing)*(*den) [n1] )
3334      {
3335        symmetric = FALSE;
3336      }
3337    }
3338
3339    if( symmetric==FALSE )
3340    {
3341            // ---------------------------------------------
3342            //  the spectrum is not symmetric => degenerate
3343            //  principal part
3344            // ---------------------------------------------
3345
3346      *L = (lists)omAllocBin( slists_bin);
3347      (*L)->Init( 1 );
3348      (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3349      (*L)->m[0].data = (void*)(long)mu;
3350
3351      return spectrumDegenerate;
3352    }
3353  }
3354
3355  *L = (lists)omAllocBin( slists_bin);
3356
3357  (*L)->Init( 6 );
3358
3359  (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3360  (*L)->m[1].rtyp = INT_CMD;    //  geometrical genus
3361  (*L)->m[2].rtyp = INT_CMD;    //  number of spectrum values
3362  (*L)->m[3].rtyp = INTVEC_CMD; //  nominators
3363  (*L)->m[4].rtyp = INTVEC_CMD; //  denomiantors
3364  (*L)->m[5].rtyp = INTVEC_CMD; //  multiplicities
3365
3366  (*L)->m[0].data = (void*)(long)mu;
3367  (*L)->m[1].data = (void*)(long)pg;
3368  (*L)->m[2].data = (void*)(long)n;
3369  (*L)->m[3].data = (void*)nom;
3370  (*L)->m[4].data = (void*)den;
3371  (*L)->m[5].data = (void*)mult;
3372
3373  return  spectrumOK;
3374}
3375
3376spectrumState   spectrumCompute( poly h,lists *L,int fast )
3377{
3378  int i;
3379
3380  #ifdef SPECTRUM_DEBUG
3381  #ifdef SPECTRUM_PRINT
3382  #ifdef SPECTRUM_IOSTREAM
3383    cout << "spectrumCompute\n";
3384    if( fast==0 ) cout << "    no optimization" << endl;
3385    if( fast==1 ) cout << "    weight optimization" << endl;
3386    if( fast==2 ) cout << "    symmetry optimization" << endl;
3387  #else
3388    fprintf( stdout,"spectrumCompute\n" );
3389    if( fast==0 ) fprintf( stdout,"    no optimization\n" );
3390    if( fast==1 ) fprintf( stdout,"    weight optimization\n" );
3391    if( fast==2 ) fprintf( stdout,"    symmetry optimization\n" );
3392  #endif
3393  #endif
3394  #endif
3395
3396  // ----------------------
3397  //  check if  h  is zero
3398  // ----------------------
3399
3400  if( h==(poly)NULL )
3401  {
3402    return  spectrumZero;
3403  }
3404
3405  // ----------------------------------
3406  //  check if  h  has a constant term
3407  // ----------------------------------
3408
3409  if( hasConstTerm( h, currRing ) )
3410  {
3411    return  spectrumBadPoly;
3412  }
3413
3414  // --------------------------------
3415  //  check if  h  has a linear term
3416  // --------------------------------
3417
3418  if( hasLinearTerm( h, currRing ) )
3419  {
3420    *L = (lists)omAllocBin( slists_bin);
3421    (*L)->Init( 1 );
3422    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3423    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3424
3425    return  spectrumNoSingularity;
3426  }
3427
3428  // ----------------------------------
3429  //  compute the jacobi ideal of  (h)
3430  // ----------------------------------
3431
3432  ideal J = NULL;
3433  J = idInit( rVar(currRing),1 );
3434
3435  #ifdef SPECTRUM_DEBUG
3436  #ifdef SPECTRUM_PRINT
3437  #ifdef SPECTRUM_IOSTREAM
3438    cout << "\n   computing the Jacobi ideal...\n";
3439  #else
3440    fprintf( stdout,"\n   computing the Jacobi ideal...\n" );
3441  #endif
3442  #endif
3443  #endif
3444
3445  for( i=0; i<rVar(currRing); i++ )
3446  {
3447    J->m[i] = pDiff( h,i+1); //j );
3448
3449    #ifdef SPECTRUM_DEBUG
3450    #ifdef SPECTRUM_PRINT
3451    #ifdef SPECTRUM_IOSTREAM
3452      cout << "        ";
3453    #else
3454      fprintf( stdout,"        " );
3455    #endif
3456      pWrite( J->m[i] );
3457    #endif
3458    #endif
3459  }
3460
3461  // --------------------------------------------
3462  //  compute a standard basis  stdJ  of  jac(h)
3463  // --------------------------------------------
3464
3465  #ifdef SPECTRUM_DEBUG
3466  #ifdef SPECTRUM_PRINT
3467  #ifdef SPECTRUM_IOSTREAM
3468    cout << endl;
3469    cout << "    computing a standard basis..." << endl;
3470  #else
3471    fprintf( stdout,"\n" );
3472    fprintf( stdout,"    computing a standard basis...\n" );
3473  #endif
3474  #endif
3475  #endif
3476
3477  ideal stdJ = kStd(J,currQuotient,isNotHomog,NULL);
3478  idSkipZeroes( stdJ );
3479
3480  #ifdef SPECTRUM_DEBUG
3481  #ifdef SPECTRUM_PRINT
3482    for( i=0; i<IDELEMS(stdJ); i++ )
3483    {
3484      #ifdef SPECTRUM_IOSTREAM
3485        cout << "        ";
3486      #else
3487        fprintf( stdout,"        " );
3488      #endif
3489
3490      pWrite( stdJ->m[i] );
3491    }
3492  #endif
3493  #endif
3494
3495  idDelete( &J );
3496
3497  // ------------------------------------------
3498  //  check if the  h  has a singularity
3499  // ------------------------------------------
3500
3501  if( hasOne( stdJ, currRing ) )
3502  {
3503    // -------------------------------
3504    //  h is smooth in the origin
3505    //  return only the Milnor number
3506    // -------------------------------
3507
3508    *L = (lists)omAllocBin( slists_bin);
3509    (*L)->Init( 1 );
3510    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
3511    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
3512
3513    return  spectrumNoSingularity;
3514  }
3515
3516  // ------------------------------------------
3517  //  check if the singularity  h  is isolated
3518  // ------------------------------------------
3519
3520  for( i=rVar(currRing); i>0; i-- )
3521  {
3522    if( hasAxis( stdJ,i, currRing )==FALSE )
3523    {
3524      return  spectrumNotIsolated;
3525    }
3526  }
3527
3528  // ------------------------------------------
3529  //  compute the highest corner  hc  of  stdJ
3530  // ------------------------------------------
3531
3532  #ifdef SPECTRUM_DEBUG
3533  #ifdef SPECTRUM_PRINT
3534  #ifdef SPECTRUM_IOSTREAM
3535    cout << "\n    computing the highest corner...\n";
3536  #else
3537    fprintf( stdout,"\n    computing the highest corner...\n" );
3538  #endif
3539  #endif
3540  #endif
3541
3542  poly hc = (poly)NULL;
3543
3544  scComputeHC( stdJ,currQuotient, 0,hc );
3545
3546  if( hc!=(poly)NULL )
3547  {
3548    pGetCoeff(hc) = nInit(1);
3549
3550    for( i=rVar(currRing); i>0; i-- )
3551    {
3552      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3553    }
3554    pSetm( hc );
3555  }
3556  else
3557  {
3558    return  spectrumNoHC;
3559  }
3560
3561  #ifdef SPECTRUM_DEBUG
3562  #ifdef SPECTRUM_PRINT
3563  #ifdef SPECTRUM_IOSTREAM
3564    cout << "       ";
3565  #else
3566    fprintf( stdout,"       " );
3567  #endif
3568    pWrite( hc );
3569  #endif
3570  #endif
3571
3572  // ----------------------------------------
3573  //  compute the Newton polygon  nph  of  h
3574  // ----------------------------------------
3575
3576  #ifdef SPECTRUM_DEBUG
3577  #ifdef SPECTRUM_PRINT
3578  #ifdef SPECTRUM_IOSTREAM
3579    cout << "\n    computing the newton polygon...\n";
3580  #else
3581    fprintf( stdout,"\n    computing the newton polygon...\n" );
3582  #endif
3583  #endif
3584  #endif
3585
3586  newtonPolygon nph( h, currRing );
3587
3588  #ifdef SPECTRUM_DEBUG
3589  #ifdef SPECTRUM_PRINT
3590    cout << nph;
3591  #endif
3592  #endif
3593
3594  // -----------------------------------------------
3595  //  compute the weight corner  wc  of  (stdj,nph)
3596  // -----------------------------------------------
3597
3598  #ifdef SPECTRUM_DEBUG
3599  #ifdef SPECTRUM_PRINT
3600  #ifdef SPECTRUM_IOSTREAM
3601    cout << "\n    computing the weight corner...\n";
3602  #else
3603    fprintf( stdout,"\n    computing the weight corner...\n" );
3604  #endif
3605  #endif
3606  #endif
3607
3608  poly    wc = ( fast==0 ? pCopy( hc ) :
3609               ( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
3610              /* fast==2 */computeWC( nph,
3611                      ((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
3612
3613  #ifdef SPECTRUM_DEBUG
3614  #ifdef SPECTRUM_PRINT
3615  #ifdef SPECTRUM_IOSTREAM
3616    cout << "        ";
3617  #else
3618    fprintf( stdout,"        " );
3619  #endif
3620    pWrite( wc );
3621  #endif
3622  #endif
3623
3624  // -------------
3625  //  compute  NF
3626  // -------------
3627
3628  #ifdef SPECTRUM_DEBUG
3629  #ifdef SPECTRUM_PRINT
3630  #ifdef SPECTRUM_IOSTREAM
3631    cout << "\n    computing NF...\n" << endl;
3632  #else
3633    fprintf( stdout,"\n    computing NF...\n" );
3634  #endif
3635  #endif
3636  #endif
3637
3638  spectrumPolyList NF( &nph );
3639
3640  computeNF( stdJ,hc,wc,&NF, currRing );
3641
3642  #ifdef SPECTRUM_DEBUG
3643  #ifdef SPECTRUM_PRINT
3644    cout << NF;
3645  #ifdef SPECTRUM_IOSTREAM
3646    cout << endl;
3647  #else
3648    fprintf( stdout,"\n" );
3649  #endif
3650  #endif
3651  #endif
3652
3653  // ----------------------------
3654  //  compute the spectrum of  h
3655  // ----------------------------
3656//  spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
3657
3658  return spectrumStateFromList(NF, L, fast );
3659}
3660
3661// ----------------------------------------------------------------------------
3662//  this procedure is called from the interpreter
3663// ----------------------------------------------------------------------------
3664//  first  = polynomial
3665//  result = list of spectrum numbers
3666// ----------------------------------------------------------------------------
3667
3668void spectrumPrintError(spectrumState state)
3669{
3670  switch( state )
3671  {
3672    case spectrumZero:
3673      WerrorS( "polynomial is zero" );
3674      break;
3675    case spectrumBadPoly:
3676      WerrorS( "polynomial has constant term" );
3677      break;
3678    case spectrumNoSingularity:
3679      WerrorS( "not a singularity" );
3680      break;
3681    case spectrumNotIsolated:
3682      WerrorS( "the singularity is not isolated" );
3683      break;
3684    case spectrumNoHC:
3685      WerrorS( "highest corner cannot be computed" );
3686      break;
3687    case spectrumDegenerate:
3688      WerrorS( "principal part is degenerate" );
3689      break;
3690    case spectrumOK:
3691      break;
3692
3693    default:
3694      WerrorS( "unknown error occurred" );
3695      break;
3696  }
3697}
3698
3699BOOLEAN spectrumProc( leftv result,leftv first )
3700{
3701  spectrumState state = spectrumOK;
3702
3703  // -------------------
3704  //  check consistency
3705  // -------------------
3706
3707  //  check for a local ring
3708
3709  if( !ringIsLocal(currRing ) )
3710  {
3711    WerrorS( "only works for local orderings" );
3712    state = spectrumWrongRing;
3713  }
3714
3715  //  no quotient rings are allowed
3716
3717  else if( currRing->qideal != NULL )
3718  {
3719    WerrorS( "does not work in quotient rings" );
3720    state = spectrumWrongRing;
3721  }
3722  else
3723  {
3724    lists   L    = (lists)NULL;
3725    int     flag = 1; // weight corner optimization is safe
3726
3727    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3728
3729    if( state==spectrumOK )
3730    {
3731      result->rtyp = LIST_CMD;
3732      result->data = (char*)L;
3733    }
3734    else
3735    {
3736      spectrumPrintError(state);
3737    }
3738  }
3739
3740  return  (state!=spectrumOK);
3741}
3742
3743// ----------------------------------------------------------------------------
3744//  this procedure is called from the interpreter
3745// ----------------------------------------------------------------------------
3746//  first  = polynomial
3747//  result = list of spectrum numbers
3748// ----------------------------------------------------------------------------
3749
3750BOOLEAN spectrumfProc( leftv result,leftv first )
3751{
3752  spectrumState state = spectrumOK;
3753
3754  // -------------------
3755  //  check consistency
3756  // -------------------
3757
3758  //  check for a local polynomial ring
3759
3760  if( currRing->OrdSgn != -1 )
3761  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
3762  // or should we use:
3763  //if( !ringIsLocal( ) )
3764  {
3765    WerrorS( "only works for local orderings" );
3766    state = spectrumWrongRing;
3767  }
3768  else if( currRing->qideal != NULL )
3769  {
3770    WerrorS( "does not work in quotient rings" );
3771    state = spectrumWrongRing;
3772  }
3773  else
3774  {
3775    lists   L    = (lists)NULL;
3776    int     flag = 2; // symmetric optimization
3777
3778    state = spectrumCompute( (poly)first->Data( ),&L,flag );
3779
3780    if( state==spectrumOK )
3781    {
3782      result->rtyp = LIST_CMD;
3783      result->data = (char*)L;
3784    }
3785    else
3786    {
3787      spectrumPrintError(state);
3788    }
3789  }
3790
3791  return  (state!=spectrumOK);
3792}
3793
3794// ----------------------------------------------------------------------------
3795//  check if a list is a spectrum
3796//  check for:
3797//      list has 6 elements
3798//      1st element is int (mu=Milnor number)
3799//      2nd element is int (pg=geometrical genus)
3800//      3rd element is int (n =number of different spectrum numbers)
3801//      4th element is intvec (num=numerators)
3802//      5th element is intvec (den=denomiantors)
3803//      6th element is intvec (mul=multiplicities)
3804//      exactly n numerators
3805//      exactly n denominators
3806//      exactly n multiplicities
3807//      mu>0
3808//      pg>=0
3809//      n>0
3810//      num>0
3811//      den>0
3812//      mul>0
3813//      symmetriy with respect to numberofvariables/2
3814//      monotony
3815//      mu = sum of all multiplicities
3816//      pg = sum of all multiplicities where num/den<=1
3817// ----------------------------------------------------------------------------
3818
3819semicState  list_is_spectrum( lists l )
3820{
3821    // -------------------
3822    //  check list length
3823    // -------------------
3824
3825    if( l->nr < 5 )
3826    {
3827        return  semicListTooShort;
3828    }
3829    else if( l->nr > 5 )
3830    {
3831        return  semicListTooLong;
3832    }
3833
3834    // -------------
3835    //  check types
3836    // -------------
3837
3838    if( l->m[0].rtyp != INT_CMD )
3839    {
3840        return  semicListFirstElementWrongType;
3841    }
3842    else if( l->m[1].rtyp != INT_CMD )
3843    {
3844        return  semicListSecondElementWrongType;
3845    }
3846    else if( l->m[2].rtyp != INT_CMD )
3847    {
3848        return  semicListThirdElementWrongType;
3849    }
3850    else if( l->m[3].rtyp != INTVEC_CMD )
3851    {
3852        return  semicListFourthElementWrongType;
3853    }
3854    else if( l->m[4].rtyp != INTVEC_CMD )
3855    {
3856        return  semicListFifthElementWrongType;
3857    }
3858    else if( l->m[5].rtyp != INTVEC_CMD )
3859    {
3860        return  semicListSixthElementWrongType;
3861    }
3862
3863    // -------------------------
3864    //  check number of entries
3865    // -------------------------
3866
3867    int     mu = (int)(long)(l->m[0].Data( ));
3868    int     pg = (int)(long)(l->m[1].Data( ));
3869    int     n  = (int)(long)(l->m[2].Data( ));
3870
3871    if( n <= 0 )
3872    {
3873        return  semicListNNegative;
3874    }
3875
3876    intvec  *num = (intvec*)l->m[3].Data( );
3877    intvec  *den = (intvec*)l->m[4].Data( );
3878    intvec  *mul = (intvec*)l->m[5].Data( );
3879
3880    if( n != num->length( ) )
3881    {
3882        return  semicListWrongNumberOfNumerators;
3883    }
3884    else if( n != den->length( ) )
3885    {
3886        return  semicListWrongNumberOfDenominators;
3887    }
3888    else if( n != mul->length( ) )
3889    {
3890        return  semicListWrongNumberOfMultiplicities;
3891    }
3892
3893    // --------
3894    //  values
3895    // --------
3896
3897    if( mu <= 0 )
3898    {
3899        return  semicListMuNegative;
3900    }
3901    if( pg < 0 )
3902    {
3903        return  semicListPgNegative;
3904    }
3905
3906    int i;
3907
3908    for( i=0; i<n; i++ )
3909    {
3910        if( (*num)[i] <= 0 )
3911        {
3912            return  semicListNumNegative;
3913        }
3914        if( (*den)[i] <= 0 )
3915        {
3916            return  semicListDenNegative;
3917        }
3918        if( (*mul)[i] <= 0 )
3919        {
3920            return  semicListMulNegative;
3921        }
3922    }
3923
3924    // ----------------
3925    //  check symmetry
3926    // ----------------
3927
3928    int     j;
3929
3930    for( i=0, j=n-1; i<=j; i++,j-- )
3931    {
3932        if( (*num)[i] != rVar(currRing)*((*den)[i]) - (*num)[j] ||
3933            (*den)[i] != (*den)[j] ||
3934            (*mul)[i] != (*mul)[j] )
3935        {
3936            return  semicListNotSymmetric;
3937        }
3938    }
3939
3940    // ----------------
3941    //  check monotony
3942    // ----------------
3943
3944    for( i=0, j=1; i<n/2; i++,j++ )
3945    {
3946        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
3947        {
3948            return  semicListNotMonotonous;
3949        }
3950    }
3951
3952    // ---------------------
3953    //  check Milnor number
3954    // ---------------------
3955
3956    for( mu=0, i=0; i<n; i++ )
3957    {
3958        mu += (*mul)[i];
3959    }
3960
3961    if( mu != (int)(long)(l->m[0].Data( )) )
3962    {
3963        return  semicListMilnorWrong;
3964    }
3965
3966    // -------------------------
3967    //  check geometrical genus
3968    // -------------------------
3969
3970    for( pg=0, i=0; i<n; i++ )
3971    {
3972        if( (*num)[i]<=(*den)[i] )
3973        {
3974            pg += (*mul)[i];
3975        }
3976    }
3977
3978    if( pg != (int)(long)(l->m[1].Data( )) )
3979    {
3980        return  semicListPGWrong;
3981    }
3982
3983    return  semicOK;
3984}
3985
3986// ----------------------------------------------------------------------------
3987//  this procedure is called from the interpreter
3988// ----------------------------------------------------------------------------
3989//  first  = list of spectrum numbers
3990//  second = list of spectrum numbers
3991//  result = sum of the two lists
3992// ----------------------------------------------------------------------------
3993
3994BOOLEAN spaddProc( leftv result,leftv first,leftv second )
3995{
3996    semicState  state;
3997
3998    // -----------------
3999    //  check arguments
4000    // -----------------
4001
4002    lists l1 = (lists)first->Data( );
4003    lists l2 = (lists)second->Data( );
4004
4005    if( (state=list_is_spectrum( l1 )) != semicOK )
4006    {
4007        WerrorS( "first argument is not a spectrum:" );
4008        list_error( state );
4009    }
4010    else if( (state=list_is_spectrum( l2 )) != semicOK )
4011    {
4012        WerrorS( "second argument is not a spectrum:" );
4013        list_error( state );
4014    }
4015    else
4016    {
4017        spectrum s1= spectrumFromList ( l1 );
4018        spectrum s2= spectrumFromList ( l2 );
4019        spectrum sum( s1+s2 );
4020
4021        result->rtyp = LIST_CMD;
4022        result->data = (char*)(getList(sum));
4023    }
4024
4025    return  (state!=semicOK);
4026}
4027
4028// ----------------------------------------------------------------------------
4029//  this procedure is called from the interpreter
4030// ----------------------------------------------------------------------------
4031//  first  = list of spectrum numbers
4032//  second = integer
4033//  result = the multiple of the first list by the second factor
4034// ----------------------------------------------------------------------------
4035
4036BOOLEAN spmulProc( leftv result,leftv first,leftv second )
4037{
4038    semicState  state;
4039
4040    // -----------------
4041    //  check arguments
4042    // -----------------
4043
4044    lists   l = (lists)first->Data( );
4045    int     k = (int)(long)second->Data( );
4046
4047    if( (state=list_is_spectrum( l ))!=semicOK )
4048    {
4049        WerrorS( "first argument is not a spectrum" );
4050        list_error( state );
4051    }
4052    else if( k < 0 )
4053    {
4054        WerrorS( "second argument should be positive" );
4055        state = semicMulNegative;
4056    }
4057    else
4058    {
4059        spectrum s= spectrumFromList( l );
4060        spectrum product( k*s );
4061
4062        result->rtyp = LIST_CMD;
4063        result->data = (char*)getList(product);
4064    }
4065
4066    return  (state!=semicOK);
4067}
4068
4069// ----------------------------------------------------------------------------
4070//  this procedure is called from the interpreter
4071// ----------------------------------------------------------------------------
4072//  first  = list of spectrum numbers
4073//  second = list of spectrum numbers
4074//  result = semicontinuity index
4075// ----------------------------------------------------------------------------
4076
4077BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
4078{
4079  semicState  state;
4080  BOOLEAN qh=(((int)(long)w->Data())==1);
4081
4082  // -----------------
4083  //  check arguments
4084  // -----------------
4085
4086  lists l1 = (lists)u->Data( );
4087  lists l2 = (lists)v->Data( );
4088
4089  if( (state=list_is_spectrum( l1 ))!=semicOK )
4090  {
4091    WerrorS( "first argument is not a spectrum" );
4092    list_error( state );
4093  }
4094  else if( (state=list_is_spectrum( l2 ))!=semicOK )
4095  {
4096    WerrorS( "second argument is not a spectrum" );
4097    list_error( state );
4098  }
4099  else
4100  {
4101    spectrum s1= spectrumFromList( l1 );
4102    spectrum s2= spectrumFromList( l2 );
4103
4104    res->rtyp = INT_CMD;
4105    if (qh)
4106      res->data = (void*)(long)(s1.mult_spectrumh( s2 ));
4107    else
4108      res->data = (void*)(long)(s1.mult_spectrum( s2 ));
4109  }
4110
4111  // -----------------
4112  //  check status
4113  // -----------------
4114
4115  return  (state!=semicOK);
4116}
4117BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
4118{
4119  sleftv tmp;
4120  memset(&tmp,0,sizeof(tmp));
4121  tmp.rtyp=INT_CMD;
4122  /* tmp.data = (void *)0;  -- done by memset */
4123
4124  return  semicProc3(res,u,v,&tmp);
4125}
4126
4127#endif
4128
4129//from mpr_inout.cc
4130extern void nPrint(number n);
4131
4132BOOLEAN loNewtonP( leftv res, leftv arg1 )
4133{
4134  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4135  return FALSE;
4136}
4137
4138BOOLEAN loSimplex( leftv res, leftv args )
4139{
4140  if ( !(rField_is_long_R(currRing)) )
4141  {
4142    WerrorS("Ground field not implemented!");
4143    return TRUE;
4144  }
4145
4146  simplex * LP;
4147  matrix m;
4148
4149  leftv v= args;
4150  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4151    return TRUE;
4152  else
4153    m= (matrix)(v->CopyD());
4154
4155  LP = new simplex(MATROWS(m),MATCOLS(m));
4156  LP->mapFromMatrix(m);
4157
4158  v= v->next;
4159  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
4160    return TRUE;
4161  else
4162    LP->m= (int)(long)(v->Data());
4163
4164  v= v->next;
4165  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
4166    return TRUE;
4167  else
4168    LP->n= (int)(long)(v->Data());
4169
4170  v= v->next;
4171  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
4172    return TRUE;
4173  else
4174    LP->m1= (int)(long)(v->Data());
4175
4176  v= v->next;
4177  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
4178    return TRUE;
4179  else
4180    LP->m2= (int)(long)(v->Data());
4181
4182  v= v->next;
4183  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
4184    return TRUE;
4185  else
4186    LP->m3= (int)(long)(v->Data());
4187
4188#ifdef mprDEBUG_PROT
4189  Print("m (constraints) %d\n",LP->m);
4190  Print("n (columns) %d\n",LP->n);
4191  Print("m1 (<=) %d\n",LP->m1);
4192  Print("m2 (>=) %d\n",LP->m2);
4193  Print("m3 (==) %d\n",LP->m3);
4194#endif
4195
4196  LP->compute();
4197
4198  lists lres= (lists)omAlloc( sizeof(slists) );
4199  lres->Init( 6 );
4200
4201  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4202  lres->m[0].data=(void*)LP->mapToMatrix(m);
4203
4204  lres->m[1].rtyp= INT_CMD;   // found a solution?
4205  lres->m[1].data=(void*)(long)LP->icase;
4206
4207  lres->m[2].rtyp= INTVEC_CMD;
4208  lres->m[2].data=(void*)LP->posvToIV();
4209
4210  lres->m[3].rtyp= INTVEC_CMD;
4211  lres->m[3].data=(void*)LP->zrovToIV();
4212
4213  lres->m[4].rtyp= INT_CMD;
4214  lres->m[4].data=(void*)(long)LP->m;
4215
4216  lres->m[5].rtyp= INT_CMD;
4217  lres->m[5].data=(void*)(long)LP->n;
4218
4219  res->data= (void*)lres;
4220
4221  return FALSE;
4222}
4223
4224BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4225{
4226  ideal gls = (ideal)(arg1->Data());
4227  int imtype= (int)(long)arg2->Data();
4228
4229  uResultant::resMatType mtype= determineMType( imtype );
4230
4231  // check input ideal ( = polynomial system )
4232  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4233  {
4234    return TRUE;
4235  }
4236
4237  uResultant *resMat= new uResultant( gls, mtype, false );
4238  if (resMat!=NULL)
4239  {
4240    res->rtyp = MODUL_CMD;
4241    res->data= (void*)resMat->accessResMat()->getMatrix();
4242    if (!errorreported) delete resMat;
4243  }
4244  return errorreported;
4245}
4246
4247BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4248{
4249
4250  poly gls;
4251  gls= (poly)(arg1->Data());
4252  int howclean= (int)(long)arg3->Data();
4253
4254  if ( !(rField_is_R(currRing) ||
4255         rField_is_Q(currRing) ||
4256         rField_is_long_R(currRing) ||
4257         rField_is_long_C(currRing)) )
4258  {
4259    WerrorS("Ground field not implemented!");
4260    return TRUE;
4261  }
4262
4263  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4264                          rField_is_long_C(currRing)) )
4265  {
4266    unsigned long int ii = (unsigned long int)arg2->Data();
4267    setGMPFloatDigits( ii, ii );
4268  }
4269
4270  if ( gls == NULL || pIsConstant( gls ) )
4271  {
4272    WerrorS("Input polynomial is constant!");
4273    return TRUE;
4274  }
4275
4276  int ldummy;
4277  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4278  //  int deg= pDeg( gls );
4279  //  int len= pLength( gls );
4280  int i,vpos=0;
4281  poly piter;
4282  lists elist;
4283  lists rlist;
4284
4285  elist= (lists)omAlloc( sizeof(slists) );
4286  elist->Init( 0 );
4287
4288  if ( rVar(currRing) > 1 )
4289  {
4290    piter= gls;
4291    for ( i= 1; i <= rVar(currRing); i++ )
4292      if ( pGetExp( piter, i ) )
4293      {
4294        vpos= i;
4295        break;
4296      }
4297    while ( piter )
4298    {
4299      for ( i= 1; i <= rVar(currRing); i++ )
4300        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4301        {
4302          WerrorS("The input polynomial must be univariate!");
4303          return TRUE;
4304        }
4305      pIter( piter );
4306    }
4307  }
4308
4309  rootContainer * roots= new rootContainer();
4310  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4311  piter= gls;
4312  for ( i= deg; i >= 0; i-- )
4313  {
4314    //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4315    if ( piter && pTotaldegree(piter) == i )
4316    {
4317      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4318      //nPrint( pcoeffs[i] );PrintS("  ");
4319      pIter( piter );
4320    }
4321    else
4322    {
4323      pcoeffs[i]= nInit(0);
4324    }
4325  }
4326
4327#ifdef mprDEBUG_PROT
4328  for (i=deg; i >= 0; i--)
4329  {
4330    nPrint( pcoeffs[i] );PrintS("  ");
4331  }
4332  PrintLn();
4333#endif
4334
4335  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4336  roots->solver( howclean );
4337
4338  int elem= roots->getAnzRoots();
4339  char *dummy;
4340  int j;
4341
4342  rlist= (lists)omAlloc( sizeof(slists) );
4343  rlist->Init( elem );
4344
4345  if (rField_is_long_C(currRing))
4346  {
4347    for ( j= 0; j < elem; j++ )
4348    {
4349      rlist->m[j].rtyp=NUMBER_CMD;
4350      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4351      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4352    }
4353  }
4354  else
4355  {
4356    for ( j= 0; j < elem; j++ )
4357    {
4358      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4359      rlist->m[j].rtyp=STRING_CMD;
4360      rlist->m[j].data=(void *)dummy;
4361    }
4362  }
4363
4364  elist->Clean();
4365  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4366
4367  // this is (via fillContainer) the same data as in root
4368  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4369  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4370
4371  delete roots;
4372
4373  res->rtyp= LIST_CMD;
4374  res->data= (void*)rlist;
4375
4376  return FALSE;
4377}
4378
4379BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4380{
4381  int i;
4382  ideal p,w;
4383  p= (ideal)arg1->Data();
4384  w= (ideal)arg2->Data();
4385
4386  // w[0] = f(p^0)
4387  // w[1] = f(p^1)
4388  // ...
4389  // p can be a vector of numbers (multivariate polynom)
4390  //   or one number (univariate polynom)
4391  // tdg = deg(f)
4392
4393  int n= IDELEMS( p );
4394  int m= IDELEMS( w );
4395  int tdg= (int)(long)arg3->Data();
4396
4397  res->data= (void*)NULL;
4398
4399  // check the input
4400  if ( tdg < 1 )
4401  {
4402    WerrorS("Last input parameter must be > 0!");
4403    return TRUE;
4404  }
4405  if ( n != rVar(currRing) )
4406  {
4407    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4408    return TRUE;
4409  }
4410  if ( m != (int)pow((double)tdg+1,(double)n) )
4411  {
4412    Werror("Size of second input ideal must be equal to %d!",
4413      (int)pow((double)tdg+1,(double)n));
4414    return TRUE;
4415  }
4416  if ( !(rField_is_Q(currRing) /* ||
4417         rField_is_R() || rField_is_long_R() ||
4418         rField_is_long_C()*/ ) )
4419         {
4420    WerrorS("Ground field not implemented!");
4421    return TRUE;
4422  }
4423
4424  number tmp;
4425  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4426  for ( i= 0; i < n; i++ )
4427  {
4428    pevpoint[i]=nInit(0);
4429    if (  (p->m)[i] )
4430    {
4431      tmp = pGetCoeff( (p->m)[i] );
4432      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4433      {
4434        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4435        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4436        return TRUE;
4437      }
4438    } else tmp= NULL;
4439    if ( !nIsZero(tmp) )
4440    {
4441      if ( !pIsConstant((p->m)[i]))
4442      {
4443        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4444        WerrorS("Elements of first input ideal must be numbers!");
4445        return TRUE;
4446      }
4447      pevpoint[i]= nCopy( tmp );
4448    }
4449  }
4450
4451  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4452  for ( i= 0; i < m; i++ )
4453  {
4454    wresults[i]= nInit(0);
4455    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4456    {
4457      if ( !pIsConstant((w->m)[i]))
4458      {
4459        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4460        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4461        WerrorS("Elements of second input ideal must be numbers!");
4462        return TRUE;
4463      }
4464      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4465    }
4466  }
4467
4468  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4469  number *ncpoly= vm.interpolateDense( wresults );
4470  // do not free ncpoly[]!!
4471  poly rpoly= vm.numvec2poly( ncpoly );
4472
4473  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4474  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4475
4476  res->data= (void*)rpoly;
4477  return FALSE;
4478}
4479
4480BOOLEAN nuUResSolve( leftv res, leftv args )
4481{
4482  leftv v= args;
4483
4484  ideal gls;
4485  int imtype;
4486  int howclean;
4487
4488  // get ideal
4489  if ( v->Typ() != IDEAL_CMD )
4490    return TRUE;
4491  else gls= (ideal)(v->Data());
4492  v= v->next;
4493
4494  // get resultant matrix type to use (0,1)
4495  if ( v->Typ() != INT_CMD )
4496    return TRUE;
4497  else imtype= (int)(long)v->Data();
4498  v= v->next;
4499
4500  if (imtype==0)
4501  {
4502    ideal test_id=idInit(1,1);
4503    int j;
4504    for(j=IDELEMS(gls)-1;j>=0;j--)
4505    {
4506      if (gls->m[j]!=NULL)
4507      {
4508        test_id->m[0]=gls->m[j];
4509        intvec *dummy_w=id_QHomWeight(test_id, currRing);
4510        if (dummy_w!=NULL)
4511        {
4512          WerrorS("Newton polytope not of expected dimension");
4513          delete dummy_w;
4514          return TRUE;
4515        }
4516      }
4517    }
4518  }
4519
4520  // get and set precision in digits ( > 0 )
4521  if ( v->Typ() != INT_CMD )
4522    return TRUE;
4523  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4524                          rField_is_long_C(currRing)) )
4525  {
4526    unsigned long int ii=(unsigned long int)v->Data();
4527    setGMPFloatDigits( ii, ii );
4528  }
4529  v= v->next;
4530
4531  // get interpolation steps (0,1,2)
4532  if ( v->Typ() != INT_CMD )
4533    return TRUE;
4534  else howclean= (int)(long)v->Data();
4535
4536  uResultant::resMatType mtype= determineMType( imtype );
4537  int i,count;
4538  lists listofroots= NULL;
4539  number smv= NULL;
4540  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4541
4542  //emptylist= (lists)omAlloc( sizeof(slists) );
4543  //emptylist->Init( 0 );
4544
4545  //res->rtyp = LIST_CMD;
4546  //res->data= (void *)emptylist;
4547
4548  // check input ideal ( = polynomial system )
4549  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4550  {
4551    return TRUE;
4552  }
4553
4554  uResultant * ures;
4555  rootContainer ** iproots;
4556  rootContainer ** muiproots;
4557  rootArranger * arranger;
4558
4559  // main task 1: setup of resultant matrix
4560  ures= new uResultant( gls, mtype );
4561  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4562  {
4563    WerrorS("Error occurred during matrix setup!");
4564    return TRUE;
4565  }
4566
4567  // if dense resultant, check if minor nonsingular
4568  if ( mtype == uResultant::denseResMat )
4569  {
4570    smv= ures->accessResMat()->getSubDet();
4571#ifdef mprDEBUG_PROT
4572    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4573#endif
4574    if ( nIsZero(smv) )
4575    {
4576      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4577      return TRUE;
4578    }
4579  }
4580
4581  // main task 2: Interpolate specialized resultant polynomials
4582  if ( interpolate_det )
4583    iproots= ures->interpolateDenseSP( false, smv );
4584  else
4585    iproots= ures->specializeInU( false, smv );
4586
4587  // main task 3: Interpolate specialized resultant polynomials
4588  if ( interpolate_det )
4589    muiproots= ures->interpolateDenseSP( true, smv );
4590  else
4591    muiproots= ures->specializeInU( true, smv );
4592
4593#ifdef mprDEBUG_PROT
4594  int c= iproots[0]->getAnzElems();
4595  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4596  c= muiproots[0]->getAnzElems();
4597  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4598#endif
4599
4600  // main task 4: Compute roots of specialized polys and match them up
4601  arranger= new rootArranger( iproots, muiproots, howclean );
4602  arranger->solve_all();
4603
4604  // get list of roots
4605  if ( arranger->success() )
4606  {
4607    arranger->arrange();
4608    listofroots= listOfRoots(arranger, gmp_output_digits );
4609  }
4610  else
4611  {
4612    WerrorS("Solver was unable to find any roots!");
4613    return TRUE;
4614  }
4615
4616  // free everything
4617  count= iproots[0]->getAnzElems();
4618  for (i=0; i < count; i++) delete iproots[i];
4619  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4620  count= muiproots[0]->getAnzElems();
4621  for (i=0; i < count; i++) delete muiproots[i];
4622  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4623
4624  delete ures;
4625  delete arranger;
4626  nDelete( &smv );
4627
4628  res->data= (void *)listofroots;
4629
4630  //emptylist->Clean();
4631  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4632
4633  return FALSE;
4634}
4635
4636// from mpr_numeric.cc
4637lists listOfRoots( rootArranger* self, const unsigned int oprec )
4638{
4639  int i,j;
4640  int count= self->roots[0]->getAnzRoots(); // number of roots
4641  int elem= self->roots[0]->getAnzElems();  // number of koordinates per root
4642
4643  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
4644
4645  if ( self->found_roots )
4646  {
4647    listofroots->Init( count );
4648
4649    for (i=0; i < count; i++)
4650    {
4651      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
4652      onepoint->Init(elem);
4653      for ( j= 0; j < elem; j++ )
4654      {
4655        if ( !rField_is_long_C(currRing) )
4656        {
4657          onepoint->m[j].rtyp=STRING_CMD;
4658          onepoint->m[j].data=(void *)complexToStr((*self->roots[j])[i],oprec, currRing->cf);
4659        }
4660        else
4661        {
4662          onepoint->m[j].rtyp=NUMBER_CMD;
4663          onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
4664        }
4665        onepoint->m[j].next= NULL;
4666        onepoint->m[j].name= NULL;
4667      }
4668      listofroots->m[i].rtyp=LIST_CMD;
4669      listofroots->m[i].data=(void *)onepoint;
4670      listofroots->m[j].next= NULL;
4671      listofroots->m[j].name= NULL;
4672    }
4673
4674  }
4675  else
4676  {
4677    listofroots->Init( 0 );
4678  }
4679
4680  return listofroots;
4681}
4682
4683// from ring.cc
4684void rSetHdl(idhdl h)
4685{
4686  ring rg = NULL;
4687  if (h!=NULL)
4688  {
4689//   Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
4690    rg = IDRING(h);
4691    if (rg==NULL) return; //id <>NULL, ring==NULL
4692    omCheckAddrSize((ADDRESS)h,sizeof(idrec));
4693    if (IDID(h))  // OB: ????
4694      omCheckAddr((ADDRESS)IDID(h));
4695    rTest(rg);
4696  }
4697
4698  // clean up history
4699  if (sLastPrinted.RingDependend())
4700  {
4701    sLastPrinted.CleanUp();
4702    memset(&sLastPrinted,0,sizeof(sleftv));
4703  }
4704
4705  // test for valid "currRing":
4706  if ((rg!=NULL) && (rg->idroot==NULL))
4707  {
4708    ring old=rg;
4709    rg=rAssure_HasComp(rg);
4710    if (old!=rg)
4711    {
4712      rKill(old);
4713      IDRING(h)=rg;
4714    }
4715  }
4716   /*------------ change the global ring -----------------------*/
4717  rChangeCurrRing(rg);
4718  currRingHdl = h;
4719}
4720
4721BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
4722{
4723  int last = 0, o=0, n = 1, i=0, typ = 1, j;
4724  sleftv *sl = ord;
4725
4726  // determine nBlocks
4727  while (sl!=NULL)
4728  {
4729    intvec *iv = (intvec *)(sl->data);
4730    if (((*iv)[1]==ringorder_c)||((*iv)[1]==ringorder_C))
4731      i++;
4732    else if ((*iv)[1]==ringorder_L)
4733    {
4734      R->bitmask=(*iv)[2];
4735      n--;
4736    }
4737    else if (((*iv)[1]!=ringorder_a)
4738    && ((*iv)[1]!=ringorder_a64)
4739    && ((*iv)[1]!=ringorder_am))
4740      o++;
4741    n++;
4742    sl=sl->next;
4743  }
4744  // check whether at least one real ordering
4745  if (o==0)
4746  {
4747    WerrorS("invalid combination of orderings");
4748    return TRUE;
4749  }
4750  // if no c/C ordering is given, increment n
4751  if (i==0) n++;
4752  else if (i != 1)
4753  {
4754    // throw error if more than one is given
4755    WerrorS("more than one ordering c/C specified");
4756    return TRUE;
4757  }
4758
4759  // initialize fields of R
4760  R->order=(int *)omAlloc0(n*sizeof(int));
4761  R->block0=(int *)omAlloc0(n*sizeof(int));
4762  R->block1=(int *)omAlloc0(n*sizeof(int));
4763  R->wvhdl=(int**)omAlloc0(n*sizeof(int_ptr));
4764
4765  int *weights=(int*)omAlloc0((R->N+1)*sizeof(int));
4766
4767  // init order, so that rBlocks works correctly
4768  for (j=0; j < n-1; j++)
4769    R->order[j] = (int) ringorder_unspec;
4770  // set last _C order, if no c/C order was given
4771  if (i == 0) R->order[n-2] = ringorder_C;
4772
4773  /* init orders */
4774  sl=ord;
4775  n=-1;
4776  while (sl!=NULL)
4777  {
4778    intvec *iv;
4779    iv = (intvec *)(sl->data);
4780    if ((*iv)[1]!=ringorder_L)
4781    {
4782      n++;
4783
4784      /* the format of an ordering:
4785       *  iv[0]: factor
4786       *  iv[1]: ordering
4787       *  iv[2..end]: weights
4788       */
4789      R->order[n] = (*iv)[1];
4790      typ=1;
4791      switch ((*iv)[1])
4792      {
4793          case ringorder_ws:
4794          case ringorder_Ws:
4795            typ=-1;
4796          case ringorder_wp:
4797          case ringorder_Wp:
4798            R->wvhdl[n]=(int*)omAlloc((iv->length()-1)*sizeof(int));
4799            R->block0[n] = last+1;
4800            for (i=2; i<iv->length(); i++)
4801            {
4802              R->wvhdl[n][i-2] = (*iv)[i];
4803              last++;
4804              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4805            }
4806            R->block1[n] = last;
4807            break;
4808          case ringorder_ls:
4809          case ringorder_ds:
4810          case ringorder_Ds:
4811          case ringorder_rs:
4812            typ=-1;
4813          case ringorder_lp:
4814          case ringorder_dp:
4815          case ringorder_Dp:
4816          case ringorder_rp:
4817            R->block0[n] = last+1;
4818            if (iv->length() == 3) last+=(*iv)[2];
4819            else last += (*iv)[0];
4820            R->block1[n] = last;
4821            //if ((R->block0[n]>R->block1[n])
4822            //|| (R->block1[n]>rVar(R)))
4823            //{
4824            //  R->block1[n]=rVar(R);
4825            //  //WerrorS("ordering larger than number of variables");
4826            //  break;
4827            //}
4828            if (rCheckIV(iv)) return TRUE;
4829            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4830            {
4831              if (weights[i]==0) weights[i]=typ;
4832            }
4833            break;
4834
4835          case ringorder_s: // no 'rank' params!
4836          {
4837
4838            if(iv->length() > 3)
4839              return TRUE;
4840
4841            if(iv->length() == 3)
4842            {
4843              const int s = (*iv)[2];
4844              R->block0[n] = s;
4845              R->block1[n] = s;
4846            }
4847            break;
4848          }
4849          case ringorder_IS:
4850          {
4851            if(iv->length() != 3) return TRUE;
4852
4853            const int s = (*iv)[2];
4854
4855            if( 1 < s || s < -1 ) return TRUE;
4856
4857            R->block0[n] = s;
4858            R->block1[n] = s;
4859            break;
4860          }
4861          case ringorder_S:
4862          case ringorder_c:
4863          case ringorder_C:
4864          {
4865            if (rCheckIV(iv)) return TRUE;
4866            break;
4867          }
4868          case ringorder_aa:
4869          case ringorder_a:
4870          {
4871            R->block0[n] = last+1;
4872            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4873            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int));
4874            for (i=2; i<iv->length(); i++)
4875            {
4876              R->wvhdl[n][i-2]=(*iv)[i];
4877              last++;
4878              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4879            }
4880            last=R->block0[n]-1;
4881            break;
4882          }
4883          case ringorder_am:
4884          {
4885            R->block0[n] = last+1;
4886            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4887            R->wvhdl[n] = (int*)omAlloc(iv->length()*sizeof(int));
4888            if (R->block1[n]- R->block0[n]+2>=iv->length())
4889               WarnS("missing module weights");
4890            for (i=2; i<=(R->block1[n]-R->block0[n]+2); i++)
4891            {
4892              R->wvhdl[n][i-2]=(*iv)[i];
4893              last++;
4894              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4895            }
4896            R->wvhdl[n][i-2]=iv->length() -3 -(R->block1[n]- R->block0[n]);
4897            for (; i<iv->length(); i++)
4898            {
4899              R->wvhdl[n][i-1]=(*iv)[i];
4900            }
4901            last=R->block0[n]-1;
4902            break;
4903          }
4904          case ringorder_a64:
4905          {
4906            R->block0[n] = last+1;
4907            R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4908            R->wvhdl[n] = (int*)omAlloc((iv->length()-1)*sizeof(int64));
4909            int64 *w=(int64 *)R->wvhdl[n];
4910            for (i=2; i<iv->length(); i++)
4911            {
4912              w[i-2]=(*iv)[i];
4913              last++;
4914              if (weights[last]==0) weights[last]=(*iv)[i]*typ;
4915            }
4916            last=R->block0[n]-1;
4917            break;
4918          }
4919          case ringorder_M:
4920          {
4921            int Mtyp=rTypeOfMatrixOrder(iv);
4922            if (Mtyp==0) return TRUE;
4923            if (Mtyp==-1) typ = -1;
4924
4925            R->wvhdl[n] =( int *)omAlloc((iv->length()-1)*sizeof(int));
4926            for (i=2; i<iv->length();i++)
4927              R->wvhdl[n][i-2]=(*iv)[i];
4928
4929            R->block0[n] = last+1;
4930            last += (int)sqrt((double)(iv->length()-2));
4931            R->block1[n] = last;
4932            for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4933            {
4934              if (weights[i]==0) weights[i]=typ;
4935            }
4936            break;
4937          }
4938
4939          case ringorder_no:
4940            R->order[n] = ringorder_unspec;
4941            return TRUE;
4942
4943          default:
4944            Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
4945            R->order[n] = ringorder_unspec;
4946            return TRUE;
4947      }
4948    }
4949    sl=sl->next;
4950  }
4951
4952  // check for complete coverage
4953  while ( n >= 0 && (
4954          (R->order[n]==ringorder_c)
4955      ||  (R->order[n]==ringorder_C)
4956      ||  (R->order[n]==ringorder_s)
4957      ||  (R->order[n]==ringorder_S)
4958      ||  (R->order[n]==ringorder_IS)
4959                    )) n--;
4960
4961  assume( n >= 0 );
4962
4963  if (R->block1[n] != R->N)
4964  {
4965    if (((R->order[n]==ringorder_dp) ||
4966         (R->order[n]==ringorder_ds) ||
4967         (R->order[n]==ringorder_Dp) ||
4968         (R->order[n]==ringorder_Ds) ||
4969         (R->order[n]==ringorder_rp) ||
4970         (R->order[n]==ringorder_rs) ||
4971         (R->order[n]==ringorder_lp) ||
4972         (R->order[n]==ringorder_ls))
4973        &&
4974        R->block0[n] <= R->N)
4975    {
4976      R->block1[n] = R->N;
4977    }
4978    else
4979    {
4980      Werror("mismatch of number of vars (%d) and ordering (%d vars)",
4981             R->N,R->block1[n]);
4982      return TRUE;
4983    }
4984  }
4985  // find OrdSgn:
4986  R->OrdSgn = 1;
4987  for(i=1;i<=R->N;i++)
4988  { if (weights[i]<0) { R->OrdSgn=-1;break; }}
4989  omFree(weights);
4990  return FALSE;
4991}
4992
4993BOOLEAN rSleftvList2StringArray(sleftv* sl, char** p)
4994{
4995
4996  while(sl!=NULL)
4997  {
4998    if (sl->Name() == sNoName)
4999    {
5000      if (sl->Typ()==POLY_CMD)
5001      {
5002        sleftv s_sl;
5003        iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
5004        if (s_sl.Name() != sNoName)
5005          *p = omStrDup(s_sl.Name());
5006        else
5007          *p = NULL;
5008        sl->next = s_sl.next;
5009        s_sl.next = NULL;
5010        s_sl.CleanUp();
5011        if (*p == NULL) return TRUE;
5012      }
5013      else
5014        return TRUE;
5015    }
5016    else
5017      *p = omStrDup(sl->Name());
5018    p++;
5019    sl=sl->next;
5020  }
5021  return FALSE;
5022}
5023
5024const short MAX_SHORT = 32767; // (1 << (sizeof(short)*8)) - 1;
5025
5026////////////////////
5027//
5028// rInit itself:
5029//
5030// INPUT:  s: name, pn: ch & parameter (names), rv: variable (names)
5031//         ord: ordering
5032// RETURN: currRingHdl on success
5033//         NULL        on error
5034// NOTE:   * makes new ring to current ring, on success
5035//         * considers input sleftv's as read-only
5036//idhdl rInit(char *s, sleftv* pn, sleftv* rv, sleftv* ord)
5037ring rInit(sleftv* pn, sleftv* rv, sleftv* ord)
5038{
5039#ifdef HAVE_RINGS
5040  //unsigned int ringtype = 0;
5041  int_number modBase = NULL;
5042  unsigned int modExponent = 1;
5043#endif
5044  int float_len=0;
5045  int float_len2=0;
5046  ring R = NULL;
5047  //BOOLEAN ffChar=FALSE;
5048
5049  /* ch -------------------------------------------------------*/
5050  // get ch of ground field
5051
5052  // allocated ring
5053  R = (ring) omAlloc0Bin(sip_sring_bin);
5054
5055  coeffs cf = NULL;
5056
5057  assume( pn != NULL );
5058  const int P = pn->listLength();
5059
5060  if (pn->Typ()==INT_CMD)
5061  {
5062    int ch = (int)(long)pn->Data();
5063
5064    /* parameter? -------------------------------------------------------*/
5065    pn = pn->next;
5066
5067    if (pn == NULL) // no params!?
5068    {
5069      if (ch!=0)
5070      {
5071        int ch2=IsPrime(ch);
5072        if ((ch<2)||(ch!=ch2))
5073        {
5074          Warn("%d is invalid as characteristic of the ground field. 32003 is used.", ch);
5075          ch=32003;
5076        }
5077        cf = nInitChar(n_Zp, (void*)(long)ch);
5078      }
5079      else
5080        cf = nInitChar(n_Q, (void*)(long)ch);
5081    }
5082    else
5083    {
5084      const int pars = pn->listLength();
5085
5086      assume( pars > 0 );
5087
5088      // predefined finite field: (p^k, a)
5089      if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
5090      {
5091        GFInfo param;
5092
5093        param.GFChar = ch;
5094        param.GFDegree = 1;
5095        param.GFPar_name = pn->name;
5096
5097        cf = nInitChar(n_GF, &param);
5098      }
5099      else // (0/p, a, b, ..., z)
5100      {
5101        if ((ch!=0) && (ch!=IsPrime(ch))) 
5102        {
5103          WerrorS("too many parameters");
5104          goto rInitError;
5105        }
5106
5107        char ** names = (char**)omAlloc0(pars * sizeof(char_ptr));
5108
5109        if (rSleftvList2StringArray(pn, names))
5110        {
5111          WerrorS("parameter expected");
5112          goto rInitError;
5113        }
5114
5115        TransExtInfo extParam;
5116
5117        extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5118        for(int i=pars-1; i>=0;i--)
5119        {
5120          omFree(names[i]);
5121        }
5122        omFree(names);
5123
5124        cf = nInitChar(n_transExt, &extParam);
5125      }
5126    }
5127
5128//    if (cf==NULL) goto rInitError;
5129    assume( cf != NULL );
5130  }
5131  else if ((pn->name != NULL)
5132  && ((strcmp(pn->name,"real")==0) || (strcmp(pn->name,"complex")==0)))
5133  {
5134    BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5135    if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5136    {
5137      float_len=(int)(long)pn->next->Data();
5138      float_len2=float_len;
5139      pn=pn->next;
5140      if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5141      {
5142        float_len2=(int)(long)pn->next->Data();
5143        pn=pn->next;
5144      }
5145    }
5146
5147    if (!complex_flag)
5148      complex_flag= pn->next != NULL;
5149    if( !complex_flag && (float_len2 <= (short)SHORT_REAL_LENGTH))
5150       cf=nInitChar(n_R, NULL);
5151    else // longR or longC?
5152    {
5153       LongComplexInfo param;
5154
5155       param.float_len = si_min (float_len, 32767);
5156       param.float_len2 = si_min (float_len2, 32767);
5157
5158       // set the parameter name
5159       if (complex_flag)
5160       {
5161         if (param.float_len < SHORT_REAL_LENGTH)
5162         {
5163           param.float_len= SHORT_REAL_LENGTH;
5164           param.float_len2= SHORT_REAL_LENGTH;
5165         }
5166         if (pn->next == NULL)
5167           param.par_name=(const char*)"i"; //default to i
5168         else
5169           param.par_name = (const char*)pn->next->name;
5170       }
5171
5172       cf = nInitChar(complex_flag ? n_long_C: n_long_R, (void*)&param);
5173    }
5174    assume( cf != NULL );
5175  }
5176#ifdef HAVE_RINGS
5177  else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5178  {
5179    modBase = (int_number) omAlloc(sizeof(mpz_t));
5180    mpz_init_set_si(modBase, 0);
5181    if (pn->next!=NULL)
5182    {
5183      if (pn->next->Typ()==INT_CMD)
5184      {
5185        mpz_set_ui(modBase, (int)(long) pn->next->Data());
5186        pn=pn->next;
5187        if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5188        {
5189          modExponent = (long) pn->next->Data();
5190          pn=pn->next;
5191        }
5192        while ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5193        {
5194          mpz_mul_ui(modBase, modBase, (int)(long) pn->next->Data());
5195          pn=pn->next;
5196        }
5197      }
5198      else if (pn->next->Typ()==BIGINT_CMD)
5199      {
5200        number p=(number)pn->next->CopyD();
5201        nlGMP(p,(number)modBase,coeffs_BIGINT);
5202        nlDelete(&p,coeffs_BIGINT);
5203      }
5204    }
5205    else
5206      cf=nInitChar(n_Z,NULL);
5207
5208    if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_cmp_ui(modBase, 0) < 0))
5209    {
5210      Werror("Wrong ground ring specification (module is 1)");
5211      goto rInitError;
5212    }
5213    if (modExponent < 1)
5214    {
5215      Werror("Wrong ground ring specification (exponent smaller than 1");
5216      goto rInitError;
5217    }
5218    // module is 0 ---> integers ringtype = 4;
5219    // we have an exponent
5220    if (modExponent > 1 && cf == NULL)
5221    {
5222      if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(NATNUMBER)))
5223      {
5224        /* this branch should be active for modExponent = 2..32 resp. 2..64,
5225           depending on the size of a long on the respective platform */
5226        //ringtype = 1;       // Use Z/2^ch
5227        cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5228        omFreeSize (modBase, sizeof (mpz_t));
5229      }
5230      else
5231      {
5232        if (mpz_cmp_ui(modBase,0)==0)
5233        {
5234          WerrorS("modulus must not be 0");
5235          goto rInitError;
5236        }
5237        //ringtype = 3;
5238        ZnmInfo info;
5239        info.base= modBase;
5240        info.exp= modExponent;
5241        cf=nInitChar(n_Znm,(void*) &info); //exponent is missing
5242      }
5243    }
5244    // just a module m > 1
5245    else if (cf == NULL)
5246    {
5247      if (mpz_cmp_ui(modBase,0)==0)
5248      {
5249        WerrorS("modulus must not be 0");
5250        goto rInitError;
5251      }
5252      //ringtype = 2;
5253      ZnmInfo info;
5254      info.base= modBase;
5255      info.exp= modExponent;
5256      cf=nInitChar(n_Zn,(void*) &info);
5257    }
5258    assume( cf != NULL );
5259  }
5260#endif
5261  // ring NEW = OLD, (), (); where OLD is a polynomial ring...
5262  else if ((pn->Typ()==RING_CMD) && (P == 1))
5263  {
5264    TransExtInfo extParam;
5265    extParam.r = (ring)pn->Data();
5266    cf = nInitChar(n_transExt, &extParam);
5267  }
5268  else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5269  {
5270    AlgExtInfo extParam;
5271    extParam.r = (ring)pn->Data();
5272
5273    cf = nInitChar(n_algExt, &extParam);   // Q[a]/<minideal>
5274  }
5275  else
5276  {
5277    Werror("Wrong or unknown ground field specification");
5278#ifndef NDEBUG
5279    sleftv* p = pn;
5280    while (p != NULL)
5281    {
5282      Print( "pn[%p]: type: %d [%s]: %p, name: %s", (void*)p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5283      PrintLn();
5284      p = p->next;
5285    }
5286#endif
5287    goto rInitError;
5288  }
5289//  pn=pn->next;
5290
5291  /*every entry in the new ring is initialized to 0*/
5292
5293  /* characteristic -----------------------------------------------*/
5294  /* input: 0 ch=0 : Q     parameter=NULL    ffChar=FALSE   float_len
5295   *         0    1 : Q(a,...)        *names         FALSE
5296   *         0   -1 : R               NULL           FALSE  0
5297   *         0   -1 : R               NULL           FALSE  prec. >6
5298   *         0   -1 : C               *names         FALSE  prec. 0..?
5299   *         p    p : Fp              NULL           FALSE
5300   *         p   -p : Fp(a)           *names         FALSE
5301   *         q    q : GF(q=p^n)       *names         TRUE
5302  */
5303  if (cf==NULL)
5304  {
5305    Werror("Invalid ground field specification");
5306    goto rInitError;
5307//    const int ch=32003;
5308//    cf=nInitChar(n_Zp, (void*)(long)ch);
5309  }
5310
5311  assume( R != NULL );
5312
5313  R->cf = cf;
5314
5315  /* names and number of variables-------------------------------------*/
5316  {
5317    int l=rv->listLength();
5318
5319    if (l>MAX_SHORT)
5320    {
5321      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5322       goto rInitError;
5323    }
5324    R->N = l; /*rv->listLength();*/
5325  }
5326  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5327  if (rSleftvList2StringArray(rv, R->names))
5328  {
5329    WerrorS("name of ring variable expected");
5330    goto rInitError;
5331  }
5332
5333  /* check names and parameters for conflicts ------------------------- */
5334  rRenameVars(R); // conflicting variables will be renamed
5335  /* ordering -------------------------------------------------------------*/
5336  if (rSleftvOrdering2Ordering(ord, R))
5337    goto rInitError;
5338
5339  // Complete the initialization
5340  if (rComplete(R,1))
5341    goto rInitError;
5342
5343/*#ifdef HAVE_RINGS
5344// currently, coefficients which are ring elements require a global ordering:
5345  if (rField_is_Ring(R) && (R->OrdSgn==-1))
5346  {
5347    WerrorS("global ordering required for these coefficients");
5348    goto rInitError;
5349  }
5350#endif*/
5351
5352  rTest(R);
5353
5354  // try to enter the ring into the name list
5355  // need to clean up sleftv here, before this ring can be set to
5356  // new currRing or currRing can be killed beacuse new ring has
5357  // same name
5358  if (pn != NULL) pn->CleanUp();
5359  if (rv != NULL) rv->CleanUp();
5360  if (ord != NULL) ord->CleanUp();
5361  //if ((tmp = enterid(s, myynest, RING_CMD, &IDROOT))==NULL)
5362  //  goto rInitError;
5363
5364  //memcpy(IDRING(tmp),R,sizeof(*R));
5365  // set current ring
5366  //omFreeBin(R,  ip_sring_bin);
5367  //return tmp;
5368  return R;
5369
5370  // error case:
5371  rInitError:
5372  if  ((R != NULL)&&(R->cf!=NULL)) rDelete(R);
5373  if (pn != NULL) pn->CleanUp();
5374  if (rv != NULL) rv->CleanUp();
5375  if (ord != NULL) ord->CleanUp();
5376  return NULL;
5377}
5378
5379ring rSubring(ring org_ring, sleftv* rv)
5380{
5381  ring R = rCopy0(org_ring);
5382  int *perm=(int *)omAlloc0((org_ring->N+1)*sizeof(int));
5383  int n = rBlocks(org_ring), i=0, j;
5384
5385  /* names and number of variables-------------------------------------*/
5386  {
5387    int l=rv->listLength();
5388    if (l>MAX_SHORT)
5389    {
5390      Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5391       goto rInitError;
5392    }
5393    R->N = l; /*rv->listLength();*/
5394  }
5395  omFree(R->names);
5396  R->names   = (char **)omAlloc0(R->N * sizeof(char_ptr));
5397  if (rSleftvList2StringArray(rv, R->names))
5398  {
5399    WerrorS("name of ring variable expected");
5400    goto rInitError;
5401  }
5402
5403  /* check names for subring in org_ring ------------------------- */
5404  {
5405    i=0;
5406
5407    for(j=0;j<R->N;j++)
5408    {
5409      for(;i<org_ring->N;i++)
5410      {
5411        if (strcmp(org_ring->names[i],R->names[j])==0)
5412        {
5413          perm[i+1]=j+1;
5414          break;
5415        }
5416      }
5417      if (i>org_ring->N)
5418      {
5419        Werror("variable %d (%s) not in basering",j+1,R->names[j]);
5420        break;
5421      }
5422    }
5423  }
5424  //Print("perm=");
5425  //for(i=1;i<org_ring->N;i++) Print("v%d -> v%d\n",i,perm[i]);
5426  /* ordering -------------------------------------------------------------*/
5427
5428  for(i=0;i<n;i++)
5429  {
5430    int min_var=-1;
5431    int max_var=-1;
5432    for(j=R->block0[i];j<=R->block1[i];j++)
5433    {
5434      if (perm[j]>0)
5435      {
5436        if (min_var==-1) min_var=perm[j];
5437        max_var=perm[j];
5438      }
5439    }
5440    if (min_var!=-1)
5441    {
5442      //Print("block %d: old %d..%d, now:%d..%d\n",
5443      //      i,R->block0[i],R->block1[i],min_var,max_var);
5444      R->block0[i]=min_var;
5445      R->block1[i]=max_var;
5446      if (R->wvhdl[i]!=NULL)
5447      {
5448        omFree(R->wvhdl[i]);
5449        R->wvhdl[i]=(int*)omAlloc0((max_var-min_var+1)*sizeof(int));
5450        for(j=org_ring->block0[i];j<=org_ring->block1[i];j++)
5451        {
5452          if (perm[j]>0)
5453          {
5454            R->wvhdl[i][perm[j]-R->block0[i]]=
5455                org_ring->wvhdl[i][j-org_ring->block0[i]];
5456            //Print("w%d=%d (orig_w%d)\n",perm[j],R->wvhdl[i][perm[j]-R->block0[i]],j);
5457          }
5458        }
5459      }
5460    }
5461    else
5462    {
5463      if(R->block0[i]>0)
5464      {
5465        //Print("skip block %d\n",i);
5466        R->order[i]=ringorder_unspec;
5467        if (R->wvhdl[i] !=NULL) omFree(R->wvhdl[i]);
5468        R->wvhdl[i]=NULL;
5469      }
5470      //else Print("keep block %d\n",i);
5471    }
5472  }
5473  i=n-1;
5474  while(i>0)
5475  {
5476    // removed unneded blocks
5477    if(R->order[i-1]==ringorder_unspec)
5478    {
5479      for(j=i;j<=n;j++)
5480      {
5481        R->order[j-1]=R->order[j];
5482        R->block0[j-1]=R->block0[j];
5483        R->block1[j-1]=R->block1[j];
5484        if (R->wvhdl[j-1] !=NULL) omFree(R->wvhdl[j-1]);
5485        R->wvhdl[j-1]=R->wvhdl[j];
5486      }
5487      R->order[n]=ringorder_unspec;
5488      n--;
5489    }
5490    i--;
5491  }
5492  n=rBlocks(org_ring)-1;
5493  while (R->order[n]==0)  n--;
5494  while (R->order[n]==ringorder_unspec)  n--;
5495  if ((R->order[n]==ringorder_c) ||  (R->order[n]==ringorder_C)) n--;
5496  if (R->block1[n] != R->N)
5497  {
5498    if (((R->order[n]==ringorder_dp) ||
5499         (R->order[n]==ringorder_ds) ||
5500         (R->order[n]==ringorder_Dp) ||
5501         (R->order[n]==ringorder_Ds) ||
5502         (R->order[n]==ringorder_rp) ||
5503         (R->order[n]==ringorder_rs) ||
5504         (R->order[n]==ringorder_lp) ||
5505         (R->order[n]==ringorder_ls))
5506        &&
5507        R->block0[n] <= R->N)
5508    {
5509      R->block1[n] = R->N;
5510    }
5511    else
5512    {
5513      Werror("mismatch of number of vars (%d) and ordering (%d vars) in block %d",
5514             R->N,R->block1[n],n);
5515      return NULL;
5516    }
5517  }
5518  omFree(perm);
5519  // find OrdSgn:
5520  R->OrdSgn = org_ring->OrdSgn; // IMPROVE!
5521  //for(i=1;i<=R->N;i++)
5522  //{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
5523  //omFree(weights);
5524  // Complete the initialization
5525  if (rComplete(R,1))
5526    goto rInitError;
5527
5528  rTest(R);
5529
5530  if (rv != NULL) rv->CleanUp();
5531
5532  return R;
5533
5534  // error case:
5535  rInitError:
5536  if  (R != NULL) rDelete(R);
5537  if (rv != NULL) rv->CleanUp();
5538  return NULL;
5539}
5540
5541void rKill(ring r)
5542{
5543  if ((r->ref<=0)&&(r->order!=NULL))
5544  {
5545#ifdef RDEBUG
5546    if (traceit &TRACE_SHOW_RINGS) Print("kill ring %lx\n",(long)r);
5547#endif
5548    if (r->qideal!=NULL)
5549    {
5550      id_Delete(&r->qideal, r);
5551      r->qideal = NULL;
5552    }
5553    int j;
5554#ifdef USE_IILOCALRING
5555    for (j=0;j<iiRETURNEXPR_len;j++)
5556    {
5557      if (iiLocalRing[j]==r)
5558      {
5559        if (j<myynest) Warn("killing the basering for level %d",j);
5560        iiLocalRing[j]=NULL;
5561      }
5562    }
5563#else /* USE_IILOCALRING */
5564//#endif /* USE_IILOCALRING */
5565    {
5566      proclevel * nshdl = procstack;
5567      int lev=myynest-1;
5568
5569      for(; nshdl != NULL; nshdl = nshdl->next)
5570      {
5571        if (nshdl->cRing==r)
5572        {
5573          Warn("killing the basering for level %d",lev);
5574          nshdl->cRing=NULL;
5575          nshdl->cRingHdl=NULL;
5576        }
5577      }
5578    }
5579#endif /* USE_IILOCALRING */
5580// any variables depending on r ?
5581    while (r->idroot!=NULL)
5582    {
5583      killhdl2(r->idroot,&(r->idroot),r);
5584    }
5585    if (r==currRing)
5586    {
5587      // all dependend stuff is done, clean global vars:
5588      if (r->qideal!=NULL)
5589      {
5590        currQuotient=NULL;
5591      }
5592      if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
5593      if (sLastPrinted.RingDependend())
5594      {
5595        sLastPrinted.CleanUp();
5596      }
5597      if ((myynest>0) && (iiRETURNEXPR.RingDependend()))
5598      {
5599        WerrorS("return value depends on local ring variable (export missing ?)");
5600        iiRETURNEXPR.CleanUp();
5601      }
5602      currRing=NULL;
5603      currRingHdl=NULL;
5604    }
5605
5606    /* nKillChar(r); will be called from inside of rDelete */
5607    rDelete(r);
5608    return;
5609  }
5610  r->ref--;
5611}
5612
5613void rKill(idhdl h)
5614{
5615  ring r = IDRING(h);
5616  int ref=0;
5617  if (r!=NULL)
5618  {
5619    ref=r->ref;
5620    rKill(r);
5621  }
5622  if (h==currRingHdl)
5623  {
5624    if (ref<=0) { currRing=NULL; currRingHdl=NULL;}
5625    else
5626    {
5627      currRingHdl=rFindHdl(r,currRingHdl,NULL);
5628    }
5629  }
5630}
5631
5632idhdl rSimpleFindHdl(ring r, idhdl root, idhdl n)
5633{
5634  //idhdl next_best=NULL;
5635  idhdl h=root;
5636  while (h!=NULL)
5637  {
5638    if (((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD))
5639    && (h!=n)
5640    && (IDRING(h)==r)
5641    )
5642    {
5643   //   if (IDLEV(h)==myynest)
5644   //     return h;
5645   //   if ((IDLEV(h)==0) || (next_best==NULL))
5646   //     next_best=h;
5647   //   else if (IDLEV(next_best)<IDLEV(h))
5648   //     next_best=h;
5649      return h;
5650    }
5651    h=IDNEXT(h);
5652  }
5653  //return next_best;
5654  return NULL;
5655}
5656
5657extern BOOLEAN jjPROC(leftv res, leftv u, leftv v);
5658ideal kGroebner(ideal F, ideal Q)
5659{
5660  //test|=Sy_bit(OPT_PROT);
5661  idhdl save_ringhdl=currRingHdl;
5662  ideal resid;
5663  idhdl new_ring=NULL;
5664  if ((currRingHdl==NULL) || (IDRING(currRingHdl)!=currRing))
5665  {
5666    currRingHdl=enterid(omStrDup(" GROEBNERring"),0,RING_CMD,&IDROOT,FALSE);
5667    new_ring=currRingHdl;
5668    IDRING(currRingHdl)=currRing;
5669  }
5670  sleftv v; memset(&v,0,sizeof(v)); v.rtyp=IDEAL_CMD; v.data=(char *) F;
5671  idhdl h=ggetid("groebner");
5672  sleftv u; memset(&u,0,sizeof(u)); u.rtyp=IDHDL; u.data=(char *) h;
5673            u.name=IDID(h);
5674
5675  sleftv res; memset(&res,0,sizeof(res));
5676  if(jjPROC(&res,&u,&v))
5677  {
5678    resid=kStd(F,Q,testHomog,NULL);
5679  }
5680  else
5681  {
5682    //printf("typ:%d\n",res.rtyp);
5683    resid=(ideal)(res.data);
5684  }
5685  // cleanup GROEBNERring, save_ringhdl, u,v,(res )
5686  if (new_ring!=NULL)
5687  {
5688    idhdl h=IDROOT;
5689    if (h==new_ring) IDROOT=h->next;
5690    else
5691    {
5692      while ((h!=NULL) &&(h->next!=new_ring)) h=h->next;
5693      if (h!=NULL) h->next=h->next->next;
5694    }
5695    if (h!=NULL) omFreeSize(h,sizeof(*h));
5696  }
5697  currRingHdl=save_ringhdl;
5698  u.CleanUp();
5699  v.CleanUp();
5700  return resid;
5701}
5702
5703static void jjINT_S_TO_ID(int n,int *e, leftv res)
5704{
5705  if (n==0) n=1;
5706  ideal l=idInit(n,1);
5707  int i;
5708  poly p;
5709  for(i=rVar(currRing);i>0;i--)
5710  {
5711    if (e[i]>0)
5712    {
5713      n--;
5714      p=pOne();
5715      pSetExp(p,i,1);
5716      pSetm(p);
5717      l->m[n]=p;
5718      if (n==0) break;
5719    }
5720  }
5721  res->data=(char*)l;
5722  setFlag(res,FLAG_STD);
5723  omFreeSize((ADDRESS)e,(rVar(currRing)+1)*sizeof(int));
5724}
5725BOOLEAN jjVARIABLES_P(leftv res, leftv u)
5726{
5727  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5728  int n=pGetVariables((poly)u->Data(),e);
5729  jjINT_S_TO_ID(n,e,res);
5730  return FALSE;
5731}
5732
5733BOOLEAN jjVARIABLES_ID(leftv res, leftv u)
5734{
5735  int *e=(int *)omAlloc0((rVar(currRing)+1)*sizeof(int));
5736  ideal I=(ideal)u->Data();
5737  int i;
5738  int n=0;
5739  for(i=I->nrows*I->ncols-1;i>=0;i--)
5740  {
5741    int n0=pGetVariables(I->m[i],e);
5742    if (n0>n) n=n0;
5743  }
5744  jjINT_S_TO_ID(n,e,res);
5745  return FALSE;
5746}
5747
5748void paPrint(const char *n,package p)
5749{
5750  Print(" %s (",n);
5751  switch (p->language)
5752  {
5753    case LANG_SINGULAR: PrintS("S"); break;
5754    case LANG_C:        PrintS("C"); break;
5755    case LANG_TOP:      PrintS("T"); break;
5756    case LANG_NONE:     PrintS("N"); break;
5757    default:            PrintS("U");
5758  }
5759  if(p->libname!=NULL)
5760  Print(",%s", p->libname);
5761  PrintS(")");
5762}
5763
5764BOOLEAN iiApplyINTVEC(leftv res, leftv a, int op, leftv proc)
5765{
5766  intvec *aa=(intvec*)a->Data();
5767  intvec *r=ivCopy(aa);
5768  sleftv tmp_out;
5769  sleftv tmp_in;
5770  BOOLEAN bo=FALSE;
5771  for(int i=0;i<aa->length(); i++)
5772  {
5773    memset(&tmp_in,0,sizeof(tmp_in));
5774    tmp_in.rtyp=INT_CMD;
5775    tmp_in.data=(void*)(long)(*aa)[i];
5776    if (proc==NULL)
5777      bo=iiExprArith1(&tmp_out,&tmp_in,op);
5778    else
5779      bo=jjPROC(&tmp_out,proc,&tmp_in);
5780    if (bo || (tmp_out.rtyp!=INT_CMD))
5781    {
5782      if (r!=NULL) delete r;
5783      Werror("apply fails at index %d",i+1);
5784      return TRUE;
5785    }
5786    (*r)[i]=(int)(long)tmp_out.data;
5787  }
5788  res->data=(void*)r;
5789  return FALSE;
5790}
5791BOOLEAN iiApplyBIGINTMAT(leftv res, leftv a, int op, leftv proc)
5792{
5793  WerrorS("not implemented");
5794  return TRUE;
5795}
5796BOOLEAN iiApplyIDEAL(leftv res, leftv a, int op, leftv proc)
5797{
5798  WerrorS("not implemented");
5799  return TRUE;
5800}
5801BOOLEAN iiApplyLIST(leftv res, leftv a, int op, leftv proc)
5802{
5803  lists aa=(lists)a->Data();
5804  lists r=(lists)omAlloc0Bin(slists_bin); r->Init(aa->nr+1);
5805  sleftv tmp_out;
5806  sleftv tmp_in;
5807  BOOLEAN bo=FALSE;
5808  for(int i=0;i<=aa->nr; i++)
5809  {
5810    memset(&tmp_in,0,sizeof(tmp_in));
5811    tmp_in.Copy(&(aa->m[i]));
5812    if (proc==NULL)
5813      bo=iiExprArith1(&tmp_out,&tmp_in,op);
5814    else
5815      bo=jjPROC(&tmp_out,proc,&tmp_in);
5816    tmp_in.CleanUp();
5817    if (bo)
5818    {
5819      if (r!=NULL) r->Clean();
5820      Werror("apply fails at index %d",i+1);
5821      return TRUE;
5822    }
5823    memcpy(&(r->m[i]),&tmp_out,sizeof(sleftv));
5824  }
5825  res->data=(void*)r;
5826  return FALSE;
5827  WerrorS("not implemented");
5828  return TRUE;
5829}
5830BOOLEAN iiApply(leftv res, leftv a, int op, leftv proc)
5831{
5832  memset(res,0,sizeof(sleftv));
5833  res->rtyp=a->Typ();
5834  switch (res->rtyp /*a->Typ()*/)
5835  {
5836    case INTVEC_CMD:
5837    case INTMAT_CMD:
5838        return iiApplyINTVEC(res,a,op,proc);
5839    case BIGINTMAT_CMD:
5840        return iiApplyBIGINTMAT(res,a,op,proc);
5841    case IDEAL_CMD:
5842    case MODUL_CMD:
5843    case MATRIX_CMD:
5844        return iiApplyIDEAL(res,a,op,proc);
5845    case LIST_CMD:
5846        return iiApplyLIST(res,a,op,proc);
5847  }
5848  WerrorS("first argument to `apply` must allow an index");
5849  return TRUE;
5850}
Note: See TracBrowser for help on using the repository browser.