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