Context Navigation

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

spielwiese

Last change on this file since 6cc7f5 was 6cc7f5, checked in by Burcin Erocal <burcin@…>, 13 years ago
Convert Singular/ipshell.cc to libpolys (not done).
Property mode set to `100644`
File size: 134.1 KB

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

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

Download in other formats: