source: git/Singular/ipshell.cc @ 6b4ff12

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