source: git/Singular/ipshell.cc @ 922a71f

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