source: git/Singular/ipshell.cc @ 3c0498

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