Context Navigation

source: git/Singular/ipshell.cc @ d3e630

spielwiese

Last change on this file since d3e630 was d3e630, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: code cleanup (see ChangeLog) git-svn-id: file:///usr/local/Singular/svn/trunk@12029 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 124.8 KB

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

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

Download in other formats: