Context Navigation

source: git/Singular/ipshell.cc @ af6ba3

spielwiese

Last change on this file since af6ba3 was af6ba3, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: listvar of package from master
Property mode set to `100644`
File size: 139.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: