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,®,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,®,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*)¶m);
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, ¶m);
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*)¶m);
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