Context Navigation

source: git/Singular/ipshell.cc @ e050c2

fieker-DuValspielwiese

Last change on this file since e050c2 was 7d58b6, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: use currRing instead of currRingHdl (from master)
Property mode set to `100644`
File size: 135.4 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	#ifdef HAVE_FACTORY
2824	BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
2825	{
2826	res->data=singclap_resultant((poly)u->CopyD(),(poly)v->CopyD(),
2827	(poly)w->CopyD(), currRing);
2828	return errorreported;
2829	}
2830	BOOLEAN jjCHARSERIES(leftv res, leftv u)
2831	{
2832	#if 1
2833	Werror("Sorry: not yet re-factored: see libpolys/polys/clapsing.cc");
2834	return FALSE;
2835	#else
2836	res->data=singclap_irrCharSeries((ideal)u->Data(), currRing);
2837	return (res->data==NULL);
2838	#endif
2839	}
2840	#endif
2841
2842	// from semic.cc
2843	#ifdef HAVE_SPECTRUM
2844
2845	// ----------------------------------------------------------------------------
2846	// Initialize a spectrum deep from a singular lists
2847	// ----------------------------------------------------------------------------
2848
2849	void copy_deep( spectrum& spec, lists l )
2850	{
2851	spec.mu = (int)(long)(l->m[0].Data( ));
2852	spec.pg = (int)(long)(l->m[1].Data( ));
2853	spec.n = (int)(long)(l->m[2].Data( ));
2854
2855	spec.copy_new( spec.n );
2856
2857	intvec num = (intvec)l->m[3].Data( );
2858	intvec den = (intvec)l->m[4].Data( );
2859	intvec mul = (intvec)l->m[5].Data( );
2860
2861	for( int i=0; i<spec.n; i++ )
2862	{
2863	spec.s[i] = (Rational)((num)[i])/(Rational)((den)[i]);
2864	spec.w[i] = (*mul)[i];
2865	}
2866	}
2867
2868	// ----------------------------------------------------------------------------
2869	// singular lists constructor for spectrum
2870	// ----------------------------------------------------------------------------
2871
2872	spectrum /former spectrum::spectrum ( lists l )/
2873	spectrumFromList( lists l )
2874	{
2875	spectrum result;
2876	copy_deep( result, l );
2877	return result;
2878	}
2879
2880	// ----------------------------------------------------------------------------
2881	// generate a Singular lists from a spectrum
2882	// ----------------------------------------------------------------------------
2883
2884	/* former spectrum::thelist ( void )*/
2885	lists getList( spectrum& spec )
2886	{
2887	lists L = (lists)omAllocBin( slists_bin);
2888
2889	L->Init( 6 );
2890
2891	intvec *num = new intvec( spec.n );
2892	intvec *den = new intvec( spec.n );
2893	intvec *mult = new intvec( spec.n );
2894
2895	for( int i=0; i<spec.n; i++ )
2896	{
2897	(*num) [i] = spec.s[i].get_num_si( );
2898	(*den) [i] = spec.s[i].get_den_si( );
2899	(*mult)[i] = spec.w[i];
2900	}
2901
2902	L->m[0].rtyp = INT_CMD; // milnor number
2903	L->m[1].rtyp = INT_CMD; // geometrical genus
2904	L->m[2].rtyp = INT_CMD; // # of spectrum numbers
2905	L->m[3].rtyp = INTVEC_CMD; // numerators
2906	L->m[4].rtyp = INTVEC_CMD; // denomiantors
2907	L->m[5].rtyp = INTVEC_CMD; // multiplicities
2908
2909	L->m[0].data = (void*)spec.mu;
2910	L->m[1].data = (void*)spec.pg;
2911	L->m[2].data = (void*)spec.n;
2912	L->m[3].data = (void*)num;
2913	L->m[4].data = (void*)den;
2914	L->m[5].data = (void*)mult;
2915
2916	return L;
2917	}
2918	// from spectrum.cc
2919	// ----------------------------------------------------------------------------
2920	// print out an error message for a spectrum list
2921	// ----------------------------------------------------------------------------
2922
2923	typedef enum
2924	{
2925	semicOK,
2926	semicMulNegative,
2927
2928	semicListTooShort,
2929	semicListTooLong,
2930
2931	semicListFirstElementWrongType,
2932	semicListSecondElementWrongType,
2933	semicListThirdElementWrongType,
2934	semicListFourthElementWrongType,
2935	semicListFifthElementWrongType,
2936	semicListSixthElementWrongType,
2937
2938	semicListNNegative,
2939	semicListWrongNumberOfNumerators,
2940	semicListWrongNumberOfDenominators,
2941	semicListWrongNumberOfMultiplicities,
2942
2943	semicListMuNegative,
2944	semicListPgNegative,
2945	semicListNumNegative,
2946	semicListDenNegative,
2947	semicListMulNegative,
2948
2949	semicListNotSymmetric,
2950	semicListNotMonotonous,
2951
2952	semicListMilnorWrong,
2953	semicListPGWrong
2954
2955	} semicState;
2956
2957	void list_error( semicState state )
2958	{
2959	switch( state )
2960	{
2961	case semicListTooShort:
2962	WerrorS( "the list is too short" );
2963	break;
2964	case semicListTooLong:
2965	WerrorS( "the list is too long" );
2966	break;
2967
2968	case semicListFirstElementWrongType:
2969	WerrorS( "first element of the list should be int" );
2970	break;
2971	case semicListSecondElementWrongType:
2972	WerrorS( "second element of the list should be int" );
2973	break;
2974	case semicListThirdElementWrongType:
2975	WerrorS( "third element of the list should be int" );
2976	break;
2977	case semicListFourthElementWrongType:
2978	WerrorS( "fourth element of the list should be intvec" );
2979	break;
2980	case semicListFifthElementWrongType:
2981	WerrorS( "fifth element of the list should be intvec" );
2982	break;
2983	case semicListSixthElementWrongType:
2984	WerrorS( "sixth element of the list should be intvec" );
2985	break;
2986
2987	case semicListNNegative:
2988	WerrorS( "first element of the list should be positive" );
2989	break;
2990	case semicListWrongNumberOfNumerators:
2991	WerrorS( "wrong number of numerators" );
2992	break;
2993	case semicListWrongNumberOfDenominators:
2994	WerrorS( "wrong number of denominators" );
2995	break;
2996	case semicListWrongNumberOfMultiplicities:
2997	WerrorS( "wrong number of multiplicities" );
2998	break;
2999
3000	case semicListMuNegative:
3001	WerrorS( "the Milnor number should be positive" );
3002	break;
3003	case semicListPgNegative:
3004	WerrorS( "the geometrical genus should be nonnegative" );
3005	break;
3006	case semicListNumNegative:
3007	WerrorS( "all numerators should be positive" );
3008	break;
3009	case semicListDenNegative:
3010	WerrorS( "all denominators should be positive" );
3011	break;
3012	case semicListMulNegative:
3013	WerrorS( "all multiplicities should be positive" );
3014	break;
3015
3016	case semicListNotSymmetric:
3017	WerrorS( "it is not symmetric" );
3018	break;
3019	case semicListNotMonotonous:
3020	WerrorS( "it is not monotonous" );
3021	break;
3022
3023	case semicListMilnorWrong:
3024	WerrorS( "the Milnor number is wrong" );
3025	break;
3026	case semicListPGWrong:
3027	WerrorS( "the geometrical genus is wrong" );
3028	break;
3029
3030	default:
3031	WerrorS( "unspecific error" );
3032	break;
3033	}
3034	}
3035	// ----------------------------------------------------------------------------
3036	// this is the main spectrum computation function
3037	// ----------------------------------------------------------------------------
3038
3039	enum spectrumState
3040	{
3041	spectrumOK,
3042	spectrumZero,
3043	spectrumBadPoly,
3044	spectrumNoSingularity,
3045	spectrumNotIsolated,
3046	spectrumDegenerate,
3047	spectrumWrongRing,
3048	spectrumNoHC,
3049	spectrumUnspecErr
3050	};
3051
3052	// from splist.cc
3053	// ----------------------------------------------------------------------------
3054	// Compute the spectrum of a spectrumPolyList
3055	// ----------------------------------------------------------------------------
3056
3057	/* former spectrumPolyList::spectrum ( lists, int) /
3058	spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L,int fast )
3059	{
3060	spectrumPolyNode **node = &speclist.root;
3061	spectrumPolyNode *search;
3062
3063	poly f,tmp;
3064	int found,cmp;
3065
3066	Rational smax( ( fast==0 ? 0 : rVar(currRing) ),
3067	( fast==2 ? 2 : 1 ) );
3068
3069	Rational weight_prev( 0,1 );
3070
3071	int mu = 0; // the milnor number
3072	int pg = 0; // the geometrical genus
3073	int n = 0; // number of different spectral numbers
3074	int z = 0; // number of spectral number equal to smax
3075
3076	int k = 0;
3077
3078	while( (node)!=(spectrumPolyNode)NULL &&
3079	( fast==0 \|\| (*node)->weight<=smax ) )
3080	{
3081	// ---------------------------------------
3082	// determine the first normal form which
3083	// contains the monomial node->mon
3084	// ---------------------------------------
3085
3086	found = FALSE;
3087	search = *node;
3088
3089	while( search!=(spectrumPolyNode*)NULL && found==FALSE )
3090	{
3091	if( search->nf!=(poly)NULL )
3092	{
3093	f = search->nf;
3094
3095	do
3096	{
3097	// --------------------------------
3098	// look for (*node)->mon in f
3099	// --------------------------------
3100
3101	cmp = pCmp( (*node)->mon,f );
3102
3103	if( cmp<0 )
3104	{
3105	f = pNext( f );
3106	}
3107	else if( cmp==0 )
3108	{
3109	// -----------------------------
3110	// we have found a normal form
3111	// -----------------------------
3112
3113	found = TRUE;
3114
3115	// normalize coefficient
3116
3117	number inv = nInvers( pGetCoeff( f ) );
3118	pMult_nn( search->nf,inv );
3119	nDelete( &inv );
3120
3121	// exchange normal forms
3122
3123	tmp = (*node)->nf;
3124	(*node)->nf = search->nf;
3125	search->nf = tmp;
3126	}
3127	}
3128	while( cmp<0 && f!=(poly)NULL );
3129	}
3130	search = search->next;
3131	}
3132
3133	if( found==FALSE )
3134	{
3135	// ------------------------------------------------
3136	// the weight of node->mon is a spectrum number
3137	// ------------------------------------------------
3138
3139	mu++;
3140
3141	if( (*node)->weight<=(Rational)1 ) pg++;
3142	if( (*node)->weight==smax ) z++;
3143	if( (*node)->weight>weight_prev ) n++;
3144
3145	weight_prev = (*node)->weight;
3146	node = &((*node)->next);
3147	}
3148	else
3149	{
3150	// -----------------------------------------------
3151	// determine all other normal form which contain
3152	// the monomial node->mon
3153	// replace for node->mon its normal form
3154	// -----------------------------------------------
3155
3156	while( search!=(spectrumPolyNode*)NULL )
3157	{
3158	if( search->nf!=(poly)NULL )
3159	{
3160	f = search->nf;
3161
3162	do
3163	{
3164	// --------------------------------
3165	// look for (*node)->mon in f
3166	// --------------------------------
3167
3168	cmp = pCmp( (*node)->mon,f );
3169
3170	if( cmp<0 )
3171	{
3172	f = pNext( f );
3173	}
3174	else if( cmp==0 )
3175	{
3176	search->nf = pSub( search->nf,
3177	ppMult_nn( (*node)->nf,pGetCoeff( f ) ) );
3178	pNorm( search->nf );
3179	}
3180	}
3181	while( cmp<0 && f!=(poly)NULL );
3182	}
3183	search = search->next;
3184	}
3185	speclist.delete_node( node );
3186	}
3187
3188	}
3189
3190	// --------------------------------------------------------
3191	// fast computation exploits the symmetry of the spectrum
3192	// --------------------------------------------------------
3193
3194	if( fast==2 )
3195	{
3196	mu = 2*mu - z;
3197	n = ( z > 0 ? 2n - 1 : 2n );
3198	}
3199
3200	// --------------------------------------------------------
3201	// compute the spectrum numbers with their multiplicities
3202	// --------------------------------------------------------
3203
3204	intvec *nom = new intvec( n );
3205	intvec *den = new intvec( n );
3206	intvec *mult = new intvec( n );
3207
3208	int count = 0;
3209	int multiplicity = 1;
3210
3211	for( search=speclist.root; search!=(spectrumPolyNode*)NULL &&
3212	( fast==0 \|\| search->weight<=smax );
3213	search=search->next )
3214	{
3215	if( search->next==(spectrumPolyNode*)NULL \|\|
3216	search->weight<search->next->weight )
3217	{
3218	(*nom) [count] = search->weight.get_num_si( );
3219	(*den) [count] = search->weight.get_den_si( );
3220	(*mult)[count] = multiplicity;
3221
3222	multiplicity=1;
3223	count++;
3224	}
3225	else
3226	{
3227	multiplicity++;
3228	}
3229	}
3230
3231	// --------------------------------------------------------
3232	// fast computation exploits the symmetry of the spectrum
3233	// --------------------------------------------------------
3234
3235	if( fast==2 )
3236	{
3237	int n1,n2;
3238	for( n1=0, n2=n-1; n1<n2; n1++, n2-- )
3239	{
3240	(nom) [n2] = rVar(currRing)(den)[n1]-(nom)[n1];
3241	(den) [n2] = (den)[n1];
3242	(mult)[n2] = (mult)[n1];
3243	}
3244	}
3245
3246	// -----------------------------------
3247	// test if the spectrum is symmetric
3248	// -----------------------------------
3249
3250	if( fast==0 \|\| fast==1 )
3251	{
3252	int symmetric=TRUE;
3253
3254	for( int n1=0, n2=n-1 ; n1<n2 && symmetric==TRUE; n1++, n2-- )
3255	{
3256	if( (mult)[n1]!=(mult)[n2] \|\|
3257	(den) [n1]!= (den)[n2] \|\|
3258	(nom)[n1]+(nom)[n2]!=rVar(currRing)(den) [n1] )
3259	{
3260	symmetric = FALSE;
3261	}
3262	}
3263
3264	if( symmetric==FALSE )
3265	{
3266	// ---------------------------------------------
3267	// the spectrum is not symmetric => degenerate
3268	// principal part
3269	// ---------------------------------------------
3270
3271	*L = (lists)omAllocBin( slists_bin);
3272	(*L)->Init( 1 );
3273	(*L)->m[0].rtyp = INT_CMD; // milnor number
3274	(L)->m[0].data = (void)mu;
3275
3276	return spectrumDegenerate;
3277	}
3278	}
3279
3280	*L = (lists)omAllocBin( slists_bin);
3281
3282	(*L)->Init( 6 );
3283
3284	(*L)->m[0].rtyp = INT_CMD; // milnor number
3285	(*L)->m[1].rtyp = INT_CMD; // geometrical genus
3286	(*L)->m[2].rtyp = INT_CMD; // number of spectrum values
3287	(*L)->m[3].rtyp = INTVEC_CMD; // nominators
3288	(*L)->m[4].rtyp = INTVEC_CMD; // denomiantors
3289	(*L)->m[5].rtyp = INTVEC_CMD; // multiplicities
3290
3291	(L)->m[0].data = (void)mu;
3292	(L)->m[1].data = (void)pg;
3293	(L)->m[2].data = (void)n;
3294	(L)->m[3].data = (void)nom;
3295	(L)->m[4].data = (void)den;
3296	(L)->m[5].data = (void)mult;
3297
3298	return spectrumOK;
3299	}
3300
3301	spectrumState spectrumCompute( poly h,lists *L,int fast )
3302	{
3303	int i,j;
3304
3305	#ifdef SPECTRUM_DEBUG
3306	#ifdef SPECTRUM_PRINT
3307	#ifdef SPECTRUM_IOSTREAM
3308	cout << "spectrumCompute\n";
3309	if( fast==0 ) cout << " no optimization" << endl;
3310	if( fast==1 ) cout << " weight optimization" << endl;
3311	if( fast==2 ) cout << " symmetry optimization" << endl;
3312	#else
3313	fprintf( stdout,"spectrumCompute\n" );
3314	if( fast==0 ) fprintf( stdout," no optimization\n" );
3315	if( fast==1 ) fprintf( stdout," weight optimization\n" );
3316	if( fast==2 ) fprintf( stdout," symmetry optimization\n" );
3317	#endif
3318	#endif
3319	#endif
3320
3321	// ----------------------
3322	// check if h is zero
3323	// ----------------------
3324
3325	if( h==(poly)NULL )
3326	{
3327	return spectrumZero;
3328	}
3329
3330	// ----------------------------------
3331	// check if h has a constant term
3332	// ----------------------------------
3333
3334	if( hasConstTerm( h, currRing ) )
3335	{
3336	return spectrumBadPoly;
3337	}
3338
3339	// --------------------------------
3340	// check if h has a linear term
3341	// --------------------------------
3342
3343	if( hasLinearTerm( h, currRing ) )
3344	{
3345	*L = (lists)omAllocBin( slists_bin);
3346	(*L)->Init( 1 );
3347	(*L)->m[0].rtyp = INT_CMD; // milnor number
3348	/* (L)->m[0].data = (void)0;a -- done by Init */
3349
3350	return spectrumNoSingularity;
3351	}
3352
3353	// ----------------------------------
3354	// compute the jacobi ideal of (h)
3355	// ----------------------------------
3356
3357	ideal J = NULL;
3358	J = idInit( rVar(currRing),1 );
3359
3360	#ifdef SPECTRUM_DEBUG
3361	#ifdef SPECTRUM_PRINT
3362	#ifdef SPECTRUM_IOSTREAM
3363	cout << "\n computing the Jacobi ideal...\n";
3364	#else
3365	fprintf( stdout,"\n computing the Jacobi ideal...\n" );
3366	#endif
3367	#endif
3368	#endif
3369
3370	for( i=0; i<rVar(currRing); i++ )
3371	{
3372	J->m[i] = pDiff( h,i+1); //j );
3373
3374	#ifdef SPECTRUM_DEBUG
3375	#ifdef SPECTRUM_PRINT
3376	#ifdef SPECTRUM_IOSTREAM
3377	cout << " ";
3378	#else
3379	fprintf( stdout," " );
3380	#endif
3381	pWrite( J->m[i] );
3382	#endif
3383	#endif
3384	}
3385
3386	// --------------------------------------------
3387	// compute a standard basis stdJ of jac(h)
3388	// --------------------------------------------
3389
3390	#ifdef SPECTRUM_DEBUG
3391	#ifdef SPECTRUM_PRINT
3392	#ifdef SPECTRUM_IOSTREAM
3393	cout << endl;
3394	cout << " computing a standard basis..." << endl;
3395	#else
3396	fprintf( stdout,"\n" );
3397	fprintf( stdout," computing a standard basis...\n" );
3398	#endif
3399	#endif
3400	#endif
3401
3402	ideal stdJ = kStd(J,currQuotient,isNotHomog,NULL);
3403	idSkipZeroes( stdJ );
3404
3405	#ifdef SPECTRUM_DEBUG
3406	#ifdef SPECTRUM_PRINT
3407	for( i=0; i<IDELEMS(stdJ); i++ )
3408	{
3409	#ifdef SPECTRUM_IOSTREAM
3410	cout << " ";
3411	#else
3412	fprintf( stdout," " );
3413	#endif
3414
3415	pWrite( stdJ->m[i] );
3416	}
3417	#endif
3418	#endif
3419
3420	idDelete( &J );
3421
3422	// ------------------------------------------
3423	// check if the h has a singularity
3424	// ------------------------------------------
3425
3426	if( hasOne( stdJ, currRing ) )
3427	{
3428	// -------------------------------
3429	// h is smooth in the origin
3430	// return only the Milnor number
3431	// -------------------------------
3432
3433	*L = (lists)omAllocBin( slists_bin);
3434	(*L)->Init( 1 );
3435	(*L)->m[0].rtyp = INT_CMD; // milnor number
3436	/* (L)->m[0].data = (void)0;a -- done by Init */
3437
3438	return spectrumNoSingularity;
3439	}
3440
3441	// ------------------------------------------
3442	// check if the singularity h is isolated
3443	// ------------------------------------------
3444
3445	for( i=rVar(currRing); i>0; i-- )
3446	{
3447	if( hasAxis( stdJ,i, currRing )==FALSE )
3448	{
3449	return spectrumNotIsolated;
3450	}
3451	}
3452
3453	// ------------------------------------------
3454	// compute the highest corner hc of stdJ
3455	// ------------------------------------------
3456
3457	#ifdef SPECTRUM_DEBUG
3458	#ifdef SPECTRUM_PRINT
3459	#ifdef SPECTRUM_IOSTREAM
3460	cout << "\n computing the highest corner...\n";
3461	#else
3462	fprintf( stdout,"\n computing the highest corner...\n" );
3463	#endif
3464	#endif
3465	#endif
3466
3467	poly hc = (poly)NULL;
3468
3469	scComputeHC( stdJ,currQuotient, 0,hc );
3470
3471	if( hc!=(poly)NULL )
3472	{
3473	pGetCoeff(hc) = nInit(1);
3474
3475	for( i=rVar(currRing); i>0; i-- )
3476	{
3477	if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
3478	}
3479	pSetm( hc );
3480	}
3481	else
3482	{
3483	return spectrumNoHC;
3484	}
3485
3486	#ifdef SPECTRUM_DEBUG
3487	#ifdef SPECTRUM_PRINT
3488	#ifdef SPECTRUM_IOSTREAM
3489	cout << " ";
3490	#else
3491	fprintf( stdout," " );
3492	#endif
3493	pWrite( hc );
3494	#endif
3495	#endif
3496
3497	// ----------------------------------------
3498	// compute the Newton polygon nph of h
3499	// ----------------------------------------
3500
3501	#ifdef SPECTRUM_DEBUG
3502	#ifdef SPECTRUM_PRINT
3503	#ifdef SPECTRUM_IOSTREAM
3504	cout << "\n computing the newton polygon...\n";
3505	#else
3506	fprintf( stdout,"\n computing the newton polygon...\n" );
3507	#endif
3508	#endif
3509	#endif
3510
3511	newtonPolygon nph( h, currRing );
3512
3513	#ifdef SPECTRUM_DEBUG
3514	#ifdef SPECTRUM_PRINT
3515	cout << nph;
3516	#endif
3517	#endif
3518
3519	// -----------------------------------------------
3520	// compute the weight corner wc of (stdj,nph)
3521	// -----------------------------------------------
3522
3523	#ifdef SPECTRUM_DEBUG
3524	#ifdef SPECTRUM_PRINT
3525	#ifdef SPECTRUM_IOSTREAM
3526	cout << "\n computing the weight corner...\n";
3527	#else
3528	fprintf( stdout,"\n computing the weight corner...\n" );
3529	#endif
3530	#endif
3531	#endif
3532
3533	poly wc = ( fast==0 ? pCopy( hc ) :
3534	( fast==1 ? computeWC( nph,(Rational)rVar(currRing), currRing ) :
3535	/* fast==2 */computeWC( nph,
3536	((Rational)rVar(currRing))/(Rational)2, currRing ) ) );
3537
3538	#ifdef SPECTRUM_DEBUG
3539	#ifdef SPECTRUM_PRINT
3540	#ifdef SPECTRUM_IOSTREAM
3541	cout << " ";
3542	#else
3543	fprintf( stdout," " );
3544	#endif
3545	pWrite( wc );
3546	#endif
3547	#endif
3548
3549	// -------------
3550	// compute NF
3551	// -------------
3552
3553	#ifdef SPECTRUM_DEBUG
3554	#ifdef SPECTRUM_PRINT
3555	#ifdef SPECTRUM_IOSTREAM
3556	cout << "\n computing NF...\n" << endl;
3557	#else
3558	fprintf( stdout,"\n computing NF...\n" );
3559	#endif
3560	#endif
3561	#endif
3562
3563	spectrumPolyList NF( &nph );
3564
3565	computeNF( stdJ,hc,wc,&NF, currRing );
3566
3567	#ifdef SPECTRUM_DEBUG
3568	#ifdef SPECTRUM_PRINT
3569	cout << NF;
3570	#ifdef SPECTRUM_IOSTREAM
3571	cout << endl;
3572	#else
3573	fprintf( stdout,"\n" );
3574	#endif
3575	#endif
3576	#endif
3577
3578	// ----------------------------
3579	// compute the spectrum of h
3580	// ----------------------------
3581	// spectrumState spectrumStateFromList( spectrumPolyList& speclist, lists *L, int fast );
3582
3583	return spectrumStateFromList(NF, L, fast );
3584	}
3585
3586	// ----------------------------------------------------------------------------
3587	// this procedure is called from the interpreter
3588	// ----------------------------------------------------------------------------
3589	// first = polynomial
3590	// result = list of spectrum numbers
3591	// ----------------------------------------------------------------------------
3592
3593	void spectrumPrintError(spectrumState state)
3594	{
3595	switch( state )
3596	{
3597	case spectrumZero:
3598	WerrorS( "polynomial is zero" );
3599	break;
3600	case spectrumBadPoly:
3601	WerrorS( "polynomial has constant term" );
3602	break;
3603	case spectrumNoSingularity:
3604	WerrorS( "not a singularity" );
3605	break;
3606	case spectrumNotIsolated:
3607	WerrorS( "the singularity is not isolated" );
3608	break;
3609	case spectrumNoHC:
3610	WerrorS( "highest corner cannot be computed" );
3611	break;
3612	case spectrumDegenerate:
3613	WerrorS( "principal part is degenerate" );
3614	break;
3615	case spectrumOK:
3616	break;
3617
3618	default:
3619	WerrorS( "unknown error occurred" );
3620	break;
3621	}
3622	}
3623
3624	BOOLEAN spectrumProc( leftv result,leftv first )
3625	{
3626	spectrumState state = spectrumOK;
3627
3628	// -------------------
3629	// check consistency
3630	// -------------------
3631
3632	// check for a local ring
3633
3634	if( !ringIsLocal(currRing ) )
3635	{
3636	WerrorS( "only works for local orderings" );
3637	state = spectrumWrongRing;
3638	}
3639
3640	// no quotient rings are allowed
3641
3642	else if( currRing->qideal != NULL )
3643	{
3644	WerrorS( "does not work in quotient rings" );
3645	state = spectrumWrongRing;
3646	}
3647	else
3648	{
3649	lists L = (lists)NULL;
3650	int flag = 1; // weight corner optimization is safe
3651
3652	state = spectrumCompute( (poly)first->Data( ),&L,flag );
3653
3654	if( state==spectrumOK )
3655	{
3656	result->rtyp = LIST_CMD;
3657	result->data = (char*)L;
3658	}
3659	else
3660	{
3661	spectrumPrintError(state);
3662	}
3663	}
3664
3665	return (state!=spectrumOK);
3666	}
3667
3668	// ----------------------------------------------------------------------------
3669	// this procedure is called from the interpreter
3670	// ----------------------------------------------------------------------------
3671	// first = polynomial
3672	// result = list of spectrum numbers
3673	// ----------------------------------------------------------------------------
3674
3675	BOOLEAN spectrumfProc( leftv result,leftv first )
3676	{
3677	spectrumState state = spectrumOK;
3678
3679	// -------------------
3680	// check consistency
3681	// -------------------
3682
3683	// check for a local polynomial ring
3684
3685	if( currRing->OrdSgn != -1 )
3686	// ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
3687	// or should we use:
3688	//if( !ringIsLocal( ) )
3689	{
3690	WerrorS( "only works for local orderings" );
3691	state = spectrumWrongRing;
3692	}
3693	else if( currRing->qideal != NULL )
3694	{
3695	WerrorS( "does not work in quotient rings" );
3696	state = spectrumWrongRing;
3697	}
3698	else
3699	{
3700	lists L = (lists)NULL;
3701	int flag = 2; // symmetric optimization
3702
3703	state = spectrumCompute( (poly)first->Data( ),&L,flag );
3704
3705	if( state==spectrumOK )
3706	{
3707	result->rtyp = LIST_CMD;
3708	result->data = (char*)L;
3709	}
3710	else
3711	{
3712	spectrumPrintError(state);
3713	}
3714	}
3715
3716	return (state!=spectrumOK);
3717	}
3718
3719	// ----------------------------------------------------------------------------
3720	// check if a list is a spectrum
3721	// check for:
3722	// list has 6 elements
3723	// 1st element is int (mu=Milnor number)
3724	// 2nd element is int (pg=geometrical genus)
3725	// 3rd element is int (n =number of different spectrum numbers)
3726	// 4th element is intvec (num=numerators)
3727	// 5th element is intvec (den=denomiantors)
3728	// 6th element is intvec (mul=multiplicities)
3729	// exactly n numerators
3730	// exactly n denominators
3731	// exactly n multiplicities
3732	// mu>0
3733	// pg>=0
3734	// n>0
3735	// num>0
3736	// den>0
3737	// mul>0
3738	// symmetriy with respect to numberofvariables/2
3739	// monotony
3740	// mu = sum of all multiplicities
3741	// pg = sum of all multiplicities where num/den<=1
3742	// ----------------------------------------------------------------------------
3743
3744	semicState list_is_spectrum( lists l )
3745	{
3746	// -------------------
3747	// check list length
3748	// -------------------
3749
3750	if( l->nr < 5 )
3751	{
3752	return semicListTooShort;
3753	}
3754	else if( l->nr > 5 )
3755	{
3756	return semicListTooLong;
3757	}
3758
3759	// -------------
3760	// check types
3761	// -------------
3762
3763	if( l->m[0].rtyp != INT_CMD )
3764	{
3765	return semicListFirstElementWrongType;
3766	}
3767	else if( l->m[1].rtyp != INT_CMD )
3768	{
3769	return semicListSecondElementWrongType;
3770	}
3771	else if( l->m[2].rtyp != INT_CMD )
3772	{
3773	return semicListThirdElementWrongType;
3774	}
3775	else if( l->m[3].rtyp != INTVEC_CMD )
3776	{
3777	return semicListFourthElementWrongType;
3778	}
3779	else if( l->m[4].rtyp != INTVEC_CMD )
3780	{
3781	return semicListFifthElementWrongType;
3782	}
3783	else if( l->m[5].rtyp != INTVEC_CMD )
3784	{
3785	return semicListSixthElementWrongType;
3786	}
3787
3788	// -------------------------
3789	// check number of entries
3790	// -------------------------
3791
3792	int mu = (int)(long)(l->m[0].Data( ));
3793	int pg = (int)(long)(l->m[1].Data( ));
3794	int n = (int)(long)(l->m[2].Data( ));
3795
3796	if( n <= 0 )
3797	{
3798	return semicListNNegative;
3799	}
3800
3801	intvec num = (intvec)l->m[3].Data( );
3802	intvec den = (intvec)l->m[4].Data( );
3803	intvec mul = (intvec)l->m[5].Data( );
3804
3805	if( n != num->length( ) )
3806	{
3807	return semicListWrongNumberOfNumerators;
3808	}
3809	else if( n != den->length( ) )
3810	{
3811	return semicListWrongNumberOfDenominators;
3812	}
3813	else if( n != mul->length( ) )
3814	{
3815	return semicListWrongNumberOfMultiplicities;
3816	}
3817
3818	// --------
3819	// values
3820	// --------
3821
3822	if( mu <= 0 )
3823	{
3824	return semicListMuNegative;
3825	}
3826	if( pg < 0 )
3827	{
3828	return semicListPgNegative;
3829	}
3830
3831	int i;
3832
3833	for( i=0; i<n; i++ )
3834	{
3835	if( (*num)[i] <= 0 )
3836	{
3837	return semicListNumNegative;
3838	}
3839	if( (*den)[i] <= 0 )
3840	{
3841	return semicListDenNegative;
3842	}
3843	if( (*mul)[i] <= 0 )
3844	{
3845	return semicListMulNegative;
3846	}
3847	}
3848
3849	// ----------------
3850	// check symmetry
3851	// ----------------
3852
3853	int j;
3854
3855	for( i=0, j=n-1; i<=j; i++,j-- )
3856	{
3857	if( (num)[i] != rVar(currRing)((den)[i]) - (num)[j] \|\|
3858	(den)[i] != (den)[j] \|\|
3859	(mul)[i] != (mul)[j] )
3860	{
3861	return semicListNotSymmetric;
3862	}
3863	}
3864
3865	// ----------------
3866	// check monotony
3867	// ----------------
3868
3869	for( i=0, j=1; i<n/2; i++,j++ )
3870	{
3871	if( (num)[i](den)[j] >= (num)[j](den)[i] )
3872	{
3873	return semicListNotMonotonous;
3874	}
3875	}
3876
3877	// ---------------------
3878	// check Milnor number
3879	// ---------------------
3880
3881	for( mu=0, i=0; i<n; i++ )
3882	{
3883	mu += (*mul)[i];
3884	}
3885
3886	if( mu != (int)(long)(l->m[0].Data( )) )
3887	{
3888	return semicListMilnorWrong;
3889	}
3890
3891	// -------------------------
3892	// check geometrical genus
3893	// -------------------------
3894
3895	for( pg=0, i=0; i<n; i++ )
3896	{
3897	if( (num)[i]<=(den)[i] )
3898	{
3899	pg += (*mul)[i];
3900	}
3901	}
3902
3903	if( pg != (int)(long)(l->m[1].Data( )) )
3904	{
3905	return semicListPGWrong;
3906	}
3907
3908	return semicOK;
3909	}
3910
3911	// ----------------------------------------------------------------------------
3912	// this procedure is called from the interpreter
3913	// ----------------------------------------------------------------------------
3914	// first = list of spectrum numbers
3915	// second = list of spectrum numbers
3916	// result = sum of the two lists
3917	// ----------------------------------------------------------------------------
3918
3919	BOOLEAN spaddProc( leftv result,leftv first,leftv second )
3920	{
3921	semicState state;
3922
3923	// -----------------
3924	// check arguments
3925	// -----------------
3926
3927	lists l1 = (lists)first->Data( );
3928	lists l2 = (lists)second->Data( );
3929
3930	if( (state=list_is_spectrum( l1 )) != semicOK )
3931	{
3932	WerrorS( "first argument is not a spectrum:" );
3933	list_error( state );
3934	}
3935	else if( (state=list_is_spectrum( l2 )) != semicOK )
3936	{
3937	WerrorS( "second argument is not a spectrum:" );
3938	list_error( state );
3939	}
3940	else
3941	{
3942	spectrum s1= spectrumFromList ( l1 );
3943	spectrum s2= spectrumFromList ( l2 );
3944	spectrum sum( s1+s2 );
3945
3946	result->rtyp = LIST_CMD;
3947	result->data = (char*)(getList(sum));
3948	}
3949
3950	return (state!=semicOK);
3951	}
3952
3953	// ----------------------------------------------------------------------------
3954	// this procedure is called from the interpreter
3955	// ----------------------------------------------------------------------------
3956	// first = list of spectrum numbers
3957	// second = integer
3958	// result = the multiple of the first list by the second factor
3959	// ----------------------------------------------------------------------------
3960
3961	BOOLEAN spmulProc( leftv result,leftv first,leftv second )
3962	{
3963	semicState state;
3964
3965	// -----------------
3966	// check arguments
3967	// -----------------
3968
3969	lists l = (lists)first->Data( );
3970	int k = (int)(long)second->Data( );
3971
3972	if( (state=list_is_spectrum( l ))!=semicOK )
3973	{
3974	WerrorS( "first argument is not a spectrum" );
3975	list_error( state );
3976	}
3977	else if( k < 0 )
3978	{
3979	WerrorS( "second argument should be positive" );
3980	state = semicMulNegative;
3981	}
3982	else
3983	{
3984	spectrum s= spectrumFromList( l );
3985	spectrum product( k*s );
3986
3987	result->rtyp = LIST_CMD;
3988	result->data = (char*)getList(product);
3989	}
3990
3991	return (state!=semicOK);
3992	}
3993
3994	// ----------------------------------------------------------------------------
3995	// this procedure is called from the interpreter
3996	// ----------------------------------------------------------------------------
3997	// first = list of spectrum numbers
3998	// second = list of spectrum numbers
3999	// result = semicontinuity index
4000	// ----------------------------------------------------------------------------
4001
4002	BOOLEAN semicProc3 ( leftv res,leftv u,leftv v,leftv w )
4003	{
4004	semicState state;
4005	BOOLEAN qh=(((int)(long)w->Data())==1);
4006
4007	// -----------------
4008	// check arguments
4009	// -----------------
4010
4011	lists l1 = (lists)u->Data( );
4012	lists l2 = (lists)v->Data( );
4013
4014	if( (state=list_is_spectrum( l1 ))!=semicOK )
4015	{
4016	WerrorS( "first argument is not a spectrum" );
4017	list_error( state );
4018	}
4019	else if( (state=list_is_spectrum( l2 ))!=semicOK )
4020	{
4021	WerrorS( "second argument is not a spectrum" );
4022	list_error( state );
4023	}
4024	else
4025	{
4026	spectrum s1= spectrumFromList( l1 );
4027	spectrum s2= spectrumFromList( l2 );
4028
4029	res->rtyp = INT_CMD;
4030	if (qh)
4031	res->data = (void*)(s1.mult_spectrumh( s2 ));
4032	else
4033	res->data = (void*)(s1.mult_spectrum( s2 ));
4034	}
4035
4036	// -----------------
4037	// check status
4038	// -----------------
4039
4040	return (state!=semicOK);
4041	}
4042	BOOLEAN semicProc ( leftv res,leftv u,leftv v )
4043	{
4044	sleftv tmp;
4045	memset(&tmp,0,sizeof(tmp));
4046	tmp.rtyp=INT_CMD;
4047	/* tmp.data = (void )0; -- done by memset /
4048
4049	return semicProc3(res,u,v,&tmp);
4050	}
4051
4052	#endif
4053
4054	//from mpr_inout.cc
4055	extern void nPrint(number n);
4056
4057	BOOLEAN loNewtonP( leftv res, leftv arg1 )
4058	{
4059	res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4060	return FALSE;
4061	}
4062
4063	BOOLEAN loSimplex( leftv res, leftv args )
4064	{
4065	if ( !(rField_is_long_R(currRing)) )
4066	{
4067	WerrorS("Ground field not implemented!");
4068	return TRUE;
4069	}
4070
4071	simplex * LP;
4072	matrix m;
4073
4074	leftv v= args;
4075	if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4076	return TRUE;
4077	else
4078	m= (matrix)(v->CopyD());
4079
4080	LP = new simplex(MATROWS(m),MATCOLS(m));
4081	LP->mapFromMatrix(m);
4082
4083	v= v->next;
4084	if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4085	return TRUE;
4086	else
4087	LP->m= (int)(long)(v->Data());
4088
4089	v= v->next;
4090	if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4091	return TRUE;
4092	else
4093	LP->n= (int)(long)(v->Data());
4094
4095	v= v->next;
4096	if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4097	return TRUE;
4098	else
4099	LP->m1= (int)(long)(v->Data());
4100
4101	v= v->next;
4102	if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4103	return TRUE;
4104	else
4105	LP->m2= (int)(long)(v->Data());
4106
4107	v= v->next;
4108	if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4109	return TRUE;
4110	else
4111	LP->m3= (int)(long)(v->Data());
4112
4113	#ifdef mprDEBUG_PROT
4114	Print("m (constraints) %d\n",LP->m);
4115	Print("n (columns) %d\n",LP->n);
4116	Print("m1 (<=) %d\n",LP->m1);
4117	Print("m2 (>=) %d\n",LP->m2);
4118	Print("m3 (==) %d\n",LP->m3);
4119	#endif
4120
4121	LP->compute();
4122
4123	lists lres= (lists)omAlloc( sizeof(slists) );
4124	lres->Init( 6 );
4125
4126	lres->m[0].rtyp= MATRIX_CMD; // output matrix
4127	lres->m[0].data=(void*)LP->mapToMatrix(m);
4128
4129	lres->m[1].rtyp= INT_CMD; // found a solution?
4130	lres->m[1].data=(void*)LP->icase;
4131
4132	lres->m[2].rtyp= INTVEC_CMD;
4133	lres->m[2].data=(void*)LP->posvToIV();
4134
4135	lres->m[3].rtyp= INTVEC_CMD;
4136	lres->m[3].data=(void*)LP->zrovToIV();
4137
4138	lres->m[4].rtyp= INT_CMD;
4139	lres->m[4].data=(void*)LP->m;
4140
4141	lres->m[5].rtyp= INT_CMD;
4142	lres->m[5].data=(void*)LP->n;
4143
4144	res->data= (void*)lres;
4145
4146	return FALSE;
4147	}
4148
4149	BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
4150	{
4151	ideal gls = (ideal)(arg1->Data());
4152	int imtype= (int)(long)arg2->Data();
4153
4154	uResultant::resMatType mtype= determineMType( imtype );
4155
4156	// check input ideal ( = polynomial system )
4157	if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4158	{
4159	return TRUE;
4160	}
4161
4162	uResultant *resMat= new uResultant( gls, mtype, false );
4163	if (resMat!=NULL)
4164	{
4165	res->rtyp = MODUL_CMD;
4166	res->data= (void*)resMat->accessResMat()->getMatrix();
4167	if (!errorreported) delete resMat;
4168	}
4169	return errorreported;
4170	}
4171
4172	BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
4173	{
4174
4175	poly gls;
4176	gls= (poly)(arg1->Data());
4177	int howclean= (int)(long)arg3->Data();
4178
4179	if ( !(rField_is_R(currRing) \|\|
4180	rField_is_Q(currRing) \|\|
4181	rField_is_long_R(currRing) \|\|
4182	rField_is_long_C(currRing)) )
4183	{
4184	WerrorS("Ground field not implemented!");
4185	return TRUE;
4186	}
4187
4188	if ( !(rField_is_R(currRing) \|\| rField_is_long_R(currRing) \|\| \
4189	rField_is_long_C(currRing)) )
4190	{
4191	unsigned long int ii = (unsigned long int)arg2->Data();
4192	setGMPFloatDigits( ii, ii );
4193	}
4194
4195	if ( gls == NULL \|\| pIsConstant( gls ) )
4196	{
4197	WerrorS("Input polynomial is constant!");
4198	return TRUE;
4199	}
4200
4201	int ldummy;
4202	int deg= currRing->pLDeg( gls, &ldummy, currRing );
4203	// int deg= pDeg( gls );
4204	int len= pLength( gls );
4205	int i,vpos=0;
4206	poly piter;
4207	lists elist;
4208	lists rlist;
4209
4210	elist= (lists)omAlloc( sizeof(slists) );
4211	elist->Init( 0 );
4212
4213	if ( rVar(currRing) > 1 )
4214	{
4215	piter= gls;
4216	for ( i= 1; i <= rVar(currRing); i++ )
4217	if ( pGetExp( piter, i ) )
4218	{
4219	vpos= i;
4220	break;
4221	}
4222	while ( piter )
4223	{
4224	for ( i= 1; i <= rVar(currRing); i++ )
4225	if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4226	{
4227	WerrorS("The input polynomial must be univariate!");
4228	return TRUE;
4229	}
4230	pIter( piter );
4231	}
4232	}
4233
4234	rootContainer * roots= new rootContainer();
4235	number * pcoeffs= (number )omAlloc( (deg+1) sizeof( number ) );
4236	piter= gls;
4237	for ( i= deg; i >= 0; i-- )
4238	{
4239	//if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4240	if ( piter && pTotaldegree(piter) == i )
4241	{
4242	pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4243	//nPrint( pcoeffs[i] );PrintS(" ");
4244	pIter( piter );
4245	}
4246	else
4247	{
4248	pcoeffs[i]= nInit(0);
4249	}
4250	}
4251
4252	#ifdef mprDEBUG_PROT
4253	for (i=deg; i >= 0; i--)
4254	{
4255	nPrint( pcoeffs[i] );PrintS(" ");
4256	}
4257	PrintLn();
4258	#endif
4259
4260	roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4261	roots->solver( howclean );
4262
4263	int elem= roots->getAnzRoots();
4264	char *out;
4265	char *dummy;
4266	int j;
4267
4268	rlist= (lists)omAlloc( sizeof(slists) );
4269	rlist->Init( elem );
4270
4271	if (rField_is_long_C(currRing))
4272	{
4273	for ( j= 0; j < elem; j++ )
4274	{
4275	rlist->m[j].rtyp=NUMBER_CMD;
4276	rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4277	//rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4278	}
4279	}
4280	else
4281	{
4282	for ( j= 0; j < elem; j++ )
4283	{
4284	dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4285	rlist->m[j].rtyp=STRING_CMD;
4286	rlist->m[j].data=(void *)dummy;
4287	}
4288	}
4289
4290	elist->Clean();
4291	//omFreeSize( (ADDRESS) elist, sizeof(slists) );
4292
4293	// this is (via fillContainer) the same data as in root
4294	//for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4295	//omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4296
4297	delete roots;
4298
4299	res->rtyp= LIST_CMD;
4300	res->data= (void*)rlist;
4301
4302	return FALSE;
4303	}
4304
4305	BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
4306	{
4307	int i;
4308	ideal p,w;
4309	p= (ideal)arg1->Data();
4310	w= (ideal)arg2->Data();
4311
4312	// w[0] = f(p^0)
4313	// w[1] = f(p^1)
4314	// ...
4315	// p can be a vector of numbers (multivariate polynom)
4316	// or one number (univariate polynom)
4317	// tdg = deg(f)
4318
4319	int n= IDELEMS( p );
4320	int m= IDELEMS( w );
4321	int tdg= (int)(long)arg3->Data();
4322
4323	res->data= (void*)NULL;
4324
4325	// check the input
4326	if ( tdg < 1 )
4327	{
4328	WerrorS("Last input parameter must be > 0!");
4329	return TRUE;
4330	}
4331	if ( n != rVar(currRing) )
4332	{
4333	Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4334	return TRUE;
4335	}
4336	if ( m != (int)pow((double)tdg+1,(double)n) )
4337	{
4338	Werror("Size of second input ideal must be equal to %d!",
4339	(int)pow((double)tdg+1,(double)n));
4340	return TRUE;
4341	}
4342	if ( !(rField_is_Q(currRing) /* \|\|
4343	rField_is_R() \|\| rField_is_long_R() \|\|
4344	rField_is_long_C()*/ ) )
4345	{
4346	WerrorS("Ground field not implemented!");
4347	return TRUE;
4348	}
4349
4350	number tmp;
4351	number pevpoint= (number )omAlloc( n * sizeof( number ) );
4352	for ( i= 0; i < n; i++ )
4353	{
4354	pevpoint[i]=nInit(0);
4355	if ( (p->m)[i] )
4356	{
4357	tmp = pGetCoeff( (p->m)[i] );
4358	if ( nIsZero(tmp) \|\| nIsOne(tmp) \|\| nIsMOne(tmp) )
4359	{
4360	omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4361	WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4362	return TRUE;
4363	}
4364	} else tmp= NULL;
4365	if ( !nIsZero(tmp) )
4366	{
4367	if ( !pIsConstant((p->m)[i]))
4368	{
4369	omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4370	WerrorS("Elements of first input ideal must be numbers!");
4371	return TRUE;
4372	}
4373	pevpoint[i]= nCopy( tmp );
4374	}
4375	}
4376
4377	number wresults= (number )omAlloc( m * sizeof( number ) );
4378	for ( i= 0; i < m; i++ )
4379	{
4380	wresults[i]= nInit(0);
4381	if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4382	{
4383	if ( !pIsConstant((w->m)[i]))
4384	{
4385	omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4386	omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4387	WerrorS("Elements of second input ideal must be numbers!");
4388	return TRUE;
4389	}
4390	wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4391	}
4392	}
4393
4394	vandermonde vm( m, n, tdg, pevpoint, FALSE );
4395	number *ncpoly= vm.interpolateDense( wresults );
4396	// do not free ncpoly[]!!
4397	poly rpoly= vm.numvec2poly( ncpoly );
4398
4399	omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4400	omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4401
4402	res->data= (void*)rpoly;
4403	return FALSE;
4404	}
4405
4406	BOOLEAN nuUResSolve( leftv res, leftv args )
4407	{
4408	leftv v= args;
4409
4410	ideal gls;
4411	int imtype;
4412	int howclean;
4413
4414	// get ideal
4415	if ( v->Typ() != IDEAL_CMD )
4416	return TRUE;
4417	else gls= (ideal)(v->Data());
4418	v= v->next;
4419
4420	// get resultant matrix type to use (0,1)
4421	if ( v->Typ() != INT_CMD )
4422	return TRUE;
4423	else imtype= (int)(long)v->Data();
4424	v= v->next;
4425
4426	if (imtype==0)
4427	{
4428	ideal test_id=idInit(1,1);
4429	int j;
4430	for(j=IDELEMS(gls)-1;j>=0;j--)
4431	{
4432	if (gls->m[j]!=NULL)
4433	{
4434	test_id->m[0]=gls->m[j];
4435	intvec *dummy_w=id_QHomWeight(test_id, currRing);
4436	if (dummy_w!=NULL)
4437	{
4438	WerrorS("Newton polytope not of expected dimension");
4439	delete dummy_w;
4440	return TRUE;
4441	}
4442	}
4443	}
4444	}
4445
4446	// get and set precision in digits ( > 0 )
4447	if ( v->Typ() != INT_CMD )
4448	return TRUE;
4449	else if ( !(rField_is_R(currRing) \|\| rField_is_long_R(currRing) \|\| \
4450	rField_is_long_C(currRing)) )
4451	{
4452	unsigned long int ii=(unsigned long int)v->Data();
4453	setGMPFloatDigits( ii, ii );
4454	}
4455	v= v->next;
4456
4457	// get interpolation steps (0,1,2)
4458	if ( v->Typ() != INT_CMD )
4459	return TRUE;
4460	else howclean= (int)(long)v->Data();
4461
4462	uResultant::resMatType mtype= determineMType( imtype );
4463	int i,c,count;
4464	lists listofroots= NULL;
4465	lists emptylist;
4466	number smv= NULL;
4467	BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4468
4469	//emptylist= (lists)omAlloc( sizeof(slists) );
4470	//emptylist->Init( 0 );
4471
4472	//res->rtyp = LIST_CMD;
4473	//res->data= (void *)emptylist;
4474
4475	// check input ideal ( = polynomial system )
4476	if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4477	{
4478	return TRUE;
4479	}
4480
4481	uResultant * ures;
4482	rootContainer ** iproots;
4483	rootContainer ** muiproots;
4484	rootArranger * arranger;
4485
4486	// main task 1: setup of resultant matrix
4487	ures= new uResultant( gls, mtype );
4488	if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4489	{
4490	WerrorS("Error occurred during matrix setup!");
4491	return TRUE;
4492	}
4493
4494	// if dense resultant, check if minor nonsingular
4495	if ( mtype == uResultant::denseResMat )
4496	{
4497	smv= ures->accessResMat()->getSubDet();
4498	#ifdef mprDEBUG_PROT
4499	PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4500	#endif
4501	if ( nIsZero(smv) )
4502	{
4503	WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4504	return TRUE;
4505	}
4506	}
4507
4508	// main task 2: Interpolate specialized resultant polynomials
4509	if ( interpolate_det )
4510	iproots= ures->interpolateDenseSP( false, smv );
4511	else
4512	iproots= ures->specializeInU( false, smv );
4513
4514	// main task 3: Interpolate specialized resultant polynomials
4515	if ( interpolate_det )
4516	muiproots= ures->interpolateDenseSP( true, smv );
4517	else
4518	muiproots= ures->specializeInU( true, smv );
4519
4520	#ifdef mprDEBUG_PROT
4521	c= iproots[0]->getAnzElems();
4522	for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4523	c= muiproots[0]->getAnzElems();
4524	for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4525	#endif
4526
4527	// main task 4: Compute roots of specialized polys and match them up
4528	arranger= new rootArranger( iproots, muiproots, howclean );
4529	arranger->solve_all();
4530
4531	// get list of roots
4532	if ( arranger->success() )
4533	{
4534	arranger->arrange();
4535	listofroots= listOfRoots(arranger, gmp_output_digits );
4536	}
4537	else
4538	{
4539	WerrorS("Solver was unable to find any roots!");
4540	return TRUE;
4541	}
4542
4543	// free everything
4544	count= iproots[0]->getAnzElems();
4545	for (i=0; i < count; i++) delete iproots[i];
4546	omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4547	count= muiproots[0]->getAnzElems();
4548	for (i=0; i < count; i++) delete muiproots[i];
4549	omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4550
4551	delete ures;
4552	delete arranger;
4553	nDelete( &smv );
4554
4555	res->data= (void *)listofroots;
4556
4557	//emptylist->Clean();
4558	// omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4559
4560	return FALSE;
4561	}
4562
4563	// from mpr_numeric.cc
4564	lists listOfRoots( rootArranger* self, const unsigned int oprec )
4565	{
4566	int i,j,tr;
4567	int count= self->roots[0]->getAnzRoots(); // number of roots
4568	int elem= self->roots[0]->getAnzElems(); // number of koordinates per root
4569
4570	lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
4571
4572	if ( self->found_roots )
4573	{
4574	listofroots->Init( count );
4575
4576	for (i=0; i < count; i++)
4577	{
4578	lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
4579	onepoint->Init(elem);
4580	for ( j= 0; j < elem; j++ )
4581	{
4582	if ( !rField_is_long_C(currRing) )
4583	{
4584	onepoint->m[j].rtyp=STRING_CMD;
4585	onepoint->m[j].data=(void )complexToStr((self->roots[j])[i],oprec, currRing->cf);
4586	}
4587	else
4588	{
4589	onepoint->m[j].rtyp=NUMBER_CMD;
4590	onepoint->m[j].data=(void *)n_Copy((number)(self->roots[j]->getRoot(i)), currRing->cf);
4591	}
4592	onepoint->m[j].next= NULL;
4593	onepoint->m[j].name= NULL;
4594	}
4595	listofroots->m[i].rtyp=LIST_CMD;
4596	listofroots->m[i].data=(void *)onepoint;
4597	listofroots->m[j].next= NULL;
4598	listofroots->m[j].name= NULL;
4599	}
4600
4601	}
4602	else
4603	{
4604	listofroots->Init( 0 );
4605	}
4606
4607	return listofroots;
4608	}
4609
4610	// from ring.cc
4611	void rSetHdl(idhdl h)
4612	{
4613	int i;
4614	ring rg = NULL;
4615	if (h!=NULL)
4616	{
4617	// Print(" new ring:%s (l:%d)\n",IDID(h),IDLEV(h));
4618	rg = IDRING(h);
4619	if (rg==NULL) return; //id <>NULL, ring==NULL
4620	omCheckAddrSize((ADDRESS)h,sizeof(idrec));
4621	if (IDID(h)) // OB: ????
4622	omCheckAddr((ADDRESS)IDID(h));
4623	rTest(rg);
4624	}
4625
4626	// clean up history
4627	if (sLastPrinted.RingDependend())
4628	{
4629	sLastPrinted.CleanUp();
4630	memset(&sLastPrinted,0,sizeof(sleftv));
4631	}
4632
4633	// test for valid "currRing":
4634	if ((rg!=NULL) && (rg->idroot==NULL))
4635	{
4636	ring old=rg;
4637	rg=rAssure_HasComp(rg);
4638	if (old!=rg)
4639	{
4640	rKill(old);
4641	IDRING(h)=rg;
4642	}
4643	}
4644	/------------ change the global ring -----------------------/
4645	rChangeCurrRing(rg);
4646	currRingHdl = h;
4647	}
4648
4649	BOOLEAN rSleftvOrdering2Ordering(sleftv *ord, ring R)
4650	{
4651	int last = 0, o=0, n = 1, i=0, typ = 1, j;
4652	sleftv *sl = ord;
4653
4654	// determine nBlocks
4655	while (sl!=NULL)
4656	{
4657	intvec iv = (intvec )(sl->data);
4658	if (((iv)[1]==ringorder_c)\|\|((iv)[1]==ringorder_C))
4659	i++;
4660	else if ((*iv)[1]==ringorder_L)
4661	{
4662	R->bitmask=(*iv)[2];
4663	n--;
4664	}
4665	else if (((*iv)[1]!=ringorder_a)
4666	&& ((*iv)[1]!=ringorder_a64))
4667	o++;
4668	n++;
4669	sl=sl->next;
4670	}
4671	// check whether at least one real ordering
4672	if (o==0)
4673	{
4674	WerrorS("invalid combination of orderings");
4675	return TRUE;
4676	}
4677	// if no c/C ordering is given, increment n
4678	if (i==0) n++;
4679	else if (i != 1)
4680	{
4681	// throw error if more than one is given
4682	WerrorS("more than one ordering c/C specified");
4683	return TRUE;
4684	}
4685
4686	// initialize fields of R
4687	R->order=(int )omAlloc0(nsizeof(int));
4688	R->block0=(int )omAlloc0(nsizeof(int));
4689	R->block1=(int )omAlloc0(nsizeof(int));
4690	R->wvhdl=(int*)omAlloc0(nsizeof(int_ptr));
4691
4692	int weights=(int)omAlloc0((R->N+1)*sizeof(int));
4693
4694	// init order, so that rBlocks works correctly
4695	for (j=0; j < n-1; j++)
4696	R->order[j] = (int) ringorder_unspec;
4697	// set last _C order, if no c/C order was given
4698	if (i == 0) R->order[n-2] = ringorder_C;
4699
4700	/* init orders */
4701	sl=ord;
4702	n=-1;
4703	while (sl!=NULL)
4704	{
4705	intvec *iv;
4706	iv = (intvec *)(sl->data);
4707	if ((*iv)[1]!=ringorder_L)
4708	{
4709	n++;
4710
4711	/* the format of an ordering:
4712	* iv[0]: factor
4713	* iv[1]: ordering
4714	* iv[2..end]: weights
4715	*/
4716	R->order[n] = (*iv)[1];
4717	typ=1;
4718	switch ((*iv)[1])
4719	{
4720	case ringorder_ws:
4721	case ringorder_Ws:
4722	typ=-1;
4723	case ringorder_wp:
4724	case ringorder_Wp:
4725	R->wvhdl[n]=(int)omAlloc((iv->length()-1)sizeof(int));
4726	R->block0[n] = last+1;
4727	for (i=2; i<iv->length(); i++)
4728	{
4729	R->wvhdl[n][i-2] = (*iv)[i];
4730	last++;
4731	if (weights[last]==0) weights[last]=(iv)[i]typ;
4732	}
4733	R->block1[n] = last;
4734	break;
4735	case ringorder_ls:
4736	case ringorder_ds:
4737	case ringorder_Ds:
4738	case ringorder_rs:
4739	typ=-1;
4740	case ringorder_lp:
4741	case ringorder_dp:
4742	case ringorder_Dp:
4743	case ringorder_rp:
4744	R->block0[n] = last+1;
4745	if (iv->length() == 3) last+=(*iv)[2];
4746	else last += (*iv)[0];
4747	R->block1[n] = last;
4748	//if ((R->block0[n]>R->block1[n])
4749	//\|\| (R->block1[n]>rVar(R)))
4750	//{
4751	// R->block1[n]=rVar(R);
4752	// //WerrorS("ordering larger than number of variables");
4753	// break;
4754	//}
4755	if (rCheckIV(iv)) return TRUE;
4756	for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4757	{
4758	if (weights[i]==0) weights[i]=typ;
4759	}
4760	break;
4761
4762	case ringorder_s: // no 'rank' params!
4763	{
4764
4765	if(iv->length() > 3)
4766	return TRUE;
4767
4768	if(iv->length() == 3)
4769	{
4770	const int s = (*iv)[2];
4771	R->block0[n] = s;
4772	R->block1[n] = s;
4773	}
4774	break;
4775	}
4776	case ringorder_IS:
4777	{
4778	if(iv->length() != 3) return TRUE;
4779
4780	const int s = (*iv)[2];
4781
4782	if( 1 < s \|\| s < -1 ) return TRUE;
4783
4784	R->block0[n] = s;
4785	R->block1[n] = s;
4786	break;
4787	}
4788	case ringorder_S:
4789	case ringorder_c:
4790	case ringorder_C:
4791	{
4792	if (rCheckIV(iv)) return TRUE;
4793	break;
4794	}
4795	case ringorder_aa:
4796	case ringorder_a:
4797	{
4798	R->block0[n] = last+1;
4799	R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4800	R->wvhdl[n] = (int)omAlloc((iv->length()-1)sizeof(int));
4801	for (i=2; i<iv->length(); i++)
4802	{
4803	R->wvhdl[n][i-2]=(*iv)[i];
4804	last++;
4805	if (weights[last]==0) weights[last]=(iv)[i]typ;
4806	}
4807	last=R->block0[n]-1;
4808	break;
4809	}
4810	case ringorder_a64:
4811	{
4812	R->block0[n] = last+1;
4813	R->block1[n] = si_min(last+iv->length()-2 , rVar(R));
4814	R->wvhdl[n] = (int)omAlloc((iv->length()-1)sizeof(int64));
4815	int64 w=(int64 )R->wvhdl[n];
4816	for (i=2; i<iv->length(); i++)
4817	{
4818	w[i-2]=(*iv)[i];
4819	last++;
4820	if (weights[last]==0) weights[last]=(iv)[i]typ;
4821	}
4822	last=R->block0[n]-1;
4823	break;
4824	}
4825	case ringorder_M:
4826	{
4827	int Mtyp=rTypeOfMatrixOrder(iv);
4828	if (Mtyp==0) return TRUE;
4829	if (Mtyp==-1) typ = -1;
4830
4831	R->wvhdl[n] =( int )omAlloc((iv->length()-1)sizeof(int));
4832	for (i=2; i<iv->length();i++)
4833	R->wvhdl[n][i-2]=(*iv)[i];
4834
4835	R->block0[n] = last+1;
4836	last += (int)sqrt((double)(iv->length()-2));
4837	R->block1[n] = last;
4838	for(i=si_min(rVar(R),R->block1[n]);i>=R->block0[n];i--)
4839	{
4840	if (weights[i]==0) weights[i]=typ;
4841	}
4842	break;
4843	}
4844
4845	case ringorder_no:
4846	R->order[n] = ringorder_unspec;
4847	return TRUE;
4848
4849	default:
4850	Werror("Internal Error: Unknown ordering %d", (*iv)[1]);
4851	R->order[n] = ringorder_unspec;
4852	return TRUE;
4853	}
4854	}
4855	sl=sl->next;
4856	}
4857
4858	// check for complete coverage
4859	while ( n >= 0 && (
4860	(R->order[n]==ringorder_c)
4861	\|\| (R->order[n]==ringorder_C)
4862	\|\| (R->order[n]==ringorder_s)
4863	\|\| (R->order[n]==ringorder_S)
4864	\|\| (R->order[n]==ringorder_IS)
4865	)) n--;
4866
4867	assume( n >= 0 );
4868
4869	if (R->block1[n] != R->N)
4870	{
4871	if (((R->order[n]==ringorder_dp) \|\|
4872	(R->order[n]==ringorder_ds) \|\|
4873	(R->order[n]==ringorder_Dp) \|\|
4874	(R->order[n]==ringorder_Ds) \|\|
4875	(R->order[n]==ringorder_rp) \|\|
4876	(R->order[n]==ringorder_rs) \|\|
4877	(R->order[n]==ringorder_lp) \|\|
4878	(R->order[n]==ringorder_ls))
4879	&&
4880	R->block0[n] <= R->N)
4881	{
4882	R->block1[n] = R->N;
4883	}
4884	else
4885	{
4886	Werror("mismatch of number of vars (%d) and ordering (%d vars)",
4887	R->N,R->block1[n]);
4888	return TRUE;
4889	}
4890	}
4891	// find OrdSgn:
4892	R->OrdSgn = 1;
4893	for(i=1;i<=R->N;i++)
4894	{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
4895	omFree(weights);
4896	return FALSE;
4897	}
4898
4899	BOOLEAN rSleftvList2StringArray(sleftv* sl, char** p)
4900	{
4901
4902	while(sl!=NULL)
4903	{
4904	if (sl->Name() == sNoName)
4905	{
4906	if (sl->Typ()==POLY_CMD)
4907	{
4908	sleftv s_sl;
4909	iiConvert(POLY_CMD,ANY_TYPE,-1,sl,&s_sl);
4910	if (s_sl.Name() != sNoName)
4911	*p = omStrDup(s_sl.Name());
4912	else
4913	*p = NULL;
4914	sl->next = s_sl.next;
4915	s_sl.next = NULL;
4916	s_sl.CleanUp();
4917	if (*p == NULL) return TRUE;
4918	}
4919	else
4920	return TRUE;
4921	}
4922	else
4923	*p = omStrDup(sl->Name());
4924	p++;
4925	sl=sl->next;
4926	}
4927	return FALSE;
4928	}
4929
4930	const short MAX_SHORT = SHRT_MAX; // (1 << (sizeof(short)*8)) - 1;
4931
4932	////////////////////
4933	//
4934	// rInit itself:
4935	//
4936	// INPUT: s: name, pn: ch & parameter (names), rv: variable (names)
4937	// ord: ordering
4938	// RETURN: currRingHdl on success
4939	// NULL on error
4940	// NOTE: * makes new ring to current ring, on success
4941	// * considers input sleftv's as read-only
4942	//idhdl rInit(char s, sleftv pn, sleftv* rv, sleftv* ord)
4943	ring rInit(sleftv* pn, sleftv* rv, sleftv* ord)
4944	{
4945	#ifdef HAVE_RINGS
4946	unsigned int ringtype = 0;
4947	int_number modBase = NULL;
4948	unsigned int modExponent = 1;
4949	#endif
4950	int float_len=0;
4951	int float_len2=0;
4952	ring R = NULL;
4953	idhdl tmp = NULL;
4954	BOOLEAN ffChar=FALSE;
4955	int typ = 1;
4956
4957	/* ch -------------------------------------------------------*/
4958	// get ch of ground field
4959	int numberOfAllocatedBlocks;
4960
4961	// allocated ring
4962	R = (ring) omAlloc0Bin(sip_sring_bin);
4963
4964	coeffs cf = NULL;
4965
4966	assume( pn != NULL );
4967	const int P = pn->listLength();
4968
4969	if (pn->Typ()==INT_CMD)
4970	{
4971	int ch = (int)(long)pn->Data();
4972
4973	/* parameter? -------------------------------------------------------*/
4974	pn = pn->next;
4975
4976	if (pn == NULL) // no params!?
4977	{
4978	if (ch!=0)
4979	{
4980	ch=IsPrime(ch);
4981	cf = nInitChar(n_Zp, (void*)(long)ch);
4982	}
4983	else
4984	cf = nInitChar(n_Q, (void*)(long)ch);
4985	}
4986	else
4987	{
4988	const int pars = pn->listLength();
4989
4990	assume( pars > 0 );
4991
4992	// predefined finite field: (p^k, a)
4993	if ((ch!=0) && (ch!=IsPrime(ch)) && (pars == 1))
4994	{
4995	GFInfo param;
4996
4997	param.GFChar = ch;
4998	param.GFDegree = 1;
4999	param.GFPar_name = pn->name;
5000
5001	cf = nInitChar(n_GF, &param);
5002	}
5003	else // (0/p, a, b, ..., z)
5004	{
5005	assume( (ch == 0) \|\| (ch==IsPrime(ch)) );
5006
5007	// if ((pars > 1) && (ffChar))
5008	// {
5009	// WerrorS("too many parameters");
5010	// goto rInitError;
5011	// }
5012
5013	char names = (char)omAlloc0(pars * sizeof(char_ptr));
5014
5015	if (rSleftvList2StringArray(pn, names))
5016	{
5017	WerrorS("parameter expected");
5018	goto rInitError;
5019	}
5020
5021	TransExtInfo extParam;
5022
5023	extParam.r = rDefault( ch, pars, names); // Q/Zp [ p_1, ... p_pars ]
5024
5025	cf = nInitChar(n_transExt, &extParam);
5026	}
5027	}
5028
5029	// if (cf==NULL) goto rInitError;
5030	assume( cf != NULL );
5031	}
5032	else if ((pn->name != NULL)
5033	&& ((strcmp(pn->name,"real")==0) \|\| (strcmp(pn->name,"complex")==0)))
5034	{
5035	BOOLEAN complex_flag=(strcmp(pn->name,"complex")==0);
5036	const int ch=0;
5037	if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5038	{
5039	float_len=(int)(long)pn->next->Data();
5040	float_len2=float_len;
5041	pn=pn->next;
5042	if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5043	{
5044	float_len2=(int)(long)pn->next->Data();
5045	pn=pn->next;
5046	}
5047	}
5048	assume( float_len <= float_len2 );
5049
5050	if( !complex_flag && (float_len2 <= (short)SHORT_REAL_LENGTH) )
5051	cf=nInitChar(n_R, NULL);
5052	else // longR or longC?
5053	{
5054	LongComplexInfo param;
5055
5056	param.float_len = float_len;
5057	param.float_len2 = float_len2;
5058
5059	// set the parameter name
5060	if (complex_flag) {
5061	if (pn->next == NULL)
5062	param.par_name=(const char*)"i"; //default to i
5063	else
5064	param.par_name = (const char*)pn->next->name;
5065	}
5066
5067	cf = nInitChar(complex_flag ? n_long_C: n_long_R, (void*)&param);
5068	}
5069	assume( cf != NULL );
5070	}
5071	#ifdef HAVE_RINGS
5072	else if ((pn->name != NULL) && (strcmp(pn->name, "integer") == 0))
5073	{
5074	modBase = (int_number) omAlloc(sizeof(mpz_t));
5075	mpz_init_set_si(modBase, 0);
5076	if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5077	{
5078	mpz_set_ui(modBase, (int)(long) pn->next->Data());
5079	pn=pn->next;
5080	if ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5081	{
5082	modExponent = (long) pn->next->Data();
5083	pn=pn->next;
5084	}
5085	while ((pn->next!=NULL) && (pn->next->Typ()==INT_CMD))
5086	{
5087	mpz_mul_ui(modBase, modBase, (int)(long) pn->next->Data());
5088	pn=pn->next;
5089	}
5090	}
5091	else
5092	cf=nInitChar(n_Z,NULL);
5093
5094	if ((mpz_cmp_ui(modBase, 1) == 0) && (mpz_cmp_ui(modBase, 0) < 0))
5095	{
5096	Werror("Wrong ground ring specification (module is 1)");
5097	goto rInitError;
5098	}
5099	if (modExponent < 1)
5100	{
5101	Werror("Wrong ground ring specification (exponent smaller than 1");
5102	goto rInitError;
5103	}
5104	// module is 0 ---> integers ringtype = 4;
5105	// we have an exponent
5106	if (modExponent > 1)
5107	{
5108	const int ch = modExponent;
5109	if ((mpz_cmp_ui(modBase, 2) == 0) && (modExponent <= 8*sizeof(NATNUMBER)))
5110	{
5111	/* this branch should be active for modExponent = 2..32 resp. 2..64,
5112	depending on the size of a long on the respective platform */
5113	ringtype = 1; // Use Z/2^ch
5114	cf=nInitChar(n_Z2m,(void*)(long)modExponent);
5115	}
5116	else
5117	{
5118	ringtype = 3;
5119	cf=nInitChar(n_Zn,(void*)(long)modBase);
5120	}
5121	}
5122	// just a module m > 1
5123	else
5124	{
5125	ringtype = 2;
5126	const int ch = mpz_get_ui(modBase);
5127	cf=nInitChar(n_Zn,(void*)(long)ch);
5128	}
5129	assume( cf != NULL );
5130	}
5131	#endif
5132	// ring NEW = OLD, (), (); where OLD is a polynomial ring...
5133	else if ((pn->Typ()==RING_CMD) && (P == 1))
5134	{
5135	TransExtInfo extParam;
5136	extParam.r = (ring)pn->Data();
5137	cf = nInitChar(n_transExt, &extParam);
5138	}
5139	else if ((pn->Typ()==QRING_CMD) && (P == 1)) // same for qrings - which should be fields!?
5140	{
5141	AlgExtInfo extParam;
5142	extParam.r = (ring)pn->Data();
5143	extParam.i = (extParam.r->qideal);
5144
5145	cf = nInitChar(n_algExt, &extParam); // Q[a]/<minideal>
5146	}
5147	else
5148	{
5149	sleftv* p = pn;
5150	Werror("Wrong or unknown ground field specification");
5151	#ifndef NDEBUG
5152	while (p != NULL)
5153	{
5154	Print( "pn[%p]: type: %d [%s]: %p, name: %s", p, p->Typ(), Tok2Cmdname(p->Typ()), p->Data(), (p->name == NULL? "NULL" : p->name) );
5155	PrintLn();
5156	p = p->next;
5157	}
5158	#endif
5159	goto rInitError;
5160	}
5161	// pn=pn->next;
5162
5163	int l, last;
5164	sleftv * sl;
5165	/every entry in the new ring is initialized to 0/
5166
5167	/* characteristic -----------------------------------------------*/
5168	/* input: 0 ch=0 : Q parameter=NULL ffChar=FALSE float_len
5169	* 0 1 : Q(a,...) *names FALSE
5170	* 0 -1 : R NULL FALSE 0
5171	* 0 -1 : R NULL FALSE prec. >6
5172	* 0 -1 : C *names FALSE prec. 0..?
5173	* p p : Fp NULL FALSE
5174	* p -p : Fp(a) *names FALSE
5175	* q q : GF(q=p^n) *names TRUE
5176	*/
5177	l = 0;
5178
5179	if (cf==NULL)
5180	{
5181	const int ch=32003;
5182	Warn("Invalid ground field specification: using the default field: Z_{%d}", ch);
5183	cf=nInitChar(n_Zp, (void*)(long)ch);
5184	}
5185
5186	assume( R != NULL );
5187
5188	R->cf = cf;
5189
5190	#ifdef HAVE_RINGS
5191	// the following should have beed set already into cf, right?!
5192	// R->cf->ringtype = ringtype;
5193	// R->cf->modBase = modBase;
5194	// R->cf->modExponent = modExponent;
5195	#endif
5196
5197	/* names and number of variables-------------------------------------*/
5198	{
5199	int l=rv->listLength();
5200
5201	if (l>MAX_SHORT)
5202	{
5203	Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5204	goto rInitError;
5205	}
5206	R->N = l; /rv->listLength();/
5207	}
5208	R->names = (char *)omAlloc0(R->N sizeof(char_ptr));
5209	if (rSleftvList2StringArray(rv, R->names))
5210	{
5211	WerrorS("name of ring variable expected");
5212	goto rInitError;
5213	}
5214
5215	/* check names and parameters for conflicts ------------------------- */
5216	rRenameVars(R); // conflicting variables will be renamed
5217	/* ordering -------------------------------------------------------------*/
5218	if (rSleftvOrdering2Ordering(ord, R))
5219	goto rInitError;
5220
5221	// Complete the initialization
5222	if (rComplete(R,1))
5223	goto rInitError;
5224
5225	#ifdef HABE_RINGS
5226	// currently, coefficients which are ring elements require a global ordering:
5227	if (rField_is_Ring(R) && (R->OrdSgn==-1))
5228	{
5229	WerrorS("global ordering required for these coefficients");
5230	goto rInitError;
5231	}
5232	#endif
5233
5234	rTest(R);
5235
5236	// try to enter the ring into the name list
5237	// need to clean up sleftv here, before this ring can be set to
5238	// new currRing or currRing can be killed beacuse new ring has
5239	// same name
5240	if (pn != NULL) pn->CleanUp();
5241	if (rv != NULL) rv->CleanUp();
5242	if (ord != NULL) ord->CleanUp();
5243	//if ((tmp = enterid(s, myynest, RING_CMD, &IDROOT))==NULL)
5244	// goto rInitError;
5245
5246	//memcpy(IDRING(tmp),R,sizeof(*R));
5247	// set current ring
5248	//omFreeBin(R, ip_sring_bin);
5249	//return tmp;
5250	return R;
5251
5252	// error case:
5253	rInitError:
5254	if ((R != NULL)&&(R->cf!=NULL)) rDelete(R);
5255	if (pn != NULL) pn->CleanUp();
5256	if (rv != NULL) rv->CleanUp();
5257	if (ord != NULL) ord->CleanUp();
5258	return NULL;
5259	}
5260
5261	ring rSubring(ring org_ring, sleftv* rv)
5262	{
5263	ring R = rCopy0(org_ring);
5264	int perm=(int )omAlloc0((org_ring->N+1)*sizeof(int));
5265	int last = 0, o=0, n = rBlocks(org_ring), i=0, typ = 1, j;
5266
5267	/* names and number of variables-------------------------------------*/
5268	{
5269	int l=rv->listLength();
5270	if (l>MAX_SHORT)
5271	{
5272	Werror("too many ring variables(%d), max is %d",l,MAX_SHORT);
5273	goto rInitError;
5274	}
5275	R->N = l; /rv->listLength();/
5276	}
5277	omFree(R->names);
5278	R->names = (char *)omAlloc0(R->N sizeof(char_ptr));
5279	if (rSleftvList2StringArray(rv, R->names))
5280	{
5281	WerrorS("name of ring variable expected");
5282	goto rInitError;
5283	}
5284
5285	/* check names for subring in org_ring ------------------------- */
5286	{
5287	i=0;
5288
5289	for(j=0;j<R->N;j++)
5290	{
5291	for(;i<org_ring->N;i++)
5292	{
5293	if (strcmp(org_ring->names[i],R->names[j])==0)
5294	{
5295	perm[i+1]=j+1;
5296	break;
5297	}
5298	}
5299	if (i>org_ring->N)
5300	{
5301	Werror("variable %d (%s) not in basering",j+1,R->names[j]);
5302	break;
5303	}
5304	}
5305	}
5306	//Print("perm=");
5307	//for(i=1;i<org_ring->N;i++) Print("v%d -> v%d\n",i,perm[i]);
5308	/* ordering -------------------------------------------------------------*/
5309
5310	for(i=0;i<n;i++)
5311	{
5312	int min_var=-1;
5313	int max_var=-1;
5314	for(j=R->block0[i];j<=R->block1[i];j++)
5315	{
5316	if (perm[j]>0)
5317	{
5318	if (min_var==-1) min_var=perm[j];
5319	max_var=perm[j];
5320	}
5321	}
5322	if (min_var!=-1)
5323	{
5324	//Print("block %d: old %d..%d, now:%d..%d\n",
5325	// i,R->block0[i],R->block1[i],min_var,max_var);
5326	R->block0[i]=min_var;
5327	R->block1[i]=max_var;
5328	if (R->wvhdl[i]!=NULL)
5329	{
5330	omFree(R->wvhdl[i]);
5331	R->wvhdl[i]=(int)omAlloc0((max_var-min_var+1)sizeof(int));
5332	for(j=org_ring->block0[i];j<=org_ring->block1[i];j++)
5333	{
5334	if (perm[j]>0)
5335	{
5336	R->wvhdl[i][perm[j]-R->block0[i]]=
5337	org_ring->wvhdl[i][j-org_ring->block0[i]];
5338	//Print("w%d=%d (orig_w%d)\n",perm[j],R->wvhdl[i][perm[j]-R->block0[i]],j);
5339	}
5340	}
5341	}
5342	}
5343	else
5344	{
5345	if(R->block0[i]>0)
5346	{
5347	//Print("skip block %d\n",i);
5348	R->order[i]=ringorder_unspec;
5349	if (R->wvhdl[i] !=NULL) omFree(R->wvhdl[i]);
5350	R->wvhdl[i]=NULL;
5351	}
5352	//else Print("keep block %d\n",i);
5353	}
5354	}
5355	i=n-1;
5356	while(i>0)
5357	{
5358	// removed unneded blocks
5359	if(R->order[i-1]==ringorder_unspec)
5360	{
5361	for(j=i;j<=n;j++)
5362	{
5363	R->order[j-1]=R->order[j];
5364	R->block0[j-1]=R->block0[j];
5365	R->block1[j-1]=R->block1[j];
5366	if (R->wvhdl[j-1] !=NULL) omFree(R->wvhdl[j-1]);
5367	R->wvhdl[j-1]=R->wvhdl[j];
5368	}
5369	R->order[n]=ringorder_unspec;
5370	n--;
5371	}
5372	i--;
5373	}
5374	n=rBlocks(org_ring)-1;
5375	while (R->order[n]==0) n--;
5376	while (R->order[n]==ringorder_unspec) n--;
5377	if ((R->order[n]==ringorder_c) \|\| (R->order[n]==ringorder_C)) n--;
5378	if (R->block1[n] != R->N)
5379	{
5380	if (((R->order[n]==ringorder_dp) \|\|
5381	(R->order[n]==ringorder_ds) \|\|
5382	(R->order[n]==ringorder_Dp) \|\|
5383	(R->order[n]==ringorder_Ds) \|\|
5384	(R->order[n]==ringorder_rp) \|\|
5385	(R->order[n]==ringorder_rs) \|\|
5386	(R->order[n]==ringorder_lp) \|\|
5387	(R->order[n]==ringorder_ls))
5388	&&
5389	R->block0[n] <= R->N)
5390	{
5391	R->block1[n] = R->N;
5392	}
5393	else
5394	{
5395	Werror("mismatch of number of vars (%d) and ordering (%d vars) in block %d",
5396	R->N,R->block1[n],n);
5397	return NULL;
5398	}
5399	}
5400	omFree(perm);
5401	// find OrdSgn:
5402	R->OrdSgn = org_ring->OrdSgn; // IMPROVE!
5403	//for(i=1;i<=R->N;i++)
5404	//{ if (weights[i]<0) { R->OrdSgn=-1;break; }}
5405	//omFree(weights);
5406	// Complete the initialization
5407	if (rComplete(R,1))
5408	goto rInitError;
5409
5410	rTest(R);
5411
5412	if (rv != NULL) rv->CleanUp();
5413
5414	return R;
5415
5416	// error case:
5417	rInitError:
5418	if (R != NULL) rDelete(R);
5419	if (rv != NULL) rv->CleanUp();
5420	return NULL;
5421	}
5422
5423	void rKill(ring r)
5424	{
5425	if ((r->ref<=0)&&(r->order!=NULL))
5426	{
5427	#ifdef RDEBUG
5428	if (traceit &TRACE_SHOW_RINGS) Print("kill ring %lx\n",(long)r);
5429	#endif
5430	if (r->qideal!=NULL)
5431	{
5432	id_Delete(&r->qideal, r);
5433	r->qideal = NULL;
5434	}
5435	int i=1;
5436	int j;
5437	int *pi=r->order;
5438	#ifdef USE_IILOCALRING
5439	for (j=0;j<iiRETURNEXPR_len;j++)
5440	{
5441	if (iiLocalRing[j]==r)
5442	{
5443	if (j<myynest) Warn("killing the basering for level %d",j);
5444	iiLocalRing[j]=NULL;
5445	}
5446	}
5447	#else /* USE_IILOCALRING */
5448	//#endif /* USE_IILOCALRING */
5449	{
5450	proclevel * nshdl = procstack;
5451	int lev=myynest-1;
5452
5453	for(; nshdl != NULL; nshdl = nshdl->next)
5454	{
5455	if (nshdl->cRing==r)
5456	{
5457	Warn("killing the basering for level %d",lev);
5458	nshdl->cRing=NULL;
5459	nshdl->cRingHdl=NULL;
5460	}
5461	}
5462	}
5463	#endif /* USE_IILOCALRING */
5464	// any variables depending on r ?
5465	while (r->idroot!=NULL)
5466	{
5467	killhdl2(r->idroot,&(r->idroot),r);
5468	}
5469	if (r==currRing)
5470	{
5471	// all dependend stuff is done, clean global vars:
5472	if (r->qideal!=NULL)
5473	{
5474	currQuotient=NULL;
5475	}
5476	if ((currRing->ppNoether)!=NULL) pDelete(&(currRing->ppNoether));
5477	if (sLastPrinted.RingDependend())
5478	{
5479	sLastPrinted.CleanUp();
5480	}
5481	if ((myynest>0) && (iiRETURNEXPR[myynest].RingDependend()))
5482	{
5483	WerrorS("return value depends on local ring variable (export missing ?)");
5484	iiRETURNEXPR[myynest].CleanUp();
5485	}
5486	currRing=NULL;
5487	currRingHdl=NULL;
5488	}
5489
5490	/* nKillChar(r); will be called from inside of rDelete */
5491	rDelete(r);
5492	return;
5493	}
5494	r->ref--;
5495	}
5496
5497	void rKill(idhdl h)
5498	{
5499	ring r = IDRING(h);
5500	int ref=0;
5501	if (r!=NULL)
5502	{
5503	ref=r->ref;
5504	rKill(r);
5505	}
5506	if (h==currRingHdl)
5507	{
5508	if (ref<=0) { currRing=NULL; currRingHdl=NULL;}
5509	else
5510	{
5511	currRingHdl=rFindHdl(r,currRingHdl,NULL);
5512	}
5513	}
5514	}
5515
5516	idhdl rSimpleFindHdl(ring r, idhdl root, idhdl n)
5517	{
5518	//idhdl next_best=NULL;
5519	idhdl h=root;
5520	while (h!=NULL)
5521	{
5522	if (((IDTYP(h)==RING_CMD)\|\|(IDTYP(h)==QRING_CMD))
5523	&& (h!=n)
5524	&& (IDRING(h)==r)
5525	)
5526	{
5527	// if (IDLEV(h)==myynest)
5528	// return h;
5529	// if ((IDLEV(h)==0) \|\| (next_best==NULL))
5530	// next_best=h;
5531	// else if (IDLEV(next_best)<IDLEV(h))
5532	// next_best=h;
5533	return h;
5534	}
5535	h=IDNEXT(h);
5536	}
5537	//return next_best;
5538	return NULL;
5539	}
5540
5541	extern BOOLEAN jjPROC(leftv res, leftv u, leftv v);
5542	ideal kGroebner(ideal F, ideal Q)
5543	{
5544	//test\|=Sy_bit(OPT_PROT);
5545	idhdl save_ringhdl=currRingHdl;
5546	ideal resid;
5547	idhdl new_ring=NULL;
5548	if ((currRingHdl==NULL) \|\| (IDRING(currRingHdl)!=currRing))
5549	{
5550	currRingHdl=enterid(omStrDup(" GROEBNERring"),0,RING_CMD,&IDROOT,FALSE);
5551	new_ring=currRingHdl;
5552	IDRING(currRingHdl)=currRing;
5553	}
5554	sleftv v; memset(&v,0,sizeof(v)); v.rtyp=IDEAL_CMD; v.data=(char *) F;
5555	idhdl h=ggetid("groebner");
5556	sleftv u; memset(&u,0,sizeof(u)); u.rtyp=IDHDL; u.data=(char *) h;
5557	u.name=IDID(h);
5558
5559	sleftv res; memset(&res,0,sizeof(res));
5560	if(jjPROC(&res,&u,&v))
5561	{
5562	resid=kStd(F,Q,testHomog,NULL);
5563	}
5564	else
5565	{
5566	//printf("typ:%d\n",res.rtyp);
5567	resid=(ideal)(res.data);
5568	}
5569	// cleanup GROEBNERring, save_ringhdl, u,v,(res )
5570	if (new_ring!=NULL)
5571	{
5572	idhdl h=IDROOT;
5573	if (h==new_ring) IDROOT=h->next;
5574	else
5575	{
5576	while ((h!=NULL) &&(h->next!=new_ring)) h=h->next;
5577	if (h!=NULL) h->next=h->next->next;
5578	}
5579	if (h!=NULL) omFreeSize(h,sizeof(*h));
5580	}
5581	currRingHdl=save_ringhdl;
5582	u.CleanUp();
5583	v.CleanUp();
5584	return resid;
5585	}
5586
5587	static void jjINT_S_TO_ID(int n,int *e, leftv res)
5588	{
5589	if (n==0) n=1;
5590	ideal l=idInit(n,1);
5591	int i;
5592	poly p;
5593	for(i=rVar(currRing);i>0;i--)
5594	{
5595	if (e[i]>0)
5596	{
5597	n--;
5598	p=pOne();
5599	pSetExp(p,i,1);
5600	pSetm(p);
5601	l->m[n]=p;
5602	if (n==0) break;
5603	}
5604	}
5605	res->data=(char*)l;
5606	setFlag(res,FLAG_STD);
5607	omFreeSize((ADDRESS)e,(rVar(currRing)+1)*sizeof(int));
5608	}
5609	BOOLEAN jjVARIABLES_P(leftv res, leftv u)
5610	{
5611	int e=(int )omAlloc0((rVar(currRing)+1)*sizeof(int));
5612	int n=pGetVariables((poly)u->Data(),e);
5613	jjINT_S_TO_ID(n,e,res);
5614	return FALSE;
5615	}
5616
5617	BOOLEAN jjVARIABLES_ID(leftv res, leftv u)
5618	{
5619	int e=(int )omAlloc0((rVar(currRing)+1)*sizeof(int));
5620	ideal I=(ideal)u->Data();
5621	int i;
5622	int n=0;
5623	for(i=I->nrows*I->ncols-1;i>=0;i--)
5624	{
5625	n=pGetVariables(I->m[i],e);
5626	}
5627	jjINT_S_TO_ID(n,e,res);
5628	return FALSE;
5629	}

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

Download in other formats: