source: git/Singular/ipshell.cc @ f9ffa1

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