Context Navigation

brnoeth.lib @ 4713c6

fieker-DuValspielwiese

Last change on this file since 4713c6 was ae0203, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: from 2-0 git-svn-id: file:///usr/local/Singular/svn/trunk@7044 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 137.5 KB

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

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

Download in other formats:

Original Format