source: git/Singular/ipshell.cc @ 85fd90

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