Context Navigation

brnoeth.lib @ 380a17b

spielwiese

Last change on this file since 380a17b was 380a17b, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: new version numbers for libs
Property mode set to `100644`
File size: 140.1 KB

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

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

Download in other formats:

Original Format