Context Navigation

brnoeth.lib @ 7db9488

spielwiese

Last change on this file since 7db9488 was 7db9488, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: headlines git-svn-id: file:///usr/local/Singular/svn/trunk@4914 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 130.4 KB

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

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

Download in other formats:

Original Format