source: git/Singular/ipshell.cc @ d8c7a6

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