source: git/Singular/ipshell.cc @ 8c1285

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