Context Navigation

source: git/Singular/ipshell.cc @ 24a9587

spielwiese

Last change on this file since 24a9587 was 06c0b3, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: newstruct changes fropm master
Property mode set to `100644`
File size: 138.1 KB

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

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

Download in other formats: