source: git/Singular/ipshell.cc @ 4b7db8

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