source: git/Singular/ipshell.cc @ 47ab5b

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