source: git/Singular/ipshell.cc @ 883eacf

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