source: git/Singular/ipshell.cc @ b9642b1

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