source: git/Singular/ipshell.cc @ c25b56c

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