Context Navigation

source: git/Singular/ipshell.cc @ 51c97f

spielwiese

Last change on this file since 51c97f was 51c97f, checked in by Hans Schoenemann <hannes@…>, 12 years ago
partially seperate feResources
Property mode set to `100644`
File size: 135.6 KB

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

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

Download in other formats: