source: git/Singular/ipshell.cc @ fba8c99

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