source: git/Singular/ipshell.cc @ f0168a5

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