Context Navigation

source: git/Singular/ipshell.cc @ 96f9be8

spielwiese

Last change on this file since 96f9be8 was 96f9be8, checked in by Hans Schönemann <hannes@…>, 15 years ago
* IsPrime tests divisibility by primes from factory: smallprimes (<=32749) If candidate is larger than 32749^2, we test now by all uneven numbers in 32749..sqrt(cand) git-svn-id: file:///usr/local/Singular/svn/trunk@12096 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 125.0 KB

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

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

Download in other formats: