source: git/Singular/ipshell.cc @ 6ce030f

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