Context Navigation

brnoeth.lib @ 667ba1

fieker-DuValspielwiese

Last change on this file since 667ba1 was e820b2, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
fixes due to minpoly' killer solution CHG: since setting minpoly kills local objects - declare them again if needed... (psi in brnoeth.lib) FIX: minpoly may have killed the map already... (pre_HN in LIB/hnoether.lib)
Property mode set to `100644`
File size: 139.9 KB

Line
1	//
2	version="$Id$";
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	kill psi;
2878	number A=a;
2879	number B=imap(base_r,B);
2880	}
2881	else
2882	{
2883	if (deg(subst(P[2],y,1))==1)
2884	{
2885	// the point is non-rational but the needed minpoly is just P[1]
2886	// careful : the parameter will be called "a" anyway
2887	poly P1=P[1];
2888	poly P2=P[2];
2889	ring ES=(char(basering),a),(x,y,t),ls;
2890	map psi=base_r,0,a;
2891	minpoly=number(psi(P1));
2892	kill psi;
2893	poly aux1=imap(base_r,P2);
2894	poly aux2=subst(aux1,y,a);
2895	number A=-number(subst(aux2,x,0));
2896	number B=a;
2897	}
2898	else
2899	{
2900	// this is the most complicated case of non-rational point
2901	poly P1=P[1];
2902	poly P2=P[2];
2903	int p=char(basering);
2904	int Q=p^ext;
2905	ring aux_r=(Q,a),(x,y,t),ls;
2906	string minpoly_string=string(minpoly);
2907	ring ES=(char(basering),a),(x,y,t),ls;
2908	execute("minpoly="+minpoly_string+";");
2909	poly P_1=imap(base_r,P1);
2910	poly P_2=imap(base_r,P2);
2911	ideal factors1=factorize(P_1,1);
2912	number B=-number(subst(factors1[1],y,0));
2913	poly P_0=subst(P_2,y,B);
2914	ideal factors2=factorize(P_0,1);
2915	number A=-number(subst(factors2[1],x,0));
2916	kill aux_r;
2917	}
2918	}
2919	list POINT=list();
2920	POINT[1]=A;
2921	POINT[2]=B;
2922	export(POINT);
2923	map phi=base_r,x+A,y+B;
2924	poly LOC_EQ=phi(CHI);
2925	kill A,B,phi;
2926	list LLL=hnexpansion(LOC_EQ,"ess");
2927	if (typeof(LLL[1])=="ring") {
2928	def altring=basering;
2929	def HNEring = LLL[1];
2930	setring HNEring;
2931	def L@HNE = hne[1];
2932	kill hne;
2933	export(L@HNE);
2934	int m=nrows(L@HNE[1]);
2935	int n=ncols(L@HNE[1]);
2936	intvec Li2=L@HNE[2];
2937	int Li3=L@HNE[3];
2938	setring ES;
2939	string newa=subfield(HNEring);
2940	poly paux=importdatum(HNEring,"L@HNE[4]",newa);
2941	matrix Maux[m][n];
2942	int i,j;
2943	for (i=1;i<=m;i=i+1)
2944	{
2945	for (j=1;j<=n;j=j+1)
2946	{
2947	Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
2948	string(j)+"]",newa);
2949	}
2950	}
2951	}
2952	else {
2953	def L@HNE=LLL[1][1];
2954	matrix Maux=L@HNE[1];
2955	intvec Li2=L@HNE[2];
2956	int Li3=L@HNE[3];
2957	poly paux=L@HNE[4];
2958	def HNEring=basering;
2959	setring HNEring;
2960	}
2961	kill LLL;
2962
2963	list BRANCH=list();
2964	BRANCH[1]=Maux;
2965	BRANCH[2]=Li2;
2966	BRANCH[3]=Li3;
2967	BRANCH[4]=paux;
2968	export(BRANCH);
2969	list PARAMETRIZATION=param(BRANCH,0);
2970	export(PARAMETRIZATION);
2971	kill HNEring;
2972	setring base_r;
2973	list answer=list();
2974	answer[1]=ES;
2975	answer[2]=ext;
2976	kill ES;
2977	return(answer);
2978	}
2979	///////////////////////////////////////////////////////////////////////////////
2980	static proc intersection_div (poly H,list CURVE)
2981	"USAGE: intersection_div(H,CURVE), where H is a homogeneous polynomial
2982	in ring Proj_R=p,(x,y,z),lp and CURVE is the list of data for
2983	the given curve
2984	CREATE: new places which had not been computed in advance if necessary
2985	those places are stored each one in a local ring where you find
2986	lists POINT,BRANCH,PARAMETRIZATION for the place living in that
2987	ring; the degree of the point/place in such a ring is stored in
2988	an intvec, and the base points in the remaining list
2989	Everything is exported in a list @EXTRA@ inside the ring
2990	aff_r=CURVE[1][1] and returned with the updated CURVE
2991	RETURN: list with the intersection divisor (intvec) between the underlying
2992	curve and H=0, and the list CURVE updated
2993	SEE ALSO: Adj_div, NSplaces, closed_points
2994	NOTE: The procedure must be called from the ring Proj_R=CURVE[1][2]
2995	(projective).
2996	If @EXTRA@ already exists, the new places are added to the
2997	previous data.
2998	"
2999	{
3000	// computes the intersection divisor of H and the curve CHI
3001	// returns a list with (possibly) "extra places" and it must be called
3002	// inside Proj_R
3003	// in case of extra places, some local rings ES(1) ... ES(m) are created
3004	// together with an integer list "extra_dgs" containing the degrees of
3005	// those places
3006	intvec opgt=option(get);
3007	option(redSB);
3008	intvec interdiv;
3009	def BRing=basering;
3010	int Tz1=deg(H);
3011	list Places=CURVE[3];
3012	int N=size(Places);
3013	def aff_r=CURVE[1][1];
3014	setring aff_r;
3015	if (defined(@EXTRA@)==0)
3016	{
3017	list @EXTRA@=list();
3018	list EP=list();
3019	list ES=list();
3020	list extra_dgs=list();
3021	}
3022	else
3023	{
3024	list EP=@EXTRA@[1];
3025	list ES=@EXTRA@[2];
3026	list extra_dgs=@EXTRA@[3];
3027	}
3028	int NN=size(extra_dgs);
3029	int counterEPl=0;
3030	setring BRing;
3031	poly h=subst(H,z,1);
3032	int Tz2=deg(h);
3033	int Tz3=Tz1-Tz2;
3034	int i,j,k,l,m,n,s,np,NP,I_N;
3035	if (Tz3==0)
3036	{
3037	// if this still does not work -> try always with ALL points in
3038	// Inf_Points !!!!
3039	poly Hinf=subst(H,z,0);
3040	setring aff_r;
3041	// compute the points at infinity of H and see which of them are in
3042	// Inf_Points
3043	poly h=imap(BRing,h);
3044	poly hinf=imap(BRing,Hinf);
3045	ideal pinf=factorize(hinf,1);
3046	list TIP=Inf_Points[1]+Inf_Points[2];
3047	s=size(TIP);
3048	NP=size(pinf);
3049	for (i=1;i<=NP;i=i+1)
3050	{
3051	for (j=1;j<=s;j=j+1)
3052	{
3053	if (pinf[i]==TIP[j][1])
3054	{
3055	np=size(TIP[j][2]);
3056	for (k=1;k<=np;k=k+1)
3057	{
3058	n=TIP[j][2][k];
3059	l=Places[n][1];
3060	m=Places[n][2];
3061	def SS=CURVE[5][l][1];
3062	setring SS;
3063	// local equation h of H
3064	if (POINTS[m][2]==number(1))
3065	{
3066	number A=POINTS[m][1];
3067	map psi=BRing,x+A,1,y;
3068	kill A;
3069	}
3070	else
3071	{
3072	map psi=BRing,1,x,y;
3073	}
3074	poly h=psi(H);
3075	I_N=local_IN(h,m)[1];
3076	interdiv[n]=I_N;
3077	kill h,psi;
3078	setring aff_r;
3079	kill SS;
3080	}
3081	break;
3082	}
3083	}
3084	}
3085	kill hinf,pinf;
3086	}
3087	else
3088	{
3089	// H is a multiple of z and hence all the points in Inf_Points intersect
3090	// with H
3091	setring aff_r;
3092	poly h=imap(BRing,h);
3093	list TIP=Inf_Points[1]+Inf_Points[2];
3094	s=size(TIP);
3095	for (j=1;j<=s;j=j+1)
3096	{
3097	np=size(TIP[j][2]);
3098	for (k=1;k<=np;k=k+1)
3099	{
3100	n=TIP[j][2][k];
3101	l=Places[n][1];
3102	m=Places[n][2];
3103	def SS=CURVE[5][l][1];
3104	setring SS;
3105	// local equation h of H
3106	if (POINTS[m][2]==number(1))
3107	{
3108	number A=POINTS[m][1];
3109	map psi=BRing,x+A,1,y;
3110	kill A;
3111	}
3112	else
3113	{
3114	map psi=BRing,1,x,y;
3115	}
3116	poly h=psi(H);
3117	I_N=local_IN(h,m)[1];
3118	interdiv[n]=I_N;
3119	kill h,psi;
3120	setring aff_r;
3121	kill SS;
3122	}
3123	}
3124	}
3125	// compute common affine points of H and CHI
3126	ideal CAL=h,CHI;
3127	CAL=std(CAL);
3128	if (CAL<>1)
3129	{
3130	list TAP=list();
3131	TAP=closed_points(CAL);
3132	NP=size(TAP);
3133	list auxP=list();
3134	int dP;
3135	for (i=1;i<=NP;i=i+1)
3136	{
3137	if (belongs(TAP[i],Aff_SLocus)==1)
3138	{
3139	// search the point in the list Aff_SPoints
3140	j=isInLP(TAP[i],Aff_SPoints);
3141	np=size(Aff_SPoints[j][2]);
3142	for (k=1;k<=np;k=k+1)
3143	{
3144	n=Aff_SPoints[j][2][k];
3145	l=Places[n][1];
3146	m=Places[n][2];
3147	def SS=CURVE[5][l][1];
3148	setring SS;
3149	// local equation h of H
3150	number A=POINTS[m][1];
3151	number B=POINTS[m][2];
3152	map psi=BRing,x+A,y+B,1;
3153	poly h=psi(H);
3154	I_N=local_IN(h,m)[1];
3155	interdiv[n]=I_N;
3156	kill A,B,h,psi;
3157	setring aff_r;
3158	kill SS;
3159	}
3160	}
3161	else
3162	{
3163	auxP=list();
3164	auxP[1]=TAP[i];
3165	dP=degree_P(auxP);
3166	if (defined(Aff_Points(dP))<>0)
3167	{
3168	// search the point in the list Aff_Points(dP)
3169	j=isInLP(TAP[i],Aff_Points(dP));
3170	n=Aff_Points(dP)[j][2][1];
3171	l=Places[n][1];
3172	m=Places[n][2];
3173	def SS=CURVE[5][l][1];
3174	setring SS;
3175	// local equation h of H
3176	number A=POINTS[m][1];
3177	number B=POINTS[m][2];
3178	map psi=BRing,x+A,y+B,1;
3179	poly h=psi(H);
3180	I_N=local_IN(h,m)[1];
3181	interdiv[n]=I_N;
3182	kill A,B,h,psi;
3183	setring aff_r;
3184	kill SS;
3185	}
3186	else
3187	{
3188	// previously check if it is an existing "extra place"
3189	j=isInLP(TAP[i],EP);
3190	if (j>0)
3191	{
3192	def SS=ES[j];
3193	setring SS;
3194	// local equation h of H
3195	number A=POINT[1];
3196	number B=POINT[2];
3197	map psi=BRing,x+A,y+B,1;
3198	poly h=psi(H);
3199	I_N=local_IN(h,0)[1];
3200	interdiv[N+j]=I_N;
3201	setring aff_r;
3202	kill SS;
3203	}
3204	else
3205	{
3206	// then we must create a new "extra place"
3207	counterEPl=counterEPl+1;
3208	list EXTRA_PLACE=extra_place(TAP[i]);
3209	def SS=EXTRA_PLACE[1];
3210	ES[NN+counterEPl]=SS;
3211	extra_dgs[NN+counterEPl]=EXTRA_PLACE[2];
3212	EP[NN+counterEPl]=list();
3213	EP[NN+counterEPl][1]=TAP[i];
3214	EP[NN+counterEPl][2]=0;
3215	setring SS;
3216	// local equation h of H
3217	number A=POINT[1];
3218	number B=POINT[2];
3219	map psi=BRing,x+A,y+B,1;
3220	poly h=psi(H);
3221	I_N=local_IN(h,0)[1];
3222	kill A,B,h,psi;
3223	interdiv[N+NN+counterEPl]=I_N;
3224	setring aff_r;
3225	kill SS;
3226	}
3227	}
3228	}
3229	}
3230	kill TAP,auxP;
3231	}
3232	kill h,CAL,TIP;
3233	@EXTRA@[1]=EP;
3234	@EXTRA@[2]=ES;
3235	@EXTRA@[3]=extra_dgs;
3236	kill EP;
3237	list update_CURVE=CURVE;
3238	if (size(extra_dgs)>0)
3239	{
3240	export(@EXTRA@);
3241	update_CURVE[1][1]=basering;
3242	}
3243	else
3244	{
3245	kill @EXTRA@;
3246	}
3247	setring BRing;
3248	kill h;
3249	kill aff_r;
3250	list answer=list();
3251	answer[1]=interdiv;
3252	answer[2]=update_CURVE;
3253	option(set,opgt);
3254	return(answer);
3255	}
3256	///////////////////////////////////////////////////////////////////////////////
3257	static proc local_eq (poly H,SS,int m)
3258	{
3259	// computes a local equation of polynomial H in the ring SS related to the place
3260	// "m"
3261	// list POINT/POINTS is searched depending on wether m=0 or m>0 respectively
3262	// warning : the procedure must be called from ring "Proj_R" and returns a
3263	// string
3264	def BRing=basering;
3265	setring SS;
3266	if (m>0)
3267	{
3268	if (POINTS[m][3]==number(1))
3269	{
3270	number A=POINTS[m][1];
3271	number B=POINTS[m][2];
3272	map psi=BRing,x+A,y+B,1;
3273	}
3274	else
3275	{
3276	if (POINTS[m][2]==number(1))
3277	{
3278	number A=POINTS[m][1];
3279	map psi=BRing,x+A,1,y;
3280	}
3281	else
3282	{
3283	map psi=BRing,1,x,y;
3284	}
3285	}
3286	}
3287	else
3288	{
3289	number A=POINT[1];
3290	number B=POINT[2];
3291	map psi=BRing,x+A,y+B,1;
3292	}
3293	poly h=psi(H);
3294	string str_h=string(h);
3295	setring BRing;
3296	return(str_h);
3297	}
3298	///////////////////////////////////////////////////////////////////////////////
3299	static proc min_wt_rmat (matrix M)
3300	{
3301	// finds the row of M with minimum non-zero entries, i.e. minimum
3302	// "Hamming-weight"
3303	int m=nrows(M);
3304	int n=ncols(M);
3305	int i,j;
3306	int Hwt=0;
3307	for (j=1;j<=n;j=j+1)
3308	{
3309	if (M[1,j]<>0)
3310	{
3311	Hwt=Hwt+1;
3312	}
3313	}
3314	int minHwt=Hwt;
3315	int k=1;
3316	for (i=2;i<=m;i=i+1)
3317	{
3318	Hwt=0;
3319	for (j=1;j<=n;j=j+1)
3320	{
3321	if (M[i,j]<>0)
3322	{
3323	Hwt=Hwt+1;
3324	}
3325	}
3326	if (Hwt<minHwt)
3327	{
3328	minHwt=Hwt;
3329	k=i;
3330	}
3331	}
3332	return(k);
3333	}
3334	///////////////////////////////////////////////////////////////////////////////
3335
3336	// ============================================================================
3337	// ***** MAIN PROCEDURE : the Brill-Noether algorithm ******
3338	// ============================================================================
3339
3340	proc BrillNoether (intvec G,list CURVE)
3341	"USAGE: BrillNoether(G,CURVE); G an intvec, CURVE a list
3342	RETURN: list of ideals (each of them with two homogeneous generators,
3343	which represent the numerator, resp. denominator, of a rational
3344	function).@*
3345	The corresponding rational functions form a vector space basis of the
3346	linear system L(G), G a rational divisor over a non-singular curve.
3347	NOTE: The procedure must be called from the ring CURVE[1][2], where
3348	CURVE is the output of the procedure @code{NSplaces}. @*
3349	The intvec G represents a rational divisor supported on the closed
3350	places of CURVE[3] (e.g. @code{G=2,0,-1;} means 2 times the closed
3351	place 1 minus 1 times the closed place 3).
3352	SEE ALSO: Adj_div, NSplaces, Weierstrass
3353	EXAMPLE: example BrillNoether; shows an example
3354	"
3355	{
3356	// computes a vector basis for the space L(G),
3357	// where G is a given rational divisor over the non-singular curve
3358	// returns : list of ideals in R each with 2 elements H,Ho such that
3359	// the set of functions {H/Ho} is the searched basis
3360	// warning : the conductor and sufficiently many points of the plane
3361	// curve should be computed in advance, in list CURVE
3362	// the algorithm of Brill-Noether is carried out in the procedure
3363	def BRing=basering;
3364	int degX=CURVE[2][1];
3365	list Places=CURVE[3];
3366	intvec Conductor=CURVE[4];
3367	if (deg_D(G,Places)<0)
3368	{
3369	return(list());
3370	}
3371	intvec nuldiv;
3372	if (G==nuldiv)
3373	{
3374	list quickL=list();
3375	ideal quickId;
3376	quickId[1]=1;
3377	quickId[2]=1;
3378	quickL[1]=quickId;
3379	return(quickL);
3380	}
3381	intvec J=max_D(G,nuldiv)+Conductor;
3382	int n=estim_n(J,degX,Places);
3383	dbprint(printlevel+1,"Forms of degree "+string(n)+" : ");
3384	matrix W=nmultiples(n,degX,CHI);
3385	kill nFORMS(n-degX);
3386	list update_CURVE=CURVE;
3387	matrix V=interpolating_forms(J,n,update_CURVE);
3388	matrix VmW=supplement(W,V);
3389	int k=min_wt_rmat(VmW);
3390	int N=size(nFORMS(n));
3391	matrix H0[1][N];
3392	int i,j;
3393	for (i=1;i<=N;i=i+1) { H0[1,i]=VmW[k,i]; }
3394	poly Ho;
3395	for (i=1;i<=N;i=i+1) { Ho=Ho+(H0[1,i]*nFORMS(n)[i]); }
3396	list INTERD=intersection_div(Ho,update_CURVE);
3397	intvec NHo=INTERD[1];
3398	update_CURVE=INTERD[2];
3399	intvec AR=NHo-G;
3400	matrix V2=interpolating_forms(AR,n,update_CURVE);
3401	def aux_RING=update_CURVE[1][1];
3402	setring aux_RING;
3403	if (defined(@EXTRA@)<>0)
3404	{
3405	kill @EXTRA@;
3406	}
3407	setring BRing;
3408	update_CURVE[1][1]=aux_RING;
3409	kill aux_RING;
3410	matrix B0=supplement(W,V2);
3411	if (Hamming_wt(B0)==0)
3412	{
3413	return(list());
3414	}
3415	int ld=nrows(B0);
3416	list Bres=list();
3417	ideal HH;
3418	poly H;
3419	for (j=1;j<=ld;j=j+1)
3420	{
3421	H=0;
3422	for (i=1;i<=N;i=i+1)
3423	{
3424	H=H+(B0[j,i]*nFORMS(n)[i]);
3425	}
3426	HH=H,Ho;
3427	Bres[j]=simplifyRF(HH);
3428	}
3429	kill nFORMS(n);
3430	dbprint(printlevel+1," ");
3431	dbprint(printlevel+2,"Vector basis successfully computed ");
3432	dbprint(printlevel+1," ");
3433	return(Bres);
3434	}
3435	example
3436	{
3437	"EXAMPLE:"; echo = 2;
3438	int plevel=printlevel;
3439	printlevel=-1;
3440	ring s=2,(x,y),lp;
3441	list C=Adj_div(x3y+y3+x);
3442	C=NSplaces(1..4,C);
3443	// the first 3 Places in C[3] are of degree 1.
3444	// we define the rational divisor G = 4C[3][1]+4C[3][3] (of degree 8):
3445	intvec G=4,0,4;
3446	def R=C[1][2];
3447	setring R;
3448	list LG=BrillNoether(G,C);
3449	// here is the vector basis of L(G):
3450	LG;
3451	printlevel=plevel;
3452	}
3453	///////////////////////////////////////////////////////////////////////////////
3454
3455	// * procedures for dealing with "RATIONAL FUNCTIONS" over a plane curve *
3456	// a rational function F may be given by (homogeneous) ideal or (affine) poly
3457	// (or number)
3458
3459	static proc polytoRF (F)
3460	{
3461	// converts a polynomial (or number) into a "rational function" of type "ideal"
3462	// warning : it must be called inside "R" and the polynomial is expected to be affine
3463	ideal RF;
3464	RF[1]=homog(F,z);
3465	RF[2]=z^(deg(F));
3466	return(RF);
3467	}
3468	///////////////////////////////////////////////////////////////////////////////
3469	static proc simplifyRF (ideal F)
3470	{
3471	// simplifies a rational function f/g extracting the gcd(f,g)
3472	// maybe add a "restriction" to the curve "CHI" but it is not easy to
3473	// programm
3474	poly auxp=gcd(F[1],F[2]);
3475	return(ideal(division(F,auxp)[1]));
3476	}
3477	///////////////////////////////////////////////////////////////////////////////
3478	static proc sumRF (F,G)
3479	{
3480	// sum of two "rational functions" F,G given by either a pair
3481	// numerator/denominator or a poly
3482	if ( typeof(F)=="ideal" && typeof(G)=="ideal" )
3483	{
3484	ideal FG;
3485	FG[1]=F[1]G[2]+F[2]G[1];
3486	FG[2]=F[2]*G[2];
3487	return(simplifyRF(FG));
3488	}
3489	else
3490	{
3491	if (typeof(F)=="ideal")
3492	{
3493	ideal GG=polytoRF(G);
3494	ideal FG;
3495	FG[1]=F[1]GG[2]+F[2]GG[1];
3496	FG[2]=F[2]*GG[2];
3497	return(simplifyRF(FG));
3498	}
3499	else
3500	{
3501	if (typeof(G)=="ideal")
3502	{
3503	ideal FF=polytoRF(F);
3504	ideal FG;
3505	FG[1]=FF[1]G[2]+FF[2]G[1];
3506	FG[2]=FF[2]*G[2];
3507	return(simplifyRF(FG));
3508	}
3509	else
3510	{
3511	return(F+G);
3512	}
3513	}
3514	}
3515	}
3516	///////////////////////////////////////////////////////////////////////////////
3517	static proc negRF (F)
3518	{
3519	// returns -F as rational function
3520	if (typeof(F)=="ideal")
3521	{
3522	ideal FF=F;
3523	FF[1]=-F[1];
3524	return(FF);
3525	}
3526	else
3527	{
3528	return(-F);
3529	}
3530	}
3531	///////////////////////////////////////////////////////////////////////////////
3532	static proc escprodRF (l,F)
3533	{
3534	// computes l*F as rational function
3535	// l should be either a number or a polynomial of degree zero
3536	if (typeof(F)=="ideal")
3537	{
3538	ideal lF=F;
3539	lF[1]=l*F[1];
3540	return(lF);
3541	}
3542	else
3543	{
3544	return(l*F);
3545	}
3546	}
3547	///////////////////////////////////////////////////////////////////////////////
3548
3549	// ****** procedures to compute Weierstrass semigroups ******
3550
3551	static proc orderRF (ideal F,SS,int m)
3552	"USAGE: orderRF(F,SS,m); F an ideal, SS a ring and m an integer
3553	RETURN: list with the order (int) and the leading coefficient (number)
3554	NOTE: F represents a rational function, thus the procedure must be
3555	called from global ring R or R(d).
3556	SS contains the name of a local ring where rational places are
3557	stored, and then we take that which is in position m in the
3558	corresponding lists of data.
3559	The order of F at the place given by SS,m is returned together
3560	with the coefficient of minimum degree in the corresponding power
3561	series.
3562	"
3563	{
3564	// computes the order of a rational function F at a RATIONAL place given by
3565	// a local ring SS and a position "m" inside SS
3566	// warning : the procedure must be called from global projective ring "R" or
3567	// "R(i)"
3568	// returns a list with the order (int) and the "leading coefficient" (number)
3569	def BR=basering;
3570	poly f=F[1];
3571	string sf=local_eq(f,SS,m);
3572	poly g=F[2];
3573	string sg=local_eq(g,SS,m);
3574	setring SS;
3575	execute("poly ff="+sf+";");
3576	execute("poly gg="+sg+";");
3577	list o1=local_IN(ff,m);
3578	list o2=local_IN(gg,m);
3579	int oo=o1[1]-o2[1];
3580	number lc=o1[2]/o2[2];
3581	setring BR;
3582	number LC=number(imap(SS,lc));
3583	return(list(oo,LC));
3584	}
3585	///////////////////////////////////////////////////////////////////////////////
3586
3587	proc Weierstrass (int P,int m,list CURVE)
3588	"USAGE: Weierstrass( i, m, CURVE ); i,m integers and CURVE a list
3589	RETURN: list WS of two lists:
3590	@format
3591	WS[1] list of integers (Weierstr. semigroup of the curve at place i up to m)
3592	WS[2] list of ideals (the associated rational functions)
3593	@end format
3594	NOTE: The procedure must be called from the ring CURVE[1][2],
3595	where CURVE is the output of the procedure @code{NSplaces}.
3596	@* i represents the place CURVE[3][i].
3597	@* Rational functions are represented by numerator/denominator
3598	in form of ideals with two homogeneous generators.
3599	WARNING: The place must be rational, i.e., necessarily CURVE[3][i][1]=1. @*
3600	SEE ALSO: Adj_div, NSplaces, BrillNoether
3601	EXAMPLE: example Weierstrass; shows an example
3602	"
3603	{
3604	// computes the Weierstrass semigroup at a RATIONAL place P up to a bound "m"
3605	// together with the functions achieving each value up to m, via
3606	// Brill-Noether
3607	// returns 2 lists : the first consists of all the poles up to m in
3608	// increasing order and the second consists of the corresponging rational
3609	// functions
3610	list Places=CURVE[3];
3611	intvec pl=Places[P];
3612	int dP=pl[1];
3613	int nP=pl[2];
3614	if (dP<>1)
3615	{
3616	ERROR("the given place is not defined over the prime field");
3617	}
3618	if (m<=0)
3619	{
3620	if (m==0)
3621	{
3622	list semig=list();
3623	int auxint=0;
3624	semig[1]=auxint;
3625	list funcs=list();
3626	ideal auxF;
3627	auxF[1]=1;
3628	auxF[2]=1;
3629	funcs[1]=auxF;
3630	return(list(semig,funcs));
3631	}
3632	else
3633	{
3634	ERROR("second argument must be non-negative");
3635	}
3636	}
3637	int auxint=0;
3638	ideal auxF;
3639	auxF[1]=1;
3640	auxF[2]=1;
3641	// Brill-Noether algorithm
3642	intvec mP;
3643	mP[P]=m;
3644	list LmP=BrillNoether(mP,CURVE);
3645	int lmP=size(LmP);
3646	if (lmP==1)
3647	{
3648	list semig=list();
3649	semig[1]=auxint;
3650	list funcs=list();
3651	funcs[1]=auxF;
3652	return(list(semig,funcs));
3653	}
3654	def SS=CURVE[5][dP][1];
3655	list ordLmP=list();
3656	list polLmP=list();
3657	list sortpol=list();
3658	int maxpol;
3659	int i,j,k;
3660	for (i=1;i<=lmP-1;i=i+1)
3661	{
3662	for (j=1;j<=lmP-i+1;j=j+1)
3663	{
3664	ordLmP[j]=orderRF(LmP[j],SS,nP);
3665	polLmP[j]=-ordLmP[j][1];
3666	}
3667	sortpol=sort(polLmP);
3668	polLmP=sortpol[1];
3669	maxpol=polLmP[lmP-i+1];
3670	LmP=permute_L(LmP,sortpol[2]);
3671	ordLmP=permute_L(ordLmP,sortpol[2]);
3672	// triangulate LmP
3673	for (k=1;k<=lmP-i;k=k+1)
3674	{
3675	if (polLmP[lmP-i+1-k]==maxpol)
3676	{
3677	LmP[lmP-i+1-k]=sumRF(LmP[lmP-i+1-k],negRF(escprodRF(
3678	ordLmP[lmP-i+1-k][2]/ordLmP[lmP-i+1][2],LmP[lmP-i+1])));
3679	}
3680	else
3681	{
3682	break;
3683	}
3684	}
3685	}
3686	polLmP[1]=auxint;
3687	LmP[1]=auxF;
3688	return(list(polLmP,LmP));
3689	}
3690	example
3691	{
3692	"EXAMPLE:"; echo = 2;
3693	int plevel=printlevel;
3694	printlevel=-1;
3695	ring s=2,(x,y),lp;
3696	list C=Adj_div(x3y+y3+x);
3697	C=NSplaces(1..4,C);
3698	def R=C[1][2];
3699	setring R;
3700	// Place C[3][1] has degree 1 (i.e it is rational);
3701	list WS=Weierstrass(1,7,C);
3702	// the first part of the list is the Weierstrass semigroup up to 7 :
3703	WS[1];
3704	// and the second part are the corresponding functions :
3705	WS[2];
3706	printlevel=plevel;
3707	}
3708	///////////////////////////////////////////////////////////////////////////////
3709
3710	// axiliary procedure for permuting a list or intvec
3711
3712	proc permute_L (L,P)
3713	"USAGE: permute_L( L, P ); L,P either intvecs or lists
3714	RETURN: list obtained from L by applying the permutation given by P.
3715	NOTE: If P is a list, all entries must be integers.
3716	SEE ALSO: sys_code, AGcode_Omega, prepSV
3717	EXAMPLE: example permute_L; shows an example
3718	"
3719	{
3720	// permutes the list L according to the permutation P (both intvecs or
3721	// lists of integers)
3722	int s=size(L);
3723	int n=size(P);
3724	int i;
3725	if (s<n)
3726	{
3727	for (i=s+1;i<=n;i=i+1)
3728	{
3729	L[i]=0;
3730	}
3731	s=size(L);
3732	}
3733	list auxL=L;
3734	for (i=1;i<=n;i=i+1)
3735	{
3736	auxL[i]=L[P[i]];
3737	}
3738	return(auxL);
3739	}
3740	example
3741	{
3742	"EXAMPLE:"; echo = 2;
3743	list L=list();
3744	L[1]="a";
3745	L[2]="b";
3746	L[3]="c";
3747	L[4]="d";
3748	intvec P=1,3,4,2;
3749	// the list L is permuted according to P :
3750	permute_L(L,P);
3751	}
3752	///////////////////////////////////////////////////////////////////////////////
3753	static proc evalRF (ideal F,SS,int m)
3754	"USAGE: evalRF(F,SS,m), where F is an ideal, SS is a ring and m is an
3755	integer
3756	RETURN: the evaluation (number) of F at the place given by SS,m if it is
3757	well-defined
3758	NOTE: F represents a rational function, thus the procedure must be
3759	called from R or R(d).
3760	SS contains the name of a local ring where rational places are
3761	stored, and then we take that which is in position m in the
3762	corresponding lists of data.
3763	"
3764	{
3765	// evaluates a rational function F at a RATIONAL place given by
3766	// a local ring SS and a position "m" inside SS
3767	list olc=orderRF(F,SS,m);
3768	int oo=olc[1];
3769	if (oo==0)
3770	{
3771	return(olc[2]);
3772	}
3773	else
3774	{
3775	if (oo>0)
3776	{
3777	return(number(0));
3778	}
3779	else
3780	{
3781	ERROR("the function is not well-defined at the given place");
3782	}
3783	}
3784	}
3785	///////////////////////////////////////////////////////////////////////////////
3786	//
3787	// ****** procedures for constructing AG codes ******
3788	//
3789	///////////////////////////////////////////////////////////////////////////////
3790
3791	static proc gen_mat (list LF,intvec LP,RP)
3792	"USAGE: gen_mat(LF,LP,RP); LF list of rational functions,
3793	LP intvec of rational places and RP a local ring
3794	RETURN: a generator matrix of the evaluation code given by LF and LP
3795	SEE ALSO: extcurve
3796	KEYWORDS: evaluation codes
3797	NOTE: Rational places are searched in local ring RP
3798	The procedure must be called from R or R(d) fromlist CURVE
3799	after having executed extcurve(d,CURVE)
3800	"
3801	{
3802	// computes a generator matrix (with numbers) of the evaluation code given
3803	// by a list of rational functions LF and a list of RATIONAL places LP
3804	int m=size(LF);
3805	int n=size(LP);
3806	matrix GM[m][n];
3807	int i,j;
3808	for (i=1;i<=m;i=i+1)
3809	{
3810	for (j=1;j<=n;j=j+1)
3811	{
3812	GM[i,j]=evalRF(LF[i],RP,LP[j]);
3813	}
3814	}
3815	return(GM);
3816	}
3817	///////////////////////////////////////////////////////////////////////////////
3818
3819	proc dual_code (matrix G)
3820	"USAGE: dual_code(G); G a matrix of numbers
3821	RETURN: a generator matrix of the dual code generated by G
3822	NOTE: The input should be a matrix G of numbers. @*
3823	The output is also a parity check matrix for the code defined by G
3824	KEYWORDS: linear code, dual
3825	EXAMPLE: example dual_code; shows an example
3826	"
3827	{
3828	// computes the dual code of C given by a generator matrix G
3829	// i.e. computes a parity-check matrix H of C
3830	// conversely : computes also G if H is given
3831	return(Ker(G));
3832	}
3833	example
3834	{
3835	"EXAMPLE:"; echo = 2;
3836	ring s=2,T,lp;
3837	// here is the Hamming code of length 7 and dimension 3
3838	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;
3839	print(G);
3840	matrix H=dual_code(G);
3841	print(H);
3842	}
3843	///////////////////////////////////////////////////////////////////////////////
3844
3845	// ======================================================================
3846	// ********* initial test for disjointness *************
3847	// ======================================================================
3848
3849	static proc disj_divs (intvec H,intvec P,list EC)
3850	{
3851	int s1=size(H);
3852	int s2=size(P);
3853	list PLACES=EC[3];
3854	def BRing=basering;
3855	def auxR=EC[1][5];
3856	setring auxR;
3857	int s=res_deg();
3858	setring BRing;
3859	kill auxR;
3860	int i,j,k,d,l,N,M;
3861	intvec auxIV,iw;
3862	for (i=1;i<=s;i=i+1)
3863	{
3864	if ( (s mod i) == 0 )
3865	{
3866	if (typeof(EC[5][i])=="list")
3867	{
3868	def auxR=EC[5][i][1];
3869	setring auxR;
3870	auxIV[i]=size(POINTS);
3871	setring BRing;
3872	kill auxR;
3873	}
3874	else
3875	{
3876	auxIV[i]=0;
3877	}
3878	}
3879	else
3880	{
3881	auxIV[i]=0;
3882	}
3883	}
3884	for (i=1;i<=s1;i=i+1)
3885	{
3886	if (H[i]<>0)
3887	{
3888	iw=PLACES[i];
3889	d=iw[1];
3890	if ( (s mod d) == 0 )
3891	{
3892	l=iw[2];
3893	// check that this place(s) is/are not in sup(D)
3894	if (d==1)
3895	{
3896	for (j=1;j<=s2;j=j+1)
3897	{
3898	if (l==P[j])
3899	{
3900	return(0);
3901	}
3902	}
3903	}
3904	else
3905	{
3906	N=0;
3907	for (j=1;j<d;j=j+1)
3908	{
3909	N=N+j*auxIV[j];
3910	}
3911	N=N+d*(l-1);
3912	M=N+d;
3913	for (k=N+1;k<=M;k=k+1)
3914	{
3915	for (j=1;j<=s2;j=j+1)
3916	{
3917	if (k==P[j])
3918	{
3919	return(0);
3920	}
3921	}
3922	}
3923	}
3924	}
3925	}
3926	}
3927	kill auxIV,iw;
3928	return(1);
3929	}
3930	///////////////////////////////////////////////////////////////////////////////
3931
3932	proc AGcode_L (intvec G,intvec D,list EC)
3933	"USAGE: AGcode_L( G, D, EC ); G,D intvec, EC a list
3934	RETURN: a generator matrix for the evaluation AG code defined by the
3935	divisors G and D.
3936	NOTE: The procedure must be called within the ring EC[1][4],
3937	where EC is the output of @code{extcurve(d)} (or within
3938	the ring EC[1][2] if d=1). @*
3939	The entry i in the intvec D refers to the i-th rational
3940	place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
3941	The intvec G represents a rational divisor (see @ref{BrillNoether}
3942	for more details).@*
3943	The code evaluates the vector space basis of L(G) at the rational
3944	places given by D.
3945	WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
3946	not checked by the algorithm.
3947	G and D should have disjoint supports (checked by the algorithm).
3948	SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_Omega
3949	EXAMPLE: example AGcode_L; shows an example
3950	"
3951	{
3952	// returns a generator matrix for the evaluation AG code given by G and D
3953	// G must be a divisor defined over the prime field and D an intvec of
3954	// "rational places"
3955	// it must be called inside R or R(d) and requires previously "extcurve(d)"
3956	def BR=basering;
3957	if (disj_divs(G,D,EC)==0)
3958	{
3959	dbprint(printlevel+3,"? the divisors do not have disjoint supports,
3960	0-matrix returned ?");
3961	matrix answer;
3962	return(answer);
3963	}
3964	if (res_deg()>1)
3965	{
3966	def R=EC[1][2];
3967	setring R;
3968	list LG=BrillNoether(G,EC);
3969	setring BR;
3970	list LG=imap(R,LG);
3971	setring R;
3972	kill LG;
3973	setring BR;
3974	kill R;
3975	}
3976	else
3977	{
3978	list LG=BrillNoether(G,EC);
3979	}
3980	def RP=EC[1][5];
3981	matrix M=gen_mat(LG,D,RP);
3982	kill LG;
3983	return(M);
3984	}
3985	example
3986	{
3987	"EXAMPLE:"; echo = 2;
3988	int plevel=printlevel;
3989	printlevel=-1;
3990	ring s=2,(x,y),lp;
3991	list HC=Adj_div(x3+y2+y);
3992	HC=NSplaces(1..2,HC);
3993	HC=extcurve(2,HC);
3994	def ER=HC[1][4];
3995	setring ER;
3996	intvec G=5; // the rational divisor G = 5*HC[3][1]
3997	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
3998	// let us construct the corresponding evaluation AG code :
3999	matrix C=AGcode_L(G,D,HC);
4000	// here is a linear code of type [8,5,>=3] over F_4
4001	print(C);
4002	printlevel=plevel;
4003	}
4004	///////////////////////////////////////////////////////////////////////////////
4005
4006	proc AGcode_Omega (intvec G,intvec D,list EC)
4007	"USAGE: AGcode_Omega( G, D, EC ); G,D intvec, EC a list
4008	RETURN: a generator matrix for the residual AG code defined by the
4009	divisors G and D.
4010	NOTE: The procedure must be called within the ring EC[1][4],
4011	where EC is the output of @code{extcurve(d)} (or within
4012	the ring EC[1][2] if d=1). @*
4013	The entry i in the intvec D refers to the i-th rational
4014	place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
4015	The intvec G represents a rational divisor (see @ref{BrillNoether}
4016	for more details).@*
4017	The code computes the residues of a vector space basis of
4018	@math{\Omega(G-D)} at the rational places given by D.
4019	WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
4020	not checked by the algorithm.
4021	G and D should have disjoint supports (checked by the algorithm).
4022	SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_L
4023	EXAMPLE: example AGcode_Omega; shows an example
4024	"
4025	{
4026	// returns a generator matrix for the residual AG code given by G and D
4027	// G must be a divisor defined over the prime field and D an intvec or
4028	// "rational places"
4029	// it must be called inside R or R(d) and requires previously "extcurve(d)"
4030	return(dual_code(AGcode_L(G,D,EC)));
4031	}
4032	example
4033	{
4034	"EXAMPLE:"; echo = 2;
4035	int plevel=printlevel;
4036	printlevel=-1;
4037	ring s=2,(x,y),lp;
4038	list HC=Adj_div(x3+y2+y);
4039	HC=NSplaces(1..2,HC);
4040	HC=extcurve(2,HC);
4041	def ER=HC[1][4];
4042	setring ER;
4043	intvec G=5; // the rational divisor G = 5*HC[3][1]
4044	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4045	// let us construct the corresponding residual AG code :
4046	matrix C=AGcode_Omega(G,D,HC);
4047	// here is a linear code of type [8,3,>=5] over F_4
4048	print(C);
4049	printlevel=plevel;
4050	}
4051	///////////////////////////////////////////////////////////////////////////////
4052
4053	// ============================================================================
4054	// ***** auxiliary procedure to define AG codes over extensions ******
4055	// ============================================================================
4056
4057	proc extcurve (int d,list CURVE)
4058	"USAGE: extcurve( d, CURVE ); d an integer, CURVE a list
4059	RETURN: list L which is the update of the list CURVE with additional entries
4060	@format
4061	L[1][3]: ring (p,a),(x,y),lp (affine),
4062	L[1][4]: ring (p,a),(x,y,z),lp (projective),
4063	L[1][5]: ring (p,a),(x,y,t),ls (local),
4064	L[2][3]: int (the number of rational places),
4065	@end format
4066	the rings being defined over a field extension of degree d. @*
4067	If d<2 then @code{extcurve(d,CURVE);} creates a list L which
4068	is the update of the list CURVE with additional entries
4069	@format
4070	L[1][5]: ring p,(x,y,t),ls,
4071	L[2][3]: int (the number of computed places over the base field).
4072	@end format
4073	In both cases, in the ring L[1][5] lists with the data for all the
4074	computed rational places (after a field extension of degree d) are
4075	created (see @ref{Adj_div}):
4076	@format
4077	lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS.
4078	@end format
4079	NOTE: The list CURVE should be the output of @code{NSplaces},
4080	and must contain (at least) one place of degree d. @*
4081	You actually need all the places with degree dividing d.
4082	Otherwise, not all the places are computed, but only part of them. @*
4083	This procedure must be executed before constructing AG codes,
4084	even if no extension is needed. The ring L[1][4] must be active
4085	when constructing codes over the field extension.@*
4086	SEE ALSO: closed_points, Adj_div, NSplaces, AGcode_L, AGcode_Omega
4087	EXAMPLE: example extcurve; shows an example
4088	"
4089	{
4090	// extends the underlying curve and all its associated objects to a larger
4091	// base field in order to evaluate points over such a extension
4092	// if d<2 then the only change is that a local ring "RatPl" (which is a
4093	// copy of "S(1)") is created in order to store the rational places
4094	// where we can do evaluations
4095	// otherwise, such a ring contains all places which are rational over the
4096	// extension
4097	// warning : list Places does not change so that only divisors which are
4098	// "rational over the prime field" are allowed; this probably will
4099	// change in the future
4100	// warning : the places in RatPl are ranged in increasing degree, respecting
4101	// the order from list Places and placing the conjugate branches all
4102	// together
4103	def BR=basering;
4104	list ext_CURVE=CURVE;
4105	if (d<2)
4106	{
4107	def SS=CURVE[5][1][1];
4108	ring RatPl=char(basering),(x,y,t),ls;
4109	list POINTS=imap(SS,POINTS);
4110	list LOC_EQS=imap(SS,LOC_EQS);
4111	list BRANCHES=imap(SS,BRANCHES);
4112	list PARAMETRIZATIONS=imap(SS,PARAMETRIZATIONS);
4113	export(POINTS);
4114	export(LOC_EQS);
4115	export(BRANCHES);
4116	export(PARAMETRIZATIONS);
4117	int NrRatPl=size(POINTS);
4118	ext_CURVE[2][3]=NrRatPl;
4119	setring BR;
4120	ext_CURVE[1][5]=RatPl;
4121	dbprint(printlevel+1,"");
4122	dbprint(printlevel+2,"Total number of rational places : "+string(NrRatPl));
4123	dbprint(printlevel+1,"");
4124	kill RatPl;
4125	return(ext_CURVE);
4126	}
4127	else
4128	{
4129	// exclude the case when no place of degree d was previously computed/found
4130	int dd=size(CURVE[5]);
4131	if (dd<d)
4132	{
4133	ERROR("you did not compute/find any place of degree "+string(d));
4134	}
4135	if (typeof(CURVE[5][d])<>"list")
4136	{
4137	ERROR("you did not compute/find any place of degree "+string(d));
4138	}
4139
4140	// Define the ring "RatPl" :
4141	def S(d)=CURVE[5][d][1];
4142	setring S(d);
4143	string smin=string(minpoly);
4144	setring BR;
4145	ring RatPl=(char(basering),a),(x,y,t),ls;
4146	execute("minpoly="+smin+";");
4147
4148	// import data from local ring S(d)
4149	list POINTS=imap(S(d),POINTS);
4150	list LOC_EQS=imap(S(d),LOC_EQS);
4151	list BRANCHES=imap(S(d),BRANCHES);
4152	list PARAMETRIZATIONS=imap(S(d),PARAMETRIZATIONS);
4153	kill S(d);
4154
4155	// conjugate data from S(d)
4156	int s=size(POINTS);
4157	int counter=0;
4158	int piv=0;
4159	int i,j,k;
4160	for (j=1;j<=s;j=j+1)
4161	{
4162	counter=counter+1;
4163	piv=counter;
4164	for (k=1;k<d;k=k+1)
4165	{
4166	POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4167	LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4168	BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4169	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4170	PARAMETRIZATIONS[piv],1),counter);
4171	counter=counter+1;
4172	piv=counter;
4173	}
4174	}
4175
4176	// now the same thing for smaller base fields
4177	string olda;
4178	poly paux;
4179	intvec iv,iw;
4180	int ii,jj,kk,m,n;
4181	for (i=d-1;i>=2;i=i-1)
4182	{
4183	if ( (d mod i)==0 )
4184	{
4185	if (typeof(CURVE[5][i])=="list")
4186	{
4187	def S(i)=CURVE[5][i][1];
4188	// embedd S(i) inside basering == RatPl " ==S(d) "
4189	olda=subfield(S(i));
4190	setring S(i);
4191	s=size(POINTS);
4192	setring RatPl;
4193	// import data from S(i)
4194	for (j=s;j>=1;j=j-1)
4195	{
4196	counter=0;
4197	POINTS=insert(POINTS,list(),0);
4198	POINTS[1][1]=number(importdatum(S(i),"POINTS["+string(j)
4199	+"][1]",olda));
4200	POINTS[1][2]=number(importdatum(S(i),"POINTS["+string(j)
4201	+"][2]",olda));
4202	POINTS[1][3]=number(importdatum(S(i),"POINTS["+string(j)
4203	+"][3]",olda));
4204	LOC_EQS=insert(LOC_EQS,importdatum(S(i),"LOC_EQS["+string(j)
4205	+"]",olda),0);
4206	BRANCHES=insert(BRANCHES,list(),0);
4207	setring S(i);
4208	m=nrows(BRANCHES[j][1]);
4209	n=ncols(BRANCHES[j][1]);
4210	iv=BRANCHES[j][2];
4211	kk=BRANCHES[j][3];
4212	poly par@1=subst(PARAMETRIZATIONS[j][1][1],t,x);
4213	poly par@2=subst(PARAMETRIZATIONS[j][1][2],t,x);
4214	export(par@1);
4215	export(par@2);
4216	iw=PARAMETRIZATIONS[j][2];
4217	setring RatPl;
4218	paux=importdatum(S(i),"BRANCHES["+string(j)+"][4]",olda);
4219	matrix Maux[m][n];
4220	for (ii=1;ii<=m;ii=ii+1)
4221	{
4222	for (jj=1;jj<=n;jj=jj+1)
4223	{
4224	Maux[ii,jj]=importdatum(S(i),"BRANCHES["+string(j)
4225	+"][1]["+string(ii)+","+string(jj)+"]",olda);
4226	}
4227	}
4228	BRANCHES[1][1]=Maux;
4229	BRANCHES[1][2]=iv;
4230	BRANCHES[1][3]=kk;
4231	BRANCHES[1][4]=paux;
4232	kill Maux;
4233	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,list(),0);
4234	PARAMETRIZATIONS[1][1]=ideal(0);
4235	PARAMETRIZATIONS[1][1][1]=importdatum(S(i),"par@1",olda);
4236	PARAMETRIZATIONS[1][1][2]=importdatum(S(i),"par@2",olda);
4237	PARAMETRIZATIONS[1][1][1]=subst(PARAMETRIZATIONS[1][1][1],x,t);
4238	PARAMETRIZATIONS[1][1][2]=subst(PARAMETRIZATIONS[1][1][2],x,t);
4239	PARAMETRIZATIONS[1][2]=iw;
4240	for (k=1;k<i;k=k+1)
4241	{
4242	counter=counter+1;
4243	piv=counter;
4244	POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4245	LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4246	BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4247	PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4248	PARAMETRIZATIONS[piv],1),counter);
4249	}
4250	setring S(i);
4251	kill par@1,par@2;
4252	setring RatPl;
4253	}
4254	kill S(i);
4255	}
4256	}
4257	}
4258	kill iw;
4259	kill paux;
4260
4261	// finally add points from S(1)
4262	if (typeof(CURVE[5][1])=="list")
4263	{
4264	def S(1)=CURVE[5][1][1];
4265	POINTS=imap(S(1),POINTS)+POINTS;
4266	LOC_EQS=imap(S(1),LOC_EQS)+LOC_EQS;
4267	BRANCHES=imap(S(1),BRANCHES)+BRANCHES;
4268	PARAMETRIZATIONS=imap(S(1),PARAMETRIZATIONS)+PARAMETRIZATIONS;
4269	}
4270
4271	// prepare data for output
4272	export(POINTS);
4273	export(LOC_EQS);
4274	export(BRANCHES);
4275	export(PARAMETRIZATIONS);
4276	int NrRatPl=size(POINTS);
4277	ext_CURVE[2][3]=NrRatPl;
4278	setring BR;
4279	ext_CURVE[1][5]=RatPl;
4280	ring r(d)=(char(basering),a),(x,y),lp;
4281	execute("minpoly="+smin+";");
4282	setring BR;
4283	ext_CURVE[1][3]=r(d);
4284	ring R(d)=(char(basering),a),(x,y,z),lp;
4285	execute("minpoly="+smin+";");
4286	setring BR;
4287	ext_CURVE[1][4]=R(d);
4288	dbprint(printlevel+1,"");
4289	dbprint(printlevel+2,"Total number of rational places : NrRatPl = "
4290	+string(NrRatPl));
4291	dbprint(printlevel+1,"");
4292	kill S(1);
4293	kill R(d);
4294	kill RatPl;
4295	return(ext_CURVE);
4296	}
4297	}
4298	example
4299	{
4300	"EXAMPLE:"; echo = 2;
4301	int plevel=printlevel;
4302	printlevel=-1;
4303	ring s=2,(x,y),lp;
4304	list C=Adj_div(x5+y2+y);
4305	C=NSplaces(1..4,C);
4306	// since we have all points up to degree 4, we can extend the curve
4307	// to that extension, in order to get rational points over F_16;
4308	C=extcurve(4,C);
4309	// e.g., display the basepoint of place no. 32:
4310	def R=C[1][5];
4311	setring R;
4312	POINTS[32];
4313	printlevel=plevel;
4314	}
4315	///////////////////////////////////////////////////////////////////////////////
4316
4317	// specific procedures for linear/AG codes
4318
4319	static proc Hamming_wt (matrix A)
4320	"USAGE: Hamming_wt(A), where A is any matrix
4321	RETURN: the Hamming weight (number of non-zero entries) of the matrix A
4322	"
4323	{
4324	// computes the Hamming weight (number of non-zero entries) of any matrix
4325	// notice that "words" are represented by matrices of size 1xn
4326	// computing the Hamming distance between two matrices can be done by
4327	// Hamming_wt(A-B)
4328	int m=nrows(A);
4329	int n=ncols(A);
4330	int i,j;
4331	int w=0;
4332	for (i=1;i<=m;i=i+1)
4333	{
4334	for (j=1;j<=n;j=j+1)
4335	{
4336	if (A[i,j]<>0)
4337	{
4338	w=w+1;
4339	}
4340	}
4341	}
4342	return(w);
4343	}
4344	///////////////////////////////////////////////////////////////////////////////
4345
4346	// Basic Algorithm of Skorobogatov and Vladut for decoding AG codes
4347	// warning : the user must choose carefully the parameters for the code and
4348	// the decoding since they will never be checked by the procedures
4349
4350	proc prepSV (intvec G,intvec D,intvec F,list EC)
4351	"USAGE: prepSV( G, D, F, EC ); G,D,F intvecs and EC a list
4352	RETURN: list E of size n+3, where n=size(D). All its entries but E[n+3]
4353	are matrices:
4354	@format
4355	E[1]: parity check matrix for the current AG code
4356	E[2] ... E[n+2]: matrices used in the procedure decodeSV
4357	E[n+3]: intvec with
4358	E[n+3][1]: correction capacity @math{epsilon} of the algorithm
4359	E[n+3][2]: designed Goppa distance @math{delta} of the current AG code
4360	@end format
4361	NOTE: Computes the preprocessing for the basic (Skorobogatov-Vladut)
4362	decoding algorithm.@*
4363	The procedure must be called within the ring EC[1][4], where EC is
4364	the output of @code{extcurve(d)} (or in the ring EC[1][2] if d=1) @*
4365	The intvec G and F represent rational divisors (see
4366	@ref{BrillNoether} for more details).@*
4367	The intvec D refers to rational places (see @ref{AGcode_Omega}
4368	for more details.).
4369	The current AG code is @code{AGcode_Omega(G,D,EC)}.@*
4370	If you know the exact minimum distance d and you want to use it in
4371	@code{decodeSV} instead of @math{delta}, you can change the value
4372	of E[n+3][2] to d before applying decodeSV.
4373	If you have a systematic encoding for the current code and want to
4374	keep it during the decoding, you must previously permute D (using
4375	@code{permute_L(D,P);}), e.g., according to the permutation
4376	P=L[3], L being the output of @code{sys_code}.
4377	WARNINGS: F must be a divisor with support disjoint from the support of D and
4378	with degree @math{epsilon + genus}, where
4379	@math{epsilon:=[(deg(G)-3genus+1)/2]}.@
4380	G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
4381	not checked by the algorithm.
4382	G and D should also have disjoint supports (checked by the
4383	algorithm).
4384	KEYWORDS: SV-decoding algorithm, preprocessing
4385	SEE ALSO: extcurve, AGcode_Omega, decodeSV, sys_code, permute_L
4386	EXAMPLE: example prepSV; shows an example
4387	"
4388	{
4389	if (disj_divs(F,D,EC)==0)
4390	{
4391	dbprint(printlevel+3,"? the divisors do not have disjoint supports,
4392	empty list returned ?");
4393	return(list());
4394	}
4395	list E=list();
4396	list Places=EC[3];
4397	int m=deg_D(G,Places);
4398	int genusX=EC[2][2];
4399	int e=(m+1-3*genusX) div 2;
4400	if (e<1)
4401	{
4402	dbprint(printlevel+3,"? the correction capacity of the basic algorithm
4403	is zero, empty list returned ?");
4404	return(list());
4405	}
4406	// deg(F)==e+genusX should be satisfied, and sup(D),sup(F) should be
4407	// disjoint !!!!
4408	int n=size(D);
4409	// 2*genusX-2<m<n should also be satisfied !!!!
4410	matrix EE=AGcode_L(G,D,EC);
4411	int l=nrows(EE);
4412	E[1]=EE;
4413	matrix GP=AGcode_L(F,D,EC);
4414	int r=nrows(GP);
4415	intvec H=G-F;
4416	matrix HP=AGcode_L(H,D,EC);
4417	int s=nrows(HP);
4418	int i,j,k;
4419	kill EE;
4420	for (k=1;k<=n;k=k+1)
4421	{
4422	E[k+1]=GP[1,k]*submat(HP,1..s,k..k);
4423	for (i=2;i<=r;i=i+1)
4424	{
4425	E[k+1]=concat(E[k+1],GP[i,k]*submat(HP,1..s,k..k));
4426	}
4427	}
4428	E[n+2]=GP;
4429	intvec IW=e,m+2-2*genusX;
4430	E[n+3]=IW;
4431	kill IW;
4432	return(E);
4433	}
4434	example
4435	{
4436	"EXAMPLE:"; echo = 2;
4437	int plevel=printlevel;
4438	printlevel=-1;
4439	ring s=2,(x,y),lp;
4440	list HC=Adj_div(x3+y2+y);
4441	HC=NSplaces(1..2,HC);
4442	HC=extcurve(2,HC);
4443	def ER=HC[1][4];
4444	setring ER;
4445	intvec G=5; // the rational divisor G = 5*HC[3][1]
4446	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4447	// construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4448	matrix C=AGcode_Omega(G,D,HC);
4449	// we can correct 1 error and the genus is 1, thus F must have degree 2
4450	// and support disjoint from that of D;
4451	intvec F=2;
4452	list SV=prepSV(G,D,F,HC);
4453	// now everything is prepared to decode with the basic algorithm;
4454	// for example, here is a parity check matrix to compute the syndrome :
4455	print(SV[1]);
4456	// and here you have the correction capacity of the algorithm :
4457	int epsilon=SV[size(D)+3][1];
4458	epsilon;
4459	printlevel=plevel;
4460	}
4461	///////////////////////////////////////////////////////////////////////////////
4462
4463	proc decodeSV (matrix y,list K)
4464	"USAGE: decodeSV( y, K ); y a row-matrix and K a list
4465	RETURN: a codeword (row-matrix) if possible, resp. the 0-matrix (of size
4466	1) if decoding is impossible.
4467	For decoding the basic (Skorobogatov-Vladut) decoding algorithm
4468	is applied.
4469	NOTE: The list_expression should be the output K of the procedure
4470	@code{prepSV}.@*
4471	The matrix_expression should be a (1 x n)-matrix, where
4472	n = ncols(K[1]).@*
4473	The decoding may fail if the number of errors is greater than
4474	the correction capacity of the algorithm.
4475	KEYWORDS: SV-decoding algorithm
4476	SEE ALSO: extcurve, AGcode_Omega, prepSV
4477	EXAMPLE: example decodeSV; shows an example
4478	"
4479	{
4480	// decodes y with the "basic decoding algorithm", if possible
4481	// requires the preprocessing done by the procedure "prepSV"
4482	// the procedure must be called from ring R or R(d)
4483	// returns either a codeword (matrix) of none (in case of too many errors)
4484	matrix syndr=K[1]*transpose(y);
4485	if (Hamming_wt(syndr)==0)
4486	{
4487	return(y);
4488	}
4489	matrix Ey=y[1,1]*K[2];
4490	int n=ncols(y);
4491	int i;
4492	for (i=2;i<=n;i=i+1)
4493	{
4494	Ey=Ey+y[1,i]*K[i+1];
4495	}
4496	matrix Ky=get_NZsol(Ey);
4497	if (Hamming_wt(Ky)==0)
4498	{
4499	dbprint(printlevel+3,"? no error-locator found ?");
4500	dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4501	matrix answer;
4502	return(answer);
4503	}
4504	int r=nrows(K[n+2]);
4505	matrix ErrLoc[1][n];
4506	list Z=list();
4507	list NZ=list();
4508	int j;
4509	for (j=1;j<=n;j=j+1)
4510	{
4511	for (i=1;i<=r;i=i+1)
4512	{
4513	ErrLoc[1,j]=ErrLoc[1,j]+K[n+2][i,j]*Ky[1,i];
4514	}
4515	if (ErrLoc[1,j]==0)
4516	{
4517	Z=insert(Z,j,size(Z));
4518	}
4519	else
4520	{
4521	NZ=insert(NZ,j,size(NZ));
4522	}
4523	}
4524	int k=size(NZ);
4525	int l=nrows(K[1]);
4526	int s=l+k;
4527	matrix A[s][n];
4528	matrix b[s][1];
4529	for (i=1;i<=l;i=i+1)
4530	{
4531	for (j=1;j<=n;j=j+1)
4532	{
4533	A[i,j]=K[1][i,j];
4534	}
4535	b[i,1]=syndr[i,1];
4536	}
4537	for (i=1;i<=k;i=i+1)
4538	{
4539	A[l+i,NZ[i]]=number(1);
4540	}
4541	intvec opgt=option(get);
4542	option(redSB);
4543	matrix L=transpose(syz(concat(A,-b)));
4544	if (nrows(L)==1)
4545	{
4546	if (L[1,n+1]<>0)
4547	{
4548	poly pivote=L[1,n+1];
4549	matrix sol=submat(L,1..1,1..n);
4550	if (pivote<>1)
4551	{
4552	sol=(number(1)/number(pivote))*sol;
4553	}
4554	// check at least that the number of comitted errors is less than half
4555	// the Goppa distance
4556	// imposing Hamming_wt(sol)<=K[n+3][1] would be more correct, but maybe
4557	// is too strong
4558	// on the other hand, if Hamming_wt(sol) is too large the decoding may
4559	// not be acceptable
4560	if ( Hamming_wt(sol) <= (K[n+3][2]-1) div 2 )
4561	{
4562	option(set,opgt);
4563	return(y-sol);
4564	}
4565	else
4566	{
4567	dbprint(printlevel+3,"? non-acceptable decoding ?");
4568	}
4569	}
4570	else
4571	{
4572	dbprint(printlevel+3,"? no solution found ?");
4573	}
4574	}
4575	else
4576	{
4577	dbprint(printlevel+3,"? non-unique solution ?");
4578	}
4579	option(set,opgt);
4580	dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4581	matrix answer;
4582	return(answer);
4583	}
4584	example
4585	{
4586	"EXAMPLE:"; echo = 2;
4587	int plevel=printlevel;
4588	printlevel=-1;
4589	ring s=2,(x,y),lp;
4590	list HC=Adj_div(x3+y2+y);
4591	HC=NSplaces(1..2,HC);
4592	HC=extcurve(2,HC);
4593	def ER=HC[1][4];
4594	setring ER;
4595	intvec G=5; // the rational divisor G = 5*HC[3][1]
4596	intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4
4597	// construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4598	matrix C=AGcode_Omega(G,D,HC);
4599	// we can correct 1 error and the genus is 1, thus F must have degree 2
4600	// and support disjoint from that of D
4601	intvec F=2;
4602	list SV=prepSV(G,D,F,HC);
4603	// now we produce 1 error on the zero-codeword :
4604	matrix y[1][8];
4605	y[1,3]=a;
4606	// and then we decode :
4607	print(decodeSV(y,SV));
4608	printlevel=plevel;
4609	}
4610	///////////////////////////////////////////////////////////////////////////////
4611
4612	proc sys_code (matrix C)
4613	"USAGE: sys_code(C); C is a matrix of constants
4614	RETURN: list L with:
4615	@format
4616	L[1] is the generator matrix in standard form of an equivalent code,
4617	L[2] is the parity check matrix in standard form of such code,
4618	L[3] is an intvec which represents the needed permutation.
4619	@end format
4620	NOTE: Computes a systematic code which is equivalent to the given one.@*
4621	The input should be a matrix of numbers.@*
4622	The output has to be interpreted as follows: if the input was
4623	the generator matrix of an AG code then one should apply the
4624	permutation L[3] to the divisor D of rational points by means
4625	of @code{permute_L(D,L[3]);} before continuing to work with the
4626	code (for instance, if you want to use the systematic encoding
4627	together with a decoding algorithm).
4628	KEYWORDS: linear code, systematic
4629	SEE ALSO: permute_L, AGcode_Omega, prepSV
4630	EXAMPLE: example sys_code; shows an example
4631	"
4632	{
4633	// computes a systematic code which is equivalent to that given by C
4634	int i,j,k,l,h,r;
4635	int m=nrows(C);
4636	int n=ncols(C);
4637	int mr=m;
4638	matrix A[m][n]=C;
4639	poly c,p;
4640	list corners=list();
4641	if(m>n)
4642	{
4643	mr=n;
4644	}
4645	// first the matrix A will be reduced with elementary operations by rows
4646	for(i=1;i<=mr;i=i+1)
4647	{
4648	if((i+l)>n)
4649	{
4650	// the process is finished
4651	break;
4652	}
4653	// look for a non-zero element
4654	if(A[i,i+l]==0)
4655	{
4656	h=i;
4657	p=0;
4658	for (j=i+1;j<=m;j=j+1)
4659	{
4660	c=A[j,i+l];
4661	if (c!=0)
4662	{
4663	p=c;
4664	h=j;
4665	break;
4666	}
4667	}
4668	if (h!=i)
4669	{
4670	// permutation of rows i and h
4671	for (j=1;j<=n;j=j+1)
4672	{
4673	c=A[i,j];
4674	A[i,j]=A[h,j];
4675	A[h,j]=c;
4676	}
4677	}
4678	if(p==0)
4679	{
4680	// non-zero element not found in the current column
4681	l=l+1;
4682	continue;
4683	}
4684	}
4685	// non-zero element was found in "strategic" position
4686	corners[i]=i+l;
4687	// make zeros below that position
4688	for(j=i+1;j<=m;j=j+1)
4689	{
4690	c=A[j,i+l]/A[i,i+l];
4691	for(k=i+l+1;k<=n;k=k+1)
4692	{
4693	A[j,k]=A[j,k]-A[i,k]*c;
4694	}
4695	A[j,i+l]=0;
4696	A[j,i]=0;
4697	}
4698	// the rank is at least r
4699	// when the process stops the last r is actually the true rank of A=a
4700	r=i;
4701	}
4702	if (r<m)
4703	{
4704	ERROR("the given matrix does not have maximum rank");
4705	}
4706	// set the corners to the beginning and construct the permutation
4707	intvec PCols=1..n;
4708	for (j=1;j<=m;j=j+1)
4709	{
4710	if (corners[j]>j)
4711	{
4712	// interchange columns j and corners[j]
4713	for (i=1;i<=m;i=i+1)
4714	{
4715	c=A[i,j];
4716	A[i,j]=A[i,corners[j]];
4717	A[i,corners[j]]=c;
4718	}
4719	k=PCols[j];
4720	PCols[j]=PCols[corners[j]];
4721	PCols[corners[j]]=k;
4722	}
4723	}
4724	// convert the diagonal into ones
4725	for (i=1;i<=m;i=i+1)
4726	{
4727	for (j=i+1;j<=n;j=j+1)
4728	{
4729	A[i,j]=A[i,j]/A[i,i];
4730	}
4731	A[i,i]=1;
4732	}
4733	// complete a block with the identity matrix
4734	for (k=1;k<m;k=k+1)
4735	{
4736	for (i=k+1;i<=m;i=i+1)
4737	{
4738	for (j=i+1;j<=n;j=j+1)
4739	{
4740	A[k,j]=A[k,j]-A[i,j]*A[k,i];
4741	}
4742	A[k,i]=0;
4743	}
4744	}
4745	// compute a parity-check matrix in standard form
4746	matrix B=concat(-transpose(submat(A,1..m,m+1..n)),diag(1,n-m));
4747	list L=list();
4748	L[1]=A;
4749	L[2]=B;
4750	L[3]=PCols;
4751	return(L);
4752	}
4753	example
4754	{
4755	"EXAMPLE:"; echo = 2;
4756	ring s=3,T,lp;
4757	matrix C[2][5]=0,1,0,1,1,0,1,0,0,1;
4758	print(C);
4759	list L=sys_code(C);
4760	L[3];
4761	// here is the generator matrix in standard form
4762	print(L[1]);
4763	// here is the control matrix in standard form
4764	print(L[2]);
4765	// we can check that both codes are dual to each other
4766	print(L[1]*transpose(L[2]));
4767	}
4768
4769	////////////// this procedure should be moved to hnoether.lib /////////
4770	static proc ratdevelop (poly chi)
4771	{
4772	int ring_is_changed;
4773	int k1=res_deg();
4774	list L=hnexpansion(chi,"ess");
4775	if (typeof(L[1])=="ring") {
4776	def altring=basering;
4777	def HNring = L[1]; setring HNring;
4778	def HND = hne;
4779	kill hne;
4780	ring_is_changed=1;
4781	}
4782	else {
4783	def HND=L[1];
4784	}
4785	int k2=res_deg();
4786	if (k1<k2)
4787	{
4788	int N=size(HND);
4789	list CB=list();
4790	int i,j,k,d;
4791	intvec dgs;
4792	list auxL=list();
4793	intvec auxiv;
4794	for (i=1;i<=N;i++)
4795	{
4796	CB[i]=conj_bs(HND[i],k1);
4797	dgs[i]=size(CB[i]);
4798	}
4799	for (i=1;i<N;i++)
4800	{
4801	for (j=i+1;j<=N;j++)
4802	{
4803	if (dgs[i]==dgs[j])
4804	{
4805	d=dgs[i];
4806	for (k=1;k<=d;k++)
4807	{
4808	if (intersection(CB[i][1],CB[j][k])==-1)
4809	{
4810	HND=delete(HND,j);
4811	CB=delete(CB,j);
4812	// dgs=delete(dgs,j);
4813	execute("auxL="+string(dgs)+";");
4814	auxL=delete(auxL,j);
4815	execute("dgs="+string(auxL)+";");
4816	N=N-1;
4817	j=j-1;
4818	break;
4819	}
4820	}
4821	}
4822	}
4823	}
4824	}
4825	if (ring_is_changed==0) { return(HND); }
4826	else { export HND; setring altring; return(HNring); }
4827	}
4828
4829
4830
4831	///////////////////////////////////////////////////////////////////////////////
4832
4833	/*
4834	// ============================================================================
4835	// ***** ADDITIONAL INFORMATION ABOUT THE LIBRARY ******
4836	// ============================================================================
4837
4838	A SINGULAR library for plane curves, Weierstrass semigroups and AG codes
4839	Also available via http://wmatem.eis.uva.es/~ignfar/singular/
4840
4841	PREVIOUS WARNINGS :
4842
4843	(1) The procedures will work only in positive characteristic
4844	(2) The base field must be prime (this may change in the future)
4845	This limitation is not too serious, since in coding theory
4846	the used curves are usually (if not always) defined over a
4847	prime field, and extensions are only considered for
4848	evaluating functions in a field with many points;
4849	by the same reason, auxiliary divisors are usually defined
4850	over the prime field,
4851	with the exception of that consisting of "rational points"
4852	(3) The curve must be absolutely irreducible (but it is not checked)
4853	(4) Only (algebraic projective) plane curves are considered
4854
4855	GENERAL CONCEPTS :
4856
4857	(1) An affine point P is represented by a std of a prime ideal,
4858	and an intvec containing the position of the places above P
4859	in the list of Places; if the point is at infinity, the ideal is
4860	changed by a homogeneous irreducible polynomial in two variables
4861	(2) A place is represented by :
4862	a base point (list of homogeneous coordinates),
4863	a local equation for the curve at the base point,
4864	a Hamburger-Noether expansion of the corresponding branch,
4865	and a local parametrization (in "t") of such branch; everything is
4866	stored in a local ring "_[5][d][1]", d being the degree of the place,
4867	by means of lists "POINTS,LOC_EQS,BRANCHES,PARAMETRIZATIONS", and
4868	the degrees of the base points corresponding to the places in the
4869	ring "_[5][d][1]" are stored in an intvec "_[5][d][2]"
4870	(3) A divisor is represented by an intvec, where the integer at the
4871	position i means the coefficient of the i-th place in the divisor
4872	(4) Rational functions are represented by numerator/denominator
4873	in form of ideals with two homogeneous generators
4874
4875	OUTLINE/EXAMPLE OF THE USE OF THE LIBRARY :
4876
4877	Plane curves :
4878
4879	(1.0) ring s=p,(x,y[,z]),lp;
4880
4881	Notice that if you use 3 variables, then the equation
4882	of the curve is assumed to be a homogeneous polynomial.
4883
4884	(1.1) list CURVE=Adj_div(equation);
4885
4886	In CURVE[3] are listed all the (singular closed) places
4887	with their corresponding degrees; thus, you can now decide
4888	how many other points you want to compute with NSplaces.
4889
4890	(1.2) CURVE=NSplaces(degrees,CURVE);
4891
4892	"degrees" is an intvec with the degrees we want to compute
4893
4894	(1.3) CURVE=extcurve(extension,CURVE);
4895
4896	The rational places over the extension are ranged in
4897	the ring CURVE[1][5] with the following rules:
4898
4899	(i) all the representatives of the same closed point
4900	are listed in consecutive positions;
4901	(ii) if deg(P)<deg(Q), then the representatives of P
4902	are listed before those of Q;
4903	(iii) if two closed points P,Q have the same degree,
4904	then the representatives of P are listed before
4905	if P appears before in the list CURVE[3].
4906
4907	Rational functions :
4908
4909	(2.0) def R=CURVE[1][2];
4910	setring R;
4911	(2.1) list LG=BrillNoether(intvec divisor,CURVE);
4912	(2.2) list WS=Weierstrass(int place,int bound,CURVE);
4913
4914	Algebraic Geometry codes :
4915
4916	(3.0) def ER=CURVE[1][4]; // if extension>1; else use R instead
4917	setring ER;
4918
4919	Now you have to decide the divisor G and the sequence of
4920	rational points D to use for constructing the codes;
4921	first notice that the syntax for G and D is different:
4922
4923	(a) G is a divisor supported on the closed places of
4924	CURVE[3], and you must say just the coefficient
4925	of each such a point; for example, G=2,0,-1 means
4926	2 times the place 1 minus 1 times the place 3.
4927
4928	(b) D is a sequence of rational points (all different
4929	and counted 1 time each one), whose data are read
4930	from the lists inside CURVE[1][5] and now you must
4931	say just the order how you range the chosen point;
4932	for example, D=2,4,1 means that you choose the
4933	rational places 1,2,4 and you range them as 2,4,1.
4934
4935	(3.1) matrix C=AGcode_L(divisor,places,CURVE);
4936
4937	(3.2) AGcode_Omega(divisor,places,CURVE);
4938
4939	In the same way as for defining the codes, the auxiliary
4940	divisor F must have disjoint support to that of D, and
4941	its degree has to be given by a formula (see help prepSV).
4942
4943	(3.3) list SV=prepSV(divisor,places,auxiliary_divisor,CURVE);
4944
4945	(3.4) decodeSV(word,SV);
4946
4947	Special Issues :
4948
4949	(I) AG codes with systematic encoding :
4950
4951	matrix C=AGcode_Omega(G,D,CURVE);
4952	list CODE=sys_code(G);
4953	C=CODE[1]; // generator matrix in standard form
4954	D=permute_L(D,CODE[3]); // suitable permutation of coordinates
4955	list SV=prepSV(G,D,F,CURVE);
4956	SV[1]=CODE[2]; // parity-check matrix in standard form
4957
4958	(II) Using the true minimum distance d for decoding :
4959
4960	matrix C=AGcode_Omega(G,D,CURVE);
4961	int n=size(D);
4962	list SV=prepSV(G,D,F,CURVE);
4963	SV[n+3][2]=d; // then use decodeSV(y,SV);
4964
4965
4966	// ============================================================================
4967	// * Some "macros" with typical examples of curves in Coding Theory **
4968	// ============================================================================
4969
4970
4971	proc Klein ()
4972	{
4973	list KLEIN=Adj_div(x3y+y3+x);
4974	KLEIN=NSplaces(1..3,KLEIN);
4975	KLEIN=extcurve(3,KLEIN);
4976	dbprint(printlevel+1,"Klein quartic over F_8 successfully constructed");
4977	return(KLEIN);
4978	}
4979
4980	proc Hermite (int m)
4981	{
4982	int p=char(basering);
4983	int r=p^m;
4984	list HERMITE=Adj_div(y^r+y-x^(r+1));
4985	HERMITE=NSplaces(1..2*m,HERMITE);
4986	HERMITE=extcurve(2*m,HERMITE);
4987	dbprint(printlevel+1,"Hermitian curve over F_("+string(r)+"^2)
4988	successfully constructed");
4989	return(HERMITE);
4990	}
4991
4992	proc Suzuki ()
4993	{
4994	list SUZUKI=Adj_div(x10+x3+y8+y);
4995	SUZUKI=NSplaces(1..3,SUZUKI);
4996	SUZUKI=extcurve(3,SUZUKI);
4997	dbprint(printlevel+1,"Suzuki curve over F_8 successfully constructed");
4998	return(SUZUKI);
4999	}
5000
5001	// **** Other interesting examples :
5002
5003	// A hyperelliptic curve with 33 rational points over F_16
5004
5005	list CURVE=Adj_div(x5+y2+y);
5006	CURVE=NSplaces(1..4,CURVE);
5007	CURVE=extcurve(4,CURVE);
5008
5009	// A nice curve with 113 rational points over F_64
5010
5011	list CURVE=Adj_div(y9+y8+xy6+x2y3+y2+x3);
5012	intvec ww=1,2,3,6;
5013	CURVE=NSplaces(ww,CURVE);
5014	CURVE=extcurve(6,CURVE);
5015
5016	*/

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

Download in other formats:

Original Format