Context Navigation

brnoeth.lib @ 224420b

spielwiese

Last change on this file since 224420b was 224420b, checked in by Frank Seelisch <seelisch@…>, 14 years ago
removed duplicate proc declarations when doing LIB "all.lib"; i.e., no more "redefine" prompts git-svn-id: file:///usr/local/Singular/svn/trunk@12529 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 139.9 KB

Line
1	version="$Id$";
2	category="Coding theory";
3	info="
4	LIBRARY: brnoeth.lib Brill-Noether Algorithm, Weierstrass-SG and AG-codes
5	AUTHORS: Jose Ignacio Farran Martin, ignfar@eis.uva.es
6	Christoph Lossen, lossen@mathematik.uni-kl.de
7
8	OVERVIEW:
9	Implementation of the Brill-Noether algorithm for solving the
10	Riemann-Roch problem and applications to Algebraic Geometry codes.
11	The computation of Weierstrass semigroups is also implemented.@*
12	The procedures are intended only for plane (singular) curves defined over
13	a prime field of positive characteristic.@*
14	For more information about the library see the end of the file brnoeth.lib.
15
16	MAIN PROCEDURES:
17	Adj_div(f) computes the conductor of a curve
18	NSplaces(h,A) computes non-singular places with given degrees
19	BrillNoether(D,C) computes a vector space basis of the linear system L(D)
20	Weierstrass(P,m,C) computes the Weierstrass semigroup of C at P up to m
21	extcurve(d,C) extends the curve C to an extension of degree d
22	AGcode_L(G,D,E) computes the evaluation AG code with divisors G and D
23	AGcode_Omega(G,D,E) computes the residual AG code with divisors G and D
24	prepSV(G,D,F,E) preprocessing for the basic decoding algorithm
25	decodeSV(y,K) decoding of a word with the basic decoding algorithm
26
27	AUXILIARY PROCEDURES:
28	closed_points(I) computes the zero-set of a zero-dim. ideal in 2 vars
29	dual_code(C) computes the dual code
30	sys_code(C) computes an equivalent systematic code
31	permute_L(L,P) applies a permutation to a list
32
33	KEYWORDS: Weierstrass semigroup; Algebraic Geometry codes;
34	Brill-Noether algorithm
35	";
36
37	LIB "matrix.lib";
38	LIB "triang.lib";
39	LIB "hnoether.lib";
40	LIB "inout.lib";
41
42	///////////////////////////////////////////////////////////////////////////////
43
44	// **********************************************************
45	// * POINTS, PLACES AND DIVISORS OF (SINGULAR) PLANE CURVES *
46	// **********************************************************
47
48	proc closed_points (ideal I)
49	"USAGE: closed_points(I); I an ideal
50	RETURN: list of prime ideals (each a Groebner basis), corresponding to
51	the (distinct affine closed) points of V(I)
52	NOTE: The ideal must have dimension 0, the basering must have 2
53	variables, the ordering must be lp, and the base field must
54	be finite and prime.@*
55	It might be convenient to set the option(redSB) in advance.
56	SEE ALSO: triang_lib
57	EXAMPLE: example closed_points; shows an example
58	"
59	{
60	ideal II=std(I);
61	if (II==1)
62	{
63	return(list());
64	}
65	list TL=triangMH(II);
66	int s=size(TL);
67	list L=list();
68	int i,j,k;
69	ideal Facts;
70	poly G2;
71	for (i=1;i<=s;i=i+1)
72	{
73	Facts=factorize(TL[i][1],1);
74	k=size(Facts);
75	G2=TL[i][2];
76	for (j=1;j<=k;j=j+1)
77	{
78	L=L+pd2(Facts[j],G2);
79	}
80	}
81	// eliminate possible repetitions
82	s=size(L);
83	list LP=list();
84	LP[1]=std(L[1]);
85	int counter=1;
86	for (i=2;i<=s;i=i+1)
87	{
88	if (isPinL(L[i],LP)==0)
89	{
90	counter=counter+1;
91	LP[counter]=std(L[i]);
92	}
93	}
94	return(LP);
95	}
96	example
97	{
98	"EXAMPLE:"; echo = 2;
99	ring s=2,(x,y),lp;
100	// this is just the affine plane over F_4 :
101	ideal I=x4+x,y4+y;
102	list L=closed_points(I);
103	// and here you have all the points :
104	L;
105	}
106	///////////////////////////////////////////////////////////////////////////////
107
108	static proc pd2 (poly g1,poly g2)
109	{
110	// If g1,g2 is a std. resp. lex. in (x,y) then the procedure
111	// factorizes g2 in the "extension given by g1"
112	// (then g1 must be irreducible) and returns a list of
113	// ideals with always g1 as first component and the
114	// distinct factors of g2 as second components
115	list L=list();
116	ideal J=g1;
117	int i,s;
118	if (deg(g1)==1)
119	{
120	poly A=-subst(g1,var(2),0);
121	poly B=subst(g2,var(2),A);
122	ideal facts=factorize(B,1);
123	s=size(facts);
124	for (i=1;i<=s;i=i+1)
125	{
126	J[2]=facts[i];
127	L[i]=J;
128	}
129	}
130	else
131	{
132	def BR=basering;
133	poly A=g1;
134	poly B=g2;
135	ring raux1=char(basering),(x,y,@a),lp;
136	poly G;
137	ring raux2=(char(basering),@a),(x,y),lp;
138	map psi=BR,x,@a;
139	minpoly=number(psi(A));
140	poly f=psi(B);
141	ideal facts=factorize(f,1);
142	s=size(facts);
143	poly g;
144	string sg;
145	for (i=1;i<=s;i=i+1)
146	{
147	g=facts[i];
148	sg=string(g);
149	setring raux1;
150	execute("G="+sg+";");
151	G=subst(G,@a,y);
152	setring BR;
153	map ppssii=raux1,var(1),var(2),0;
154	J[2]=ppssii(G);
155	L[i]=J;
156	kill ppssii;
157	setring raux2;
158	}
159	setring BR;
160	kill raux1,raux2;
161	}
162	return(L);
163	}
164	///////////////////////////////////////////////////////////////////////////////
165	static proc isPinL (ideal P,list L)
166	{
167	// checks if a (plane) point P is in a list of (plane) points L
168	// by just comparing generators
169	// it works only if all (prime) ideals are given in a "canonical way",
170	// namely:
171	// the first generator is monic and irreducible,
172	// and depends only on the second variable,
173	// and the second one is monic in the first variable
174	// and irreducible over the field extension determined by
175	// the second variable and the first generator as minpoly
176	int s=size(L);
177	int i;
178	for (i=1;i<=s;i=i+1)
179	{
180	if ( P[1]==L[i][1] && P[2]==L[i][2] )
181	{
182	return(1);
183	}
184	}
185	return(0);
186	}
187	///////////////////////////////////////////////////////////////////////////////
188	static proc s_locus (poly f)
189	{
190	// computes : ideal of affine singular locus
191	// the equation f must be affine
192	// warning : if there is an error message then the output is "none"
193	// option(redSB) is convenient to be set in advance
194	ideal I=f,jacob(f);
195	I=std(I);
196	if (dim(I)>0)
197	{
198	// dimension check (has to be 0)
199	ERROR("something was wrong; possibly non-reduced curve");
200	}
201	else
202	{
203	return(I);
204	}
205	}
206	///////////////////////////////////////////////////////////////////////////////
207	static proc curve (poly f)
208	"USAGE: curve(f), where f is a polynomial (affine or projective)
209	CREATE: poly CHI in both rings aff_r=p,(x,y),lp and Proj_R=p,(x,y,z),lp
210	also ideal (std) Aff_SLocus of affine singular locus in the ring
211	aff_r
212	RETURN: list (size 3) with two rings aff_r,Proj_R and an integer deg(f)
213	NOTE: f must be absolutely irreducible, but this is not checked
214	it is not implemented yet for extensions of prime fields
215	"
216	{
217	def base_r=basering;
218	ring aff_r=char(basering),(x,y),lp;
219	ring Proj_R=char(basering),(x,y,z),lp;
220	setring base_r;
221	int degX=deg(f);
222	if (nvars(basering)==2)
223	{
224	setring aff_r;
225	map embpol=base_r,x,y;
226	poly CHI=embpol(f);
227	export(CHI);
228	kill embpol;
229	ideal Aff_SLocus=s_locus(CHI);
230	export(Aff_SLocus);
231	setring Proj_R;
232	poly CHI=homog(imap(aff_r,CHI),z);
233	export(CHI);
234	setring base_r;
235	list L=list();
236	L[1]=aff_r;
237	L[2]=Proj_R;
238	L[3]=degX;
239	kill aff_r,Proj_R;
240	return(L);
241	}
242	if (nvars(basering)==3)
243	{
244	setring Proj_R;
245	map embpol=base_r,x,y,z;
246	poly CHI=embpol(f);
247	export(CHI);
248	kill embpol;
249	string s=string(subst(CHI,z,1));
250	setring aff_r;
251	execute("poly CHI="+s+";");
252	export(CHI);
253	ideal Aff_SLocus=s_locus(CHI);
254	export(Aff_SLocus);
255	setring base_r;
256	list L=list();
257	L[1]=aff_r;
258	L[2]=Proj_R;
259	L[3]=degX;
260	kill aff_r,Proj_R;
261	return(L);
262	}
263	ERROR("basering must have 2 or 3 variables");
264	}
265	///////////////////////////////////////////////////////////////////////////////
266	static proc Aff_SL (ideal ISL)
267	{
268	// computes : affine singular (closed) points as a list of lists of
269	// prime ideals and intvec (for storing the places over each point)
270	// the ideal ISL=s_locus(CHI) is assumed to be computed in advance for
271	// a plane curve CHI, and it must be given by a standard basis
272	// for our purpose the function must called with the "global" ideal
273	// "Aff_SLocus"
274	list SL=list();
275	ideal I=ISL;
276	if ( I[1] != 1 )
277	{
278	list L=list();
279	ideal aux;
280	intvec iv;
281	int i,s;
282	L=closed_points(I);
283	s=size(L);
284	for (i=1;i<=s;i=i+1)
285	{
286	aux=std(L[i]);
287	SL[i]=list(aux,iv);
288	}
289	}
290	return(SL);
291	}
292	///////////////////////////////////////////////////////////////////////////////
293	static proc inf_P (poly f)
294	{
295	// computes : all (closed) points at infinity as homogeneous polynomials
296	// output : two lists with respectively singular and non-singular points
297	intvec iv;
298	def base_r=basering;
299	ring r_auxz=char(basering),(x,y,z),lp;
300	poly f=imap(base_r,f);
301	poly F=homog(f,z); // equation of projective curve
302	poly f_inf=subst(F,z,0);
303	setring base_r;
304	poly f_inf=imap(r_auxz,f_inf);
305	ideal I=factorize(f_inf,1); // points at infinity as homogeneous
306	// polynomials
307	int s=size(I);
308	int i;
309	list IP_S=list(); // for singular points at infinity
310	list IP_NS=list(); // for non-singular points at infinity
311	int counter_S;
312	int counter_NS;
313	poly aux;
314	for (i=1;i<=s;i=i+1)
315	{
316	aux=subst(I[i],y,1);
317	if (aux==1)
318	{
319	// the point is (1:0:0)
320	setring r_auxz;
321	poly f_yz=subst(F,x,1);
322	if ( subst(subst(diff(f_yz,y),y,0),z,0)==0 &&
323	subst(subst(diff(f_yz,z),y,0),z,0)==0 )
324	{
325	// the point is singular
326	counter_S=counter_S+1;
327	kill f_yz;
328	setring base_r;
329	IP_S[counter_S]=list(I[i],iv);
330	}
331	else
332	{
333	// the point is non-singular
334	counter_NS=counter_NS+1;
335	kill f_yz;
336	setring base_r;
337	IP_NS[counter_NS]=list(I[i],iv);
338	}
339	}
340	else
341	{
342	// the point is (a:1:0) \| a is root of aux
343	if (deg(aux)==1)
344	{
345	// the point is rational and no field extension is needed
346	setring r_auxz;
347	poly f_xz=subst(F,y,1);
348	poly aux=imap(base_r,aux);
349	number A=-number(subst(aux,x,0));
350	map phi=r_auxz,x+A,0,z;
351	poly f_origin=phi(f_xz);
352	if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
353	subst(subst(diff(f_origin,z),x,0),z,0)==0 )
354	{
355	// the point is singular
356	counter_S=counter_S+1;
357	kill f_xz,aux,A,phi,f_origin;
358	setring base_r;
359	IP_S[counter_S]=list(I[i],iv);
360	}
361	else
362	{
363	// the point is non-singular
364	counter_NS=counter_NS+1;
365	kill f_xz,aux,A,phi,f_origin;
366	setring base_r;
367	IP_NS[counter_NS]=list(I[i],iv);
368	}
369	}
370	else
371	{
372	// the point is non-rational and a field extension with minpoly=aux
373	// is needed
374	ring r_ext=(char(basering),@a),(x,y,z),lp;
375	poly F=imap(r_auxz,F);
376	poly f_xz=subst(F,y,1);
377	poly aux=imap(base_r,aux);
378	minpoly=number(subst(aux,x,@a));
379	map phi=r_ext,x+@a,0,z;
380	poly f_origin=phi(f_xz);
381	if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
382	subst(subst(diff(f_origin,z),x,0),z,0)==0 )
383	{
384	// the point is singular
385	counter_S=counter_S+1;
386	setring base_r;
387	kill r_ext;
388	IP_S[counter_S]=list(I[i],iv);
389	}
390	else
391	{
392	// the point is non-singular
393	counter_NS=counter_NS+1;
394	setring base_r;
395	kill r_ext;
396	IP_NS[counter_NS]=list(I[i],iv);
397	}
398	}
399	}
400	}
401	kill r_auxz;
402	return(list(IP_S,IP_NS));
403	}
404	///////////////////////////////////////////////////////////////////////////////
405	static proc closed_points_ext (poly f,int d,ideal SL)
406	{
407	// computes : (closed affine non-singular) points over an extension of
408	// degree d
409	// remark(1) : singular points are supposed to be listed appart
410	// remark(2) : std SL=s_locus(f) is supposed to be computed in advance
411	// remark(3) : ideal SL is used to remove those points which are singular
412	// output : list of list of prime ideals with an intvec for storing the
413	// places
414	int Q=char(basering)^d; // cardinality of the extension field
415	ideal I=f,x^Q-x,y^Q-y; // ideal of the searched points
416	I=std(I);
417	if (I==1)
418	{
419	return(list());
420	}
421	list LP=list();
422	int m=size(SL);
423	list L=list();
424	ideal aux;
425	intvec iv;
426	int i,s,j,counter;
427	L=closed_points(I);
428	s=size(L);
429	for (i=1;i<=s;i=i+1)
430	{
431	aux=std(L[i]);
432	for (j=1;j<=m;j=j+1)
433	{
434	// check if singular i.e. if SL is contained in aux
435	if ( NF(SL[j],aux) != 0 )
436	{
437	counter=counter+1;
438	LP[counter]=list(aux,iv);
439	break;
440	}
441	}
442	}
443	return(LP);
444	}
445	///////////////////////////////////////////////////////////////////////////////
446	static proc degree_P (list P)
447	"USAGE: degree_P(P), where P is either a polynomial or an ideal
448	RETURN: integer with the degree of the closed point given by P
449	SEE ALSO: closed_points
450	NOTE: If P is a (homogeneous irreducible) polynomial the point is at
451	infinity, and if P is a (prime) ideal the points is affine, and
452	the ideal must be given by 2 generators: the first one irreducible
453	and depending only on y, and the second one irreducible over the
454	extension given by y with the first generator as minimal polynomial
455	"
456	{
457	// computes : the degree of a given point
458	// remark(1) : if the input is (irreducible homogeneous) polynomial => the point
459	// is at infinity
460	// remark(2) : it the input is (std. resp. lp. prime) ideal => the point is
461	// affine
462	if (typeof(P[1])=="ideal")
463	{
464	if (size(P[1])==2)
465	{
466	int d=deg(P[1][1]);
467	poly aux=subst(P[1][2],y,1);
468	d=d*deg(aux);
469	return(d);
470	}
471	else
472	{
473	// this should not happen in principle
474	ERROR("non-valid parameter");
475	}
476	}
477	else
478	{
479	if (typeof(P[1])=="poly")
480	{
481	return(deg(P[1]));
482	}
483	else
484	{
485	ERROR("parameter must have a polynomial or ideal in the first component");
486	}
487	}
488	}
489	///////////////////////////////////////////////////////////////////////////////
490	static proc closed_points_deg (poly f,int d,ideal SL)
491	{
492	// computes : (closed affine non-singular) points of degree d
493	// remark(1) : singular points are supposed to be listed appart
494	// remark(2) : std SL=s_locus(f) is supposed to be computed in advance
495	list L=closed_points_ext(f,d,SL);
496	int s=size(L);
497	int i,counter;
498	list LP=list();
499	for (i=1;i<=s;i=i+1)
500	{
501	if (degree_P(L[i])==d)
502	{
503	counter=counter+1;
504	LP[counter]=L[i];
505	}
506	}
507	return(LP);
508	}
509	///////////////////////////////////////////////////////////////////////////////
510	static proc idealSubset (ideal I,ideal J)
511	{
512	// checks wether I is contained in J and returns a boolean
513	// remark : J is assumed to be given by a standard basis
514	int s=ncols(I);
515	int i;
516	for (i=1;i<=s;i=i+1)
517	{
518	if ( NF(I[i],std(J)) != 0 )
519	{
520	return(0);
521	}
522	}
523	return(1);
524	}
525	///////////////////////////////////////////////////////////////////////////////
526	static proc belongs (list P,ideal I)
527	{
528	// checks if affine point P is contained in V(I) and returns a boolean
529	// remark : P[1] is assumed to be an ideal given by a standard basis
530	if (typeof(P[1])=="ideal")
531	{
532	return(idealSubset(I,P[1]));
533	}
534	else
535	{
536	ERROR("first argument must be an affine point");
537	}
538	}
539	///////////////////////////////////////////////////////////////////////////////
540	static proc equals (ideal I,ideal J)
541	{
542	// checks if I is equal to J and returns a boolean
543	// remark : I and J are assumed to be given by a standard basis
544	int answer=0;
545	if (idealSubset(I,J)==1)
546	{
547	if (idealSubset(J,I)==1)
548	{
549	answer=1;
550	}
551	}
552	return(answer);
553	}
554	///////////////////////////////////////////////////////////////////////////////
555	static proc isInLP (ideal P,list LP)
556	{
557	// checks if affine point P is a list LP and returns either its position or
558	// zero
559	// remark : all points in LP and P itself are assumed to be given by a
560	// standard basis
561	// warning : the procedure does not check whether the points are affine or
562	// not
563	int s=size(LP);
564	if (s==0)
565	{
566	return(0);
567	}
568	int i;
569	for (i=1;i<=s;i=i+1)
570	{
571	if (equals(P,LP[i][1])==1)
572	{
573	return(i);
574	}
575	}
576	return(0);
577	}
578	///////////////////////////////////////////////////////////////////////////////
579	static proc res_deg ()
580	{
581	// computes the residual degree of the basering with respect to its prime
582	// field
583	// warning : minpoly must depend on a parameter called "a"
584	int ext;
585	string s_m=string(minpoly);
586	if (s_m=="0")
587	{
588	ext=1;
589	}
590	else
591	{
592	ring auxr=char(basering),a,lp;
593	execute("poly minp="+s_m+";");
594	ext=deg(minp);
595	}
596	return(ext);
597	}
598	///////////////////////////////////////////////////////////////////////////////
599	static proc Frobenius (etwas,int r)
600	{
601	// applies the Frobenius map over F_{p^r} to an object defined over an
602	// extension of such field
603	// usually it is called with r=1, i.e. the Frobenius map over the prime
604	// field F_p
605	// returns always an object of the same type, and works correctly on
606	// numbers, polynomials, ideals, matrices or lists of the above types
607	// maybe : types vector and module should be added in the future, but they
608	// are not needed now
609	int q=char(basering)^r;
610	if (typeof(etwas)=="number")
611	{
612	return(etwas^q);
613	}
614	if (typeof(etwas)=="poly")
615	{
616	int s=size(etwas);
617	poly f;
618	int i;
619	for (i=1;i<=s;i=i+1)
620	{
621	f=f+(leadcoef(etwas[i])^q)*leadmonom(etwas[i]);
622	}
623	return(f);
624	}
625	if (typeof(etwas)=="ideal")
626	{
627	int s=ncols(etwas);
628	ideal I;
629	int i;
630	for (i=1;i<=s;i=i+1)
631	{
632	I[i]=Frobenius(etwas[i],r);
633	}
634	return(I);
635	}
636	if (typeof(etwas)=="matrix")
637	{
638	int m=nrows(etwas);
639	int n=ncols(etwas);
640	matrix A[m][n];
641	int i,j;
642	for (i=1;i<=m;i=i+1)
643	{
644	for (j=1;j<=n;j=j+1)
645	{
646	A[i,j]=Frobenius(etwas[i,j],r);
647	}
648	}
649	return(A);
650	}
651	if (typeof(etwas)=="list")
652	{
653	int s=size(etwas);
654	list L=list();
655	int i;
656	for (i=1;i<=s;i=i+1)
657	{
658	if (typeof(etwas[i])<>"none")
659	{
660	L[i]=Frobenius(etwas[i],r);
661	}
662	}
663	return(L);
664	}
665	return(etwas);
666	}
667	///////////////////////////////////////////////////////////////////////////////
668	static proc conj_b (list L,int r)
669	{
670	// applies the Frobenius map over F_{p^r} to a list of type HNE defined over
671	// a larger extension
672	// when r=1 it turns to be the Frobenius map over the prime field F_{p}
673	// returns : a list of type HNE which is either conjugate of the input or
674	// the same list in case of L being actually defined over the base field
675	// F_{p^r}
676	int p=char(basering);
677	int Q=p^r;
678	list LL=list();
679	int m=nrows(L[1]);
680	int n=ncols(L[1]);
681	matrix A[m][n];
682	poly f;
683	poly aux;
684	int i,j;
685	for (i=1;i<=m;i=i+1)
686	{
687	for (j=1;j<=n;j=j+1)
688	{
689	aux=L[1][i,j];
690	if (aux<>x)
691	{
692	A[i,j]=aux^Q;
693	}
694	else
695	{
696	A[i,j]=aux;
697	break;
698	}
699	}
700	}
701	m=size(L[4]);
702	for (i=1;i<=m;i=i+1)
703	{
704	f=f+(leadcoef(L[4][i])^Q)*leadmonom(L[4][i]);
705	}
706	LL[1]=A;
707	LL[2]=L[2];
708	LL[3]=L[3];
709	LL[4]=f;
710	return(LL);
711	}
712	///////////////////////////////////////////////////////////////////////////////
713	static proc grad_b (list L,int r)
714	{
715	// computes the degree of a list of type HNE which is actually defined over
716	// F_{p^r} eventhough it is given in an extension of such field
717	int gr=1;
718	int rd=res_deg() div r;
719	list LL=L;
720	int i;
721	for (i=1;i<=rd;i=i+1)
722	{
723	LL=conj_b(LL,r);
724	if ( LL[1]==L[1] && LL[4]==L[4] )
725	{
726	break;
727	}
728	else
729	{
730	gr=gr+1;
731	}
732	}
733	return(gr);
734	}
735	///////////////////////////////////////////////////////////////////////////////
736	static proc conj_bs (list L,int r)
737	{
738	// computes all the conjugates over F_{p^r} of a list of type HNE defined
739	// over an extension
740	// returns : a list of lists of type HNE, where the first one is the input
741	// list
742	// remark : notice that the degree of the branch is then the size of the
743	// output
744	list branches=list();
745	int gr=1;
746	branches[1]=L;
747	int rd=res_deg() div r;
748	list LL=L;
749	int i;
750	for (i=1;i<=rd;i=i+1)
751	{
752	LL=conj_b(LL,r);
753	if ( LL[1]==L[1] && LL[4]==L[4] )
754	{
755	break;
756	}
757	else
758	{
759	gr=gr+1;
760	branches[gr]=LL;
761	}
762	}
763	return(branches);
764	}
765	///////////////////////////////////////////////////////////////////////////////
766	static proc subfield (sf)
767	{
768	// writes the generator "a" of a subfield of the coefficients field of
769	// basering in terms of the current generator (also called "a") as a
770	// string sf is an existing ring whose coefficient field is such a subfield
771	// warning : in basering there must be a variable called "x" and subfield
772	// must not be prime
773	def base_r=basering;
774	string new_m=string(minpoly);
775	setring sf;
776	string old_m=string(minpoly);
777	if (old_m==new_m)
778	{
779	setring base_r;
780	return("a");
781	}
782	else
783	{
784	if (old_m<>string(0))
785	{
786	ring auxring=char(basering),(a,x),lp;
787	execute("poly mpol="+old_m+";");
788	mpol=subst(mpol,a,x);
789	setring base_r;
790	poly mpol=imap(auxring,mpol);
791	kill auxring;
792	string answer="? error : non-primitive element";
793	int r=res_deg();
794	int q=char(basering)^r;
795	int i;
796	number b;
797	for (i=1;i<=q-2;i=i+1)
798	{
799	b=a^i;
800	if (subst(mpol,x,b)==0)
801	{
802	answer=string(b);
803	break;
804	}
805	}
806	if (answer<>"? error : non-primitive element")
807	{
808	return(answer);
809	}
810	else
811	{
812	// list all the elements of the finite field F_q
813	int p=char(basering);
814	list FF1,FF2;
815	for (i=0;i<p;i=i+1)
816	{
817	FF1[i+1]=number(i);
818	}
819	int s,j,k;
820	for (i=1;i<r;i=i+1)
821	{
822	s=size(FF1);
823	for (j=1;j<=s;j=j+1)
824	{
825	FF1[j]=FF1[j]*a;
826	}
827	FF2=FF1;
828	for (k=1;k<p;k=k+1)
829	{
830	for (j=1;j<=s;j=j+1)
831	{
832	FF1[j]=FF1[j]+number(1);
833	}
834	FF2=FF2+FF1;
835	}
836	FF1=FF2;
837	FF2=list();
838	}
839	kill FF2;
840	for (i=1;i<=q-1;i=i+1)
841	{
842	b=FF1[i];
843	if (subst(mpol,x,b)==0)
844	{
845	answer=string(b);
846	break;
847	}
848	}
849	}
850	kill FF1;
851	return(answer);
852	}
853	else
854	{
855	dbprint(printlevel+1,"warning : minpoly=0 in the subfield;
856	you should check that nothing is wrong");
857	return(string(1));
858	}
859	}
860	}
861	///////////////////////////////////////////////////////////////////////////////
862	static proc importdatum (sf,string datum,string rel)
863	{
864	// fetchs a poly with name "datum" to the current basering from the ring sf
865	// such that the generator is given by string "rel"
866	// warning : ring sf must have only variables (x,y) and basering must have
867	// at least (x,y)
868	// warning : the case of minpoly=0 is not regarded; there you can use "imap"
869	// instead
870	def base_r=basering;
871	if (rel=="a")
872	{
873	setring sf;
874	execute("poly pdatum="+datum+";");
875	setring base_r;
876	poly pdatum=imap(sf,pdatum);
877	return(pdatum);
878	}
879	else
880	{
881	setring sf;
882	execute("string sdatum=string("+datum+");");
883	ring auxring=char(basering),(a,x,y),lp;
884	execute("poly pdatum="+sdatum+";");
885	execute("map phi=basering,"+rel+",x,y;");
886	pdatum=phi(pdatum);
887	string snewdatum=string(pdatum);
888	setring base_r;
889	execute("poly pdatum="+snewdatum+";");
890	kill auxring;
891	return(pdatum);
892	}
893	}
894	///////////////////////////////////////////////////////////////////////////////
895	static proc rationalize (lf,string datum,string rel)
896	{
897	// fetchs a poly with name "datum" to the current basering from the ring lf
898	// and larger coefficients field, where the generator of current ring is
899	// given by string "rel" and "datum" is actually defined over the small
900	// field
901	// warning : ring lf must have only variables (x,y) and basering must have
902	// at least (x,y)
903	// warning : the case of minpoly=0 is supposed unnecessary, since then
904	// "datum" should be already written in the right way, i.e. in
905	// terms of the prime field
906	def base_r=basering;
907	if (rel=="a")
908	{
909	setring lf;
910	execute("poly pdatum="+datum+";");
911	setring base_r;
912	poly pdatum=imap(lf,pdatum);
913	return(pdatum);
914	}
915	else
916	{
917	setring lf;
918	execute("string sdatum=string("+datum+");");
919	ring auxring=char(basering),(a,b,x,y,t),lp;
920	// execute("poly pdatum="+datum+";");
921	execute("poly pdatum="+sdatum+";");
922	execute("poly prel=b-("+rel+");");
923	// was before:
924	// ideal I=pdatum,prel;
925	// I=eliminate(I,a);
926	// if (size(I)<=1){ poly newdatum=I[1]; }
927	// else { ERROR("Elimination error"); }
928	// CL 09/05 new:
929	poly newdatum=reduce(pdatum,std(prel));
930	newdatum;
931	//
932	newdatum=subst(newdatum,b,a);
933	string snewdatum=string(newdatum);
934	setring base_r;
935	execute("poly newdatum="+snewdatum+";");
936	kill auxring;
937	return(newdatum);
938	}
939	}
940	///////////////////////////////////////////////////////////////////////////////
941	static proc place (intvec Pp,int sing,list CURVE)
942	{
943	// computes the "rational" places which are defined over a (closed) point
944	// Pp points to an appropriate point of the given curve
945	// creates : local rings (if they do not exist yet) and then add the
946	// places to a list
947	// each place is given basically by the coordinates of the point and a
948	// list HNdevelop
949	// returns : list with all updated data of the curve
950	// if the places already exist they are not computed again
951	// if sing==1 the point is assumed singular and computes the local conductor
952	// for all places using the local invariants of the branches
953	// if sing==2 the point is assumed singular and computes the local conductor
954	// for all places using the Dedekind formula and local parametrizations
955	// of the branches
956	// if sing<>1&2 the point is assumed non-singular and the local conductor
957	// should be zero
958	list PP=list();
959	if (Pp[1]==0)
960	{
961	if (Pp[2]==0)
962	{
963	PP=Aff_SPoints[Pp[3]];
964	}
965	if (Pp[2]==1)
966	{
967	PP=Inf_Points[1][Pp[3]];
968	}
969	if (Pp[2]==2)
970	{
971	PP=Inf_Points[2][Pp[3]];
972	}
973	}
974	else
975	{
976	PP=Aff_Points(Pp[2])[Pp[3]];
977	}
978	if (PP[2]<>0)
979	{
980	return(CURVE);
981	}
982	intvec PtoPl;
983	def base_r=basering;
984	int ext1;
985	list Places=CURVE[3];
986	intvec Conductor=CURVE[4];
987	list update_CURVE=CURVE;
988	if (typeof(PP[1])=="ideal")
989	{
990	ideal P=PP[1];
991	if (size(P)==2)
992	{
993	int d=deg(P[1]);
994	poly aux=subst(P[2],y,1);
995	d=d*deg(aux);
996	ext1=d;
997	// the point is (A:B:1) but one must distinguish several cases
998	// P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends
999	// only on "y"
1000	if (d==1)
1001	{
1002	// the point is rational
1003	number B=-number(subst(P[1],y,0));
1004	poly aux2=subst(P[2],y,B);
1005	number A=-number(subst(aux2,x,0));
1006	// the point is (A:B:1)
1007	ring local_aux=char(basering),(x,y),ls;
1008	number coord@1=imap(base_r,A);
1009	number coord@2=imap(base_r,B);
1010	number coord@3=number(1);
1011	map phi=base_r,x+coord@1,y+coord@2;
1012	poly CHI=phi(CHI);
1013	}
1014	else
1015	{
1016	if (deg(P[1])==1)
1017	{
1018	// the point is non-rational but the second component needs no
1019	// field extension
1020	number B=-number(subst(P[1],y,0));
1021	poly aux2=subst(P[2],y,B);
1022	// the point has degree d>1
1023	// careful : the parameter will be called "a" anyway
1024	ring local_aux=(char(basering),a),(x,y),ls;
1025	map psi=base_r,a,0;
1026	minpoly=number(psi(aux2));
1027	number coord@1=a;
1028	number coord@2=imap(base_r,B);
1029	number coord@3=number(1);
1030	// the point is (a:B:1)
1031	map phi=base_r,x+a,y+coord@2;
1032	poly CHI=phi(CHI);
1033	}
1034	else
1035	{
1036	if (deg(subst(P[2],y,1))==1)
1037	{
1038	// the point is non-rational but the needed minpoly is just P[1]
1039	// careful : the parameter will be called "a" anyway
1040	poly P1=P[1];
1041	poly P2=P[2];
1042	ring local_aux=(char(basering),a),(x,y),ls;
1043	map psi=base_r,0,a;
1044	minpoly=number(psi(P1));
1045	// the point looks like (A:a:1)
1046	// A is computed by substituting y=a in P[2]
1047	poly aux1=imap(base_r,P2);
1048	poly aux2=subst(aux1,y,a);
1049	number coord@1=-number(subst(aux2,x,0));
1050	number coord@2=a;
1051	number coord@3=number(1);
1052	map phi=base_r,x+coord@1,y+a;
1053	poly CHI=phi(CHI);
1054	}
1055	else
1056	{
1057	// this is the most complicated case of non-rational point
1058	// firstly : construct an extension of degree d and guess the
1059	// minpoly
1060	poly P1=P[1];
1061	poly P2=P[2];
1062	int p=char(basering);
1063	int Q=p^d;
1064	ring aux_r=(Q,a),(x,y,t),ls;
1065	string minpoly_string=string(minpoly);
1066	ring local_aux=(char(basering),a),(x,y),ls;
1067	execute("minpoly="+minpoly_string+";");
1068	// secondly : compute one root of P[1]
1069	poly P_1=imap(base_r,P1);
1070	poly P_2=imap(base_r,P2);
1071	ideal factors1=factorize(P_1,1); // hopefully this works !!!!
1072	number coord@2=-number(subst(factors1[1],y,0));
1073	// thirdly : compute one of the first components for the above root
1074	poly P_0=subst(P_2,y,coord@2);
1075	ideal factors2=factorize(P_0,1); // hopefully this works !!!!
1076	number coord@1=-number(subst(factors2[1],x,0));
1077	number coord@3=number(1);
1078	map phi=base_r,x+coord@1,y+coord@2;
1079	poly CHI=phi(CHI);
1080	kill aux_r;
1081	}
1082	}
1083	}
1084	}
1085	else
1086	{
1087	// this should not happen in principle
1088	ERROR("non-valid parameter");
1089	}
1090	}
1091	else
1092	{
1093	if (typeof(PP[1])=="poly")
1094	{
1095	poly P=PP[1];
1096	ring r_auxz=char(basering),(x,y,z),lp;
1097	poly CHI=imap(base_r,CHI);
1098	CHI=homog(CHI,z);
1099	setring base_r;
1100	poly aux=subst(P,y,1);
1101	if (aux==1)
1102	{
1103	// the point is (1:0:0)
1104	ring local_aux=char(basering),(x,y),ls;
1105	number coord@1=number(1);
1106	number coord@2=number(0);
1107	number coord@3=number(0);
1108	map Phi=r_auxz,1,x,y;
1109	poly CHI=Phi(CHI);
1110	ext1=1;
1111	}
1112	else
1113	{
1114	// the point is (A:1:0) where A is a root of aux
1115	int d=deg(aux);
1116	ext1=d;
1117	if (d==1)
1118	{
1119	// the point is rational
1120	number A=-number(subst(aux,x,0));
1121	ring local_aux=char(basering),(x,y),ls;
1122	number coord@1=imap(base_r,A);
1123	number coord@2=number(1);
1124	number coord@3=number(0);
1125	map Phi=r_auxz,x+coord@1,1,y;
1126	poly CHI=Phi(CHI);
1127	}
1128	else
1129	{
1130	// the point has degree d>1
1131	// careful : the parameter will be called "a" anyway
1132	ring local_aux=(char(basering),a),(x,y),ls;
1133	map psi=base_r,a,1;
1134	minpoly=number(psi(P));
1135	number coord@1=a;
1136	number coord@2=number(1);
1137	number coord@3=number(0);
1138	map Phi=r_auxz,x+a,1,y;
1139	poly CHI=Phi(CHI);
1140	}
1141	}
1142	kill r_auxz;
1143	}
1144	else
1145	{
1146	ERROR("a point must have a polynomial or ideal in the first component");
1147	}
1148	}
1149	export(coord@1);
1150	export(coord@2);
1151	export(coord@3);
1152	export(CHI);
1153	int i,j,k;
1154	int m,n;
1155	list LLL=ratdevelop(CHI);
1156	// list LLL=hnexpansion(CHI,"ess");
1157	if (typeof(LLL[1])=="ring") {
1158	def altring=basering;
1159	def HNEring = LLL[1]; setring HNEring;
1160	def L@HNE = HND;
1161	kill HND;
1162	// def L@HNE = hne;
1163	// kill hne;
1164
1165	}
1166	else {
1167	def L@HNE=LLL;
1168	// def L@HNE=LLL[1];
1169	def HNEring=basering;
1170	setring HNEring;
1171	}
1172	kill LLL;
1173	export(L@HNE);
1174	int n_branches=size(L@HNE);
1175	list Li_aux=list();
1176	int N_branches;
1177	int N=size(Places);
1178	if (sing==1)
1179	{
1180	list delta2=list();
1181	for (i=1;i<=n_branches;i=i+1)
1182	{
1183	delta2[i]=invariants(L@HNE[i])[5];
1184	}
1185	int dq;
1186	}
1187	int ext2=res_deg();
1188	list dgs=list();
1189	int ext_0;
1190	int check;
1191	string sss,olda,newa;
1192	if (defined(Q)==0)
1193	{
1194	int Q;
1195	}
1196	if (ext1==1)
1197	{
1198	if (ext2==1)
1199	{
1200	if (sing==1)
1201	{
1202	intmat I_mult[n_branches][n_branches];
1203	if (n_branches>1)
1204	{
1205	for (i=1;i<=n_branches-1;i=i+1)
1206	{
1207	for (j=i+1;j<=n_branches;j=j+1)
1208	{
1209	I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
1210	I_mult[j,i]=I_mult[i,j];
1211	}
1212	}
1213	}
1214	}
1215	////////
1216	check=0;
1217	////////
1218	if (size(update_CURVE[5])>0)
1219	{
1220	if (typeof(update_CURVE[5][1])=="list")
1221	{
1222	check=1;
1223	}
1224	}
1225	if (check==0)
1226	{
1227	intvec dgs_points(1);
1228	ring S(1)=char(basering),(x,y,t),ls;
1229	list BRANCHES=list();
1230	list POINTS=list();
1231	list LOC_EQS=list();
1232	list PARAMETRIZATIONS=list();
1233	export(BRANCHES);
1234	export(POINTS);
1235	export(LOC_EQS);
1236	export(PARAMETRIZATIONS);
1237	}
1238	else
1239	{
1240	intvec dgs_points(1)=update_CURVE[5][1][2];
1241	def S1=update_CURVE[5][1][1];
1242	execute("ring S(1)="+string(update_CURVE[5][1][1])+";");
1243	fetchall(S1);
1244	kill S1;
1245	}
1246	N_branches=size(BRANCHES);
1247	for (i=1;i<=n_branches;i=i+1)
1248	{
1249	dgs_points(1)[N_branches+i]=1;
1250	POINTS[N_branches+i]=list();
1251	POINTS[N_branches+i][1]=imap(local_aux,coord@1);
1252	POINTS[N_branches+i][2]=imap(local_aux,coord@2);
1253	POINTS[N_branches+i][3]=imap(local_aux,coord@3);
1254	LOC_EQS[N_branches+i]=imap(local_aux,CHI);
1255	setring HNEring;
1256	Li_aux=L@HNE[i];
1257	setring S(1);
1258	BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches+i-1);
1259	PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
1260	N=N+1;
1261	intvec iw=1,N_branches+i;
1262	Places[N]=iw;
1263	if (sing==1)
1264	{
1265	dq=delta2[i];
1266	for (j=1;j<=n_branches;j=j+1)
1267	{
1268	dq=dq+I_mult[i,j];
1269	}
1270	Conductor[N]=dq;
1271	}
1272	if (sing==2)
1273	{
1274	Conductor[N]=local_conductor(iw[2],S(1));
1275	}
1276	kill iw;
1277	PtoPl[i]=N;
1278	}
1279	setring base_r;
1280	update_CURVE[5][1]=list();
1281	update_CURVE[5][1][1]=S(1);
1282	update_CURVE[5][1][2]=dgs_points(1);
1283	}
1284	else
1285	{
1286	// we start with a rational point but we get non-rational branches
1287	// they may have different degrees and then we may need to reduce the
1288	// field extensions for each one, and finally check if the minpoly
1289	// fetches with S(i) or not
1290	// if one of the branches is rational, we may trust that it is written
1291	// correctly
1292	if (sing==1)
1293	{
1294	int n_geobrs;
1295	int counter_c;
1296	list auxgb=list();
1297	list geobrs=list();
1298	for (i=1;i<=n_branches;i=i+1)
1299	{
1300	auxgb=conj_bs(L@HNE[i],1);
1301	dgs[i]=size(auxgb);
1302	n_geobrs=n_geobrs+dgs[i];
1303	for (j=1;j<=dgs[i];j=j+1)
1304	{
1305	counter_c=counter_c+1;
1306	geobrs[counter_c]=auxgb[j];
1307	}
1308	}
1309	intmat I_mult[n_geobrs][n_geobrs];
1310	for (i=1;i<n_geobrs;i=i+1)
1311	{
1312	for (j=i+1;j<=n_geobrs;j=j+1)
1313	{
1314	I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
1315	I_mult[j,i]=I_mult[i,j];
1316	}
1317	}
1318	kill auxgb,geobrs;
1319	}
1320	else
1321	{
1322	for (i=1;i<=n_branches;i=i+1)
1323	{
1324	dgs[i]=grad_b(L[i],1);
1325	}
1326	}
1327	// the actual degree of each branch is computed and now check if the
1328	// local ring exists
1329	for (i=1;i<=n_branches;i=i+1)
1330	{
1331	ext_0=dgs[i];
1332	////////
1333	check=0;
1334	////////
1335	if (size(update_CURVE[5])>=ext_0)
1336	{
1337	if (typeof(update_CURVE[5][ext_0])=="list")
1338	{
1339	check=1;
1340	}
1341	}
1342	if (check==0)
1343	{
1344	if (ext_0>1)
1345	{
1346	if (ext_0==ext2)
1347	{
1348	sss=string(minpoly);
1349	}
1350	else
1351	{
1352	Q=char(basering)^ext_0;
1353	ring auxxx=(Q,a),z,lp;
1354	sss=string(minpoly);
1355	setring base_r;
1356	kill auxxx;
1357	}
1358	ring S(ext_0)=(char(basering),a),(x,y,t),ls;
1359	execute("minpoly="+sss+";");
1360	}
1361	else
1362	{
1363	ring S(ext_0)=char(basering),(x,y,t),ls;
1364	}
1365	intvec dgs_points(ext_0);
1366	list BRANCHES=list();
1367	list POINTS=list();
1368	list LOC_EQS=list();
1369	list PARAMETRIZATIONS=list();
1370	export(BRANCHES);
1371	export(POINTS);
1372	export(LOC_EQS);
1373	export(PARAMETRIZATIONS);
1374	}
1375	else
1376	{
1377	intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
1378	def Sext_0=update_CURVE[5][ext_0][1];
1379	setring Sext_0;
1380	string SM=string(minpoly);
1381	string SR=string(update_CURVE[5][ext_0][1]);
1382	execute("ring S("+string(ext_0)+")="+SR+";");
1383	execute("minpoly="+SM+";");
1384	kill SM,SR;
1385	fetchall(Sext_0);
1386	kill Sext_0;
1387	}
1388	N_branches=size(BRANCHES);
1389	dgs_points(ext_0)[N_branches+1]=1;
1390	POINTS[N_branches+1]=list();
1391	POINTS[N_branches+1][1]=imap(local_aux,coord@1);
1392	POINTS[N_branches+1][2]=imap(local_aux,coord@2);
1393	POINTS[N_branches+1][3]=imap(local_aux,coord@3);
1394	LOC_EQS[N_branches+1]=imap(local_aux,CHI);
1395	// now fetch the branches into the new local ring
1396	if (ext_0==1)
1397	{
1398	setring HNEring;
1399	Li_aux=L@HNE[i];
1400	setring S(1);
1401	BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches);
1402	}
1403	else
1404	{
1405	// rationalize branch
1406	setring HNEring;
1407	newa=subfield(S(ext_0));
1408	m=nrows(L@HNE[i][1]);
1409	n=ncols(L@HNE[i][1]);
1410	setring S(ext_0);
1411	list Laux=list();
1412	poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
1413	matrix Maux[m][n];
1414	for (j=1;j<=m;j=j+1)
1415	{
1416	for (k=1;k<=n;k=k+1)
1417	{
1418	Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
1419	string(j)+","+string(k)+"]",newa);
1420	}
1421	}
1422	setring HNEring;
1423	intvec Li2=L@HNE[i][2];
1424	int Li3=L@HNE[i][3];
1425	setring S(ext_0);
1426	Laux[1]=Maux;
1427	Laux[2]=Li2;
1428	Laux[3]=Li3;
1429	Laux[4]=paux;
1430	BRANCHES=insert(BRANCHES,Laux,N_branches);
1431	kill Laux,Maux,paux,Li2,Li3;
1432	}
1433	PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
1434	N=N+1;
1435	intvec iw=ext_0,N_branches+1;
1436	Places[N]=iw;
1437	if (sing==2)
1438	{
1439	Conductor[N]=local_conductor(iw[2],S(ext_0));
1440	}
1441	kill iw;
1442	PtoPl[i]=N;
1443	setring HNEring;
1444	update_CURVE[5][ext_0]=list();
1445	update_CURVE[5][ext_0][1]=S(ext_0);
1446	update_CURVE[5][ext_0][2]=dgs_points(ext_0);
1447	}
1448	if (sing==1)
1449	{
1450	int N_ini=N-n_branches;
1451	counter_c=1;
1452	for (i=1;i<=n_branches;i=i+1)
1453	{
1454	dq=delta2[i];
1455	for (j=1;j<=n_geobrs;j=j+1)
1456	{
1457	dq=dq+I_mult[counter_c,j];
1458	}
1459	Conductor[N_ini+i]=dq;
1460	counter_c=counter_c+dgs[i];
1461	}
1462	}
1463	setring base_r;
1464	}
1465	}
1466	else
1467	{
1468	if (ext1==ext2)
1469	{
1470	// the degree of the point equals to the degree of all branches
1471	// one must just fetch the minpoly's of local_aux, HNEring and S(ext2)
1472	if (sing==1)
1473	{
1474	intmat I_mult[n_branches][n_branches];
1475	if (n_branches>1)
1476	{
1477	for (i=1;i<=n_branches-1;i=i+1)
1478	{
1479	for (j=i+1;j<=n_branches;j=j+1)
1480	{
1481	I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
1482	I_mult[j,i]=I_mult[i,j];
1483	}
1484	}
1485	}
1486	}
1487	////////
1488	check=0;
1489	////////
1490	if (size(update_CURVE[5])>=ext2)
1491	{
1492	if (typeof(update_CURVE[5][ext2])=="list")
1493	{
1494	check=1;
1495	}
1496	}
1497	if (check==0)
1498	{
1499	sss=string(minpoly);
1500	ring S(ext2)=(char(basering),a),(x,y,t),ls;
1501	execute("minpoly="+sss+";");
1502	intvec dgs_points(ext2);
1503	list BRANCHES=list();
1504	list POINTS=list();
1505	list LOC_EQS=list();
1506	list PARAMETRIZATIONS=list();
1507	export(BRANCHES);
1508	export(POINTS);
1509	export(LOC_EQS);
1510	export(PARAMETRIZATIONS);
1511	}
1512	else
1513	{
1514	intvec dgs_points(ext2)=update_CURVE[5][ext2][2];
1515	def Sext2=update_CURVE[5][ext2][1];
1516	setring Sext2;
1517	string SM=string(minpoly);
1518	string SR=string(update_CURVE[5][ext2][1]);
1519	execute("ring S("+string(ext2)+")="+SR+";");
1520	execute("minpoly="+SM+";");
1521	kill SM,SR;
1522	fetchall(Sext2);
1523	kill Sext2;
1524	}
1525	N_branches=size(BRANCHES);
1526	for (i=1;i<=n_branches;i=i+1)
1527	{
1528	// fetch all the data into the new local ring
1529	olda=subfield(local_aux);
1530	dgs_points(ext2)[N_branches+i]=ext1;
1531	POINTS[N_branches+i]=list();
1532	POINTS[N_branches+i][1]=number(importdatum(local_aux,"coord@1",olda));
1533	POINTS[N_branches+i][2]=number(importdatum(local_aux,"coord@2",olda));
1534	POINTS[N_branches+i][3]=number(importdatum(local_aux,"coord@3",olda));
1535	LOC_EQS[N_branches+i]=importdatum(local_aux,"CHI",olda);
1536	newa=subfield(HNEring);
1537	setring HNEring;
1538	m=nrows(L@HNE[i][1]);
1539	n=ncols(L@HNE[i][1]);
1540	setring S(ext2);
1541	list Laux=list();
1542	poly paux=importdatum(HNEring,"L@HNE["+string(i)+"][4]",newa);
1543	matrix Maux[m][n];
1544	for (j=1;j<=m;j=j+1)
1545	{
1546	for (k=1;k<=n;k=k+1)
1547	{
1548	Maux[j,k]=importdatum(HNEring,"L@HNE["+string(i)+"][1]["+
1549	string(j)+","+string(k)+"]",newa);
1550	}
1551	}
1552	setring HNEring;
1553	intvec Li2=L@HNE[i][2];
1554	int Li3=L@HNE[i][3];
1555	setring S(ext2);
1556	Laux[1]=Maux;
1557	Laux[2]=Li2;
1558	Laux[3]=Li3;
1559	Laux[4]=paux;
1560	BRANCHES=insert(BRANCHES,Laux,N_branches+i-1);
1561	kill Laux,Maux,paux,Li2,Li3;
1562	PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
1563	N=N+1;
1564	intvec iw=ext2,N_branches+i;
1565	Places[N]=iw;
1566	if (sing==1)
1567	{
1568	dq=delta2[i];
1569	for (j=1;j<=n_branches;j=j+1)
1570	{
1571	dq=dq+I_mult[i,j];
1572	}
1573	Conductor[N]=dq;
1574	}
1575	if (sing==2)
1576	{
1577	Conductor[N]=local_conductor(iw[2],S(ext2));
1578	}
1579	kill iw;
1580	PtoPl[i]=N;
1581	}
1582	setring base_r;
1583	update_CURVE[5][ext2]=list();
1584	update_CURVE[5][ext2][1]=S(ext2);
1585	update_CURVE[5][ext2][2]=dgs_points(ext2);
1586	}
1587	else
1588	{
1589	// this is the most complicated case
1590	if (sing==1)
1591	{
1592	int n_geobrs;
1593	int counter_c;
1594	list auxgb=list();
1595	list geobrs=list();
1596	for (i=1;i<=n_branches;i=i+1)
1597	{
1598	auxgb=conj_bs(L@HNE[i],ext1);
1599	dgs[i]=size(auxgb);
1600	n_geobrs=n_geobrs+dgs[i];
1601	for (j=1;j<=dgs[i];j=j+1)
1602	{
1603	counter_c=counter_c+1;
1604	geobrs[counter_c]=auxgb[j];
1605	}
1606	}
1607	intmat I_mult[n_geobrs][n_geobrs];
1608	for (i=1;i<n_geobrs;i=i+1)
1609	{
1610	for (j=i+1;j<=n_geobrs;j=j+1)
1611	{
1612	I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
1613	I_mult[j,i]=I_mult[i,j];
1614	}
1615	}
1616	kill auxgb,geobrs;
1617	}
1618	else
1619	{
1620	for (i=1;i<=n_branches;i=i+1)
1621	{
1622	dgs[i]=grad_b(L@HNE[i],ext1);
1623	}
1624	}
1625	for (i=1;i<=n_branches;i=i+1)
1626	{
1627	// first compute the actual degree of each branch and check if the
1628	// local ring exists
1629	ext_0=ext1*dgs[i];
1630	////////
1631	check=0;
1632	////////
1633	if (size(update_CURVE[5])>=ext_0)
1634	{
1635	if (typeof(update_CURVE[5][ext_0])=="list")
1636	{
1637	check=1;
1638	}
1639	}
1640	if (check==0)
1641	{
1642	if (ext_0>ext1)
1643	{
1644	if (ext_0==ext2)
1645	{
1646	sss=string(minpoly);
1647	}
1648	else
1649	{
1650	Q=char(basering)^ext_0;
1651	ring auxxx=(Q,a),z,lp;
1652	sss=string(minpoly);
1653	setring base_r;
1654	kill auxxx;
1655	}
1656	}
1657	else
1658	{
1659	setring local_aux;
1660	sss=string(minpoly);
1661	}
1662	ring S(ext_0)=(char(basering),a),(x,y,t),ls;
1663	execute("minpoly="+sss+";");
1664	intvec dgs_points(ext_0);
1665	list BRANCHES=list();
1666	list POINTS=list();
1667	list LOC_EQS=list();
1668	list PARAMETRIZATIONS=list();
1669	export(BRANCHES);
1670	export(POINTS);
1671	export(LOC_EQS);
1672	export(PARAMETRIZATIONS);
1673	}
1674	else
1675	{
1676	intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
1677	def Sext_0=update_CURVE[5][ext_0][1];
1678	setring Sext_0;
1679	string SM=string(minpoly);
1680	string SR=string(update_CURVE[5][ext_0][1]);
1681	execute("ring S("+string(ext_0)+")="+SR+";");
1682	execute("minpoly="+SM+";");
1683	kill SM,SR;
1684	fetchall(badring);
1685	kill badring;
1686	}
1687	N_branches=size(BRANCHES);
1688	// now fetch all the data into the new local ring
1689	olda=subfield(local_aux);
1690	dgs_points(ext_0)[N_branches+1]=ext1;
1691	POINTS[N_branches+1]=list();
1692	POINTS[N_branches+1][1]=number(importdatum(local_aux,"coord@1",olda));
1693	POINTS[N_branches+1][2]=number(importdatum(local_aux,"coord@2",olda));
1694	POINTS[N_branches+1][3]=number(importdatum(local_aux,"coord@3",olda));
1695	LOC_EQS[N_branches+1]=importdatum(local_aux,"CHI",olda);
1696	setring HNEring;
1697	newa=subfield(S(ext_0));
1698	m=nrows(L@HNE[i][1]);
1699	n=ncols(L@HNE[i][1]);
1700	setring S(ext_0);
1701	list Laux=list();
1702	poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
1703	matrix Maux[m][n];
1704	for (j=1;j<=m;j=j+1)
1705	{
1706	for (k=1;k<=n;k=k+1)
1707	{
1708	Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
1709	string(j)+","+string(k)+"]",newa);
1710	}
1711	}
1712	setring HNEring;
1713	intvec Li2=L@HNE[i][2];
1714	int Li3=L@HNE[i][3];
1715	setring S(ext_0);
1716	Laux[1]=Maux;
1717	Laux[2]=Li2;
1718	Laux[3]=Li3;
1719	Laux[4]=paux;
1720	BRANCHES=insert(BRANCHES,Laux,N_branches);
1721	kill Laux,Maux,paux,Li2,Li3;
1722	PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
1723	N=N+1;
1724	intvec iw=ext_0,N_branches+1;
1725	Places[N]=iw;
1726	if (sing==2)
1727	{
1728	Conductor[N]=local_conductor(iw[2],S(ext_0));
1729	}
1730	kill iw;
1731	PtoPl[i]=N;
1732	setring HNEring;
1733	update_CURVE[5][ext_0]=list();
1734	update_CURVE[5][ext_0][1]=S(ext_0);
1735	update_CURVE[5][ext_0][2]=dgs_points(ext_0);
1736	}
1737	if (sing==1)
1738	{
1739	int N_ini=N-n_branches;
1740	counter_c=1;
1741	for (i=1;i<=n_branches;i=i+1)
1742	{
1743	dq=delta2[i];
1744	for (j=1;j<=n_geobrs;j=j+1)
1745	{
1746	dq=dq+I_mult[counter_c,j];
1747	}
1748	Conductor[N_ini+i]=dq;
1749	counter_c=counter_c+dgs[i];
1750	}
1751	}
1752	setring base_r;
1753	}
1754	}
1755	update_CURVE[3]=Places;
1756	update_CURVE[4]=Conductor;
1757	PP[2]=PtoPl;
1758	if (Pp[1]==0)
1759	{
1760	if (Pp[2]==0)
1761	{
1762	Aff_SPoints[Pp[3]]=PP;
1763	}
1764	if (Pp[2]==1)
1765	{
1766	Inf_Points[1][Pp[3]]=PP;
1767	}
1768	if (Pp[2]==2)
1769	{
1770	Inf_Points[2][Pp[3]]=PP;
1771	}
1772	}
1773	else
1774	{
1775	Aff_Points(Pp[2])[Pp[3]]=PP;
1776	}
1777	update_CURVE[1][1]=base_r;
1778	kill HNEring;
1779	kill local_aux;
1780	return(update_CURVE);
1781	}
1782	///////////////////////////////////////////////////////////////////////////////
1783	static proc local_conductor (int k,SS)
1784	{
1785	// computes the degree of the local conductor at a place of a plane curve
1786	// if the point is non-singular the result will be zero
1787	// the computation is carried out with the "Dedekind formula" via
1788	// parametrizations
1789	int a,b,Cq;
1790	def b_ring=basering;
1791	setring SS;
1792	poly fx=diff(LOC_EQS[k],x);
1793	poly fy=diff(LOC_EQS[k],y);
1794	int nr=ncols(BRANCHES[k][1]);
1795	poly xt=PARAMETRIZATIONS[k][1][1];
1796	poly yt=PARAMETRIZATIONS[k][1][2];
1797	int ordx=PARAMETRIZATIONS[k][2][1];
1798	int ordy=PARAMETRIZATIONS[k][2][2];
1799	map phi_t=basering,xt,yt,1;
1800	poly derf;
1801	if (fx<>0)
1802	{
1803	derf=fx;
1804	poly tt=diff(yt,t);
1805	b=mindeg(tt);
1806	if (ordy>-1)
1807	{
1808	while (b>=ordy)
1809	{
1810	BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1811	nr=ncols(BRANCHES[k][1]);
1812	PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1813	ordy=PARAMETRIZATIONS[k][2][2];
1814	yt=PARAMETRIZATIONS[k][1][2];
1815	tt=diff(yt,t);
1816	b=mindeg(tt);
1817	}
1818	xt=PARAMETRIZATIONS[k][1][1];
1819	ordx=PARAMETRIZATIONS[k][2][1];
1820	}
1821	poly ft=phi_t(derf);
1822	}
1823	else
1824	{
1825	derf=fy;
1826	poly tt=diff(xt,t);
1827	b=mindeg(tt);
1828	if (ordx>-1)
1829	{
1830	while (b>=ordx)
1831	{
1832	BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1833	nr=ncols(BRANCHES[k][1]);
1834	PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1835	ordx=PARAMETRIZATIONS[k][2][1];
1836	xt=PARAMETRIZATIONS[k][1][1];
1837	tt=diff(xt,t);
1838	b=mindeg(tt);
1839	}
1840	yt=PARAMETRIZATIONS[k][1][2];
1841	ordy=PARAMETRIZATIONS[k][2][2];
1842	}
1843	poly ft=phi_t(derf);
1844	}
1845	a=mindeg(ft);
1846	if ( ordx>-1 \|\| ordy>-1 )
1847	{
1848	if (ordy==-1)
1849	{
1850	while (a>ordx)
1851	{
1852	BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1853	nr=ncols(BRANCHES[k][1]);
1854	PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1855	ordx=PARAMETRIZATIONS[k][2][1];
1856	xt=PARAMETRIZATIONS[k][1][1];
1857	ft=phi_t(derf);
1858	a=mindeg(ft);
1859	}
1860	}
1861	else
1862	{
1863	if (ordx==-1)
1864	{
1865	while (a>ordy)
1866	{
1867	BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1868	nr=ncols(BRANCHES[k][1]);
1869	PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1870	ordy=PARAMETRIZATIONS[k][2][2];
1871	yt=PARAMETRIZATIONS[k][1][2];
1872	ft=phi_t(derf);
1873	a=mindeg(ft);
1874	}
1875	}
1876	else
1877	{
1878	int ordf=ordx;
1879	if (ordx>ordy)
1880	{
1881	ordf=ordy;
1882	}
1883	while (a>ordf)
1884	{
1885	BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1886	nr=ncols(BRANCHES[k][1]);
1887	PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1888	ordx=PARAMETRIZATIONS[k][2][1];
1889	ordy=PARAMETRIZATIONS[k][2][2];
1890	ordf=ordx;
1891	if (ordx>ordy)
1892	{
1893	ordf=ordy;
1894	}
1895	xt=PARAMETRIZATIONS[k][1][1];
1896	yt=PARAMETRIZATIONS[k][1][2];
1897	ft=phi_t(derf);
1898	a=mindeg(ft);
1899	}
1900	}
1901	}
1902	}
1903	Cq=a-b;
1904	setring b_ring;
1905	return(Cq);
1906	}
1907	///////////////////////////////////////////////////////////////////////////////
1908	static proc max_D (intvec D1,intvec D2)
1909	{
1910	// computes the maximum of two divisors (intvec)
1911	int s1=size(D1);
1912	int s2=size(D2);
1913	int i;
1914	if (s1>s2)
1915	{
1916	for (i=1;i<=s2;i=i+1)
1917	{
1918	if (D2[i]>D1[i])
1919	{
1920	D1[i]=D2[i];
1921	}
1922	}
1923	for (i=s2+1;i<=s1;i=i+1)
1924	{
1925	if (D1[i]<0)
1926	{
1927	D1[i]=0;
1928	}
1929	}
1930	return(D1);
1931	}
1932	else
1933	{
1934	for (i=1;i<=s1;i=i+1)
1935	{
1936	if (D1[i]>D2[i])
1937	{
1938	D2[i]=D1[i];
1939	}
1940	}
1941	for (i=s1+1;i<=s2;i=i+1)
1942	{
1943	if (D2[i]<0)
1944	{
1945	D2[i]=0;
1946	}
1947	}
1948	return(D2);
1949	}
1950	}
1951	///////////////////////////////////////////////////////////////////////////////
1952	static proc deg_D (intvec D,list PP)
1953	{
1954	// computes the degree of a divisor (intvec or list of integers)
1955	int i;
1956	int d=0;
1957	int s=size(D);
1958	for (i=1;i<=s;i=i+1)
1959	{
1960	d=d+D[i]*PP[i][1];
1961	}
1962	return(d);
1963	}
1964	///////////////////////////////////////////////////////////////////////////////
1965
1966	// ============================================================================
1967	// ***** MAIN PROCEDURES for the "preprocessing" of Brill-Noether ******
1968	// ============================================================================
1969
1970	proc Adj_div (poly f,list #)
1971	"USAGE: Adj_div( f [,l] ); f a poly, [l a list]
1972	RETURN: list L with the computed data:
1973	@format
1974	L[1] a list of rings: L[1][1]=aff_r (affine), L[1][2]=Proj_R (projective),
1975	L[2] an intvec with 2 entries (degree, genus),
1976	L[3] a list of intvec (closed places),
1977	L[4] an intvec (conductor),
1978	L[5] a list of lists:
1979	L[5][d][1] a (local) ring over an extension of degree d,
1980	L[5][d][2] an intvec (degrees of base points of places of degree d)
1981	@end format
1982	NOTE: @code{Adj_div(f);} computes and stores the fundamental data of the
1983	plane curve defined by f as needed for AG codes.
1984	In the affine ring you can find the following data:
1985	@format
1986	poly CHI: affine equation of the curve,
1987	ideal Aff_SLocus: affine singular locus (std),
1988	list Inf_Points: points at infinity
1989	Inf_Points[1]: singular points
1990	Inf_Points[2]: non-singular points,
1991	list Aff_SPoints: affine singular points (if not empty).
1992	@end format
1993	In the projective ring you can find the projective equation
1994	CHI of the curve (poly).
1995	In the local rings L[5][d][1] you find:
1996	@format
1997	list POINTS: base points of the places of degree d,
1998	list LOC_EQS: local equations of the curve at the base points,
1999	list BRANCHES: Hamburger-Noether developments of the places,
2000	list PARAMETRIZATIONS: local parametrizations of the places,
2001	@end format
2002	Each entry of the list L[3] corresponds to one closed place (i.e.,
2003	a place and all its conjugates) which is represented by an intvec
2004	of size two, the first entry is the degree of the place (in
2005	particular, it tells the local ring where to find the data
2006	describing one representative of the closed place), and the
2007	second one is the position of those data in the lists POINTS, etc.,
2008	inside this local ring.@*
2009	In the intvec L[4] (conductor) the i-th entry corresponds to the
2010	i-th entry in the list of places L[3].
2011
2012	With no optional arguments, the conductor is computed by
2013	local invariants of the singularities; otherwise it is computed
2014	by the Dedekind formula. @*
2015	An affine point is represented by a list P where P[1] is std
2016	of a prime ideal and P[2] is an intvec containing the position
2017	of the places above P in the list of closed places L[3]. @*
2018	If the point is at infinity, P[1] is a homogeneous irreducible
2019	polynomial in two variables.
2020
2021	If @code{printlevel>=0} additional comments are displayed (default:
2022	@code{printlevel=0}).
2023	WARNING: The parameter of the needed field extensions is called 'a'. Thus,
2024	there should be no global object named 'a' when executing Adj_div.
2025	KEYWORDS: Hamburger-Noether expansions; adjunction divisor
2026	SEE ALSO: closed_points, NSplaces
2027	EXAMPLE: example Adj_div; shows an example
2028	"
2029	{
2030	// computes the adjunction divisor and the genus of a (singular) plane curve
2031	// as a byproduct, it computes all the singular points with the corresponding
2032	// places and the genus of the curve
2033	// the adjunction divisor is stored in an intvec
2034	// also the non-singular places at infinity are computed
2035	// returns a list with all the computed data
2036	if (char(basering)==0)
2037	{
2038	ERROR("Base field not implemented");
2039	}
2040	if (npars(basering)>0)
2041	{
2042	ERROR("Base field not implemented");
2043	}
2044	if (defined(a)==1)
2045	{
2046	if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
2047	(typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
2048	{
2049	print("WARNING: There is an object named 'a' of the global type "
2050	+ typeof(a)+".");
2051	print("This will cause problems when "
2052	+ "executing Adj_div, as 'a' is used for "
2053	+ "the name of the parameter of the field extensions needed.");
2054	ERROR("Please rename or kill the object named 'a'");
2055	}
2056	}
2057	intvec opgt=option(get);
2058	option(redSB);
2059	def Base_R=basering;
2060	list CURVE=curve(f);
2061	def aff_r=CURVE[1];
2062	def Proj_R=CURVE[2];
2063	int degX=CURVE[3];
2064	int genusX=(degX-1)*(degX-2);
2065	genusX = genusX div 2;
2066	intvec iivv=degX,genusX;
2067	intvec Conductor;
2068	setring aff_r;
2069	dbprint(printlevel+1,"Computing affine singular points ... ");
2070	list Aff_SPoints=Aff_SL(Aff_SLocus);
2071	int s=size(Aff_SPoints);
2072	if (s>0)
2073	{
2074	export(Aff_SPoints);
2075	}
2076	dbprint(printlevel+1,"Computing all points at infinity ... ");
2077	list Inf_Points=inf_P(CHI);
2078	export(Inf_Points);
2079	list update_CURVE=list();
2080	update_CURVE[1]=list();
2081	update_CURVE[1][1]=aff_r;
2082	update_CURVE[1][2]=Proj_R;
2083	update_CURVE[2]=iivv;
2084	update_CURVE[3]=list();
2085	update_CURVE[4]=Conductor;
2086	update_CURVE[5]=list();
2087	int i;
2088	intvec pP=0,0,0;
2089	if (size(#)==0)
2090	{
2091	dbprint(printlevel+1,"Computing affine singular places ... ");
2092	if (s>0)
2093	{
2094	for (i=1;i<=s;i=i+1)
2095	{
2096	pP[3]=i;
2097	update_CURVE=place(pP,1,update_CURVE);
2098	}
2099	}
2100	dbprint(printlevel+1,"Computing singular places at infinity ... ");
2101	s=size(Inf_Points[1]);
2102	if (s>0)
2103	{
2104	pP[2]=1;
2105	for (i=1;i<=s;i=i+1)
2106	{
2107	pP[3]=i;
2108	update_CURVE=place(pP,1,update_CURVE);
2109	}
2110	}
2111	}
2112	else
2113	{
2114	dbprint(printlevel+1,"Computing affine singular places ... ");
2115	s=size(Aff_SPoints);
2116	if (s>0)
2117	{
2118	for (i=1;i<=s;i=i+1)
2119	{
2120	pP[3]=i;
2121	update_CURVE=place(pP,2,update_CURVE);
2122	}
2123	}
2124	dbprint(printlevel+1,"Computing singular places at infinity ... ");
2125	s=size(Inf_Points[1]);
2126	if (s>0)
2127	{
2128	pP[2]=1;
2129	for (i=1;i<=s;i=i+1)
2130	{
2131	pP[3]=i;
2132	update_CURVE=place(pP,2,update_CURVE);
2133	}
2134	}
2135	}
2136	dbprint(printlevel+1,"Computing non-singular places at infinity ... ");
2137	s=size(Inf_Points[2]);
2138	if (s>0)
2139	{
2140	pP[2]=2;
2141	for (i=1;i<=s;i=i+1)
2142	{
2143	pP[3]=i;
2144	update_CURVE=place(pP,0,update_CURVE);
2145	}
2146	}
2147	dbprint(printlevel+1,"Adjunction divisor computed successfully");
2148	list Places=update_CURVE[3];
2149	Conductor=update_CURVE[4];
2150	genusX = genusX - (deg_D(Conductor,Places) div 2);
2151	update_CURVE[2][2]=genusX;
2152	setring Base_R;
2153	dbprint(printlevel+1," ");
2154	dbprint(printlevel+2,"The genus of the curve is "+string(genusX));
2155	option(set,opgt);
2156	return(update_CURVE);
2157	}
2158	example
2159	{
2160	"EXAMPLE:"; echo = 2;
2161	int plevel=printlevel;
2162	printlevel=-1;
2163	ring s=2,(x,y),lp;
2164	list C=Adj_div(y9+y8+xy6+x2y3+y2+x3);
2165	def aff_R=C[1][1]; // the affine ring
2166	setring aff_R;
2167	listvar(aff_R); // data in the affine ring
2168	CHI; // affine equation of the curve
2169	Aff_SLocus; // ideal of the affine singular locus
2170	Aff_SPoints[1]; // 1st affine singular point: (1:1:1), no.1
2171	Inf_Points[1]; // singular point(s) at infinity: (1:0:0), no.4
2172	Inf_Points[2]; // list of non-singular points at infinity
2173	//
2174	pause("press RETURN");
2175	def proj_R=C[1][2]; // the projective ring
2176	setring proj_R;
2177	CHI; // projective equation of the curve
2178	C[2][1]; // degree of the curve
2179	C[2][2]; // genus of the curve
2180	C[3]; // list of computed places
2181	C[4]; // adjunction divisor (all points are singular!)
2182	//
2183	pause("press RETURN");
2184	// we look at the place(s) of degree 2 by changing to the ring
2185	C[5][2][1];
2186	def S(2)=C[5][2][1];
2187	setring S(2);
2188	POINTS; // base point(s) of place(s) of degree 2: (1:a:1)
2189	LOC_EQS; // local equation(s)
2190	PARAMETRIZATIONS; // parametrization(s) and exactness
2191	BRANCHES; // Hamburger-Noether development
2192	printlevel=plevel;
2193	}
2194	///////////////////////////////////////////////////////////////////////////////
2195
2196	proc NSplaces (intvec h,list CURVE)
2197	"USAGE: NSplaces( h, CURVE ), where h is an intvec and CURVE is a list
2198	RETURN: list L with updated data of CURVE after computing all non-singular
2199	affine closed places whose degrees are in the intvec h: @*
2200	@format
2201	in L[1][1]: (affine ring) lists Aff_Points(d) with affine non-singular
2202	(closed) points of degree d (if non-empty),
2203	in L[3]: the newly computed closed places are added,
2204	in L[5]: local rings created/updated to store (repres. of) new places.
2205	@end format
2206	See @ref{Adj_div} for a description of the entries in L.
2207	NOTE: The list_expression should be the output of the procedure Adj_div.@*
2208	Raising @code{printlevel}, additional comments are displayed
2209	(default: @code{printlevel=0}).
2210	WARNING: The parameter of the needed field extensions is called 'a'. Thus,
2211	there should be no global object named 'a' when executing NSplaces.
2212	SEE ALSO: closed_points, Adj_div
2213	EXAMPLE: example NSplaces; shows an example
2214	"
2215	{
2216	// computes all the non-singular affine (closed) places with fixed degrees
2217	// creates lists of points and the corresponding places
2218	// list CURVE must be the output of the procedure "Adj_div"
2219	// notice : if h==0 then nothing is done
2220	// warning : non-positive entries are forgotten
2221	if (defined(a)==1)
2222	{
2223	if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
2224	(typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
2225	{
2226	print("WARNING: There is an object named 'a' of the global type "
2227	+ typeof(a)+".");
2228	print("This will cause problems when "
2229	+ "executing Adj_div, as 'a' is used for "
2230	+ "the name of the parameter of the field extensions needed.");
2231	ERROR("Please rename or kill the object named 'a'");
2232	}
2233	}
2234	if (h==0)
2235	{
2236	return(CURVE);
2237	}
2238	intvec opgt=option(get);
2239	option(redSB);
2240	def Base_R=basering;
2241	def aff_r=CURVE[1][1];
2242	int M=size(h);
2243	list update_CURVE=CURVE;
2244	int i,j,l,s;
2245	setring aff_r;
2246	intvec pP=1,0,0;
2247	for (i=1;i<=M;i=i+1)
2248	{
2249	l=h[i];
2250	if (l>0)
2251	{
2252	dbprint(printlevel+1,"Computing non-singular affine places of degree "
2253	+string(l)+" ... ");
2254	if (defined(Aff_Points(l))==0)
2255	{
2256	list Aff_Points(l)=closed_points_deg(CHI,l,Aff_SLocus);
2257	s=size(Aff_Points(l));
2258	if (s>0)
2259	{
2260	export(Aff_Points(l));
2261	pP[2]=l;
2262	for (j=1;j<=s;j=j+1)
2263	{
2264	dbprint(printlevel-voice-2,"Place No.\ "+string(j));
2265	pP[3]=j;
2266	update_CURVE=place(pP,0,update_CURVE);
2267	}
2268	}
2269	}
2270	}
2271	}
2272	setring Base_R;
2273	option(set,opgt);
2274	return(update_CURVE);
2275	}
2276	example
2277	{
2278	"EXAMPLE:"; echo = 2;
2279	int plevel=printlevel;
2280	printlevel=-1;
2281	ring s=2,(x,y),lp;
2282	list C=Adj_div(x3y+y3+x);
2283	// The list of computed places:
2284	C[3];
2285	// create places up to degree 4
2286	list L=NSplaces(1..4,C);
2287	// The list of computed places is now:
2288	L[3];
2289	// e.g., affine non-singular points of degree 4 :
2290	def aff_r=L[1][1];
2291	setring aff_r;
2292	Aff_Points(4);
2293	// e.g., base point of the 1st place of degree 4 :
2294	def S(4)=L[5][4][1];
2295	setring S(4);
2296	POINTS[1];
2297	printlevel=plevel;
2298	}
2299	///////////////////////////////////////////////////////////////////////////////
2300
2301	// SPECIAL PROCEDURES FOR LINEAR ALGEBRA
2302
2303	static proc Ker (matrix A)
2304	{
2305	// warning : "lp" ordering is necessary
2306	intvec opgt=option(get);
2307	option(redSB);
2308	matrix M=transpose(syz(A));
2309	option(set,opgt);
2310	return(M);
2311	}
2312	///////////////////////////////////////////////////////////////////////////////
2313	static proc get_NZsol (matrix A)
2314	{
2315	matrix sol=Ker(A);
2316	return(submat(sol,1..1,1..ncols(sol)));
2317	}
2318	///////////////////////////////////////////////////////////////////////////////
2319	static proc supplement (matrix W,matrix V)
2320	"USAGE: supplement(W,V), where W,V are matrices of numbers such that the
2321	vector space generated by the rows of W is contained in that
2322	generated by the rows of V
2323	RETURN: matrix whose rows generate a supplementary vector space of W in V,
2324	or a zero row-matrix if <W>==<V>
2325	NOTE: W,V must be given with maximal rank
2326	"
2327	{
2328	// W and V represent independent sets of vectors and <W> is assumed to be
2329	// contained in <V>
2330	// computes matrix S whose rows are l.i. vectors s.t. <W> union <S> is a
2331	// basis of <V>
2332	// careful : the size of all vectors is assumed to be the same but it is
2333	// not checked and neither the linear independence of the vectors is checked
2334	// the trivial case W=0 is not covered by this procedure (and thus V<>0)
2335	// if <W>=<V> then a zero row-matrix is returned
2336	// warning : option(redSB) must be set in advance
2337	int n1=nrows(W);
2338	int n2=nrows(V);
2339	int s=n2-n1;
2340	if (s==0)
2341	{
2342	int n=ncols(W);
2343	matrix HH[1][n];
2344	return(HH);
2345	}
2346	matrix H=transpose(lift(transpose(V),transpose(W)));
2347	H=supplem(H);
2348	return(H*V);
2349	}
2350	///////////////////////////////////////////////////////////////////////////////
2351	static proc supplem (matrix M)
2352	"USAGE: supplem(M), where M is a matrix of numbers with maximal rank
2353	RETURN: matrix whose rows generate a supplementary vector space of <M> in
2354	k^n, where k is the base field and n is the number of columns
2355	SEE ALSO: supplement
2356	NOTE: The rank r is assumed to be 1<r<n.
2357	"
2358	{
2359	// warning : the linear independence of the rows is not checked
2360	int r=nrows(M);
2361	int n=ncols(M);
2362	int s=n-r;
2363	matrix A=M;
2364	matrix supl[s][n];
2365	int counter=0;
2366	int h=r+1;
2367	int i;
2368	for (i=1;i<=n;i=i+1)
2369	{
2370	matrix TT[1][n];
2371	TT[1,i]=1;
2372	A=transpose(concat(transpose(A),transpose(TT)));
2373	r=mat_rank(A);
2374	if (r==h)
2375	{
2376	h=h+1;
2377	counter=counter+1;
2378	supl=transpose(concat(transpose(supl),transpose(TT)));
2379	if (counter==s)
2380	{
2381	break;
2382	}
2383	}
2384	kill TT;
2385	}
2386	supl=transpose(compress(transpose(supl)));
2387	return(supl);
2388	}
2389	///////////////////////////////////////////////////////////////////////////////
2390	static proc mat_rank (matrix A)
2391	{
2392	// warning : "lp" ordering is necessary
2393	intvec opgt=option(get);
2394	option(redSB);
2395	int r=size(std(module(transpose(A))));
2396	option(set,opgt);
2397	return(r);
2398	}
2399	///////////////////////////////////////////////////////////////////////////////
2400
2401	// ***************************************************************
2402	// * PROCEDURES FOR INTERPOLATION, INTERSECTION AND EXTRA PLACES *
2403	// ***************************************************************
2404
2405	static proc estim_n (intvec Dplus,int dgX,list PL)
2406	{
2407	// computes an estimate for the degree n in the Brill-Noether algorithm
2408	int estim=2deg_D(Dplus,PL)+dgX(dgX-3);
2409	estim=estim div (2*dgX);
2410	estim=estim+1;
2411	if (estim<dgX)
2412	{
2413	estim=dgX;
2414	}
2415	return(estim);
2416	}
2417	///////////////////////////////////////////////////////////////////////////////
2418	static proc nforms (int n)
2419	{
2420	// computes the list of all homogeneous monomials of degree n>=0
2421	// exports ideal nFORMS(n) whose generators are ranged with lp order
2422	// in Proj_R and returns size(nFORMS(n))
2423	// warning : it is supposed to be called inside Proj_R
2424	// if n<=0 then nFORMS(0) is "computed fast"
2425	ideal nFORMS(n);
2426	int N;
2427	if (n>0)
2428	{
2429	N=(n+1)*(n+2);
2430	N=N div 2;
2431	N=N+1;
2432	int i,j,k;
2433	for (i=0;i<=n;i=i+1)
2434	{
2435	for (j=0;j<=n-i;j=j+1)
2436	{
2437	k=k+1;
2438	nFORMS(n)[N-k]=x^iy^jz^(n-i-j);
2439	}
2440	}
2441	export(nFORMS(n));
2442	}
2443	else
2444	{
2445	N=2;
2446	nFORMS(0)=1;
2447	export(nFORMS(0));
2448	}
2449	return(N-1);
2450	}
2451	///////////////////////////////////////////////////////////////////////////////
2452	static proc nmultiples (int n,int dgX,poly f)
2453	{
2454	// computes a basis of the space of forms of degree n which are multiple of
2455	// CHI
2456	// returns a matrix whose rows are the coordinates (related to nFORMS(n))
2457	// of such a basis
2458	// warning : it is supposed to be called inside Proj_R
2459	// warning : nFORMS(n) is created in the way, together with nFORMS(n-degX)
2460	// warning : n must be greater or equal than the degree of the curve
2461	if (defined(nFORMS(n))==0)
2462	{
2463	dbprint(printlevel+1,string(nforms(n)));
2464	}
2465	int m=n-dgX;
2466	if (defined(nFORMS(m))==0)
2467	{
2468	int k=nforms(m);
2469	}
2470	else
2471	{
2472	int k=size(nFORMS(m));
2473	}
2474	ideal nmults;
2475	int i;
2476	for (i=1;i<=k;i=i+1)
2477	{
2478	nmults[i]=f*nFORMS(m)[i];
2479	}
2480	return(transpose(lift(nFORMS(n),nmults)));
2481	}
2482	///////////////////////////////////////////////////////////////////////////////
2483	static proc interpolating_forms (intvec D,int n,list CURVE)
2484	{
2485	// computes a vector basis of the space of forms of degree n whose
2486	// intersection divisor with the curve is greater or equal than D>=0
2487	// the procedure is supposed to be called inside the ring Proj_R and
2488	// assumes that the forms nFORMS(n) are previously computed
2489	// the output is a matrix whose rows are the coordinates in nFORMS(n) of
2490	// such a basis
2491	// remark : the support of D may contain "extra" places
2492	def BR=basering;
2493	def aff_r=CURVE[1][1];
2494	int N=size(nFORMS(n));
2495	matrix totalM[1][N];
2496	int s=size(D);
2497	list Places=CURVE[3];
2498	int NPls=size(Places);
2499	int i,j,k,kk,id,ip,RR,ordx,ordy,nr,NR;
2500	if (s<=NPls)
2501	{
2502	for (i=1;i<=s;i=i+1)
2503	{
2504	if (D[i]>0)
2505	{
2506	id=Places[i][1];
2507	ip=Places[i][2];
2508	RR=D[i];
2509	def SS=CURVE[5][id][1];
2510	setring SS;
2511	poly xt=PARAMETRIZATIONS[ip][1][1];
2512	poly yt=PARAMETRIZATIONS[ip][1][2];
2513	ordx=PARAMETRIZATIONS[ip][2][1];
2514	ordy=PARAMETRIZATIONS[ip][2][2];
2515	nr=ncols(BRANCHES[ip][1]);
2516	if ( ordx>-1 \|\| ordy>-1 )
2517	{
2518	while ( ( RR>=ordx && ordx>-1 ) \|\| ( RR>=ordy && ordy>-1 ) )
2519	{
2520	BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
2521	nr=ncols(BRANCHES[ip][1]);
2522	PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
2523	xt=PARAMETRIZATIONS[ip][1][1];
2524	yt=PARAMETRIZATIONS[ip][1][2];
2525	ordx=PARAMETRIZATIONS[ip][2][1];
2526	ordy=PARAMETRIZATIONS[ip][2][2];
2527	}
2528	}
2529	if (POINTS[ip][3]==number(1))
2530	{
2531	number A=POINTS[ip][1];
2532	number B=POINTS[ip][2];
2533	map Mt=BR,A+xt,B+yt,1;
2534	kill A,B;
2535	}
2536	else
2537	{
2538	if (POINTS[ip][2]==number(1))
2539	{
2540	number A=POINTS[ip][1];
2541	map Mt=BR,A+xt,1,yt;
2542	kill A;
2543	}
2544	else
2545	{
2546	map Mt=BR,1,xt,yt;
2547	}
2548	}
2549	ideal nFORMS(n)=Mt(nFORMS(n));
2550	// rewrite properly the above conditions to obtain the local equations
2551	matrix partM[RR][N];
2552	matrix auxMC=coeffs(nFORMS(n),t);
2553	NR=nrows(auxMC);
2554	if (RR<=NR)
2555	{
2556	for (j=1;j<=RR;j=j+1)
2557	{
2558	for (k=1;k<=N;k=k+1)
2559	{
2560	partM[j,k]=number(auxMC[j,k]);
2561	}
2562	}
2563	}
2564	else
2565	{
2566	for (j=1;j<=NR;j=j+1)
2567	{
2568	for (k=1;k<=N;k=k+1)
2569	{
2570	partM[j,k]=number(auxMC[j,k]);
2571	}
2572	}
2573	for (j=NR+1;j<=RR;j=j+1)
2574	{
2575	for (k=1;k<=N;k=k+1)
2576	{
2577	partM[j,k]=number(0);
2578	}
2579	}
2580	}
2581	matrix localM=partM;
2582	matrix conjM=partM;
2583	for (j=2;j<=id;j=j+1)
2584	{
2585	conjM=Frobenius(conjM,1);
2586	localM=transpose(concat(transpose(localM),transpose(conjM)));
2587	}
2588	localM=gauss_row(localM);
2589	setring BR;
2590	totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2591	totalM=transpose(compress(transpose(totalM)));
2592	setring SS;
2593	kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
2594	setring BR;
2595	kill SS;
2596	}
2597	}
2598	}
2599	else
2600	{
2601	// distinguish between "standard" places and "extra" places
2602	for (i=1;i<=NPls;i=i+1)
2603	{
2604	if (D[i]>0)
2605	{
2606	id=Places[i][1];
2607	ip=Places[i][2];
2608	RR=D[i];
2609	def SS=CURVE[5][id][1];
2610	setring SS;
2611	poly xt=PARAMETRIZATIONS[ip][1][1];
2612	poly yt=PARAMETRIZATIONS[ip][1][2];
2613	ordx=PARAMETRIZATIONS[ip][2][1];
2614	ordy=PARAMETRIZATIONS[ip][2][2];
2615	nr=ncols(BRANCHES[ip][1]);
2616	if ( ordx>-1 \|\| ordy>-1 )
2617	{
2618	while ( ( RR>=ordx && ordx>-1 ) \|\| ( RR>=ordy && ordy>-1 ) )
2619	{
2620	BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
2621	nr=ncols(BRANCHES[ip][1]);
2622	PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
2623	xt=PARAMETRIZATIONS[ip][1][1];
2624	yt=PARAMETRIZATIONS[ip][1][2];
2625	ordx=PARAMETRIZATIONS[ip][2][1];
2626	ordy=PARAMETRIZATIONS[ip][2][2];
2627	}
2628	}
2629	if (POINTS[ip][3]==number(1))
2630	{
2631	number A=POINTS[ip][1];
2632	number B=POINTS[ip][2];
2633	map Mt=BR,A+xt,B+yt,1;
2634	kill A,B;
2635	}
2636	else
2637	{
2638	if (POINTS[ip][2]==number(1))
2639	{
2640	number A=POINTS[ip][1];
2641	map Mt=BR,A+xt,1,yt;
2642	kill A;
2643	}
2644	else
2645	{
2646	map Mt=BR,1,xt,yt;
2647	}
2648	}
2649	ideal nFORMS(n)=Mt(nFORMS(n));
2650	// rewrite properly the above conditions to obtain the local equations
2651	matrix partM[RR][N];
2652	matrix auxMC=coeffs(nFORMS(n),t);
2653	NR=nrows(auxMC);
2654	if (RR<=NR)
2655	{
2656	for (j=1;j<=RR;j=j+1)
2657	{
2658	for (k=1;k<=N;k=k+1)
2659	{
2660	partM[j,k]=number(auxMC[j,k]);
2661	}
2662	}
2663	}
2664	else
2665	{
2666	for (j=1;j<=NR;j=j+1)
2667	{
2668	for (k=1;k<=N;k=k+1)
2669	{
2670	partM[j,k]=number(auxMC[j,k]);
2671	}
2672	}
2673	for (j=NR+1;j<=RR;j=j+1)
2674	{
2675	for (k=1;k<=N;k=k+1)
2676	{
2677	partM[j,k]=number(0);
2678	}
2679	}
2680	}
2681	matrix localM=partM;
2682	matrix conjM=partM;
2683	for (j=2;j<=id;j=j+1)
2684	{
2685	conjM=Frobenius(conjM,1);
2686	localM=transpose(concat(transpose(localM),transpose(conjM)));
2687	}
2688	localM=gauss_row(localM);
2689	setring BR;
2690	totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2691	totalM=transpose(compress(transpose(totalM)));
2692	setring SS;
2693	kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
2694	setring BR;
2695	kill SS;
2696	}
2697	}
2698	k=s-NPls;
2699	int l;
2700	for (i=1;i<=k;i=i+1)
2701	{
2702	// in this case D[NPls+i]>0 is assumed to be true during the
2703	// Brill-Noether algorithm
2704	RR=D[NPls+i];
2705	setring aff_r;
2706	def SS=@EXTRA@[2][i];
2707	def extra_dgs=@EXTRA@[3];
2708	setring SS;
2709	poly xt=PARAMETRIZATION[1][1];
2710	poly yt=PARAMETRIZATION[1][2];
2711	ordx=PARAMETRIZATION[2][1];
2712	ordy=PARAMETRIZATION[2][2];
2713	nr=ncols(BRANCH[1]);
2714	if ( ordx>-1 \|\| ordy>-1 )
2715	{
2716	while ( ( RR>=ordx && ordx>-1 ) \|\| ( RR>=ordy && ordy>-1 ) )
2717	{
2718	BRANCH=extdevelop(BRANCH,2*nr);
2719	nr=ncols(BRANCH[1]);
2720	PARAMETRIZATION=param(BRANCH,0);
2721	xt=PARAMETRIZATION[1][1];
2722	yt=PARAMETRIZATION[1][2];
2723	ordx=PARAMETRIZATION[2][1];
2724	ordy=PARAMETRIZATION[2][2];
2725	}
2726	}
2727	number A=POINT[1];
2728	number B=POINT[2];
2729	map Mt=BR,A+xt,B+yt,1;
2730	kill A,B;
2731	ideal nFORMS(n)=Mt(nFORMS(n));
2732	// rewrite properly the above conditions to obtain the local equations
2733	matrix partM[RR][N];
2734	matrix auxMC=coeffs(nFORMS(n),t);
2735	NR=nrows(auxMC);
2736	if (RR<=NR)
2737	{
2738	for (j=1;j<=RR;j=j+1)
2739	{
2740	for (kk=1;kk<=N;kk=kk+1)
2741	{
2742	partM[j,kk]=number(auxMC[j,kk]);
2743	}
2744	}
2745	}
2746	else
2747	{
2748	for (j=1;j<=NR;j=j+1)
2749	{
2750	for (kk=1;kk<=N;kk=kk+1)
2751	{
2752	partM[j,kk]=number(auxMC[j,kk]);
2753	}
2754	}
2755	for (j=NR+1;j<=RR;j=j+1)
2756	{
2757	for (kk=1;kk<=N;kk=kk+1)
2758	{
2759	partM[j,kk]=number(0);
2760	}
2761	}
2762	}
2763	matrix localM=partM;
2764	matrix conjM=partM;
2765	l=extra_dgs[i];
2766	for (j=2;j<=l;j=j+1)
2767	{
2768	conjM=Frobenius(conjM,1);
2769	localM=transpose(concat(transpose(localM),transpose(conjM)));
2770	}
2771	localM=gauss_row(localM);
2772	setring BR;
2773	totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2774	totalM=transpose(compress(transpose(totalM)));
2775	setring SS;
2776	kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
2777	setring BR;
2778	kill SS;
2779	}
2780	}
2781	return(Ker(totalM));
2782	}
2783	///////////////////////////////////////////////////////////////////////////////
2784	static proc local_IN (poly h,int m)
2785	{
2786	// computes the intersection number of h and the curve CHI at a certain place
2787	// returns a list with the intersection number and the "leading coefficient"
2788	// the procedure must be called inside a local ring, h must be a local
2789	// equation with respect to the desired place, and m indicates the
2790	// number of place inside that local ring, containing lists
2791	// POINT(S)/BRANCH(ES)/PARAMETRIZATION(S) when m=0 an "extra place" is
2792	// considered
2793	def BR=basering;
2794	if (m>0)
2795	{
2796	int nr=ncols(BRANCHES[m][1]);
2797	poly xt=PARAMETRIZATIONS[m][1][1];
2798	poly yt=PARAMETRIZATIONS[m][1][2];
2799	int ordx=PARAMETRIZATIONS[m][2][1];
2800	int ordy=PARAMETRIZATIONS[m][2][2];
2801	map phi=BR,xt,yt,1;
2802	poly ht=phi(h);
2803	int inum=mindeg(ht);
2804	if ( ordx>-1 \|\| ordy>-1 )
2805	{
2806	while (((inum>ordx\|\|inum==-1) && ordx>-1) \|\| ((inum>ordy\|\|inum==-1) && ordy>-1))
2807	{
2808	BRANCHES[m]=extdevelop(BRANCHES[m],2*nr);
2809	nr=ncols(BRANCHES[m][1]);
2810	PARAMETRIZATIONS[m]=param(BRANCHES[m],0);
2811	xt=PARAMETRIZATIONS[m][1][1];
2812	yt=PARAMETRIZATIONS[m][1][2];
2813	ordx=PARAMETRIZATIONS[m][2][1];
2814	ordy=PARAMETRIZATIONS[m][2][2];
2815	phi=BR,xt,yt,1;
2816	ht=phi(h);
2817	inum=mindeg(ht);
2818	}
2819	}
2820	}
2821	else
2822	{
2823	int nr=ncols(BRANCH[1]);
2824	poly xt=PARAMETRIZATION[1][1];
2825	poly yt=PARAMETRIZATION[1][2];
2826	int ordx=PARAMETRIZATION[2][1];
2827	int ordy=PARAMETRIZATION[2][2];
2828	map phi=basering,xt,yt,1;
2829	poly ht=phi(h);
2830	int inum=mindeg(ht);
2831	if ( ordx>-1 \|\| ordy>-1 )
2832	{
2833	while (((inum>ordx\|\|inum==-1) && ordx>-1) \|\| ((inum>ordy\|\|inum==-1) && ordy>-1))
2834	{
2835	BRANCH=extdevelop(BRANCH,2*nr);
2836	nr=ncols(BRANCH[1]);
2837	PARAMETRIZATION=param(BRANCH,0);
2838	xt=PARAMETRIZATION[1][1];
2839	yt=PARAMETRIZATION[1][2];
2840	ordx=PARAMETRIZATION[2][1];
2841	ordy=PARAMETRIZATION[2][2];
2842	phi=BR,xt,yt,1;
2843	ht=phi(h);
2844	inum=mindeg(ht);
2845	}
2846	}
2847	}
2848	list answer=list();
2849	answer[1]=inum;
2850	number AA=number(coeffs(ht,t)[inum+1,1]);
2851	answer[2]=AA;
2852	return(answer);
2853	}
2854	///////////////////////////////////////////////////////////////////////////////
2855	static proc extra_place (ideal P)
2856	{
2857	// computes the "rational" place which is defined over a (closed) "extra"
2858	// point
2859	// an "extra" point will be necessarily affine, non-singular and non-rational
2860	// creates : a specific local ring to deal with the (unique) place above it
2861	// returns : list with the above local ring and the degree of the point/place
2862	// warning : the procedure must be called inside the affine ring aff_r
2863	def base_r=basering;
2864	int ext=deg(P[1]);
2865	poly aux=subst(P[2],y,1);
2866	ext=ext*deg(aux);
2867	// P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends only
2868	// on "y"
2869	if (deg(P[1])==1)
2870	{
2871	// the point is non-rational but the second component needs no field
2872	// extension
2873	number B=-number(subst(P[1],y,0));
2874	poly aux2=subst(P[2],y,B);
2875	// careful : the parameter will be called "a" anyway
2876	ring ES=(char(basering),a),(x,y,t),ls;
2877	map psi=base_r,a,0;
2878	minpoly=number(psi(aux2));
2879	kill psi;
2880	number A=a;
2881	number B=imap(base_r,B);
2882	}
2883	else
2884	{
2885	if (deg(subst(P[2],y,1))==1)
2886	{
2887	// the point is non-rational but the needed minpoly is just P[1]
2888	// careful : the parameter will be called "a" anyway
2889	poly P1=P[1];
2890	poly P2=P[2];
2891	ring ES=(char(basering),a),(x,y,t),ls;
2892	map psi=base_r,0,a;
2893	minpoly=number(psi(P1));
2894	kill psi;
2895	poly aux1=imap(base_r,P2);
2896	poly aux2=subst(aux1,y,a);
2897	number A=-number(subst(aux2,x,0));
2898	number B=a;
2899	}
2900	else
2901	{
2902	// this is the most complicated case of non-rational point
2903	poly P1=P[1];
2904	poly P2=P[2];
2905	int p=char(basering);
2906	int Q=p^ext;
2907	ring aux_r=(Q,a),(x,y,t),ls;
2908	string minpoly_string=string(minpoly);
2909	ring ES=(char(basering),a),(x,y,t),ls;
2910	execute("minpoly="+minpoly_string+";");
2911	poly P_1=imap(base_r,P1);
2912	poly P_2=imap(base_r,P2);
2913	ideal factors1=factorize(P_1,1);
2914	number B=-number(subst(factors1[1],y,0));
2915	poly P_0=subst(P_2,y,B);
2916	ideal factors2=factorize(P_0,1);
2917	number A=-number(subst(factors2[1],x,0));
2918	kill aux_r;
2919	}
2920	}
2921	list POINT=list();
2922	POINT[1]=A;
2923	POINT[2]=B;
2924	export(POINT);
2925	map phi=base_r,x+A,y+B;
2926	poly LOC_EQ=phi(CHI);
2927	kill A,B,phi;
2928	list LLL=hnexpansion(LOC_EQ,"ess");
2929	if (typeof(LLL[1])=="ring") {
2930	def altring=basering;
2931	def HNEring = LLL[1];
2932	setring HNEring;
2933	def L@HNE = hne[1];
2934	kill hne;
2935	export(L@HNE);
2936	int m=nrows(L@HNE[1]);
2937	int n=ncols(L@HNE[1]);
2938	intvec Li2=L@HNE[2];
2939	int Li3=L@HNE[3];
2940	setring ES;
2941	string newa=subfield(HNEring);
2942	poly paux=importdatum(HNEring,"L@HNE[4]",newa);
2943	matrix Maux[m][n];
2944	int i,j;
2945	for (i=1;i<=m;i=i+1)
2946	{
2947	for (j=1;j<=n;j=j+1)
2948	{
2949	Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
2950	string(j)+"]",newa);
2951	}
2952	}
2953	}
2954	else {
2955	def L@HNE=LLL[1][1];
2956	matrix Maux=L@HNE[1];
2957	intvec Li2=L@HNE[2];
2958	int Li3=L@HNE[3];
2959	poly paux=L@HNE[4];
2960	def HNEring=basering;
2961	setring HNEring;
2962	}
2963	kill LLL;
2964
2965	list BRANCH=list();
2966	BRANCH[1]=Maux;
2967	BRANCH[2]=Li2;
2968	BRANCH[3]=Li3;
2969	BRANCH[4]=paux;
2970	export(BRANCH);
2971	list PARAMETRIZATION=param(BRANCH,0);
2972	export(PARAMETRIZATION);
2973	kill HNEring;
2974	setring base_r;
2975	list answer=list();
2976	answer[1]=ES;
2977	answer[2]=ext;
2978	kill ES;
2979	return(answer);
2980	}
2981	///////////////////////////////////////////////////////////////////////////////
2982	static proc intersection_div (poly H,list CURVE)
2983	"USAGE: intersection_div(H,CURVE), where H is a homogeneous polynomial
2984	in ring Proj_R=p,(x,y,z),lp and CURVE is the list of data for
2985	the given curve
2986	CREATE: new places which had not been computed in advance if necessary
2987	those places are stored each one in a local ring where you find
2988	lists POINT,BRANCH,PARAMETRIZATION for the place living in that
2989	ring; the degree of the point/place in such a ring is stored in
2990	an intvec, and the base points in the remaining list
2991	Everything is exported in a list @EXTRA@ inside the ring
2992	aff_r=CURVE[1][1] and returned with the updated CURVE
2993	RETURN: list with the intersection divisor (intvec) between the underlying
2994	curve and H=0, and the list CURVE updated
2995	SEE ALSO: Adj_div, NSplaces, closed_points
2996	NOTE: The procedure must be called from the ring Proj_R=CURVE[1][2]
2997	(projective).
2998	If @EXTRA@ already exists, the new places are added to the
2999	previous data.
3000	"
3001	{
3002	// computes the intersection divisor of H and the curve CHI
3003	// returns a list with (possibly) "extra places" and it must be called
3004	// inside Proj_R
3005	// in case of extra places, some local rings ES(1) ... ES(m) are created
3006	// together with an integer list "extra_dgs" containing the degrees of
3007	// those places
3008	intvec opgt=option(get);
3009	option(redSB);
3010	intvec interdiv;
3011	def BRing=basering;
3012	int Tz1=deg(H);
3013	list Places=CURVE[3];
3014	int N=size(Places);
3015	def aff_r=CURVE[1][1];
3016	setring aff_r;
3017	if (defined(@EXTRA@)==0)
3018	{
3019	list @EXTRA@=list();
3020	list EP=list();
3021	list ES=list();
3022	list extra_dgs=list();
3023	}
3024	else
3025	{
3026	list EP=@EXTRA@[1];
3027	list ES=@EXTRA@[2];
3028	list extra_dgs=@EXTRA@[3];
3029	}
3030	int NN=size(extra_dgs);
3031	int counterEPl=0;
3032	setring BRing;
3033	poly h=subst(H,z,1);
3034	int Tz2=deg(h);
3035	int Tz3=Tz1-Tz2;
3036	int i,j,k,l,m,n,s,np,NP,I_N;
3037	if (Tz3==0)
3038	{
3039	// if this still does not work -> try always with ALL points in
3040	// Inf_Points !!!!
3041	poly Hinf=subst(H,z,0);
3042	setring aff_r;
3043	// compute the points at infinity of H and see which of them are in
3044	// Inf_Points
3045	poly h=imap(BRing,h);
3046	poly hinf=imap(BRing,Hinf);
3047	ideal pinf=factorize(hinf,1);
3048	list TIP=Inf_Points[1]+Inf_Points[2];
3049	s=size(TIP);
3050	NP=size(pinf);
3051	for (i=1;i<=NP;i=i+1)
3052	{
3053	for (j=1;j<=s;j=j+1)
3054	{
3055	if (pinf[i]==TIP[j][1])
3056	{
3057	np=size(TIP[j][2]);
3058	for (k=1;k<=np;k=k+1)
3059	{
3060	n=TIP[j][2][k];
3061	l=Places[n][1];
3062	m=Places[n][2];
3063	def SS=CURVE[5][l][1];
3064	setring SS;
3065	// local equation h of H
3066	if (POINTS[m][2]==number(1))
3067	{
3068	number A=POINTS[m][1];
3069	map psi=BRing,x+A,1,y;
3070	kill A;
3071	}
3072	else
3073	{
3074	map psi=BRing,1,x,y;
3075	}
3076	poly h=psi(H);
3077	I_N=local_IN(h,m)[1];
3078	interdiv[n]=I_N;
3079	kill h,psi;
3080	setring aff_r;
3081	kill SS;
3082	}
3083	break;
3084	}
3085	}
3086	}
3087	kill hinf,pinf;
3088	}
3089	else
3090	{
3091	// H is a multiple of z and hence all the points in Inf_Points intersect
3092	// with H
3093	setring aff_r;
3094	poly h=imap(BRing,h);
3095	list TIP=Inf_Points[1]+Inf_Points[2];
3096	s=size(TIP);
3097	for (j=1;j<=s;j=j+1)
3098	{
3099	np=size(TIP[j][2]);
3100	for (k=1;k<=np;k=k+1)
3101	{
3102	n=TIP[j][2][k];
3103	l=Places[n][1];
3104	m=Places[n][2];
3105	def SS=CURVE[5][l][1];
3106	setring SS;
3107	// local equation h of H
3108	if (POINTS[m][2]==number(1))
3109	{
3110	number A=POINTS[m][1];
3111	map psi=BRing,x+A,1,y;
3112	kill A;
3113	}
3114	else
3115	{
3116	map psi=BRing,1,x,y;
3117	}
3118	poly h=psi(H);
3119	I_N=local_IN(h,m)[1];
3120	interdiv[n]=I_N;
3121	kill h,psi;
3122	setring aff_r;
3123	kill SS;
3124	}
3125	}
3126	}
3127	// compute common affine points of H and CHI
3128	ideal CAL=h,CHI;
3129	CAL=std(CAL);
3130	if (CAL<>1)
3131	{
3132	list TAP=list();
3133	TAP=closed_points(CAL);
3134	NP=size(TAP);
3135	list auxP=list();
3136	int dP;
3137	for (i=1;i<=NP;i=i+1)
3138	{
3139	if (belongs(TAP[i],Aff_SLocus)==1)
3140	{
3141	// search the point in the list Aff_SPoints
3142	j=isInLP(TAP[i],Aff_SPoints);
3143	np=size(Aff_SPoints[j][2]);
3144	for (k=1;k<=np;k=k+1)
3145	{
3146	n=Aff_SPoints[j][2][k];
3147	l=Places[n][1];
3148	m=Places[n][2];
3149	def SS=CURVE[5][l][1];
3150	setring SS;
3151	// local equation h of H
3152	number A=POINTS[m][1];
3153	number B=POINTS[m][2];
3154	map psi=BRing,x+A,y+B,1;
3155	poly h=psi(H);
3156	I_N=local_IN(h,m)[1];
3157	interdiv[n]=I_N;
3158	kill A,B,h,psi;
3159	setring aff_r;
3160	kill SS;
3161	}
3162	}
3163	else
3164	{
3165	auxP=list();
3166	auxP[1]=TAP[i];
3167	dP=degree_P(auxP);
3168	if (defined(Aff_Points(dP))<>0)
3169	{
3170	// search the point in the list Aff_Points(dP)
3171	j=isInLP(TAP[i],Aff_Points(dP));
3172	n=Aff_Points(dP)[j][2][1];
3173	l=Places[n][1];
3174	m=Places[n][2];
3175	def SS=CURVE[5][l][1];
3176	setring SS;
3177	// local equation h of H
3178	number A=POINTS[m][1];
3179	number B=POINTS[m][2];
3180	map psi=BRing,x+A,y+B,1;
3181	poly h=psi(H);
3182	I_N=local_IN(h,m)[1];
3183	interdiv[n]=I_N;
3184	kill A,B,h,psi;
3185	setring aff_r;
3186	kill SS;
3187	}
3188	else
3189	{
3190	// previously check if it is an existing "extra place"
3191	j=isInLP(TAP[i],EP);
3192	if (j>0)
3193	{
3194	def SS=ES[j];
3195	setring SS;
3196	// local equation h of H
3197	number A=POINT[1];
3198	number B=POINT[2];
3199	map psi=BRing,x+A,y+B,1;
3200	poly h=psi(H);
3201	I_N=local_IN(h,0)[1];
3202	interdiv[N+j]=I_N;
3203	setring aff_r;
3204	kill SS;
3205	}
3206	else
3207	{
3208	// then we must create a new "extra place"
3209	counterEPl=counterEPl+1;
3210	list EXTRA_PLACE=extra_place(TAP[i]);
3211	def SS=EXTRA_PLACE[1];
3212	ES[NN+counterEPl]=SS;
3213	extra_dgs[NN+counterEPl]=EXTRA_PLACE[2];
3214	EP[NN+counterEPl]=list();
3215	EP[NN+counterEPl][1]=TAP[i];
3216	EP[NN+counterEPl][2]=0;
3217	setring SS;
3218	// local equation h of H
3219	number A=POINT[1];
3220	number B=POINT[2];
3221	map psi=BRing,x+A,y+B,1;
3222	poly h=psi(H);
3223	I_N=local_IN(h,0)[1];
3224	kill A,B,h,psi;
3225	interdiv[N+NN+counterEPl]=I_N;
3226	setring aff_r;
3227	kill SS;
3228	}
3229	}
3230	}
3231	}
3232	kill TAP,auxP;
3233	}
3234	kill h,CAL,TIP;
3235	@EXTRA@[1]=EP;
3236	@EXTRA@[2]=ES;
3237	@EXTRA@[3]=extra_dgs;
3238	kill EP;
3239	list update_CURVE=CURVE;
3240	if (size(extra_dgs)>0)
3241	{
3242	export(@EXTRA@);
3243	update_CURVE[1][1]=basering;
3244	}
3245	else
3246	{
3247	kill @EXTRA@;
3248	}
3249	setring BRing;
3250	kill h;
3251	kill aff_r;
3252	list answer=list();
3253	answer[1]=interdiv;
3254	answer[2]=update_CURVE;
3255	option(set,opgt);
3256	return(answer);
3257	}
3258	///////////////////////////////////////////////////////////////////////////////
3259	static proc local_eq (poly H,SS,int m)
3260	{
3261	// computes a local equation of polynomial H in the ring SS related to the place
3262	// "m"
3263	// list POINT/POINTS is searched depending on wether m=0 or m>0 respectively
3264	// warning : the procedure must be called from ring "Proj_R" and returns a
3265	// string
3266	def BRing=basering;
3267	setring SS;
3268	if (m>0)
3269	{
3270	if (POINTS[m][3]==number(1))
3271	{
3272	number A=POINTS[m][1];
3273	number B=POINTS[m][2];
3274	map psi=BRing,x+A,y+B,1;
3275	}
3276	else
3277	{
3278	if (POINTS[m][2]==number(1))
3279	{
3280	number A=POINTS[m][1];
3281	map psi=BRing,x+A,1,y;
3282	}
3283	else
3284	{
3285	map psi=BRing,1,x,y;
3286	}
3287	}
3288	}
3289	else
3290	{
3291	number A=POINT[1];
3292	number B=POINT[2];
3293	map psi=BRing,x+A,y+B,1;
3294	}
3295	poly h=psi(H);
3296	string str_h=string(h);
3297	setring BRing;
3298	return(str_h);
3299	}
3300	///////////////////////////////////////////////////////////////////////////////
3301	static proc min_wt_rmat (matrix M)
3302	{
3303	// finds the row of M with minimum non-zero entries, i.e. minimum
3304	// "Hamming-weight"
3305	int m=nrows(M);
3306	int n=ncols(M);
3307	int i,j;
3308	int Hwt=0;
3309	for (j=1;j<=n;j=j+1)
3310	{
3311	if (M[1,j]<>0)
3312	{
3313	Hwt=Hwt+1;
3314	}
3315	}
3316	int minHwt=Hwt;
3317	int k=1;
3318	for (i=2;i<=m;i=i+1)
3319	{
3320	Hwt=0;
3321	for (j=1;j<=n;j=j+1)
3322	{
3323	if (M[i,j]<>0)
3324	{
3325	Hwt=Hwt+1;
3326	}
3327	}
3328	if (Hwt<minHwt)
3329	{
3330	minHwt=Hwt;
3331	k=i;
3332	}
3333	}
3334	return(k);
3335	}
3336	///////////////////////////////////////////////////////////////////////////////
3337
3338	// ============================================================================
3339	// ***** MAIN PROCEDURE : the Brill-Noether algorithm ******
3340	// ============================================================================
3341
3342	proc BrillNoether (intvec G,list CURVE)
3343	"USAGE: BrillNoether(G,CURVE); G an intvec, CURVE a list
3344	RETURN: list of ideals (each of them with two homogeneous generators,
3345	which represent the numerator, resp. denominator, of a rational
3346	function).@*
3347	The corresponding rational functions form a vector space basis of the
3348	linear system L(G), G a rational divisor over a non-singular curve.
3349	NOTE: The procedure must be called from the ring CURVE[1][2], where
3350	CURVE is the output of the procedure @code{NSplaces}. @*
3351	The intvec G represents a rational divisor supported on the closed
3352	places of CURVE[3] (e.g. @code{G=2,0,-1;} means 2 times the closed
3353	place 1 minus 1 times the closed place 3).
3354	SEE ALSO: Adj_div, NSplaces, Weierstrass
3355	EXAMPLE: example BrillNoether; shows an example
3356	"
3357	{
3358	// computes a vector basis for the space L(G),
3359	// where G is a given rational divisor over the non-singular curve
3360	// returns : list of ideals in R each with 2 elements H,Ho such that
3361	// the set of functions {H/Ho} is the searched basis
3362	// warning : the conductor and sufficiently many points of the plane
3363	// curve should be computed in advance, in list CURVE
3364	// the algorithm of Brill-Noether is carried out in the procedure
3365	def BRing=basering;
3366	int degX=CURVE[2][1];
3367	list Places=CURVE[3];
3368	intvec Conductor=CURVE[4];
3369	if (deg_D(G,Places)<0)
3370	{
3371	return(list());
3372	}
3373	intvec nuldiv;
3374	if (G==nuldiv)
3375	{
3376	list quickL=list();
3377	ideal quickId;
3378	quickId[1]=1;
3379	quickId[2]=1;
3380	quickL[1]=quickId;
3381	return(quickL);
3382	}
3383	intvec J=max_D(G,nuldiv)+Conductor;
3384	int n=estim_n(J,degX,Places);
3385	dbprint(printlevel+1,"Forms of degree "+string(n)+" : ");
3386	matrix W=nmultiples(n,degX,CHI);
3387	kill nFORMS(n-degX);
3388	list update_CURVE=CURVE;
3389	matrix V=interpolating_forms(J,n,update_CURVE);
3390	matrix VmW=supplement(W,V);
3391	int k=min_wt_rmat(VmW);
3392	int N=size(nFORMS(n));
3393	matrix H0[1][N];
3394	int i,j;
3395	for (i=1;i<=N;i=i+1) { H0[1,i]=VmW[k,i]; }
3396	poly Ho;
3397	for (i=1;i<=N;i=i+1) { Ho=Ho+(H0[1,i]*nFORMS(n)[i]); }
3398	list INTERD=intersection_div(Ho,update_CURVE);
3399	intvec NHo=INTERD[1];
3400	update_CURVE=INTERD[2];
3401	intvec AR=NHo-G;
3402	matrix V2=interpolating_forms(AR,n,update_CURVE);
3403	def aux_RING=update_CURVE[1][1];
3404	setring aux_RING;
3405	if (defined(@EXTRA@)<>0)
3406	{
3407	kill @EXTRA@;
3408	}
3409	setring BRing;
3410	update_CURVE[1][1]=aux_RING;
3411	kill aux_RING;
3412	matrix B0=supplement(W,V2);
3413	if (Hamming_wt(B0)==0)
3414	{
3415	return(list());
3416	}
3417	int ld=nrows(B0);
3418	list Bres=list();
3419	ideal HH;
3420	poly H;
3421	for (j=1;j<=ld;j=j+1)
3422	{
3423	H=0;
3424	for (i=1;i<=N;i=i+1)
3425	{
3426	H=H+(B0[j,i]*nFORMS(n)[i]);
3427	}
3428	HH=H,Ho;
3429	Bres[j]=simplifyRF(HH);
3430	}
3431	kill nFORMS(n);
3432	dbprint(printlevel+1," ");
3433	dbprint(printlevel+2,"Vector basis successfully computed ");
3434	dbprint(printlevel+1," ");
3435	return(Bres);
3436	}
3437	example
3438	{
3439	"EXAMPLE:"; echo = 2;
3440	int plevel=printlevel;
3441	printlevel=-1;
3442	ring s=2,(x,y),lp;
3443	list C=Adj_div(x3y+y3+x);
3444	C=NSplaces(1..4,C);
3445	// the first 3 Places in C[3] are of degree 1.
3446	// we define the rational divisor G = 4C[3][1]+4C[3][3] (of degree 8):
3447	intvec G=4,0,4;
3448	def R=C[1][2];
3449	setring R;
3450	list LG=BrillNoether(G,C);
3451	// here is the vector basis of L(G):
3452	LG;
3453	printlevel=plevel;
3454	}
3455	///////////////////////////////////////////////////////////////////////////////
3456
3457	// * procedures for dealing with "RATIONAL FUNCTIONS" over a plane curve *
3458	// a rational function F may be given by (homogeneous) ideal or (affine) poly
3459	// (or number)
3460
3461	static proc polytoRF (F)
3462	{
3463	// converts a polynomial (or number) into a "rational function" of type "ideal"
3464	// warning : it must be called inside "R" and the polynomial is expected to be affine
3465	ideal RF;
3466	RF[1]=homog(F,z);
3467	RF[2]=z^(deg(F));
3468	return(RF);
3469	}
3470	///////////////////////////////////////////////////////////////////////////////
3471	static proc simplifyRF (ideal F)
3472	{
3473	// simplifies a rational function f/g extracting the gcd(f,g)
3474	// maybe add a "restriction" to the curve "CHI" but it is not easy to
3475	// programm
3476	poly auxp=gcd(F[1],F[2]);
3477	return(ideal(division(F,auxp)[1]));
3478	}
3479	///////////////////////////////////////////////////////////////////////////////
3480	static proc sumRF (F,G)
3481	{
3482	// sum of two "rational functions" F,G given by either a pair
3483	// numerator/denominator or a poly
3484	if ( typeof(F)=="ideal" && typeof(G)=="ideal" )
3485	{
3486	ideal FG;
3487	FG[1]=F[1]G[2]+F[2]G[1];
3488	FG[2]=F[2]*G[2];
3489	return(simplifyRF(FG));
3490	}
3491	else
3492	{
3493	if (typeof(F)=="ideal")
3494	{
3495	ideal GG=polytoRF(G);
3496	ideal FG;
3497	FG[1]=F[1]GG[2]+F[2]GG[1];
3498	FG[2]=F[2]*GG[2];
3499	return(simplifyRF(FG));
3500	}
3501	else
3502	{
3503	if (typeof(G)=="ideal")
3504	{
3505	ideal FF=polytoRF(F);
3506	ideal FG;
3507	FG[1]=FF[1]G[2]+FF[2]G[1];
3508	FG[2]=FF[2]*G[2];
3509	return(simplifyRF(FG));
3510	}
3511	else
3512	{
3513	return(F+G);
3514	}
3515	}
3516	}
3517	}
3518	///////////////////////////////////////////////////////////////////////////////
3519	static proc negRF (F)
3520	{
3521	// returns -F as rational function
3522	if (typeof(F)=="ideal")
3523	{
3524	ideal FF=F;
3525	FF[1]=-F[1];
3526	return(FF);
3527	}
3528	else
3529	{
3530	return(-F);
3531	}
3532	}
3533	///////////////////////////////////////////////////////////////////////////////
3534	static proc escprodRF (l,F)
3535	{
3536	// computes l*F as rational function
3537	// l should be either a number or a polynomial of degree zero
3538	if (typeof(F)=="ideal")
3539	{
3540	ideal lF=F;
3541	lF[1]=l*F[1];
3542	return(lF);
3543	}
3544	else
3545	{
3546	return(l*F);
3547	}
3548	}
3549	///////////////////////////////////////////////////////////////////////////////
3550
3551	// ****** procedures to compute Weierstrass semigroups ******
3552
3553	static proc orderRF (ideal F,SS,int m)
3554	"USAGE: orderRF(F,SS,m); F an ideal, SS a ring and m an integer
3555	RETURN: list with the order (int) and the leading coefficient (number)
3556	NOTE: F represents a rational function, thus the procedure must be
3557	called from global ring R or R(d).
3558	SS contains the name of a local ring where rational places are
3559	stored, and then we take that which is in position m in the
3560	corresponding lists of data.
3561	The order of F at the place given by SS,m is returned together
3562	with the coefficient of minimum degree in the corresponding power
3563	series.
3564	"
3565	{
3566	// computes the order of a rational function F at a RATIONAL place given by
3567	// a local ring SS and a position "m" inside SS
3568	// warning : the procedure must be called from global projective ring "R" or
3569	// "R(i)"
3570	// returns a list with the order (int) and the "leading coefficient" (number)
3571	def BR=basering;
3572	poly f=F[1];
3573	string sf=local_eq(f,SS,m);
3574	poly g=F[2];
3575	string sg=local_eq(g,SS,m);
3576	setring SS;
3577	execute("poly ff="+sf+";");
3578	execute("poly gg="+sg+";");
3579	list o1=local_IN(ff,m);
3580	list o2=local_IN(gg,m);
3581	int oo=o1[1]-o2[1];
3582	number lc=o1[2]/o2[2];
3583	setring BR;
3584	number LC=number(imap(SS,lc));
3585	return(list(oo,LC));
3586	}
3587	///////////////////////////////////////////////////////////////////////////////
3588
3589	proc Weierstrass (int P,int m,list CURVE)
3590	"USAGE: Weierstrass( i, m, CURVE ); i,m integers and CURVE a list
3591	RETURN: list WS of two lists:
3592	@format
3593	WS[1] list of integers (Weierstr. semigroup of the curve at place i up to m)
3594	WS[2] list of ideals (the associated rational functions)
3595	@end format
3596	NOTE: The procedure must be called from the ring CURVE[1][2],
3597	where CURVE is the output of the procedure @code{NSplaces}.
3598	@* i represents the place CURVE[3][i].
3599	@* Rational functions are represented by numerator/denominator
3600	in form of ideals with two homogeneous generators.
3601	WARNING: The place must be rational, i.e., necessarily CURVE[3][i][1]=1. @*
3602	SEE ALSO: Adj_div, NSplaces, BrillNoether
3603	EXAMPLE: example Weierstrass; shows an example
3604	"
3605	{
3606	// computes the Weierstrass semigroup at a RATIONAL place P up to a bound "m"
3607	// together with the functions achieving each value up to m, via
3608	// Brill-Noether
3609	// returns 2 lists : the first consists of all the poles up to m in
3610	// increasing order and the second consists of the corresponging rational
3611	// functions
3612	list Places=CURVE[3];
3613	intvec pl=Places[P];
3614	int dP=pl[1];
3615	int nP=pl[2];
3616	if (dP<>1)
3617	{
3618	ERROR("the given place is not defined over the prime field");
3619	}
3620	if (m<=0)
3621	{
3622	if (m==0)
3623	{
3624	list semig=list();
3625	int auxint=0;
3626	semig[1]=auxint;
3627	list funcs=list();
3628	ideal auxF;
3629	auxF[1]=1;
3630	auxF[2]=1;
3631	funcs[1]=auxF;
3632	return(list(semig,funcs));
3633	}
3634	else
3635	{
3636	ERROR("second argument must be non-negative");
3637	}
3638	}
3639	int auxint=0;
3640	ideal auxF;
3641	auxF[1]=1;
3642	auxF[2]=1;
3643	// Brill-Noether algorithm
3644	intvec mP;
3645	mP[P]=m;
3646	list LmP=BrillNoether(mP,CURVE);
3647	int lmP=size(LmP);
3648	if (lmP==1)
3649	{
3650	list semig=list();
3651	semig[1]=auxint;
3652	list funcs=list();
3653	funcs[1]=auxF;
3654	return(list(semig,funcs));
3655	}
3656	def SS=CURVE[5][dP][1];
3657	list ordLmP=list();
3658	list polLmP=list();
3659	list sortpol=list();
3660	int maxpol;
3661	int i,j,k;
3662	for (i=1;i<=lmP-1;i=i+1)
3663	{
3664	for (j=1;j<=lmP-i+1;j=j+1)
3665	{
3666	ordLmP[j]=orderRF(LmP[j],SS,nP);
3667	polLmP[j]=-ordLmP[j][1];
3668	}
3669	sortpol=sort(polLmP);
3670	polLmP=sortpol[1];
3671	maxpol=polLmP[lmP-i+1];
3672	LmP=permute_L(LmP,sortpol[2]);
3673	ordLmP=permute_L(ordLmP,sortpol[2]);
3674	// triangulate LmP
3675	for (k=1;k<=lmP-i;k=k+1)
3676	{
3677	if (polLmP[lmP-i+1-k]==maxpol)
3678	{
3679	LmP[lmP-i+1-k]=sumRF(LmP[lmP-i+1-k],negRF(escprodRF(
3680	ordLmP[lmP-i+1-k][2]/ordLmP[lmP-i+1][2],LmP[lmP-i+1])));
3681	}
3682	else
3683	{
3684	break;
3685	}
3686	}
3687	}
3688	polLmP[1]=auxint;
3689	LmP[1]=auxF;
3690	return(list(polLmP,LmP));
3691	}
3692	example
3693	{
3694	"EXAMPLE:"; echo = 2;
3695	int plevel=printlevel;
3696	printlevel=-1;
3697	ring s=2,(x,y),lp;
3698	list C=Adj_div(x3y+y3+x);
3699	C=NSplaces(1..4,C);
3700	def R=C[1][2];
3701	setring R;
3702	// Place C[3][1] has degree 1 (i.e it is rational);
3703	list WS=Weierstrass(1,7,C);
3704	// the first part of the list is the Weierstrass semigroup up to 7 :
3705	WS[1];
3706	// and the second part are the corresponding functions :
3707	WS[2];
3708	printlevel=plevel;
3709	}
3710	///////////////////////////////////////////////////////////////////////////////
3711
3712	// axiliary procedure for permuting a list or intvec
3713
3714	proc permute_L (L,P)
3715	"USAGE: permute_L( L, P ); L,P either intvecs or lists
3716	RETURN: list obtained from L by applying the permutation given by P.
3717	NOTE: If P is a list, all entries must be integers.
3718	SEE ALSO: sys_code, AGcode_Omega, prepSV
3719	EXAMPLE: example permute_L; shows an example
3720	"
3721	{
3722	// permutes the list L according to the permutation P (both intvecs or
3723	// lists of integers)
3724	int s=size(L);
3725	int n=size(P);
3726	int i;
3727	if (s<n)
3728	{
3729	for (i=s+1;i<=n;i=i+1)
3730	{
3731	L[i]=0;
3732	}
3733	s=size(L);
3734	}
3735	list auxL=L;
3736	for (i=1;i<=n;i=i+1)
3737	{
3738	auxL[i]=L[P[i]];
3739	}
3740	return(auxL);
3741	}
3742	example
3743	{
3744	"EXAMPLE:"; echo = 2;
3745	list L=list();
3746	L[1]="a";
3747	L[2]="b";
3748	L[3]="c";
3749	L[4]="d";
3750	intvec P=1,3,4,2;
3751	// the list L is permuted according to P :
3752	permute_L(L,P);
3753	}
3754	///////////////////////////////////////////////////////////////////////////////
3755	static proc evalRF (ideal F,SS,int m)
3756	"USAGE: evalRF(F,SS,m), where F is an ideal, SS is a ring and m is an
3757	integer
3758	RETURN: the evaluation (number) of F at the place given by SS,m if it is
3759	well-defined
3760	NOTE: F represents a rational function, thus the procedure must be
3761	called from R or R(d).
3762	SS contains the name of a local ring where rational places are
3763	stored, and then we take that which is in position m in the
3764	corresponding lists of data.
3765	"
3766	{
3767	// evaluates a rational function F at a RATIONAL place given by
3768	// a local ring SS and a position "m" inside SS
3769	list olc=orderRF(F,SS,m);
3770	int oo=olc[1];
3771	if (oo==0)
3772	{
3773	return(olc[2]);
3774	}
3775	else
3776	{
3777	if (oo>0)
3778	{
3779	return(number(0));
3780	}
3781	else
3782	{
3783	ERROR("the function is not well-defined at the given place");
3784	}
3785	}
3786	}
3787	///////////////////////////////////////////////////////////////////////////////
3788	//
3789	// ****** procedures for constructing AG codes ******
3790	//
3791	///////////////////////////////////////////////////////////////////////////////
3792
3793	static proc gen_mat (list LF,intvec LP,RP)
3794	"USAGE: gen_mat(LF,LP,RP); LF list of rational functions,
3795	LP intvec of rational places and RP a local ring
3796	RETURN: a generator matrix of the evaluation code given by LF and LP
3797	SEE ALSO: extcurve
3798	KEYWORDS: evaluation codes
3799	NOTE: Rational places are searched in local ring RP
3800	The procedure must be called from R or R(d) fromlist CURVE
3801	after having executed extcurve(d,CURVE)
3802	"
3803	{
3804	// computes a generator matrix (with numbers) of the evaluation code given
3805	// by a list of rational functions LF and a list of RATIONAL places LP
3806	int m=size(LF);
3807	int n=size(LP);
3808	matrix GM[m][n];
3809	int i,j;
3810	for (i=1;i<=m;i=i+1)
3811	{
3812	for (j=1;j<=n;j=j+1)
3813	{
3814	GM[i,j]=evalRF(LF[i],RP,LP[j]);
3815	}
3816	}
3817	return(GM);
3818	}
3819	///////////////////////////////////////////////////////////////////////////////
3820
3821	proc dual_code (matrix G)
3822	"USAGE: dual_code(G); G a matrix of numbers
3823	RETURN: a generator matrix of the dual code generated by G
3824	NOTE: The input should be a matrix G of numbers. @*
3825	The output is also a parity check matrix for the code defined by G
3826	KEYWORDS: linear code, dual
3827	EXAMPLE: example dual_code; shows an example
3828	"
3829	{
3830	// computes the dual code of C given by a generator matrix G
3831	// i.e. computes a parity-check matrix H of C
3832	// conversely : computes also G if H is given
3833	return(Ker(G));
3834	}
3835	example
3836	{
3837	"EXAMPLE:"; echo = 2;
3838	ring s=2,T,lp;
3839	// here is the Hamming code of length 7 and dimension 3
3840	matrix G[3][7]=1,0,1,0,1,0,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1;
3841	print(G);
3842	matrix H=dual_code(G);
3843	print(H);
3844	}
3845	///////////////////////////////////////////////////////////////////////////////
3846
3847	// ======================================================================
3848	// ********* initial test for disjointness *************
3849	// ======================================================================
3850
3851	static proc disj_divs (intvec H,intvec P,list EC)
3852	{
3853	int s1=size(H);
3854	int s2=size(P);
3855	list PLACES=EC[3];
3856	def BRing=basering;
3857	def auxR=EC[1][5];
3858	setring auxR;
3859	int s=res_deg();
3860	setring BRing;
3861	kill auxR;
3862	int i,j,k,d,l,N,M;
3863	intvec auxIV,iw;
3864	for (i=1;i<=s;i=i+1)
3865	{
3866	if ( (s mod i) == 0 )
3867	{
3868	if (typeof(EC[5][i])=="list")
3869	{
3870	def auxR=EC[5][i][1];
3871	setring auxR;
3872	auxIV[i]=size(POINTS);
3873	setring BRing;
3874	kill auxR;
3875	}
3876	else
3877	{
3878	auxIV[i]=0;
3879	}
3880	}
3881	else
3882	{
3883	auxIV[i]=0;
3884	}
3885	}
3886	for (i=1;i<=s1;i=i+1)
3887	{
3888	if (H[i]<>0)
3889	{
3890	iw=PLACES[i];
3891	d=iw[1];
3892	if ( (s mod d) == 0 )
3893	{
3894	l=iw[2];
3895	// check that this place(s) is/are not in sup(D)
3896	if (d==1)
3897	{
3898	for (j=1;j<=s2;j=j+1)
3899	{
3900	if (l==P[j])
3901	{
3902	return(0);
3903	}
3904	}
3905	}
3906	else
3907	{
3908	N=0;
3909	for (j=1;j<d;j=j+1)
3910	{
3911	N=N+j*auxIV[j];
3912	}
3913	N=N+d*(l-1);
3914	M=N+d;
3915	for (k=N+1;k<=M;k=k+1)
3916	{
3917	for (j=1;j<=s2;j=j+1)
3918	{
3919	if (k==P[j])
3920	{
3921	return(0);
3922	}
3923	}
3924	}
3925	}
3926	}
3927	}
3928	}
3929	kill auxIV,iw;
3930	return(1);
3931	}
3932	///////////////////////////////////////////////////////////////////////////////
3933
3934	proc AGcode_L (intvec G,intvec D,list EC)
3935	"USAGE: AGcode_L( G, D, EC ); G,D intvec, EC a list
3936	RETURN: a generator matrix for the evaluation AG code defined by the
3937	divisors G and D.
3938	NOTE: The procedure must be called within the ring EC[1][4],
3939	where EC is the output of @code{extcurve(d)} (or within
3940	the ring EC[1][2] if d=1). @*
3941	The entry i in the intvec D refers to the i-th rational
3942	place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
3943	The intvec G represents a rational divisor (see @ref{BrillNoether}
3944	for more details).@*
3945	The code evaluates the vector space basis of L(G) at the rational
3946	places given by D.
3947	WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
3948	not checked by the algorithm.
3949	G and D should have disjoint supports (checked by the algorithm).
3950	SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_Omega
3951	EXAMPLE: example AGcode_L; shows an example
3952	"
3953	{
3954	// returns a generator matrix for the evaluation AG code given by G and D
3955	// G must be a divisor defined over the prime field and D an intvec of
3956	// "rational places"
3957	// it must be called inside R or R(d) and requires previously "extcurve(d)"
3958	def BR=basering;
3959	if (disj_divs(G,D,EC)==0)
3960	{
3961	dbprint(printlevel+3,"? the divisors do not have disjoint supports,
3962	0-matrix returned ?");
3963	matrix answer;
3964	return(answer);
3965	}
3966	if (res_deg()>1)
3967	{
3968	def R=EC[1][2];
3969	setring R;
3970	list LG=BrillNoether(G,EC);
3971	setring BR;
3972	list LG=imap(R,LG);
3973	setring R;
3974	kill LG;
3975	setring BR;
3976	kill R;
3977	}
3978	else
3979	{
3980	list LG=BrillNoether(G,EC);
3981	}
3982	def RP=EC[1][5];
3983	matrix M=gen_mat(LG,D,RP);
3984	kill LG;
3985	return(M);
3986	}
3987	example
3988	{
3989	"EXAMPLE:"; echo = 2;
3990	int plevel=printlevel;
3991	printlevel=-1;
3992	ring s=2,(x,y),lp;
3993	list HC=Adj_div(x3+y2+y);
3994	HC=NSplaces(1..2,HC);
3995	HC=extcurve(2,HC);
3996	def ER=HC[1][4];
3997	setring ER;
3998	intvec G=5; // the rational divisor G = 5*HC[3][1]
3999	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4000	// let us construct the corresponding evaluation AG code :
4001	matrix C=AGcode_L(G,D,HC);
4002	// here is a linear code of type [8,5,>=3] over F_4
4003	print(C);
4004	printlevel=plevel;
4005	}
4006	///////////////////////////////////////////////////////////////////////////////
4007
4008	proc AGcode_Omega (intvec G,intvec D,list EC)
4009	"USAGE: AGcode_Omega( G, D, EC ); G,D intvec, EC a list
4010	RETURN: a generator matrix for the residual AG code defined by the
4011	divisors G and D.
4012	NOTE: The procedure must be called within the ring EC[1][4],
4013	where EC is the output of @code{extcurve(d)} (or within
4014	the ring EC[1][2] if d=1). @*
4015	The entry i in the intvec D refers to the i-th rational
4016	place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
4017	The intvec G represents a rational divisor (see @ref{BrillNoether}
4018	for more details).@*
4019	The code computes the residues of a vector space basis of
4020	@math{\Omega(G-D)} at the rational places given by D.
4021	WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
4022	not checked by the algorithm.
4023	G and D should have disjoint supports (checked by the algorithm).
4024	SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_L
4025	EXAMPLE: example AGcode_Omega; shows an example
4026	"
4027	{
4028	// returns a generator matrix for the residual AG code given by G and D
4029	// G must be a divisor defined over the prime field and D an intvec or
4030	// "rational places"
4031	// it must be called inside R or R(d) and requires previously "extcurve(d)"
4032	return(dual_code(AGcode_L(G,D,EC)));
4033	}
4034	example
4035	{
4036	"EXAMPLE:"; echo = 2;
4037	int plevel=printlevel;
4038	printlevel=-1;
4039	ring s=2,(x,y),lp;
4040	list HC=Adj_div(x3+y2+y);
4041	HC=NSplaces(1..2,HC);
4042	HC=extcurve(2,HC);
4043	def ER=HC[1][4];
4044	setring ER;
4045	intvec G=5; // the rational divisor G = 5*HC[3][1]
4046	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4047	// let us construct the corresponding residual AG code :
4048	matrix C=AGcode_Omega(G,D,HC);
4049	// here is a linear code of type [8,3,>=5] over F_4
4050	print(C);
4051	printlevel=plevel;
4052	}
4053	///////////////////////////////////////////////////////////////////////////////
4054
4055	// ============================================================================
4056	// ***** auxiliary procedure to define AG codes over extensions ******
4057	// ============================================================================
4058
4059	proc extcurve (int d,list CURVE)
4060	"USAGE: extcurve( d, CURVE ); d an integer, CURVE a list
4061	RETURN: list L which is the update of the list CURVE with additional entries
4062	@format
4063	L[1][3]: ring (p,a),(x,y),lp (affine),
4064	L[1][4]: ring (p,a),(x,y,z),lp (projective),
4065	L[1][5]: ring (p,a),(x,y,t),ls (local),
4066	L[2][3]: int (the number of rational places),
4067	@end format
4068	the rings being defined over a field extension of degree d. @*
4069	If d<2 then @code{extcurve(d,CURVE);} creates a list L which
4070	is the update of the list CURVE with additional entries
4071	@format
4072	L[1][5]: ring p,(x,y,t),ls,
4073	L[2][3]: int (the number of computed places over the base field).
4074	@end format
4075	In both cases, in the ring L[1][5] lists with the data for all the
4076	computed rational places (after a field extension of degree d) are
4077	created (see @ref{Adj_div}):
4078	@format
4079	lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS.
4080	@end format
4081	NOTE: The list CURVE should be the output of @code{NSplaces},
4082	and must contain (at least) one place of degree d. @*
4083	You actually need all the places with degree dividing d.
4084	Otherwise, not all the places are computed, but only part of them. @*
4085	This procedure must be executed before constructing AG codes,
4086	even if no extension is needed. The ring L[1][4] must be active
4087	when constructing codes over the field extension.@*
4088	SEE ALSO: closed_points, Adj_div, NSplaces, AGcode_L, AGcode_Omega
4089	EXAMPLE: example extcurve; shows an example
4090	"
4091	{
4092	// extends the underlying curve and all its associated objects to a larger
4093	// base field in order to evaluate points over such a extension
4094	// if d<2 then the only change is that a local ring "RatPl" (which is a
4095	// copy of "S(1)") is created in order to store the rational places
4096	// where we can do evaluations
4097	// otherwise, such a ring contains all places which are rational over the
4098	// extension
4099	// warning : list Places does not change so that only divisors which are
4100	// "rational over the prime field" are allowed; this probably will
4101	// change in the future
4102	// warning : the places in RatPl are ranged in increasing degree, respecting
4103	// the order from list Places and placing the conjugate branches all
4104	// together
4105	def BR=basering;
4106	list ext_CURVE=CURVE;
4107	if (d<2)
4108	{
4109	def SS=CURVE[5][1][1];
4110	ring RatPl=char(basering),(x,y,t),ls;
4111	list POINTS=imap(SS,POINTS);
4112	list LOC_EQS=imap(SS,LOC_EQS);
4113	list BRANCHES=imap(SS,BRANCHES);
4114	list PARAMETRIZATIONS=imap(SS,PARAMETRIZATIONS);
4115	export(POINTS);
4116	export(LOC_EQS);
4117	export(BRANCHES);
4118	export(PARAMETRIZATIONS);
4119	int NrRatPl=size(POINTS);
4120	ext_CURVE[2][3]=NrRatPl;
4121	setring BR;
4122	ext_CURVE[1][5]=RatPl;
4123	dbprint(printlevel+1,"");
4124	dbprint(printlevel+2,"Total number of rational places : "+string(NrRatPl));
4125	dbprint(printlevel+1,"");
4126	kill RatPl;
4127	return(ext_CURVE);
4128	}
4129	else
4130	{
4131	// exclude the case when no place of degree d was previously computed/found
4132	int dd=size(CURVE[5]);
4133	if (dd<d)
4134	{
4135	ERROR("you did not compute/find any place of degree "+string(d));
4136	}
4137	if (typeof(CURVE[5][d])<>"list")
4138	{
4139	ERROR("you did not compute/find any place of degree "+string(d));
4140	}
4141
4142	// Define the ring "RatPl" :
4143	def S(d)=CURVE[5][d][1];
4144	setring S(d);
4145	string smin=string(minpoly);
4146	setring BR;
4147	ring RatPl=(char(basering),a),(x,y,t),ls;
4148	execute("minpoly="+smin+";");
4149
4150	// import data from local ring S(d)
4151	list POINTS=imap(S(d),POINTS);
4152	list LOC_EQS=imap(S(d),LOC_EQS);
4153	list BRANCHES=imap(S(d),BRANCHES);
4154	list PARAMETRIZATIONS=imap(S(d),PARAMETRIZATIONS);
4155	kill S(d);
4156
4157	// conjugate data from S(d)
4158	int s=size(POINTS);
4159	int counter=0;
4160	int piv=0;
4161	int i,j,k;
4162	for (j=1;j<=s;j=j+1)
4163	{
4164	counter=counter+1;
4165	piv=counter;
4166	for (k=1;k<d;k=k+1)
4167	{
4168	POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4169	LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4170	BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4171	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4172	PARAMETRIZATIONS[piv],1),counter);
4173	counter=counter+1;
4174	piv=counter;
4175	}
4176	}
4177
4178	// now the same thing for smaller base fields
4179	string olda;
4180	poly paux;
4181	intvec iv,iw;
4182	int ii,jj,kk,m,n;
4183	for (i=d-1;i>=2;i=i-1)
4184	{
4185	if ( (d mod i)==0 )
4186	{
4187	if (typeof(CURVE[5][i])=="list")
4188	{
4189	def S(i)=CURVE[5][i][1];
4190	// embedd S(i) inside basering == RatPl " ==S(d) "
4191	olda=subfield(S(i));
4192	setring S(i);
4193	s=size(POINTS);
4194	setring RatPl;
4195	// import data from S(i)
4196	for (j=s;j>=1;j=j-1)
4197	{
4198	counter=0;
4199	POINTS=insert(POINTS,list(),0);
4200	POINTS[1][1]=number(importdatum(S(i),"POINTS["+string(j)
4201	+"][1]",olda));
4202	POINTS[1][2]=number(importdatum(S(i),"POINTS["+string(j)
4203	+"][2]",olda));
4204	POINTS[1][3]=number(importdatum(S(i),"POINTS["+string(j)
4205	+"][3]",olda));
4206	LOC_EQS=insert(LOC_EQS,importdatum(S(i),"LOC_EQS["+string(j)
4207	+"]",olda),0);
4208	BRANCHES=insert(BRANCHES,list(),0);
4209	setring S(i);
4210	m=nrows(BRANCHES[j][1]);
4211	n=ncols(BRANCHES[j][1]);
4212	iv=BRANCHES[j][2];
4213	kk=BRANCHES[j][3];
4214	poly par@1=subst(PARAMETRIZATIONS[j][1][1],t,x);
4215	poly par@2=subst(PARAMETRIZATIONS[j][1][2],t,x);
4216	export(par@1);
4217	export(par@2);
4218	iw=PARAMETRIZATIONS[j][2];
4219	setring RatPl;
4220	paux=importdatum(S(i),"BRANCHES["+string(j)+"][4]",olda);
4221	matrix Maux[m][n];
4222	for (ii=1;ii<=m;ii=ii+1)
4223	{
4224	for (jj=1;jj<=n;jj=jj+1)
4225	{
4226	Maux[ii,jj]=importdatum(S(i),"BRANCHES["+string(j)
4227	+"][1]["+string(ii)+","+string(jj)+"]",olda);
4228	}
4229	}
4230	BRANCHES[1][1]=Maux;
4231	BRANCHES[1][2]=iv;
4232	BRANCHES[1][3]=kk;
4233	BRANCHES[1][4]=paux;
4234	kill Maux;
4235	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,list(),0);
4236	PARAMETRIZATIONS[1][1]=ideal(0);
4237	PARAMETRIZATIONS[1][1][1]=importdatum(S(i),"par@1",olda);
4238	PARAMETRIZATIONS[1][1][2]=importdatum(S(i),"par@2",olda);
4239	PARAMETRIZATIONS[1][1][1]=subst(PARAMETRIZATIONS[1][1][1],x,t);
4240	PARAMETRIZATIONS[1][1][2]=subst(PARAMETRIZATIONS[1][1][2],x,t);
4241	PARAMETRIZATIONS[1][2]=iw;
4242	for (k=1;k<i;k=k+1)
4243	{
4244	counter=counter+1;
4245	piv=counter;
4246	POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4247	LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4248	BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4249	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4250	PARAMETRIZATIONS[piv],1),counter);
4251	}
4252	setring S(i);
4253	kill par@1,par@2;
4254	setring RatPl;
4255	}
4256	kill S(i);
4257	}
4258	}
4259	}
4260	kill iw;
4261	kill paux;
4262
4263	// finally add points from S(1)
4264	if (typeof(CURVE[5][1])=="list")
4265	{
4266	def S(1)=CURVE[5][1][1];
4267	POINTS=imap(S(1),POINTS)+POINTS;
4268	LOC_EQS=imap(S(1),LOC_EQS)+LOC_EQS;
4269	BRANCHES=imap(S(1),BRANCHES)+BRANCHES;
4270	PARAMETRIZATIONS=imap(S(1),PARAMETRIZATIONS)+PARAMETRIZATIONS;
4271	}
4272
4273	// prepare data for output
4274	export(POINTS);
4275	export(LOC_EQS);
4276	export(BRANCHES);
4277	export(PARAMETRIZATIONS);
4278	int NrRatPl=size(POINTS);
4279	ext_CURVE[2][3]=NrRatPl;
4280	setring BR;
4281	ext_CURVE[1][5]=RatPl;
4282	ring r(d)=(char(basering),a),(x,y),lp;
4283	execute("minpoly="+smin+";");
4284	setring BR;
4285	ext_CURVE[1][3]=r(d);
4286	ring R(d)=(char(basering),a),(x,y,z),lp;
4287	execute("minpoly="+smin+";");
4288	setring BR;
4289	ext_CURVE[1][4]=R(d);
4290	dbprint(printlevel+1,"");
4291	dbprint(printlevel+2,"Total number of rational places : NrRatPl = "
4292	+string(NrRatPl));
4293	dbprint(printlevel+1,"");
4294	kill S(1);
4295	kill R(d);
4296	kill RatPl;
4297	return(ext_CURVE);
4298	}
4299	}
4300	example
4301	{
4302	"EXAMPLE:"; echo = 2;
4303	int plevel=printlevel;
4304	printlevel=-1;
4305	ring s=2,(x,y),lp;
4306	list C=Adj_div(x5+y2+y);
4307	C=NSplaces(1..4,C);
4308	// since we have all points up to degree 4, we can extend the curve
4309	// to that extension, in order to get rational points over F_16;
4310	C=extcurve(4,C);
4311	// e.g., display the basepoint of place no. 32:
4312	def R=C[1][5];
4313	setring R;
4314	POINTS[32];
4315	printlevel=plevel;
4316	}
4317	///////////////////////////////////////////////////////////////////////////////
4318
4319	// specific procedures for linear/AG codes
4320
4321	static proc Hamming_wt (matrix A)
4322	"USAGE: Hamming_wt(A), where A is any matrix
4323	RETURN: the Hamming weight (number of non-zero entries) of the matrix A
4324	"
4325	{
4326	// computes the Hamming weight (number of non-zero entries) of any matrix
4327	// notice that "words" are represented by matrices of size 1xn
4328	// computing the Hamming distance between two matrices can be done by
4329	// Hamming_wt(A-B)
4330	int m=nrows(A);
4331	int n=ncols(A);
4332	int i,j;
4333	int w=0;
4334	for (i=1;i<=m;i=i+1)
4335	{
4336	for (j=1;j<=n;j=j+1)
4337	{
4338	if (A[i,j]<>0)
4339	{
4340	w=w+1;
4341	}
4342	}
4343	}
4344	return(w);
4345	}
4346	///////////////////////////////////////////////////////////////////////////////
4347
4348	// Basic Algorithm of Skorobogatov and Vladut for decoding AG codes
4349	// warning : the user must choose carefully the parameters for the code and
4350	// the decoding since they will never be checked by the procedures
4351
4352	proc prepSV (intvec G,intvec D,intvec F,list EC)
4353	"USAGE: prepSV( G, D, F, EC ); G,D,F intvecs and EC a list
4354	RETURN: list E of size n+3, where n=size(D). All its entries but E[n+3]
4355	are matrices:
4356	@format
4357	E[1]: parity check matrix for the current AG code
4358	E[2] ... E[n+2]: matrices used in the procedure decodeSV
4359	E[n+3]: intvec with
4360	E[n+3][1]: correction capacity @math{epsilon} of the algorithm
4361	E[n+3][2]: designed Goppa distance @math{delta} of the current AG code
4362	@end format
4363	NOTE: Computes the preprocessing for the basic (Skorobogatov-Vladut)
4364	decoding algorithm.@*
4365	The procedure must be called within the ring EC[1][4], where EC is
4366	the output of @code{extcurve(d)} (or in the ring EC[1][2] if d=1) @*
4367	The intvec G and F represent rational divisors (see
4368	@ref{BrillNoether} for more details).@*
4369	The intvec D refers to rational places (see @ref{AGcode_Omega}
4370	for more details.).
4371	The current AG code is @code{AGcode_Omega(G,D,EC)}.@*
4372	If you know the exact minimum distance d and you want to use it in
4373	@code{decodeSV} instead of @math{delta}, you can change the value
4374	of E[n+3][2] to d before applying decodeSV.
4375	If you have a systematic encoding for the current code and want to
4376	keep it during the decoding, you must previously permute D (using
4377	@code{permute_L(D,P);}), e.g., according to the permutation
4378	P=L[3], L being the output of @code{sys_code}.
4379	WARNINGS: F must be a divisor with support disjoint from the support of D and
4380	with degree @math{epsilon + genus}, where
4381	@math{epsilon:=[(deg(G)-3genus+1)/2]}.@
4382	G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
4383	not checked by the algorithm.
4384	G and D should also have disjoint supports (checked by the
4385	algorithm).
4386	KEYWORDS: SV-decoding algorithm, preprocessing
4387	SEE ALSO: extcurve, AGcode_Omega, decodeSV, sys_code, permute_L
4388	EXAMPLE: example prepSV; shows an example
4389	"
4390	{
4391	if (disj_divs(F,D,EC)==0)
4392	{
4393	dbprint(printlevel+3,"? the divisors do not have disjoint supports,
4394	empty list returned ?");
4395	return(list());
4396	}
4397	list E=list();
4398	list Places=EC[3];
4399	int m=deg_D(G,Places);
4400	int genusX=EC[2][2];
4401	int e=(m+1-3*genusX)/2;
4402	if (e<1)
4403	{
4404	dbprint(printlevel+3,"? the correction capacity of the basic algorithm
4405	is zero, empty list returned ?");
4406	return(list());
4407	}
4408	// deg(F)==e+genusX should be satisfied, and sup(D),sup(F) should be
4409	// disjoint !!!!
4410	int n=size(D);
4411	// 2*genusX-2<m<n should also be satisfied !!!!
4412	matrix EE=AGcode_L(G,D,EC);
4413	int l=nrows(EE);
4414	E[1]=EE;
4415	matrix GP=AGcode_L(F,D,EC);
4416	int r=nrows(GP);
4417	intvec H=G-F;
4418	matrix HP=AGcode_L(H,D,EC);
4419	int s=nrows(HP);
4420	int i,j,k;
4421	kill EE;
4422	for (k=1;k<=n;k=k+1)
4423	{
4424	E[k+1]=GP[1,k]*submat(HP,1..s,k..k);
4425	for (i=2;i<=r;i=i+1)
4426	{
4427	E[k+1]=concat(E[k+1],GP[i,k]*submat(HP,1..s,k..k));
4428	}
4429	}
4430	E[n+2]=GP;
4431	intvec IW=e,m+2-2*genusX;
4432	E[n+3]=IW;
4433	kill IW;
4434	return(E);
4435	}
4436	example
4437	{
4438	"EXAMPLE:"; echo = 2;
4439	int plevel=printlevel;
4440	printlevel=-1;
4441	ring s=2,(x,y),lp;
4442	list HC=Adj_div(x3+y2+y);
4443	HC=NSplaces(1..2,HC);
4444	HC=extcurve(2,HC);
4445	def ER=HC[1][4];
4446	setring ER;
4447	intvec G=5; // the rational divisor G = 5*HC[3][1]
4448	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4449	// construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4450	matrix C=AGcode_Omega(G,D,HC);
4451	// we can correct 1 error and the genus is 1, thus F must have degree 2
4452	// and support disjoint from that of D;
4453	intvec F=2;
4454	list SV=prepSV(G,D,F,HC);
4455	// now everything is prepared to decode with the basic algorithm;
4456	// for example, here is a parity check matrix to compute the syndrome :
4457	print(SV[1]);
4458	// and here you have the correction capacity of the algorithm :
4459	int epsilon=SV[size(D)+3][1];
4460	epsilon;
4461	printlevel=plevel;
4462	}
4463	///////////////////////////////////////////////////////////////////////////////
4464
4465	proc decodeSV (matrix y,list K)
4466	"USAGE: decodeSV( y, K ); y a row-matrix and K a list
4467	RETURN: a codeword (row-matrix) if possible, resp. the 0-matrix (of size
4468	1) if decoding is impossible.
4469	For decoding the basic (Skorobogatov-Vladut) decoding algorithm
4470	is applied.
4471	NOTE: The list_expression should be the output K of the procedure
4472	@code{prepSV}.@*
4473	The matrix_expression should be a (1 x n)-matrix, where
4474	n = ncols(K[1]).@*
4475	The decoding may fail if the number of errors is greater than
4476	the correction capacity of the algorithm.
4477	KEYWORDS: SV-decoding algorithm
4478	SEE ALSO: extcurve, AGcode_Omega, prepSV
4479	EXAMPLE: example decodeSV; shows an example
4480	"
4481	{
4482	// decodes y with the "basic decoding algorithm", if possible
4483	// requires the preprocessing done by the procedure "prepSV"
4484	// the procedure must be called from ring R or R(d)
4485	// returns either a codeword (matrix) of none (in case of too many errors)
4486	matrix syndr=K[1]*transpose(y);
4487	if (Hamming_wt(syndr)==0)
4488	{
4489	return(y);
4490	}
4491	matrix Ey=y[1,1]*K[2];
4492	int n=ncols(y);
4493	int i;
4494	for (i=2;i<=n;i=i+1)
4495	{
4496	Ey=Ey+y[1,i]*K[i+1];
4497	}
4498	matrix Ky=get_NZsol(Ey);
4499	if (Hamming_wt(Ky)==0)
4500	{
4501	dbprint(printlevel+3,"? no error-locator found ?");
4502	dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4503	matrix answer;
4504	return(answer);
4505	}
4506	int r=nrows(K[n+2]);
4507	matrix ErrLoc[1][n];
4508	list Z=list();
4509	list NZ=list();
4510	int j;
4511	for (j=1;j<=n;j=j+1)
4512	{
4513	for (i=1;i<=r;i=i+1)
4514	{
4515	ErrLoc[1,j]=ErrLoc[1,j]+K[n+2][i,j]*Ky[1,i];
4516	}
4517	if (ErrLoc[1,j]==0)
4518	{
4519	Z=insert(Z,j,size(Z));
4520	}
4521	else
4522	{
4523	NZ=insert(NZ,j,size(NZ));
4524	}
4525	}
4526	int k=size(NZ);
4527	int l=nrows(K[1]);
4528	int s=l+k;
4529	matrix A[s][n];
4530	matrix b[s][1];
4531	for (i=1;i<=l;i=i+1)
4532	{
4533	for (j=1;j<=n;j=j+1)
4534	{
4535	A[i,j]=K[1][i,j];
4536	}
4537	b[i,1]=syndr[i,1];
4538	}
4539	for (i=1;i<=k;i=i+1)
4540	{
4541	A[l+i,NZ[i]]=number(1);
4542	}
4543	intvec opgt=option(get);
4544	option(redSB);
4545	matrix L=transpose(syz(concat(A,-b)));
4546	if (nrows(L)==1)
4547	{
4548	if (L[1,n+1]<>0)
4549	{
4550	poly pivote=L[1,n+1];
4551	matrix sol=submat(L,1..1,1..n);
4552	if (pivote<>1)
4553	{
4554	sol=(number(1)/number(pivote))*sol;
4555	}
4556	// check at least that the number of comitted errors is less than half
4557	// the Goppa distance
4558	// imposing Hamming_wt(sol)<=K[n+3][1] would be more correct, but maybe
4559	// is too strong
4560	// on the other hand, if Hamming_wt(sol) is too large the decoding may
4561	// not be acceptable
4562	if ( Hamming_wt(sol) <= (K[n+3][2]-1)/2 )
4563	{
4564	option(set,opgt);
4565	return(y-sol);
4566	}
4567	else
4568	{
4569	dbprint(printlevel+3,"? non-acceptable decoding ?");
4570	}
4571	}
4572	else
4573	{
4574	dbprint(printlevel+3,"? no solution found ?");
4575	}
4576	}
4577	else
4578	{
4579	dbprint(printlevel+3,"? non-unique solution ?");
4580	}
4581	option(set,opgt);
4582	dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4583	matrix answer;
4584	return(answer);
4585	}
4586	example
4587	{
4588	"EXAMPLE:"; echo = 2;
4589	int plevel=printlevel;
4590	printlevel=-1;
4591	ring s=2,(x,y),lp;
4592	list HC=Adj_div(x3+y2+y);
4593	HC=NSplaces(1..2,HC);
4594	HC=extcurve(2,HC);
4595	def ER=HC[1][4];
4596	setring ER;
4597	intvec G=5; // the rational divisor G = 5*HC[3][1]
4598	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4599	// construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4600	matrix C=AGcode_Omega(G,D,HC);
4601	// we can correct 1 error and the genus is 1, thus F must have degree 2
4602	// and support disjoint from that of D
4603	intvec F=2;
4604	list SV=prepSV(G,D,F,HC);
4605	// now we produce 1 error on the zero-codeword :
4606	matrix y[1][8];
4607	y[1,3]=a;
4608	// and then we decode :
4609	print(decodeSV(y,SV));
4610	printlevel=plevel;
4611	}
4612	///////////////////////////////////////////////////////////////////////////////
4613
4614	proc sys_code (matrix C)
4615	"USAGE: sys_code(C); C is a matrix of constants
4616	RETURN: list L with:
4617	@format
4618	L[1] is the generator matrix in standard form of an equivalent code,
4619	L[2] is the parity check matrix in standard form of such code,
4620	L[3] is an intvec which represents the needed permutation.
4621	@end format
4622	NOTE: Computes a systematic code which is equivalent to the given one.@*
4623	The input should be a matrix of numbers.@*
4624	The output has to be interpreted as follows: if the input was
4625	the generator matrix of an AG code then one should apply the
4626	permutation L[3] to the divisor D of rational points by means
4627	of @code{permute_L(D,L[3]);} before continuing to work with the
4628	code (for instance, if you want to use the systematic encoding
4629	together with a decoding algorithm).
4630	KEYWORDS: linear code, systematic
4631	SEE ALSO: permute_L, AGcode_Omega, prepSV
4632	EXAMPLE: example sys_code; shows an example
4633	"
4634	{
4635	// computes a systematic code which is equivalent to that given by C
4636	int i,j,k,l,h,r;
4637	int m=nrows(C);
4638	int n=ncols(C);
4639	int mr=m;
4640	matrix A[m][n]=C;
4641	poly c,p;
4642	list corners=list();
4643	if(m>n)
4644	{
4645	mr=n;
4646	}
4647	// first the matrix A will be reduced with elementary operations by rows
4648	for(i=1;i<=mr;i=i+1)
4649	{
4650	if((i+l)>n)
4651	{
4652	// the process is finished
4653	break;
4654	}
4655	// look for a non-zero element
4656	if(A[i,i+l]==0)
4657	{
4658	h=i;
4659	p=0;
4660	for (j=i+1;j<=m;j=j+1)
4661	{
4662	c=A[j,i+l];
4663	if (c!=0)
4664	{
4665	p=c;
4666	h=j;
4667	break;
4668	}
4669	}
4670	if (h!=i)
4671	{
4672	// permutation of rows i and h
4673	for (j=1;j<=n;j=j+1)
4674	{
4675	c=A[i,j];
4676	A[i,j]=A[h,j];
4677	A[h,j]=c;
4678	}
4679	}
4680	if(p==0)
4681	{
4682	// non-zero element not found in the current column
4683	l=l+1;
4684	continue;
4685	}
4686	}
4687	// non-zero element was found in "strategic" position
4688	corners[i]=i+l;
4689	// make zeros below that position
4690	for(j=i+1;j<=m;j=j+1)
4691	{
4692	c=A[j,i+l]/A[i,i+l];
4693	for(k=i+l+1;k<=n;k=k+1)
4694	{
4695	A[j,k]=A[j,k]-A[i,k]*c;
4696	}
4697	A[j,i+l]=0;
4698	A[j,i]=0;
4699	}
4700	// the rank is at least r
4701	// when the process stops the last r is actually the true rank of A=a
4702	r=i;
4703	}
4704	if (r<m)
4705	{
4706	ERROR("the given matrix does not have maximum rank");
4707	}
4708	// set the corners to the beginning and construct the permutation
4709	intvec PCols=1..n;
4710	for (j=1;j<=m;j=j+1)
4711	{
4712	if (corners[j]>j)
4713	{
4714	// interchange columns j and corners[j]
4715	for (i=1;i<=m;i=i+1)
4716	{
4717	c=A[i,j];
4718	A[i,j]=A[i,corners[j]];
4719	A[i,corners[j]]=c;
4720	}
4721	k=PCols[j];
4722	PCols[j]=PCols[corners[j]];
4723	PCols[corners[j]]=k;
4724	}
4725	}
4726	// convert the diagonal into ones
4727	for (i=1;i<=m;i=i+1)
4728	{
4729	for (j=i+1;j<=n;j=j+1)
4730	{
4731	A[i,j]=A[i,j]/A[i,i];
4732	}
4733	A[i,i]=1;
4734	}
4735	// complete a block with the identity matrix
4736	for (k=1;k<m;k=k+1)
4737	{
4738	for (i=k+1;i<=m;i=i+1)
4739	{
4740	for (j=i+1;j<=n;j=j+1)
4741	{
4742	A[k,j]=A[k,j]-A[i,j]*A[k,i];
4743	}
4744	A[k,i]=0;
4745	}
4746	}
4747	// compute a parity-check matrix in standard form
4748	matrix B=concat(-transpose(submat(A,1..m,m+1..n)),diag(1,n-m));
4749	list L=list();
4750	L[1]=A;
4751	L[2]=B;
4752	L[3]=PCols;
4753	return(L);
4754	}
4755	example
4756	{
4757	"EXAMPLE:"; echo = 2;
4758	ring s=3,T,lp;
4759	matrix C[2][5]=0,1,0,1,1,0,1,0,0,1;
4760	print(C);
4761	list L=sys_code(C);
4762	L[3];
4763	// here is the generator matrix in standard form
4764	print(L[1]);
4765	// here is the control matrix in standard form
4766	print(L[2]);
4767	// we can check that both codes are dual to each other
4768	print(L[1]*transpose(L[2]));
4769	}
4770
4771	////////////// this procedure should be moved to hnoether.lib /////////
4772	static proc ratdevelop (poly chi)
4773	{
4774	int ring_is_changed;
4775	int k1=res_deg();
4776	list L=hnexpansion(chi,"ess");
4777	if (typeof(L[1])=="ring") {
4778	def altring=basering;
4779	def HNring = L[1]; setring HNring;
4780	def HND = hne;
4781	kill hne;
4782	ring_is_changed=1;
4783	}
4784	else {
4785	def HND=L[1];
4786	}
4787	int k2=res_deg();
4788	if (k1<k2)
4789	{
4790	int N=size(HND);
4791	list CB=list();
4792	int i,j,k,d;
4793	intvec dgs;
4794	list auxL=list();
4795	intvec auxiv;
4796	for (i=1;i<=N;i++)
4797	{
4798	CB[i]=conj_bs(HND[i],k1);
4799	dgs[i]=size(CB[i]);
4800	}
4801	for (i=1;i<N;i++)
4802	{
4803	for (j=i+1;j<=N;j++)
4804	{
4805	if (dgs[i]==dgs[j])
4806	{
4807	d=dgs[i];
4808	for (k=1;k<=d;k++)
4809	{
4810	if (intersection(CB[i][1],CB[j][k])==-1)
4811	{
4812	HND=delete(HND,j);
4813	CB=delete(CB,j);
4814	// dgs=delete(dgs,j);
4815	execute("auxL="+string(dgs)+";");
4816	auxL=delete(auxL,j);
4817	execute("dgs="+string(auxL)+";");
4818	N=N-1;
4819	j=j-1;
4820	break;
4821	}
4822	}
4823	}
4824	}
4825	}
4826	}
4827	if (ring_is_changed==0) { return(HND); }
4828	else { export HND; setring altring; return(HNring); }
4829	}
4830
4831
4832
4833	///////////////////////////////////////////////////////////////////////////////
4834
4835	/*
4836	// ============================================================================
4837	// ***** ADDITIONAL INFORMATION ABOUT THE LIBRARY ******
4838	// ============================================================================
4839
4840	A SINGULAR library for plane curves, Weierstrass semigroups and AG codes
4841	Also available via http://wmatem.eis.uva.es/~ignfar/singular/
4842
4843	PREVIOUS WARNINGS :
4844
4845	(1) The procedures will work only in positive characteristic
4846	(2) The base field must be prime (this may change in the future)
4847	This limitation is not too serious, since in coding theory
4848	the used curves are usually (if not always) defined over a
4849	prime field, and extensions are only considered for
4850	evaluating functions in a field with many points;
4851	by the same reason, auxiliary divisors are usually defined
4852	over the prime field,
4853	with the exception of that consisting of "rational points"
4854	(3) The curve must be absolutely irreducible (but it is not checked)
4855	(4) Only (algebraic projective) plane curves are considered
4856
4857	GENERAL CONCEPTS :
4858
4859	(1) An affine point P is represented by a std of a prime ideal,
4860	and an intvec containing the position of the places above P
4861	in the list of Places; if the point is at infinity, the ideal is
4862	changed by a homogeneous irreducible polynomial in two variables
4863	(2) A place is represented by :
4864	a base point (list of homogeneous coordinates),
4865	a local equation for the curve at the base point,
4866	a Hamburger-Noether expansion of the corresponding branch,
4867	and a local parametrization (in "t") of such branch; everything is
4868	stored in a local ring "_[5][d][1]", d being the degree of the place,
4869	by means of lists "POINTS,LOC_EQS,BRANCHES,PARAMETRIZATIONS", and
4870	the degrees of the base points corresponding to the places in the
4871	ring "_[5][d][1]" are stored in an intvec "_[5][d][2]"
4872	(3) A divisor is represented by an intvec, where the integer at the
4873	position i means the coefficient of the i-th place in the divisor
4874	(4) Rational functions are represented by numerator/denominator
4875	in form of ideals with two homogeneous generators
4876
4877	OUTLINE/EXAMPLE OF THE USE OF THE LIBRARY :
4878
4879	Plane curves :
4880
4881	(1.0) ring s=p,(x,y[,z]),lp;
4882
4883	Notice that if you use 3 variables, then the equation
4884	of the curve is assumed to be a homogeneous polynomial.
4885
4886	(1.1) list CURVE=Adj_div(equation);
4887
4888	In CURVE[3] are listed all the (singular closed) places
4889	with their corresponding degrees; thus, you can now decide
4890	how many other points you want to compute with NSplaces.
4891
4892	(1.2) CURVE=NSplaces(degrees,CURVE);
4893
4894	"degrees" is an intvec with the degrees we want to compute
4895
4896	(1.3) CURVE=extcurve(extension,CURVE);
4897
4898	The rational places over the extension are ranged in
4899	the ring CURVE[1][5] with the following rules:
4900
4901	(i) all the representatives of the same closed point
4902	are listed in consecutive positions;
4903	(ii) if deg(P)<deg(Q), then the representatives of P
4904	are listed before those of Q;
4905	(iii) if two closed points P,Q have the same degree,
4906	then the representatives of P are listed before
4907	if P appears before in the list CURVE[3].
4908
4909	Rational functions :
4910
4911	(2.0) def R=CURVE[1][2];
4912	setring R;
4913	(2.1) list LG=BrillNoether(intvec divisor,CURVE);
4914	(2.2) list WS=Weierstrass(int place,int bound,CURVE);
4915
4916	Algebraic Geometry codes :
4917
4918	(3.0) def ER=CURVE[1][4]; // if extension>1; else use R instead
4919	setring ER;
4920
4921	Now you have to decide the divisor G and the sequence of
4922	rational points D to use for constructing the codes;
4923	first notice that the syntax for G and D is different:
4924
4925	(a) G is a divisor supported on the closed places of
4926	CURVE[3], and you must say just the coefficient
4927	of each such a point; for example, G=2,0,-1 means
4928	2 times the place 1 minus 1 times the place 3.
4929
4930	(b) D is a sequence of rational points (all different
4931	and counted 1 time each one), whose data are read
4932	from the lists inside CURVE[1][5] and now you must
4933	say just the order how you range the chosen point;
4934	for example, D=2,4,1 means that you choose the
4935	rational places 1,2,4 and you range them as 2,4,1.
4936
4937	(3.1) matrix C=AGcode_L(divisor,places,CURVE);
4938
4939	(3.2) AGcode_Omega(divisor,places,CURVE);
4940
4941	In the same way as for defining the codes, the auxiliary
4942	divisor F must have disjoint support to that of D, and
4943	its degree has to be given by a formula (see help prepSV).
4944
4945	(3.3) list SV=prepSV(divisor,places,auxiliary_divisor,CURVE);
4946
4947	(3.4) decodeSV(word,SV);
4948
4949	Special Issues :
4950
4951	(I) AG codes with systematic encoding :
4952
4953	matrix C=AGcode_Omega(G,D,CURVE);
4954	list CODE=sys_code(G);
4955	C=CODE[1]; // generator matrix in standard form
4956	D=permute_L(D,CODE[3]); // suitable permutation of coordinates
4957	list SV=prepSV(G,D,F,CURVE);
4958	SV[1]=CODE[2]; // parity-check matrix in standard form
4959
4960	(II) Using the true minimum distance d for decoding :
4961
4962	matrix C=AGcode_Omega(G,D,CURVE);
4963	int n=size(D);
4964	list SV=prepSV(G,D,F,CURVE);
4965	SV[n+3][2]=d; // then use decodeSV(y,SV);
4966
4967
4968	// ============================================================================
4969	// * Some "macros" with typical examples of curves in Coding Theory **
4970	// ============================================================================
4971
4972
4973	proc Klein ()
4974	{
4975	list KLEIN=Adj_div(x3y+y3+x);
4976	KLEIN=NSplaces(1..3,KLEIN);
4977	KLEIN=extcurve(3,KLEIN);
4978	dbprint(printlevel+1,"Klein quartic over F_8 successfully constructed");
4979	return(KLEIN);
4980	}
4981
4982	proc Hermite (int m)
4983	{
4984	int p=char(basering);
4985	int r=p^m;
4986	list HERMITE=Adj_div(y^r+y-x^(r+1));
4987	HERMITE=NSplaces(1..2*m,HERMITE);
4988	HERMITE=extcurve(2*m,HERMITE);
4989	dbprint(printlevel+1,"Hermitian curve over F_("+string(r)+"^2)
4990	successfully constructed");
4991	return(HERMITE);
4992	}
4993
4994	proc Suzuki ()
4995	{
4996	list SUZUKI=Adj_div(x10+x3+y8+y);
4997	SUZUKI=NSplaces(1..3,SUZUKI);
4998	SUZUKI=extcurve(3,SUZUKI);
4999	dbprint(printlevel+1,"Suzuki curve over F_8 successfully constructed");
5000	return(SUZUKI);
5001	}
5002
5003	// **** Other interesting examples :
5004
5005	// A hyperelliptic curve with 33 rational points over F_16
5006
5007	list CURVE=Adj_div(x5+y2+y);
5008	CURVE=NSplaces(1..4,CURVE);
5009	CURVE=extcurve(4,CURVE);
5010
5011	// A nice curve with 113 rational points over F_64
5012
5013	list CURVE=Adj_div(y9+y8+xy6+x2y3+y2+x3);
5014	intvec ww=1,2,3,6;
5015	CURVE=NSplaces(ww,CURVE);
5016	CURVE=extcurve(6,CURVE);
5017
5018	*/

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

Download in other formats:

Original Format