source: git/Singular/ipshell.cc @ 9f7665

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