source: git/Singular/ipshell.cc @ 9752db

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