source: git/Singular/ipshell.cc @ 69672d

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