Context Navigation

source: git/Singular/ipshell.cc @ 47a4c2

spielwiese

Last change on this file since 47a4c2 was 47a4c2, checked in by Martin Lee <martinlee84@…>, 12 years ago
moved semicState,spectrumState and spectrumPrintError to ipshell.cc adapted member functions involving lists from semic and splist in ipshell.cc moved function declarations from spectrum.h involving leftv to ipshell.h
Property mode set to `100644`
File size: 132.9 KB

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

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

Download in other formats: