source: git/Singular/ipshell.cc @ 737a68

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