1 | ///////////////////////////////////////////////////////////////////// |
---|
2 | version="on hess.lib 4.0.1.0 Oct_2014 $"; //$Id$ |
---|
3 | category="Hess"; |
---|
4 | info=" |
---|
5 | LIBRARY: hess.lib Riemann-Roch space of divisors |
---|
6 | on function fields and curves |
---|
7 | |
---|
8 | AUTHORS: I. Stenger: stenger@mathematik.uni-kl.de |
---|
9 | |
---|
10 | OVERVIEW: |
---|
11 | Let f be an absolutely irreducible polynomial in two variables x,y. |
---|
12 | Assume that f is monic as a polynomial in y. Let F = Quot(k[x,y]/f) |
---|
13 | be the function field defined by f. |
---|
14 | Define O_F = IntCl(k[x],F) and O_(F,inf) = IntCl(k[1/x],F). |
---|
15 | We represent a divisor D on F by two fractional ideals |
---|
16 | I and J of O_F and O_(F,inf), respectively. The Riemann-Roch |
---|
17 | space L(D) is then the intersection of I^(-1) and J^(-1). |
---|
18 | |
---|
19 | PROCEDURES: |
---|
20 | RiemannRochHess() Computes a vector space basis of the |
---|
21 | Riemann-Roch space of a divisor |
---|
22 | "; |
---|
23 | |
---|
24 | LIB "integralbasis.lib"; |
---|
25 | LIB "qhmoduli.lib"; |
---|
26 | LIB "dmodloc.lib"; |
---|
27 | LIB "modstd.lib"; |
---|
28 | LIB "matrix.lib"; |
---|
29 | LIB "absfact.lib"; |
---|
30 | |
---|
31 | |
---|
32 | ///////////////////////////////////////////////////////////////////// |
---|
33 | ///////////////////// Main Procedure //////////////////////////////// |
---|
34 | |
---|
35 | //__________________________________________________________________ |
---|
36 | proc RiemannRochHess(poly f,list divisorD,string s) |
---|
37 | "USAGE: RiemannRochHess(f,divisorD,s); f polynomial, divisorD list, |
---|
38 | s string |
---|
39 | NOTE: All fractional ideals must represented by a list of size two. |
---|
40 | The first element is an ideal of k[x,y] and the second element |
---|
41 | the common denominator, i.e, a polynomial of k[x]. |
---|
42 | ASSUME: The base ring R must be a ring in two variables, say x,y, |
---|
43 | or three variables, say x,y,z. |
---|
44 | If nvars(R) = 2: |
---|
45 | - f is an absolutely irreducible polynomial, monic as a |
---|
46 | polynomial in y. |
---|
47 | - List divisorD describes a divisor D of F = Quot(k[x,y]/f). |
---|
48 | If (s = "ideals") |
---|
49 | D is given in ideal representation, i.e., divisorD is a |
---|
50 | list of size 2. |
---|
51 | divisorD[1] is the finite ideal of D, i.e., the |
---|
52 | fractional ideal of D of IntCl(k[x],F). |
---|
53 | divisorD[2] is the infinite ideal of D, i.e, |
---|
54 | the fractional ideal of D of IntCl(k[1/x],F). |
---|
55 | If (s = "free") |
---|
56 | D is given in free representation, i.e., divisorD is a list |
---|
57 | of size 2, containing the finite and infinite places of D |
---|
58 | with exponents. |
---|
59 | divisorD[i], i = 1,2, is a list. Each element of the list |
---|
60 | is again a list. The first entry is a fractional ideal, |
---|
61 | and the second an integer, the exponent of the place. |
---|
62 | If nvars(R) = 3: |
---|
63 | - f is an absolutely irreducible homogeneous polynomial |
---|
64 | describing the projective plane curve corresponding to |
---|
65 | the function field F. We assume that the dehomogenization |
---|
66 | of f w.r.t. z is monic as a polynomial in y. |
---|
67 | List divisorD describes a divisor D of F. |
---|
68 | If (s = "ideals") |
---|
69 | D is given in ideal representation, i.e., divisorD is |
---|
70 | a list of size 2. divisorD[1] is an ideal of the base |
---|
71 | ring representing the positive divisor of D and |
---|
72 | divisorD[2] is an ideal of the base ring representing the |
---|
73 | negative divisor. (i.e. D = (I) -(J)). |
---|
74 | If (s = "free") |
---|
75 | D is given in free representation, i.e., divisorD is a |
---|
76 | list of places of D. D[i][1] is an prime ideal and D[i][2] |
---|
77 | an integer, the exponent of the place. |
---|
78 | RETURN: A vector space basis of the Riemann-Roch space of D, |
---|
79 | stored in a list RRBasis. The list RRBasis contains a |
---|
80 | list, say rbasis, and a polynomial, say den. The basis of |
---|
81 | L(D) consists of all rational functions g/den, where g is |
---|
82 | an element of rbasis. |
---|
83 | EXAMPLE: exampleRiemannRochHess; shows an example |
---|
84 | " |
---|
85 | { |
---|
86 | int i,j; |
---|
87 | |
---|
88 | if ( (nvars(basering) != 2) && (nvars(basering) != 3) ) |
---|
89 | { |
---|
90 | ERROR("Basering must have two or three variables."); |
---|
91 | } |
---|
92 | int pr; |
---|
93 | |
---|
94 | // the affine case - no transformation needed |
---|
95 | if (nvars(basering) == 2) |
---|
96 | { |
---|
97 | def R_aff = basering; |
---|
98 | poly f_aff = f; |
---|
99 | } |
---|
100 | |
---|
101 | //the projective case, dehomogenize the input polynomial |
---|
102 | //and change to the affine ring |
---|
103 | if (nvars(basering) == 3) |
---|
104 | { |
---|
105 | if (homog(f)==0) |
---|
106 | { |
---|
107 | ERROR("Input polynomial must be homogeneous."); |
---|
108 | } |
---|
109 | pr = 1; |
---|
110 | def R = basering; |
---|
111 | list rlpr = ringlist(R); |
---|
112 | rlpr[2] = list(var(1),var(2)); |
---|
113 | rlpr[3] = list(list("dp",1:2),list("C",0)); |
---|
114 | def R_aff = ring(rlpr); // corresponding affine ring |
---|
115 | poly f_aff = subst(f,var(3),1); |
---|
116 | setring R_aff; |
---|
117 | poly f_aff = imap(R,f_aff); |
---|
118 | } |
---|
119 | |
---|
120 | //Checking assumptions |
---|
121 | ////////////////////////////////////////////////////////////// |
---|
122 | |
---|
123 | // No parameters or algebraic numbers are allowed. |
---|
124 | if(npars(R_aff) >0) |
---|
125 | { |
---|
126 | ERROR("No parameters or algebraic extensions are |
---|
127 | allowed in the base ring."); |
---|
128 | } |
---|
129 | |
---|
130 | // The polynomial f must be monic in the 2-th variable. |
---|
131 | matrix cs = coeffs(f_aff, var(2)); |
---|
132 | if(cs[size(cs),1] != 1) |
---|
133 | { |
---|
134 | ERROR("The input polynomial must be monic as a polynomial |
---|
135 | in the second variable."); |
---|
136 | } |
---|
137 | |
---|
138 | // The polynomial f must be irreducible. |
---|
139 | if (char(R_aff) == 0) |
---|
140 | { |
---|
141 | if(isIrreducible(f_aff) == 0) |
---|
142 | { |
---|
143 | ERROR("The input polynomial must be irreducible."); |
---|
144 | } |
---|
145 | } |
---|
146 | |
---|
147 | //Defining ring for the reduction |
---|
148 | //////////////////////////////////////////////////////////// |
---|
149 | |
---|
150 | list rl=ringlist(R_aff); |
---|
151 | rl[2] = list(var(1)); |
---|
152 | rl[3] = list(list("lp",1),list("C",0)); |
---|
153 | def Rone = ring(rl); |
---|
154 | |
---|
155 | //Compute the integral bases |
---|
156 | //////////////////////////////////////////////////////////// |
---|
157 | |
---|
158 | //integral basis of IntCL(k[x],F); |
---|
159 | list FinBasis = integralBasis(f_aff,2); |
---|
160 | |
---|
161 | // integral basis of IntCL(k[1/x],F) |
---|
162 | list Infin; |
---|
163 | list InfBasis; |
---|
164 | if(pr) |
---|
165 | // computes the geometric places of infinity as well |
---|
166 | { |
---|
167 | Infin = maxorderInfinite(f_aff,1); |
---|
168 | InfBasis = Infin[1]; |
---|
169 | } |
---|
170 | else |
---|
171 | { |
---|
172 | list Infin = maxorderInfinite(f_aff); |
---|
173 | InfBasis = Infin[1]; |
---|
174 | } |
---|
175 | |
---|
176 | if (printlevel >= 3) |
---|
177 | { |
---|
178 | "Integral basis of the maximal finite order:","";FinBasis;pause();""; |
---|
179 | "Integral basis of the maximal infinite order:";"";InfBasis;pause();""; |
---|
180 | } |
---|
181 | |
---|
182 | //Compute the Ideal representation |
---|
183 | //////////////////////////////////////////////////////////// |
---|
184 | |
---|
185 | list Fin, Inf; |
---|
186 | if ( s == "ideals") |
---|
187 | // the divisor is given by two ideals |
---|
188 | { |
---|
189 | if(pr == 1) |
---|
190 | // geometric case: the divisor is given by two homogeneous |
---|
191 | // ideals I,J representing two effective divisors, |
---|
192 | // transform the divisor in two fractional ideals |
---|
193 | // Fin and Inf |
---|
194 | { |
---|
195 | setring R; |
---|
196 | list FinBasis = imap(R_aff,FinBasis); |
---|
197 | list Infin = imap(R_aff,Infin); |
---|
198 | |
---|
199 | // the transformed divisor |
---|
200 | divisorD = divisorTrans(f,divisorD,FinBasis,Infin,"ideals"); |
---|
201 | setring R_aff; |
---|
202 | list divisorD = imap(R,divisorD); |
---|
203 | } |
---|
204 | |
---|
205 | Fin = divisorD[1]; |
---|
206 | Inf = divisorD[2]; |
---|
207 | } |
---|
208 | |
---|
209 | if (s == "free") |
---|
210 | // the divisor is given by a list of places |
---|
211 | { |
---|
212 | if (pr) |
---|
213 | { |
---|
214 | setring R; |
---|
215 | list FinBasis = imap(R_aff,FinBasis); |
---|
216 | list Infin = imap(R_aff,Infin); |
---|
217 | divisorD = divisorTrans(f,divisorD,FinBasis,Infin,"free"); |
---|
218 | setring R_aff; |
---|
219 | list divisorD = imap(R,divisorD); |
---|
220 | } |
---|
221 | // computes the ideal representation from the free rep. |
---|
222 | divisorD = Free2IdealRepresentation(f_aff,FinBasis,InfBasis,divisorD); |
---|
223 | Fin = divisorD[1]; |
---|
224 | Inf = divisorD[2]; |
---|
225 | } |
---|
226 | |
---|
227 | //Compute free bases for I and J |
---|
228 | //////////////////////////////////////////////////////////// |
---|
229 | |
---|
230 | Fin = freeGenerators(f_aff,FinBasis,Fin,1); |
---|
231 | Inf = freeGenerators(f_aff,InfBasis,Inf,2); |
---|
232 | |
---|
233 | if (printlevel >= 3) |
---|
234 | { |
---|
235 | "Integral basis of the finite ideal I:";"";Fin;pause();""; |
---|
236 | "Integral basis of the infinite ideal J:";"";Inf;pause();""; |
---|
237 | } |
---|
238 | |
---|
239 | //Compute free bases for I^(-1) and J^(-1) |
---|
240 | //////////////////////////////////////////////////////////// |
---|
241 | |
---|
242 | Fin = inverseIdeal(f_aff,FinBasis,Fin,1); |
---|
243 | Inf = inverseIdeal(f_aff,InfBasis,Inf,2); |
---|
244 | |
---|
245 | if (printlevel >= 3) |
---|
246 | { |
---|
247 | "Integral basis of the inverse ideal of I:";"";Fin;pause();""; |
---|
248 | "Integral basis of the inverse ideal of J:";"";Inf;pause();""; |
---|
249 | } |
---|
250 | |
---|
251 | //Compute the transition matrix |
---|
252 | //////////////////////////////////////////////////////////// |
---|
253 | list T = transmatrix(f_aff,Inf,Fin); |
---|
254 | |
---|
255 | |
---|
256 | //Reduce the transition matrix |
---|
257 | //////////////////////////////////////////////////////////// |
---|
258 | setring Rone; |
---|
259 | list T = imap(R_aff,T); |
---|
260 | poly denT = T[2]; |
---|
261 | list reducedT = redmatrix(T[1]); |
---|
262 | matrix redT = reducedT[1]; //reduced matrix |
---|
263 | matrix trafoT = reducedT[2]; // transformation matrix |
---|
264 | |
---|
265 | |
---|
266 | //Compute the k[x]-invariants d_j |
---|
267 | //////////////////////////////////////////////////////////// |
---|
268 | |
---|
269 | list Dcoef = computeInvariants(redT,denT); |
---|
270 | |
---|
271 | if (printlevel >= 3) |
---|
272 | { |
---|
273 | "List of k[x]-invariants:";"";Dcoef;pause();""; |
---|
274 | } |
---|
275 | |
---|
276 | //Compute the transformed basis of I^(-1) |
---|
277 | //////////////////////////////////////////////////////////// |
---|
278 | |
---|
279 | setring R_aff; |
---|
280 | matrix trafoT = imap(Rone,trafoT); |
---|
281 | |
---|
282 | |
---|
283 | // compute the basis |
---|
284 | // (v_1,...,v_n) = ((v_1)',...,(v_n)')*trafoT, |
---|
285 | // where {v_i} basis of I^(-1) |
---|
286 | |
---|
287 | list final = multiplyMatrixList(trafoT,Fin[1]); |
---|
288 | |
---|
289 | if (printlevel >= 3) |
---|
290 | { |
---|
291 | "Transformed basis of I^(-1):";"";list(final,Fin[2]);pause();""; |
---|
292 | } |
---|
293 | |
---|
294 | //Compute a k-basis of L(D) |
---|
295 | //////////////////////////////////////////////////////////// |
---|
296 | |
---|
297 | number lcDen = leadcoef(Fin[2]); |
---|
298 | list rrbasis; |
---|
299 | int counter; |
---|
300 | poly ftemp; |
---|
301 | for (i=1;i<=size(Dcoef);i++) |
---|
302 | { |
---|
303 | for (j=0; j<=(-Dcoef[i]); j++) |
---|
304 | { |
---|
305 | counter++; |
---|
306 | ftemp = (var(1)^j*final[i])/lcDen; |
---|
307 | ftemp = ftemp*(1/content(ftemp)); |
---|
308 | rrbasis[counter]= ftemp; |
---|
309 | } |
---|
310 | } |
---|
311 | list RRbasis = rrbasis, Fin[2]/lcDen; |
---|
312 | |
---|
313 | if (pr) |
---|
314 | { |
---|
315 | setring R; |
---|
316 | list RRbasis = imap(R_aff,RRbasis); |
---|
317 | return(RRbasis); |
---|
318 | } |
---|
319 | return(RRbasis); |
---|
320 | } |
---|
321 | |
---|
322 | example |
---|
323 | { "EXAMPLE:"; echo=2; |
---|
324 | ring R = 0,(x,y),dp; |
---|
325 | poly f = y^2*(y-1)^3-x^5; |
---|
326 | list A1 = list(ideal(x,y-1),1),2; |
---|
327 | list A2 = list(ideal(y^2-y+1,x),1),3; |
---|
328 | list A3 = list(ideal(1,y-x),x),-2; |
---|
329 | list D = A1,A2; |
---|
330 | list E = list(A3); |
---|
331 | RiemannRochHess(f,list(D,E),"free"); |
---|
332 | ring S = 0,(x,y,z),dp; |
---|
333 | poly f = y^2*(y-1)^3-x^5; |
---|
334 | f = homog(f,z); |
---|
335 | ideal P1 = x,y-z; |
---|
336 | ideal P2 = y^2-yz+z^2,x; |
---|
337 | ideal P3 = x-y,z; |
---|
338 | list B1 = P1,2; |
---|
339 | list B2 = P2,3; |
---|
340 | list B3 = P3,-2; |
---|
341 | list Ddivisor = B1,B2,B3; |
---|
342 | RiemannRochHess(f,Ddivisor,"free"); |
---|
343 | ideal I = intersect(P1^2,P2^3); |
---|
344 | ideal J = P3^2; |
---|
345 | RiemannRochHess(f,list(I,J),"ideals"); |
---|
346 | } |
---|
347 | |
---|
348 | |
---|
349 | //////////////////////////////////////////////////////////////////// |
---|
350 | ///////////////////// Essential Static Procedures ////////////////// |
---|
351 | |
---|
352 | //_________________________________________________________________ |
---|
353 | static proc matrixrep(poly f,list intbasis,list A) |
---|
354 | "USAGE: matrixrep(f,intbasis,A);f poly,intbasis list,A list |
---|
355 | ASSUME: The base ring R must be a ring in two variables, say x,y, and |
---|
356 | f must be monic as polynomial in the second variable. |
---|
357 | List intbasis contains an integral basis (w1,...,wn) |
---|
358 | of IntCl(k[x],F) or of IntCl(k[1/x],F). |
---|
359 | List A contains two polynomials which represent an element of |
---|
360 | Quot(R[x]). |
---|
361 | RETURN: A list, say repmatrix, containing two elements. repmatrix[1] is |
---|
362 | a matrix, say M, with entries in k[x], and repmatrix[2] a |
---|
363 | polynomial, say d, such that M_a := M/d is a representation |
---|
364 | matrix of a, i.e. a*(w1,..,wn) = (w1,..,wn)*M_a |
---|
365 | " |
---|
366 | { |
---|
367 | def Rbase = basering; |
---|
368 | list rl = ringlist(Rbase); |
---|
369 | rl[2] = list(var(2), var(1)); |
---|
370 | rl[3] = list(list("lp",1:2),list("C",0)); |
---|
371 | def Rreduction = ring(rl); // make var(2) > var(1) |
---|
372 | |
---|
373 | rl = ringlist(Rbase); |
---|
374 | rl[1] = list(char(Rbase),list(var(1)),list(list("dp",1)),ideal(0)); |
---|
375 | rl[2] = list(var(2)); |
---|
376 | rl[3] = list(list("dp",1),list("C",0)); |
---|
377 | def QF = ring(rl); // make var(1) transcendental |
---|
378 | |
---|
379 | int i,j,k,l; |
---|
380 | ideal IntBasis = intbasis[1]; |
---|
381 | poly d = intbasis[2]; |
---|
382 | int sizeIntBasis = size(IntBasis); |
---|
383 | |
---|
384 | setring Rreduction; //var(2) > var(1) |
---|
385 | list A = imap(Rbase,A); |
---|
386 | ideal IntBasis = imap(Rbase,IntBasis); |
---|
387 | ideal fred = std(imap(Rbase,f)); |
---|
388 | IntBasis = reduce(IntBasis,fred); |
---|
389 | //coeffiecient matrix w. r. t. the integral basis |
---|
390 | matrix M = coeffs(IntBasis,var(1)); |
---|
391 | |
---|
392 | setring QF; |
---|
393 | matrix M = imap(Rreduction,M); |
---|
394 | poly d = imap(Rbase,d); |
---|
395 | // coefficient matrix with denominator |
---|
396 | M = 1/d*M; |
---|
397 | // inverse of the coefficient matrix |
---|
398 | list L = luinverse(M); |
---|
399 | |
---|
400 | setring Rreduction; |
---|
401 | list multiA; |
---|
402 | for(i = 1;i<=sizeIntBasis;i++) |
---|
403 | { |
---|
404 | multiA[i] = reduce(A[1]*IntBasis[i],fred); |
---|
405 | } |
---|
406 | // coefficient matrix w.r.t. A[1]*IntBasis |
---|
407 | matrix multiM = coeffs(ideal(multiA[1..sizeIntBasis]),var(1)); |
---|
408 | |
---|
409 | // for positive characteristic - necessary if all coefficients |
---|
410 | // reduce to zero |
---|
411 | if (nrows(multiM)!= sizeIntBasis || ncols(multiM)!= sizeIntBasis) |
---|
412 | { |
---|
413 | multiM = zerosM(sizeIntBasis); |
---|
414 | } |
---|
415 | |
---|
416 | setring QF; |
---|
417 | list A = imap(Rbase,A); |
---|
418 | matrix multiM = imap(Rreduction,multiM); |
---|
419 | // multiply the both coefficient matrices |
---|
420 | matrix D = L[2]*multiM; |
---|
421 | D = 1/(d*A[2])*D; |
---|
422 | number a = contentMatrix(D); |
---|
423 | number numera = numerator(a); |
---|
424 | number denoma = denominator(a); |
---|
425 | D = (1/a)*D; |
---|
426 | |
---|
427 | setring Rbase; |
---|
428 | matrix D = imap(QF,D); |
---|
429 | poly numera = imap(QF,numera); |
---|
430 | poly denoma = imap(QF,denoma); |
---|
431 | number lcd= leadcoef(denoma); |
---|
432 | list repmatrix; |
---|
433 | repmatrix = list((numera/lcd)*D,denoma/lcd); |
---|
434 | return(repmatrix); |
---|
435 | } |
---|
436 | |
---|
437 | //__________________________________________________________________ |
---|
438 | static proc freeGenerators(poly f,list intbasis,list fracIdeal,int opt) |
---|
439 | "USAGE: freeGenerators(f,intbasis,fracIdeal,opt); f polynomial, |
---|
440 | intbasis list, fracIdeal list, opt integer |
---|
441 | ASSUME: The base ring R must be a ring in two variables, say x,y, and |
---|
442 | f must be monic as polynomial in the second variable. |
---|
443 | List intbasis contains an integral basis of: |
---|
444 | - IntCl(k[x],F), if opt = 1; |
---|
445 | - IntCl(k[1/x],F), if opt = 2; |
---|
446 | List fracIdeals contains two elements, an ideal I and a |
---|
447 | polynomial d representing the fractional ideal I/d |
---|
448 | RETURN: List of size 2, say freeGens, where freeGens[1] is an |
---|
449 | ideal J and freeGens[2] a polynomial D such that |
---|
450 | J[1]/D, ..., J[n]/D is a basis of I/d as free k[x]-(opt = 1) |
---|
451 | resp. k[1/x]-module (opt = 2) |
---|
452 | " |
---|
453 | { |
---|
454 | def Rbase = basering; |
---|
455 | list rl = ringlist(Rbase); |
---|
456 | rl[3] = list(list("C",1:2),list("dp",0)); |
---|
457 | def Rhermite = ring(rl); |
---|
458 | rl = ringlist(Rbase); |
---|
459 | rl[2] = list(var(1)); |
---|
460 | rl[3] = list(list("lp",1),list("C",0)); |
---|
461 | def Ronevar = ring(rl); // only one variable |
---|
462 | |
---|
463 | int i,j,k; |
---|
464 | ideal I = fracIdeal[1]; |
---|
465 | list T,polynomialT,denominatorT; |
---|
466 | |
---|
467 | for(i=1; i<=size(I);i++) |
---|
468 | { |
---|
469 | // compute representation matrices for every generator of I |
---|
470 | T[i] = matrixrep(f,intbasis,list(I[i],fracIdeal[2])); |
---|
471 | // list containing only the polynomial matrices |
---|
472 | polynomialT[i] = T[i][1]; |
---|
473 | // list containing the corresponding denominators |
---|
474 | denominatorT[i] = poly(T[i][2]); |
---|
475 | } |
---|
476 | poly commonden = lcm(denominatorT[1..size(denominatorT)]); |
---|
477 | for(i=1; i<=size(polynomialT);i++) |
---|
478 | { |
---|
479 | // multiply with common denominator |
---|
480 | polynomialT[i] = (commonden/denominatorT[i])*polynomialT[i]; |
---|
481 | } |
---|
482 | matrix M = concat(polynomialT); |
---|
483 | |
---|
484 | if (opt == 1) // compute a k[x]-basis |
---|
485 | { |
---|
486 | setring Rhermite; |
---|
487 | matrix M = imap(Rbase,M); |
---|
488 | |
---|
489 | // compute n generators of the big rep. matrix via Groebner basis |
---|
490 | module Mfinal = std(module(M)); |
---|
491 | |
---|
492 | setring Rbase; |
---|
493 | module Mfinal = imap(Rhermite,Mfinal); |
---|
494 | matrix P = Mfinal; |
---|
495 | } |
---|
496 | else // compute a k[1/x]-basis |
---|
497 | { |
---|
498 | list P1 = normalFormInf(list(M,commonden),"free"); |
---|
499 | matrix P = P1[1]; |
---|
500 | commonden = P1[2]; |
---|
501 | } |
---|
502 | |
---|
503 | ideal IB = intbasis[1]; |
---|
504 | list Z = multiplyMatrixList(P,IB); |
---|
505 | ideal J = ideal(Z[1..size(IB)]); |
---|
506 | poly gcdJ = idealGcd(J); |
---|
507 | poly D = intbasis[2]*commonden; |
---|
508 | poly comgcd = gcd(gcdJ,D); |
---|
509 | // cancel out common divisor |
---|
510 | J = J/comgcd; |
---|
511 | D = D/comgcd; |
---|
512 | list freeGens = J,D; |
---|
513 | return(freeGens); |
---|
514 | } |
---|
515 | |
---|
516 | //___________________________________________________________________ |
---|
517 | static proc inverseIdeal(poly f,list intbasis,list fracIdeal,int opt) |
---|
518 | "USAGE: inverseIdeal(f,intbasis,fracIdeal,opt); f polynomial, |
---|
519 | intbasis list, A list, opt integer |
---|
520 | ASSUME: The base ring R must be a ring in two variables, say x,y, |
---|
521 | and f must be monic as polynomial in the second variable. |
---|
522 | List intbasis contains an integral basis of: |
---|
523 | - IntCl(k[x],F), if opt = 1; |
---|
524 | - IntCl(k[1/x],F), if opt = 2; |
---|
525 | List fracIdeal contains two elements, an ideal I and a |
---|
526 | polynomial d representing the fractional ideal I/d |
---|
527 | RETURN: A list inverseId of size 2. inverseId[1] contains an ideal, |
---|
528 | say J, inverseId[2] a polynomial, say D. |
---|
529 | Then J/D is a free k[x]- (opt = 1) resp. k[1/x]-basis |
---|
530 | (opt = 2) for the inverse ideal of I/d |
---|
531 | " |
---|
532 | { |
---|
533 | def Rbase = basering; |
---|
534 | list rl = ringlist(Rbase); |
---|
535 | rl[2] = list(var(2), var(1)); |
---|
536 | rl[3] = list(list("lp",1:2),list("C",0)); |
---|
537 | def Rreduction = ring(rl); // make var(2) > var(1) |
---|
538 | |
---|
539 | rl = ringlist(Rbase); |
---|
540 | rl[2] = list(var(1)); |
---|
541 | rl[3] = list(list("c",0),list("lp",1)); |
---|
542 | def Rhelp = ring(rl); |
---|
543 | // ring with one variable, need for Hermite Normal Form |
---|
544 | |
---|
545 | rl = ringlist(Rbase); |
---|
546 | rl[1] = list(char(Rbase),list(var(1)),list(list("dp",1)),ideal(0)); |
---|
547 | rl[2] = list(var(2)); |
---|
548 | rl[3] = list(list("dp",1),list("C",0)); |
---|
549 | def QF = ring(rl); // make var(1) transcendental |
---|
550 | |
---|
551 | int i,j; |
---|
552 | ideal I = fracIdeal[1]; |
---|
553 | list T,polynomialT,denominatorT; |
---|
554 | for(i=1; i<=size(I);i++) |
---|
555 | { |
---|
556 | T[i] = matrixrep(f,intbasis,list(I[i],fracIdeal[2])); |
---|
557 | polynomialT[i] = transpose(T[i][1]); |
---|
558 | denominatorT[i] = poly(T[i][2]); |
---|
559 | } |
---|
560 | poly commonden = lcm(denominatorT[1..size(denominatorT)]); |
---|
561 | |
---|
562 | for(i=1; i<=size(polynomialT);i++) |
---|
563 | { |
---|
564 | polynomialT[i] = (commonden/denominatorT[i])*polynomialT[i]; |
---|
565 | } |
---|
566 | |
---|
567 | // as in freeGenerators: compute big representation matrix |
---|
568 | matrix M = concat(polynomialT); |
---|
569 | |
---|
570 | if (opt == 1) // inverse ideal in ICl(k[x],F) |
---|
571 | { |
---|
572 | setring Rhelp; |
---|
573 | matrix M = imap(Rbase,M); |
---|
574 | // computes Hermite Normal Form of the matrix via Groebner basis |
---|
575 | module TM = std(module(M)); |
---|
576 | setring Rbase; |
---|
577 | matrix TM = imap(Rhelp,TM); |
---|
578 | M = transpose(TM); |
---|
579 | } |
---|
580 | else // inverse idal in ICl(K[1/x],F) |
---|
581 | { |
---|
582 | list P1 = normalFormInf(list(M,commonden),"inverse"); |
---|
583 | M = P1[1]; |
---|
584 | commonden = P1[2]; |
---|
585 | } |
---|
586 | ideal IB = intbasis[1]; |
---|
587 | list Minverse = inverse_B(M); |
---|
588 | |
---|
589 | setring QF; |
---|
590 | list Minverse = imap(Rbase,Minverse); |
---|
591 | poly commonden = imap(Rbase,commonden); |
---|
592 | ideal IB = imap(Rbase,IB); |
---|
593 | matrix invM = (commonden/Minverse[2])*Minverse[1]; |
---|
594 | ideal Z; |
---|
595 | poly contZ; |
---|
596 | int n = size(IB); |
---|
597 | for (i=1; i<=n; i++) |
---|
598 | { |
---|
599 | Z[i]=0; |
---|
600 | for(j=1; j<=n; j++) |
---|
601 | { |
---|
602 | Z[i] = Z[i] + IB[j]*invM[j,i]; |
---|
603 | } |
---|
604 | } |
---|
605 | matrix Z1[1][n] = Z; |
---|
606 | number a = contentMatrix(Z1); |
---|
607 | number numera = numerator(a); |
---|
608 | number denoma = denominator(a); |
---|
609 | Z = (1/a)*Z; |
---|
610 | |
---|
611 | setring Rbase; |
---|
612 | ideal Z = imap(QF,Z); |
---|
613 | poly numera = imap(QF,numera); |
---|
614 | poly denoma = imap(QF,denoma); |
---|
615 | list inverseId; |
---|
616 | poly D = gcd(numera,poly(intbasis[2])); |
---|
617 | inverseId = list(Z*(numera/D),(intbasis[2]/D)*denoma); |
---|
618 | return(inverseId); |
---|
619 | } |
---|
620 | |
---|
621 | //__________________________________________________________________ |
---|
622 | static proc transmatrix(poly f,list A,list B) |
---|
623 | "USAGE: transmatrix(f,A,B); f polynomial, A list, B list |
---|
624 | ASSUME: The base ring R must be a ring in two variables, say x,y, |
---|
625 | and f must be monic as polynomial in the second variable. |
---|
626 | List A (resp. B) consists of two elements. |
---|
627 | A[1] (resp. B[1]) a list of n elements of k[x], |
---|
628 | A[2] (resp. B[2]) common denominators such that the |
---|
629 | corresponding fractions are k(x)-linear independent |
---|
630 | RETURN: List transm of size 2. transm[1] is a matrix with |
---|
631 | entries in k[x], say M, and transm[2] a polynomial in k[x], |
---|
632 | say D, such that (M/D) is the transition matrix from A to B, |
---|
633 | i.e, A[1]/A[2]*(M/D) = B[1]/B[2] |
---|
634 | " |
---|
635 | { |
---|
636 | def Rbase = basering; |
---|
637 | list rl = ringlist(Rbase); |
---|
638 | rl[2] = list(var(2), var(1)); |
---|
639 | rl[3] = list(list("lp",1:2),list("C",0)); |
---|
640 | def Rreduction = ring(rl); // make var(2) > var(1) |
---|
641 | rl = ringlist(Rbase); |
---|
642 | rl[1] = list(char(Rbase),list(var(1)),list(list("dp",1)),ideal(0)); |
---|
643 | rl[2] = list(var(2)); |
---|
644 | rl[3] = list(list("dp",1),list("C",0)); |
---|
645 | def QF = ring(rl); // make var(1) transcendental |
---|
646 | |
---|
647 | ideal a = A[1]; |
---|
648 | ideal b = B[1]; |
---|
649 | |
---|
650 | setring Rreduction; |
---|
651 | ideal a = imap(Rbase,a); |
---|
652 | ideal b = imap(Rbase,b); |
---|
653 | ideal fred = std(imap(Rbase,f)); |
---|
654 | a = reduce(a,fred); |
---|
655 | b = reduce(b,fred); |
---|
656 | // coefficient matrices w.r.t. 1,...,y^(n-1) |
---|
657 | matrix M_a = coeffs(a,var(1)); |
---|
658 | matrix M_b = coeffs(b,var(1)); |
---|
659 | |
---|
660 | setring QF; |
---|
661 | list A = imap(Rbase,A); |
---|
662 | list B = imap(Rbase,B); |
---|
663 | matrix M_a = imap(Rreduction,M_a); |
---|
664 | matrix M_b = imap(Rreduction,M_b); |
---|
665 | poly d_a = A[2]; |
---|
666 | poly d_b = B[2]; |
---|
667 | matrix N_a = 1/d_a*M_a; |
---|
668 | matrix N_b = 1/d_b*M_b; |
---|
669 | list L = luinverse(N_a); |
---|
670 | matrix M = L[2]*N_b; |
---|
671 | number cont = contentMatrix(M); |
---|
672 | number numer = numerator(cont); |
---|
673 | number denom = denominator(cont); |
---|
674 | M = (1/cont)*M; |
---|
675 | |
---|
676 | setring Rbase; |
---|
677 | matrix M = imap(QF,M); |
---|
678 | poly numer = imap(QF,numer); |
---|
679 | poly denom = imap(QF,denom); |
---|
680 | list transm; |
---|
681 | transm = numer*M,denom; |
---|
682 | return(transm); |
---|
683 | } |
---|
684 | |
---|
685 | //___________________________________________________________________ |
---|
686 | static proc maxorderInfinite(poly f, int #) |
---|
687 | "USAGE: maxorderInfinite(f[,#]); f polynomial, # integer |
---|
688 | ASSUME: The base ring must be a ring in two variables, say x,y. |
---|
689 | - f is an irreducible polynomial and the |
---|
690 | second variable describes the integral element, say y. |
---|
691 | The degree of f in y is n. |
---|
692 | - optional integer #. If # is given, then # = 1 and the |
---|
693 | infinite places of the function field F are computed if |
---|
694 | Kummer's theorem applies. |
---|
695 | RETURN: - If # is not given: |
---|
696 | A list InfBasis of size 1. InfBasis[1] contains an |
---|
697 | integral basis of IntCl(k[1/x],F), stored in a list with |
---|
698 | n = [F:k(x)] elements and a common denominator. |
---|
699 | - If # = 1: |
---|
700 | A list InfBasis of size 3. |
---|
701 | InfBasis[1] contains an integral basis of IntCl(k[1/x],F). |
---|
702 | InfBasis[2] is an integer k. |
---|
703 | If we set f1 = t^(kn)*f(1/t,y) and define u = y*t^k, |
---|
704 | then f1 is a monic polynomial in u with coefficients in |
---|
705 | k[t]. |
---|
706 | InfBasis[3] a list containing the infinite places of |
---|
707 | the function field, if the compuation was possible. Else |
---|
708 | the list is empty. |
---|
709 | " |
---|
710 | { |
---|
711 | def Rbase = basering; |
---|
712 | int k,l,i; |
---|
713 | |
---|
714 | //Necessary rings for substitution |
---|
715 | ////////////////////////////////////////////////////////////////// |
---|
716 | |
---|
717 | list rl = ringlist(Rbase); |
---|
718 | rl[1] = list(char(Rbase),list("t"),list(list("dp",1)),ideal(0)); |
---|
719 | rl[2] = list(var(1),var(2),"u"); |
---|
720 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
721 | def Rnew = ring(rl); |
---|
722 | |
---|
723 | rl = ringlist(Rbase); |
---|
724 | rl[2] = list("t","u"); |
---|
725 | def Rintegral = ring(rl); |
---|
726 | // ring with the new two variables, for the computation of the integral basis |
---|
727 | |
---|
728 | //two additional rings needed for changing t,u back to x,y |
---|
729 | rl = ringlist(Rbase); |
---|
730 | rl[2] = list(var(1),var(2),"t","u"); |
---|
731 | rl[3] = list(list("dp",1:4),list("C",0)); |
---|
732 | def Rhelp = ring(rl); |
---|
733 | |
---|
734 | rl = ringlist(Rbase); |
---|
735 | rl[1] = list(char(Rbase),list(var(1)),list(list("dp",1)),ideal(0)); |
---|
736 | rl[2] = list(var(2),"t","u"); |
---|
737 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
738 | def Rhelp1 = ring(rl); |
---|
739 | |
---|
740 | //Substitute the variables |
---|
741 | ////////////////////////////////////////////////////////////////// |
---|
742 | |
---|
743 | setring Rnew; //1 additional parameter, 1 add. variable |
---|
744 | poly f = imap(Rbase,f); |
---|
745 | poly f1 = subst(f,var(1),(1/t)); // subst. x by 1/t |
---|
746 | |
---|
747 | matrix Cf = coeffs(f,var(2)); |
---|
748 | int n = nrows(Cf)-1; // degree of f in y |
---|
749 | list degs; |
---|
750 | for(i=1; i<=n; i++) |
---|
751 | { |
---|
752 | // f = y^n + a_(n-1)(x)*y^(n-1)+.... |
---|
753 | degs[i] = list(i-1,deg(Cf[i,1])); // degs[i] =[i,deg(a_i(x))] |
---|
754 | } |
---|
755 | |
---|
756 | // find the minimum multiple of n such that t^(k*n)*f(1/t,y) |
---|
757 | // is monic after replacing y*t^k by u |
---|
758 | while (l==0) |
---|
759 | { |
---|
760 | l=1; |
---|
761 | k=k+1; |
---|
762 | for(i=1; i<=size(degs); i++) |
---|
763 | { |
---|
764 | if(degs[i][2] != -1) // if (a_i(x) != 0) |
---|
765 | { |
---|
766 | if((k*n-degs[i][2])<(k*degs[i][1])) |
---|
767 | { |
---|
768 | l=0; break; |
---|
769 | } |
---|
770 | } |
---|
771 | } |
---|
772 | } |
---|
773 | |
---|
774 | f1 = f1*t^(k*n); |
---|
775 | poly v = u/(t^k); |
---|
776 | poly fintegral = subst(f1,var(2),v); // replace y by u/t^k |
---|
777 | |
---|
778 | // compute integral basis of the new polynomial in t,u |
---|
779 | setring Rintegral; |
---|
780 | poly fintegral = imap(Rnew,fintegral); |
---|
781 | list Z = integralBasis(fintegral,2); |
---|
782 | |
---|
783 | //Check if Kummer's Theorem apply |
---|
784 | ////////////////////////////////////////////////////////////////// |
---|
785 | if (# == 1) |
---|
786 | { |
---|
787 | poly fInfinity = substitute(fintegral,t,0); |
---|
788 | list factors = factorize(fInfinity,2); |
---|
789 | list degsinf = factors[2]; |
---|
790 | int kum = 1; |
---|
791 | for (i=1; i<=size(degsinf[1]);i++) |
---|
792 | { |
---|
793 | // check if all exponents are 1 |
---|
794 | if (degsinf[1][i] != 1) |
---|
795 | { |
---|
796 | kum = 0; |
---|
797 | break; |
---|
798 | } |
---|
799 | } |
---|
800 | if (kum == 0) |
---|
801 | // check if 1,...,y^{n-1} is an integral basis |
---|
802 | { |
---|
803 | if (Z[2] == 1) |
---|
804 | { |
---|
805 | kum = 1; |
---|
806 | for (i = 1; i<=size(Z[1]); i++) |
---|
807 | { |
---|
808 | if (Z[1][i] != u^(i-1)) |
---|
809 | { |
---|
810 | kum = 0; |
---|
811 | break; |
---|
812 | } |
---|
813 | } |
---|
814 | } |
---|
815 | } |
---|
816 | } |
---|
817 | |
---|
818 | //Reverse the substitution |
---|
819 | ////////////////////////////////////////////////////////////////// |
---|
820 | |
---|
821 | setring Rnew; // t as parameter, u as variable |
---|
822 | list Z = imap(Rintegral,Z); |
---|
823 | ideal IB = Z[1]; |
---|
824 | poly d = Z[2]; |
---|
825 | IB = subst(IB,u,var(2)*(t^k)); |
---|
826 | d = subst(d,u,var(2)*(t^k)); // replace u by y*t^k |
---|
827 | |
---|
828 | if (# == 1) |
---|
829 | { |
---|
830 | if (kum == 1) |
---|
831 | { |
---|
832 | list factors = imap(Rintegral,factors); |
---|
833 | list Zfacs = factors[1]; |
---|
834 | ideal I = substitute(ideal(Zfacs[1..size(Zfacs)]),u,var(2)*(t^k)); |
---|
835 | Zfacs = I[1..size(I)]; |
---|
836 | } |
---|
837 | } |
---|
838 | |
---|
839 | setring Rhelp; // t as variable |
---|
840 | ideal IB = imap(Rnew,IB); |
---|
841 | poly d = imap(Rnew,d); |
---|
842 | matrix C = coeffs(IB,t); |
---|
843 | matrix Cden = coeffs(d,t); |
---|
844 | int p = nrows(C)-1; // maximum degree of IB in t |
---|
845 | int pden = nrows(Cden)-1; |
---|
846 | |
---|
847 | setring Rhelp1; // var(1) as parameter, t as variable |
---|
848 | ideal IB = imap(Rhelp,IB); |
---|
849 | poly d = imap(Rhelp,d); |
---|
850 | matrix Cden = imap(Rhelp,Cden); |
---|
851 | IB = subst(IB,t,(1/par(1))); // replace t by 1/x |
---|
852 | d = subst(d,t,(1/par(1))); |
---|
853 | IB = IB*par(1)^p; |
---|
854 | d = d*par(1)^pden; |
---|
855 | |
---|
856 | //Compute the infinite places |
---|
857 | ////////////////////////////////////////////////////////////////// |
---|
858 | if (# == 1) |
---|
859 | { |
---|
860 | if (kum == 1) |
---|
861 | { |
---|
862 | list Zfacs = imap(Rnew,Zfacs); |
---|
863 | list pfacs; |
---|
864 | for (i=1;i<=size(Zfacs);i++) |
---|
865 | { |
---|
866 | pfacs[i]=deg(substitute(Zfacs[i],var(1),1)); |
---|
867 | Zfacs[i]=ideal(Zfacs[i],t); |
---|
868 | Zfacs[i]=substitute(Zfacs[i],t,(1/par(1))); |
---|
869 | Zfacs[i]=list(Zfacs[i]*par(1)^pfacs[i],par(1)^pfacs[i]); |
---|
870 | } |
---|
871 | } |
---|
872 | } |
---|
873 | |
---|
874 | setring Rbase; |
---|
875 | ideal IB = imap(Rhelp1,IB); |
---|
876 | poly d = imap(Rhelp1,d); |
---|
877 | poly den; |
---|
878 | if(p >= pden) |
---|
879 | { |
---|
880 | den = d*var(1)^(p-pden); |
---|
881 | } |
---|
882 | else |
---|
883 | { |
---|
884 | IB = IB*var(1)*(pden-p); |
---|
885 | den = d; |
---|
886 | } |
---|
887 | list final = IB,den; |
---|
888 | |
---|
889 | if ( # == 1) |
---|
890 | { |
---|
891 | list InfPlaces; |
---|
892 | if (kum == 1) |
---|
893 | { |
---|
894 | list Zfacs = imap(Rhelp1,Zfacs); |
---|
895 | InfPlaces = Zfacs; |
---|
896 | list InfBasis = final,k,InfPlaces; |
---|
897 | return(InfBasis); |
---|
898 | } |
---|
899 | else |
---|
900 | { |
---|
901 | return(final,k,InfPlaces); |
---|
902 | } |
---|
903 | } |
---|
904 | |
---|
905 | list InfBasis = list(final); |
---|
906 | return(InfBasis); |
---|
907 | } |
---|
908 | |
---|
909 | //__________________________________________________________________ |
---|
910 | static proc redmatrix(matrix A) |
---|
911 | "USAGE: redmatrix(A); A matrix |
---|
912 | ASSUME: The basering must be a polynomial ring in one variable, |
---|
913 | say k[x], A a mxn-matrix with entries in k[x] |
---|
914 | RETURN: A list containing two matrices, say A' and T such that |
---|
915 | A' is the reduced matrix of A and T the |
---|
916 | transformation matrix, i.e., A' = A*T. |
---|
917 | THEORY: Consider the columns of the matrix A as a basis |
---|
918 | of a lattice. Denote the column degree of the j-th column |
---|
919 | by deg(j). Create vectors hc(j) in k^m containing the |
---|
920 | coefficient of the deg(j)-th power of x of column j. |
---|
921 | A basis is called reduced if the {hc(j)} are linearly |
---|
922 | independent. |
---|
923 | If the basis is non-reduced we proceed reduction steps |
---|
924 | until the above criterion is satisfied. |
---|
925 | In a reduction step of column j we add a certain k[x]- |
---|
926 | combination of the other columns to column j, such that |
---|
927 | the column-degree decreases. |
---|
928 | " |
---|
929 | { |
---|
930 | if (nvars(basering) != 1) |
---|
931 | { |
---|
932 | ERROR("The base ring must be a ring in one variables."); |
---|
933 | } |
---|
934 | int m = nrows(A); |
---|
935 | int n = ncols(A); |
---|
936 | matrix T = unitmat(m); |
---|
937 | list d; |
---|
938 | int i,j,counter; |
---|
939 | int h = 1; |
---|
940 | int pos, maxd; |
---|
941 | while(h) |
---|
942 | { |
---|
943 | for (j=1; j<=n; j++) |
---|
944 | { |
---|
945 | d[j]=colDegree(A,j); |
---|
946 | } |
---|
947 | matrix M[m][n]; |
---|
948 | for (j=1; j<=n; j++) |
---|
949 | { |
---|
950 | for (i=1; i<=m; i++) |
---|
951 | { |
---|
952 | M[i,j]=coefficientAt(A,j,d[j])[i]; // computes the hc(j) |
---|
953 | } |
---|
954 | } |
---|
955 | module S = M; |
---|
956 | matrix N = matrix(syz(M)); |
---|
957 | if(rank(N)==0) // hc(j) are linearly independent |
---|
958 | { |
---|
959 | h=0; |
---|
960 | } |
---|
961 | else |
---|
962 | { |
---|
963 | h = ncols(N); |
---|
964 | } |
---|
965 | list ind; |
---|
966 | list degind; |
---|
967 | |
---|
968 | for (j=1; j<=Min(intvec(1,h));j++) |
---|
969 | // minimum needed to avoid unnecessary first reduction |
---|
970 | // step in the reduced case |
---|
971 | { |
---|
972 | for (i=1; i<= ncols(M); i++) |
---|
973 | { |
---|
974 | if(N[i,j]!=0) |
---|
975 | { |
---|
976 | counter++; |
---|
977 | // position of non-zero entry |
---|
978 | ind[counter]=i; |
---|
979 | // degree of non-zero entry |
---|
980 | degind[counter]=d[i]; |
---|
981 | } |
---|
982 | } |
---|
983 | counter = 0; |
---|
984 | // determine the maximal column degree |
---|
985 | maxd = Max(degind); |
---|
986 | i=0; |
---|
987 | // find the column, where the maximum is attained |
---|
988 | while(i<=size(degind)) |
---|
989 | { |
---|
990 | i=i+1; |
---|
991 | if (degind[i] == maxd) |
---|
992 | { |
---|
993 | pos = i; |
---|
994 | break; |
---|
995 | } |
---|
996 | } |
---|
997 | pos = ind[pos]; |
---|
998 | for (i=1;i<=size(ind);i++) |
---|
999 | { |
---|
1000 | if (ind[i] != pos) |
---|
1001 | { |
---|
1002 | matrix A1=addcol(A,ind[i],N[ind[i],j]/N[pos,j]*var(1)^(maxd-d[ind[i]]),pos); |
---|
1003 | matrix T1=addcol(T,ind[i],N[ind[i],j]/N[pos,j]*var(1)^(maxd-d[ind[i]]),pos); |
---|
1004 | A=A1; |
---|
1005 | T=T1; |
---|
1006 | kill A1; |
---|
1007 | kill T1; |
---|
1008 | } |
---|
1009 | } |
---|
1010 | kill degind,ind; |
---|
1011 | } |
---|
1012 | kill M,N,S; |
---|
1013 | } |
---|
1014 | list final = A,T; |
---|
1015 | return(final); |
---|
1016 | } |
---|
1017 | |
---|
1018 | //__________________________________________________________________ |
---|
1019 | static proc normalFormInf(list K,string sopt) |
---|
1020 | "USAGE: normalFormInf(K,sopt); K list, sopt string |
---|
1021 | ASSUME: The basering must be a polynomial ring in one variable, |
---|
1022 | say k[x]. |
---|
1023 | K is a list of size 2. K[1] is a matrix, say A, with entries in |
---|
1024 | k[x] and K[2] a polynomial, say den. |
---|
1025 | A/den is a matrix with entries in k[1/x]. |
---|
1026 | sopt is either "free" or "inverse" |
---|
1027 | RETURN: A list, say norFormInf, of size 2. norFormInf[1] is a matrix, |
---|
1028 | say M, and norFormInf[2] a polynomial, say D. |
---|
1029 | The matrix M/D is a normal form of the matrix A/den which is |
---|
1030 | needed for the |
---|
1031 | - if (sopt = "free") |
---|
1032 | computation of an int. basis of an ideal of IntCl[k[1/x],F) |
---|
1033 | - if (sopt = "inverse") |
---|
1034 | computation of a int. basis of the inverse of an ideal |
---|
1035 | of IntCl[k[1/x],F) |
---|
1036 | THEORY For computing a normal form with coefficients. in k[1/x] |
---|
1037 | we have to replace 1/x by another variable, compute the normal |
---|
1038 | with the aid of Groebner bases and reverse the substitution. |
---|
1039 | " |
---|
1040 | { |
---|
1041 | def Rbase = basering; |
---|
1042 | |
---|
1043 | //_____________Necessary rings for substitution___________________ |
---|
1044 | |
---|
1045 | list rl = ringlist(Rbase); |
---|
1046 | rl[1] =list(char(Rbase),list("helppar"),list(list("dp",1)),ideal(0)); |
---|
1047 | def QF = ring(rl); // basering with additional parameter helppar |
---|
1048 | |
---|
1049 | rl = ringlist(Rbase); |
---|
1050 | rl[2] = list(var(1),var(2),"helppar"); |
---|
1051 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
1052 | def Rhelp = ring(rl); // basering with additional variable helppar |
---|
1053 | |
---|
1054 | rl = ringlist(Rbase); |
---|
1055 | rl[2] = list("helppar"); |
---|
1056 | rl[3] = list(list("c",0),list("lp",1)); |
---|
1057 | def Rhelpinv = ring(rl); // ring with one variable, ordering c |
---|
1058 | |
---|
1059 | rl = ringlist(Rbase); |
---|
1060 | rl[2] = list("helppar"); |
---|
1061 | rl[3] = list(list("C",0),list("lp",1)); |
---|
1062 | def Rhelpfree = ring(rl); // ring with one variable, ordering C |
---|
1063 | |
---|
1064 | setring QF; |
---|
1065 | list K = imap(Rbase, K); |
---|
1066 | matrix A = K[1]; |
---|
1067 | poly den = K[2]; |
---|
1068 | matrix A1 = subst(A,var(1),1/helppar); |
---|
1069 | poly den1 = subst(den, var(1),1/helppar); |
---|
1070 | A1 = (1/den1)*A1; |
---|
1071 | |
---|
1072 | // transform into a polynomial matrix with entries |
---|
1073 | // in k[helppar] |
---|
1074 | number a = contentMatrix(A1); |
---|
1075 | number numera = numerator(a); |
---|
1076 | number denoma = denominator(a); |
---|
1077 | A1 = (1/a)*A1; // A1 matrix with entries in k[helppar] |
---|
1078 | |
---|
1079 | setring Rhelp; |
---|
1080 | matrix A1 = imap(QF,A1); |
---|
1081 | poly numera = imap(QF,numera); |
---|
1082 | poly denoma = imap(QF,denoma); |
---|
1083 | number lcd= leadcoef(denoma); |
---|
1084 | A1 = (numera/lcd)*A1; |
---|
1085 | poly den = denoma/lcd; |
---|
1086 | |
---|
1087 | //_____________Compute a normal form___________________ |
---|
1088 | if (sopt == "inverse") |
---|
1089 | { |
---|
1090 | setring Rhelpinv; |
---|
1091 | matrix A1 = imap(Rhelp,A1); |
---|
1092 | module S = std(module(A1)); |
---|
1093 | A1 = S; |
---|
1094 | list mat; |
---|
1095 | int i; |
---|
1096 | int nc = ncols(A1); |
---|
1097 | // change order of columns to obtain a lower triangular |
---|
1098 | // matrix |
---|
1099 | for(i=1; i<=nc; i++) |
---|
1100 | { |
---|
1101 | mat[i]=A1[nc-i+1]; |
---|
1102 | } |
---|
1103 | A1= concat(mat); |
---|
1104 | setring Rhelp; |
---|
1105 | A1 = imap(Rhelpinv,A1); |
---|
1106 | matrix A = transpose(A1); |
---|
1107 | } |
---|
1108 | if (sopt == "free") |
---|
1109 | { |
---|
1110 | setring Rhelpfree; |
---|
1111 | matrix A1 = imap(Rhelp,A1); |
---|
1112 | module S = std(module(A1)); |
---|
1113 | A1 = S; |
---|
1114 | setring Rhelp; |
---|
1115 | matrix A = imap(Rhelpfree,A1); |
---|
1116 | } |
---|
1117 | |
---|
1118 | //_____________Reverse substitution__________________ |
---|
1119 | intvec v = 1; |
---|
1120 | // change variables back from helppar to x |
---|
1121 | def Rpar = vars2pars(v); |
---|
1122 | setring Rpar; |
---|
1123 | poly den = imap(Rhelp,den); |
---|
1124 | matrix A = imap(Rhelp,A); |
---|
1125 | A = subst(A,helppar,1/par(1)); |
---|
1126 | den = subst(den, helppar,1/par(1)); |
---|
1127 | A = (1/den)*A; |
---|
1128 | number a = contentMatrix(A); |
---|
1129 | number numera = numerator(a); |
---|
1130 | number denoma = denominator(a); |
---|
1131 | A = (1/a)*A; |
---|
1132 | |
---|
1133 | setring Rbase; |
---|
1134 | matrix A = imap(Rpar,A); |
---|
1135 | poly numera = imap(Rpar,numera); |
---|
1136 | poly denoma = imap(Rpar,denoma); |
---|
1137 | number lcd= leadcoef(denoma); |
---|
1138 | list final; |
---|
1139 | final = list((numera/lcd)*A,denoma/lcd); |
---|
1140 | |
---|
1141 | return(final); |
---|
1142 | } |
---|
1143 | |
---|
1144 | ////////////////////////////////////////////////////////////////// |
---|
1145 | ///////////////////// AUXILARY STATIC PROCEDURES ///////////////// |
---|
1146 | |
---|
1147 | //_________________________________________________________________ |
---|
1148 | static proc colDegree(matrix A, int j) |
---|
1149 | "USAGE: colDegree(A,j); A matrix, j integer |
---|
1150 | ASSUME: Base ring must be a polynomial ring in one variable, |
---|
1151 | integer j describes a column of A |
---|
1152 | RETURN: A list, containing the maximum degree occurring in the |
---|
1153 | entries of the j-th column of A and the corresponding row |
---|
1154 | " |
---|
1155 | { |
---|
1156 | if(nvars(basering) != 1) |
---|
1157 | { |
---|
1158 | ERROR("The base ring must be a ring in one variables."); |
---|
1159 | } |
---|
1160 | int d; |
---|
1161 | int nc = ncols(A); |
---|
1162 | int nr = nrows(A); |
---|
1163 | if( j< 1 || j>nc) |
---|
1164 | { |
---|
1165 | ERROR("j doesn't describe a column of the input matrix"); |
---|
1166 | } |
---|
1167 | |
---|
1168 | // ideal of the entries of the j-th column |
---|
1169 | ideal I = A[1..nr,j]; |
---|
1170 | matrix C = coeffs(I,var(1)); |
---|
1171 | // maximal degree occurring in the j-th column |
---|
1172 | d = nrows(C)-1; |
---|
1173 | return(d); |
---|
1174 | } |
---|
1175 | |
---|
1176 | //_________________________________________________________________ |
---|
1177 | static proc coefficientAt(matrix A, int j, int e); |
---|
1178 | "USAGE: coefficientAt(A,j,e); A matrix, j,e intger |
---|
1179 | ASSUME: Base ring is polynomial ring in one variable, |
---|
1180 | j describes a column of A, e non-negative integer |
---|
1181 | smaller or equal to the column degree of |
---|
1182 | the j-th column. |
---|
1183 | RETURN: A list Z of the size of the j-th column: |
---|
1184 | - if A[i,j] has a monomial with degree e, the |
---|
1185 | corresponding coefficient is stored in Z[i] |
---|
1186 | - else Z[i] = 0 |
---|
1187 | " |
---|
1188 | { |
---|
1189 | if(nvars(basering) != 1) |
---|
1190 | { |
---|
1191 | ERROR("The base ring must be a ring in one variables."); |
---|
1192 | } |
---|
1193 | int d,i; |
---|
1194 | int nc = ncols(A); |
---|
1195 | int nr = nrows(A); |
---|
1196 | if (j<1 || j>nc) |
---|
1197 | { |
---|
1198 | ERROR("j doesn't describe a column of the input matrix"); |
---|
1199 | } |
---|
1200 | ideal I = A[1..nr,j]; |
---|
1201 | matrix C = coeffs(I,var(1)); |
---|
1202 | d = nrows(C)-1; |
---|
1203 | list Z; |
---|
1204 | if (e>d || e<0) |
---|
1205 | { |
---|
1206 | ERROR("e negative or larger than the maximum degree"); |
---|
1207 | } |
---|
1208 | for (i=1; i<=ncols(C); i++) |
---|
1209 | { |
---|
1210 | Z[i]=C[e+1,i]; |
---|
1211 | } |
---|
1212 | return(Z); |
---|
1213 | } |
---|
1214 | |
---|
1215 | //__________________________________________________________________ |
---|
1216 | static proc computeInvariants(matrix T,poly denT) |
---|
1217 | "USAGE: computeInvariants(T,denT); T matrix, denT polynomial |
---|
1218 | ASSUME: Base ring has one variable,say x, T/denT is square matrix |
---|
1219 | with entries in k(x) |
---|
1220 | RETURN: The k[x]-invariants of T/denT stored in a list. |
---|
1221 | " |
---|
1222 | { |
---|
1223 | list Dcoef; |
---|
1224 | int i,j; |
---|
1225 | int counter_inv; |
---|
1226 | list S; |
---|
1227 | for (j=1; j<=ncols(T); j++) |
---|
1228 | { |
---|
1229 | S=list(); |
---|
1230 | for(i=1; i<=nrows(T); i++) |
---|
1231 | { |
---|
1232 | if (T[i,j] != 0) |
---|
1233 | { |
---|
1234 | counter_inv++; |
---|
1235 | S[counter_inv]= deg(T[i,j])-deg(denT); |
---|
1236 | } |
---|
1237 | } |
---|
1238 | Dcoef[j] = Max(S); |
---|
1239 | counter_inv = 0; |
---|
1240 | } |
---|
1241 | return(Dcoef); |
---|
1242 | } |
---|
1243 | |
---|
1244 | //__________________________________________________________________ |
---|
1245 | static proc contentMatrix(matrix A) |
---|
1246 | "USAGE: contentMatrix(A); A matrix |
---|
1247 | NOTE: Base ring is allowed to have parameters |
---|
1248 | RETURN: The content of all entries of the matrix A |
---|
1249 | " |
---|
1250 | { |
---|
1251 | def Rbase = basering; |
---|
1252 | int nv = nvars(Rbase); |
---|
1253 | list rl = ringlist(Rbase); |
---|
1254 | rl[2][nv + 1] = "helpv"; |
---|
1255 | rl[3] = list(list("dp",1:(nv+1)),list("C",0)); |
---|
1256 | def Rhelp = ring(rl); |
---|
1257 | |
---|
1258 | setring Rhelp; |
---|
1259 | matrix A = imap(Rbase,A); |
---|
1260 | poly contCol; |
---|
1261 | int i; |
---|
1262 | for (i=1; i<=ncols(A); i++) |
---|
1263 | { |
---|
1264 | contCol = contCol+ content(A[i])*(helpv^i); |
---|
1265 | } |
---|
1266 | number contmat = content(contCol); |
---|
1267 | |
---|
1268 | setring Rbase; |
---|
1269 | number contmat = imap(Rhelp,contmat); |
---|
1270 | return(contmat); |
---|
1271 | } |
---|
1272 | |
---|
1273 | //__________________________________________________________________ |
---|
1274 | static proc idealGcd(ideal A) |
---|
1275 | "USAGE: idealGcd(A); A ideal |
---|
1276 | ASSUME: Ideal in the base ring |
---|
1277 | RETURN: The greatest common divisor of the given |
---|
1278 | generators of I |
---|
1279 | { |
---|
1280 | poly med; |
---|
1281 | poly intermed; |
---|
1282 | int i; |
---|
1283 | med = A[1]; |
---|
1284 | for (i=1; i<size(A); i++){ |
---|
1285 | intermed = gcd(med,A[i+1]); |
---|
1286 | med = intermed; |
---|
1287 | } |
---|
1288 | return(med); |
---|
1289 | } |
---|
1290 | |
---|
1291 | //__________________________________________________________________ |
---|
1292 | static proc ideal2list(ideal I) |
---|
1293 | { |
---|
1294 | list D = list(I[1..size(I)]); |
---|
1295 | return(D); |
---|
1296 | } |
---|
1297 | |
---|
1298 | //__________________________________________________________________ |
---|
1299 | static proc zerosM(int n) |
---|
1300 | { |
---|
1301 | matrix A = unitmat(n); |
---|
1302 | int i; |
---|
1303 | for (i=1;i<=n;i++) |
---|
1304 | { |
---|
1305 | A[i,i]=0; |
---|
1306 | } |
---|
1307 | return(A); |
---|
1308 | } |
---|
1309 | |
---|
1310 | //__________________________________________________________________ |
---|
1311 | static proc multiplyMatrixList(matrix T,ideal L) |
---|
1312 | { |
---|
1313 | list final; |
---|
1314 | int i,j; |
---|
1315 | for (j=1; j<=ncols(T);j++) |
---|
1316 | { |
---|
1317 | final[j]=0; |
---|
1318 | for (i=1; i<=nrows(T);i++) |
---|
1319 | { |
---|
1320 | final[j] = final[j] + L[i]*T[i,j]; |
---|
1321 | } |
---|
1322 | } |
---|
1323 | return(final); |
---|
1324 | } |
---|
1325 | |
---|
1326 | //__________________________________________________________________ |
---|
1327 | static proc isIrreducible(poly f) |
---|
1328 | "USAGE: isIrreducible(f); f poly |
---|
1329 | RETURN: 1 iff the given polynomial f is irreducible; 0 otherwise. |
---|
1330 | THEORY: This test is performed by computing the absolute |
---|
1331 | factorization of f. |
---|
1332 | KEYWORDS: irreducible. |
---|
1333 | " |
---|
1334 | { |
---|
1335 | def r = basering; |
---|
1336 | def s = absFactorize(f); |
---|
1337 | setring s; |
---|
1338 | list L = absolute_factors; |
---|
1339 | int result = 0; |
---|
1340 | if (L[4] == 1){result = 1;} |
---|
1341 | setring r; |
---|
1342 | kill s; |
---|
1343 | return (result); |
---|
1344 | } |
---|
1345 | |
---|
1346 | /////////////////////////////////////////////////////////////////// |
---|
1347 | /////////// STATIC PROCEDURES FOR FRACTIONAL IDEALS /////////////// |
---|
1348 | |
---|
1349 | //_________________________________________________________________ |
---|
1350 | static proc sumFracIdeals(list A,list B) |
---|
1351 | "USAGE: sumFracIdeals(A,B); A list, B list |
---|
1352 | ASSUME: List A (resp. B) represent fractional ideals, |
---|
1353 | A[1] is an integral ideal and A[2] a common |
---|
1354 | denominator |
---|
1355 | RETURN: The sum of the two fractional ideals, again |
---|
1356 | stored in a list. |
---|
1357 | " |
---|
1358 | { |
---|
1359 | ideal K = A[1]*B[2]+A[2]*B[1]; |
---|
1360 | poly d = A[2]*B[2]; |
---|
1361 | K = simplify(K,8); |
---|
1362 | poly gcdK = idealGcd(K); |
---|
1363 | poly commonfac = gcd(gcdK,d); |
---|
1364 | list final = K/commonfac,d/commonfac; |
---|
1365 | return(final); |
---|
1366 | } |
---|
1367 | |
---|
1368 | //__________________________________________________________________ |
---|
1369 | |
---|
1370 | static proc multiplyFracIdeals(list A,list B) |
---|
1371 | "USAGE: multiplyFracIdeals(A,B); A list, B list |
---|
1372 | ASSUME: List A (resp. B) represents a fractional ideals, |
---|
1373 | A[1] is an integral ideal and A[2] a common |
---|
1374 | denominator |
---|
1375 | RETURN: The product of the two fractional ideals, |
---|
1376 | again stored in a list |
---|
1377 | NOTE: We cannot multiply the integral ideals directly, since the |
---|
1378 | integral ideals may contain 1 and are simplified |
---|
1379 | during the multiplication |
---|
1380 | Consequently, we multiply the ideals by multiplying each |
---|
1381 | element |
---|
1382 | " |
---|
1383 | { |
---|
1384 | list A1 = ideal2list(A[1]); // transform ideals to lists |
---|
1385 | list B1 = ideal2list(B[1]); |
---|
1386 | list multiplyAB = multiplyList(A1,B1); |
---|
1387 | ideal multAB = ideal(multiplyAB[1..size(multiplyAB)]); |
---|
1388 | multAB = simplify(multAB,8); |
---|
1389 | poly d = A[2]*B[2]; |
---|
1390 | return (list(multAB,d)); |
---|
1391 | } |
---|
1392 | |
---|
1393 | //_________________________________________________________________ |
---|
1394 | static proc powerFracIdeal(list A, int d) |
---|
1395 | "USAGE: powerFracIdeal(A,d); A list, d integer |
---|
1396 | ASSUME: List A represents a fractional ideal: |
---|
1397 | A[1] is an integral ideal and A[2] a common |
---|
1398 | denominator |
---|
1399 | RETURN: The d-th power of the fractional ideal, |
---|
1400 | again stored in a list. |
---|
1401 | " |
---|
1402 | { |
---|
1403 | if (d<0) |
---|
1404 | { |
---|
1405 | ERROR("only non-negative integers allowed") |
---|
1406 | } |
---|
1407 | int i; |
---|
1408 | list final = A; |
---|
1409 | for (i=1;i<d;i++) |
---|
1410 | { |
---|
1411 | final = multiplyFracIdeals(final,A); |
---|
1412 | } |
---|
1413 | return(final); |
---|
1414 | } |
---|
1415 | |
---|
1416 | //__________________________________________________________________ |
---|
1417 | static proc isEqualFracIdeal(list A,list B) |
---|
1418 | "USAGE: isEqualFracIdeal(A,B); A list, B list |
---|
1419 | ASSUME: List A (resp. B) represents a fractional ideal: |
---|
1420 | A[1] is an integral ideal and A[2] a common |
---|
1421 | denominator |
---|
1422 | RETURN: 1 if the fractional ideals are equal, 0 else |
---|
1423 | " |
---|
1424 | { |
---|
1425 | int d; |
---|
1426 | if ( isEqualId(B[2]*A[1],A[2]*B[1]) ) |
---|
1427 | { |
---|
1428 | d = 1; |
---|
1429 | } |
---|
1430 | return(d); |
---|
1431 | } |
---|
1432 | |
---|
1433 | //__________________________________________________________________ |
---|
1434 | static proc isEqualId(ideal I,ideal J) |
---|
1435 | { |
---|
1436 | int d; |
---|
1437 | if (isIncluded(I,J) && isIncluded(J,I)) |
---|
1438 | { |
---|
1439 | d = 1; |
---|
1440 | } |
---|
1441 | return(d); |
---|
1442 | } |
---|
1443 | |
---|
1444 | //__________________________________________________________________ |
---|
1445 | static proc multiplyList(list A,list B) |
---|
1446 | "USAGE: multiplyList(A,B); A list, B list |
---|
1447 | ASSUME: list A and list B contain polynomials of the |
---|
1448 | base ring |
---|
1449 | RETURN: list obtained by multiplying every element of |
---|
1450 | list A with every element of list B |
---|
1451 | " |
---|
1452 | { |
---|
1453 | int nA = size(A); |
---|
1454 | int nB = size(B); |
---|
1455 | list multiL; |
---|
1456 | int i,j; |
---|
1457 | for(i=1;i<=nA;i++) |
---|
1458 | { |
---|
1459 | for(j=1;j<=nB;j++) |
---|
1460 | { |
---|
1461 | multiL[(i-1)*nB+j]=A[i]*B[j]; |
---|
1462 | } |
---|
1463 | } |
---|
1464 | return(multiL); |
---|
1465 | } |
---|
1466 | |
---|
1467 | //////////////////////////////////////////////////////////////////// |
---|
1468 | /////// STATIC PROCEDURES FOR TRANSFORMATION OF THE INPUT ////////// |
---|
1469 | |
---|
1470 | //__________________________________________________________________ |
---|
1471 | static proc Free2IdealRepresentation(poly f,list FinBasis, |
---|
1472 | list InfBasis,list Dformal) |
---|
1473 | "USAGE: Free2IdealRepresentation(f,FinBasis,InfBasis,Dformal); f |
---|
1474 | polynomial, FinBasis list, InfBasis list, Dformal list |
---|
1475 | ASSUME: The base ring must have two variables. |
---|
1476 | - f is an absolutely irreducible polynomial defining the |
---|
1477 | function field F/k |
---|
1478 | - list FinBasis obtains an integral basis of IntCL[k[x],F) |
---|
1479 | - list InfBasis obtains an integral basis of IntCl(k[1/x],F) |
---|
1480 | - list Dformal of size 2 representing a divisor D |
---|
1481 | Dformal[1] is a list containing the finite places of D |
---|
1482 | Dformal[2] is a list containing the infinite places of D |
---|
1483 | RETURN: A list DIdeals of size 2, the ideal representation |
---|
1484 | of D. DIdeals[1] is the finite (fractional) ideal of D, |
---|
1485 | DIdeals[2] is the infinite (fractional) ideal of D. |
---|
1486 | " |
---|
1487 | { |
---|
1488 | def Rbase = basering; |
---|
1489 | |
---|
1490 | int i; |
---|
1491 | list DFinite = Dformal[1]; |
---|
1492 | list DInfinite = Dformal[2]; |
---|
1493 | list DFinitePos = ideal(1),1; |
---|
1494 | list DFiniteNeg = ideal(1),1; |
---|
1495 | list DInfinitePos = ideal(1),1; |
---|
1496 | list DInfiniteNeg = ideal(1),1; |
---|
1497 | ideal Itmp; |
---|
1498 | poly dtmp; |
---|
1499 | list Ltmp; |
---|
1500 | for (i=1;i<=size(DFinite);i++) |
---|
1501 | { |
---|
1502 | if (DFinite[i][2]>= 0) |
---|
1503 | { |
---|
1504 | if (size(DFinite[i][1][1]) == 2) |
---|
1505 | { |
---|
1506 | Itmp = DFinite[i][1][1][1]^DFinite[i][2], DFinite[i][1][1][2]^DFinite[i][2]; |
---|
1507 | dtmp = DFinite[i][1][2]^DFinite[i][2]; |
---|
1508 | Ltmp = Itmp,dtmp; |
---|
1509 | } |
---|
1510 | else |
---|
1511 | { |
---|
1512 | Ltmp = powerFracIdeal(DFinite[i][1],DFinite[i][2]); |
---|
1513 | } |
---|
1514 | DFinitePos = multiplyFracIdeals(Ltmp,DFinitePos); |
---|
1515 | } |
---|
1516 | else |
---|
1517 | { |
---|
1518 | if (size(DFinite[i][1][1]) == 2) |
---|
1519 | { |
---|
1520 | Itmp = DFinite[i][1][1][1]^(-DFinite[i][2]), DFinite[i][1][1][2]^(-DFinite[i][2]); |
---|
1521 | dtmp = DFinite[i][1][2]^(-DFinite[i][2]); |
---|
1522 | Ltmp = Itmp,dtmp; |
---|
1523 | } |
---|
1524 | else |
---|
1525 | { |
---|
1526 | Ltmp = powerFracIdeal(DFinite[i][1],-DFinite[i][2]); |
---|
1527 | } |
---|
1528 | DFiniteNeg = multiplyFracIdeals(Ltmp,DFiniteNeg); |
---|
1529 | } |
---|
1530 | } |
---|
1531 | for (i=1;i<=size(DInfinite);i++) |
---|
1532 | { |
---|
1533 | if (DInfinite[i][2]>= 0) |
---|
1534 | { |
---|
1535 | if (size(DInfinite[i][1][1]) == 2) |
---|
1536 | { |
---|
1537 | Itmp = DInfinite[i][1][1][1]^DInfinite[i][2], DInfinite[i][1][1][2]^DInfinite[i][2]; |
---|
1538 | dtmp = DInfinite[i][1][2]^DInfinite[i][2]; |
---|
1539 | Ltmp = Itmp,dtmp; |
---|
1540 | } |
---|
1541 | else |
---|
1542 | { |
---|
1543 | Ltmp = powerFracIdeal(DInfinite[i][1],DInfinite[i][2]); |
---|
1544 | } |
---|
1545 | DInfinitePos = multiplyFracIdeals(Ltmp,DInfinitePos); |
---|
1546 | } |
---|
1547 | else |
---|
1548 | { |
---|
1549 | if (size(DInfinite[i][1][1]) == 2) |
---|
1550 | { |
---|
1551 | Itmp = DInfinite[i][1][1][1]^(-DInfinite[i][2]), DInfinite[i][1][1][2]^(-DInfinite[i][2]); |
---|
1552 | dtmp = DInfinite[i][1][2]^(-DInfinite[i][2]); |
---|
1553 | Ltmp = Itmp,dtmp; |
---|
1554 | } |
---|
1555 | else |
---|
1556 | { |
---|
1557 | Ltmp = powerFracIdeal(DInfinite[i][1],-DInfinite[i][2]); |
---|
1558 | } |
---|
1559 | DInfiniteNeg = multiplyFracIdeals(Ltmp,DInfinitePos); |
---|
1560 | } |
---|
1561 | } |
---|
1562 | |
---|
1563 | DFiniteNeg = freeGenerators(f,FinBasis,DFiniteNeg,1); |
---|
1564 | DFiniteNeg = inverseIdeal(f,FinBasis,DFiniteNeg,1); |
---|
1565 | list Fin = multiplyFracIdeals(DFinitePos,DFiniteNeg); |
---|
1566 | |
---|
1567 | DInfiniteNeg = freeGenerators(f,InfBasis,DInfiniteNeg,2); |
---|
1568 | DInfiniteNeg = inverseIdeal(f,InfBasis,DInfiniteNeg,2); |
---|
1569 | list Inf = multiplyFracIdeals(DInfinitePos,DInfiniteNeg); |
---|
1570 | |
---|
1571 | return(list(Fin,Inf)); |
---|
1572 | } |
---|
1573 | |
---|
1574 | //__________________________________________________________________ |
---|
1575 | |
---|
1576 | static proc divisorTrans(poly f,list D,list FinBasis,list Infin, |
---|
1577 | string s) |
---|
1578 | "USAGE: divisorTrans(f,D,FinBasis,Infin,s); f polynomial, D list |
---|
1579 | FinBasis list, Infin list, s string |
---|
1580 | ASSUME: - f is an absolutely irreducible homogeneous polynomial in |
---|
1581 | three variables, say x,y,z, describing a projective curve |
---|
1582 | - string s is either "ideals" or "free" |
---|
1583 | "ideals": list D contains two ideals, say I and J |
---|
1584 | representing a divisor (I) - (J) |
---|
1585 | "free": list D contains lists, each list a place |
---|
1586 | and a exponent |
---|
1587 | - list FinBasis contains an integral basis of IntCl(k[x],F) |
---|
1588 | - list Infin contains an integral basis of IntCl(k[1/x],F), |
---|
1589 | an integer k and the infinite places of F if their |
---|
1590 | computation was possible |
---|
1591 | RETUN: The corresponding divisor on F in ideal representation |
---|
1592 | (needed for the procedure RiemannRochHess). |
---|
1593 | A list of size 2, say divisortrans. |
---|
1594 | If (s = "ideals"), then the transformed divisor is given in |
---|
1595 | ideal representation. |
---|
1596 | If (s = "free"), then the transformed divisor is given in |
---|
1597 | free representation. |
---|
1598 | THEORY: We compute first the places corresponding to affine |
---|
1599 | points on the curve, i.e. the places in the chart (z=1) |
---|
1600 | For the points at infinity (z=0) we allow only |
---|
1601 | non-singular points. |
---|
1602 | " |
---|
1603 | { |
---|
1604 | |
---|
1605 | def Rbase = basering; |
---|
1606 | int i,j; |
---|
1607 | poly f_aff = subst(f,var(3),1); |
---|
1608 | list rl = ringlist(Rbase); |
---|
1609 | rl[2] = list(var(1),var(2)); |
---|
1610 | rl[3] = list(list("dp",1:2),list("C",0)); |
---|
1611 | def R_aff = ring(rl); // corresponding affine ring |
---|
1612 | |
---|
1613 | // list of all infinite places of F |
---|
1614 | list InfPlaces = Infin[3]; |
---|
1615 | |
---|
1616 | InfPlaces = transformPlacesAtInfinity(f,Infin[2],InfPlaces); |
---|
1617 | //list containing the infinite places of the function field |
---|
1618 | // and the corresponding places in the chart z = 0 of |
---|
1619 | // the projective curve if the computation is possible |
---|
1620 | |
---|
1621 | if (s == "ideals") |
---|
1622 | { |
---|
1623 | ideal I = D[1]; |
---|
1624 | ideal J = D[2]; |
---|
1625 | ideal I1,J1; |
---|
1626 | |
---|
1627 | // ideals corresponding to affine points on the curve (z=1) |
---|
1628 | I1 = sat(I,var(3))[1]; |
---|
1629 | J1 = sat(J,var(3))[1]; |
---|
1630 | ideal Ifin = std(subst(I1,var(3),1)); |
---|
1631 | ideal Jfin = std(subst(J1,var(3),1)); |
---|
1632 | |
---|
1633 | // ideals corresponding to points at infinity (z=0) |
---|
1634 | list InfIdeals; |
---|
1635 | list InfPos = ideal(1),1; |
---|
1636 | list InfNeg = ideal(1),1; |
---|
1637 | list InfTmp; |
---|
1638 | int d; |
---|
1639 | ideal Iinf, Jinf; |
---|
1640 | Iinf = sat(I,I1)[1]; |
---|
1641 | Jinf = sat(J,J1)[1]; |
---|
1642 | |
---|
1643 | if (size(InfPlaces) == 0) |
---|
1644 | { |
---|
1645 | "WARNING: Infinite places have not been computed!" |
---|
1646 | } |
---|
1647 | for (i=1;i<=size(InfPlaces);i++) |
---|
1648 | { |
---|
1649 | if(isIncluded(Iinf,InfPlaces[i][1])) |
---|
1650 | { |
---|
1651 | d = sat(Iinf,InfPlaces[i][1])[2]; |
---|
1652 | InfTmp = powerFracIdeal(InfPlaces[i][2],d); |
---|
1653 | InfPos = multiplyFracIdeals(InfPos,InfTmp); |
---|
1654 | |
---|
1655 | } |
---|
1656 | if (isIncluded(Jinf,InfPlaces[i][1])) |
---|
1657 | { |
---|
1658 | d = sat(Jinf,InfPlaces[i][1])[2]; |
---|
1659 | InfTmp = powerFracIdeal(InfPlaces[i][2],d); |
---|
1660 | InfNeg = multiplyFracIdeals(InfNeg,InfTmp); |
---|
1661 | } |
---|
1662 | } |
---|
1663 | |
---|
1664 | list FinPos = Ifin,1; |
---|
1665 | list FinNeg = Jfin,1; |
---|
1666 | FinNeg = freeGenerators(f_aff,FinBasis,FinNeg,1); |
---|
1667 | FinNeg = inverseIdeal(f_aff,FinBasis,FinNeg,1); |
---|
1668 | list Fin = multiplyFracIdeals(FinPos,FinNeg); |
---|
1669 | InfNeg = freeGenerators(f_aff,Infin[1],InfNeg,2); |
---|
1670 | InfNeg = inverseIdeal(f_aff,Infin[1],InfNeg,2); |
---|
1671 | list Inf = multiplyFracIdeals(InfPos,InfNeg); |
---|
1672 | |
---|
1673 | return(list(Fin,Inf)); |
---|
1674 | } |
---|
1675 | |
---|
1676 | if ( s == "free") |
---|
1677 | { |
---|
1678 | list PlacesFin; |
---|
1679 | list PlacesInf; |
---|
1680 | ideal Itemp; |
---|
1681 | int counter_fin,counter_inf; |
---|
1682 | for (i=1;i<=size(D);i++) |
---|
1683 | { |
---|
1684 | if (reduce(var(3),std(D[i][1])) != 0) |
---|
1685 | { |
---|
1686 | counter_fin++; |
---|
1687 | Itemp = subst(D[i][1],var(3),1); |
---|
1688 | PlacesFin[counter_fin] = list(list(Itemp,1),D[i][2]); |
---|
1689 | D = delete(D,i); |
---|
1690 | i=i-1; |
---|
1691 | } |
---|
1692 | } |
---|
1693 | if (size(InfPlaces) == 0) |
---|
1694 | { |
---|
1695 | "WARNING: Infinite places have not been computed!" |
---|
1696 | PlacesInf[1] = list(list(ideal(1),1),1); |
---|
1697 | } |
---|
1698 | else |
---|
1699 | { |
---|
1700 | for (i=1;i<=size(D);i++) |
---|
1701 | { |
---|
1702 | for(j=1;j<=size(InfPlaces);j++) |
---|
1703 | { |
---|
1704 | if (isEqualId(D[i][1],InfPlaces[j][1])) |
---|
1705 | { |
---|
1706 | counter_inf++; |
---|
1707 | PlacesInf[counter_inf] = list(InfPlaces[i][2],D[i][2]); |
---|
1708 | } |
---|
1709 | } |
---|
1710 | } |
---|
1711 | } |
---|
1712 | return(list(PlacesFin,PlacesInf)); |
---|
1713 | } |
---|
1714 | } |
---|
1715 | |
---|
1716 | //___________________________________________________________________ |
---|
1717 | static proc infPlaces (poly f) |
---|
1718 | // from brnoeth.lib |
---|
1719 | { |
---|
1720 | intvec iv; |
---|
1721 | def base_r=basering; |
---|
1722 | ring r_auxz=char(basering),(x,y,z),lp; |
---|
1723 | poly F=imap(base_r,f); |
---|
1724 | poly f_inf=subst(F,z,0); |
---|
1725 | setring base_r; |
---|
1726 | poly f_inf=imap(r_auxz,f_inf); |
---|
1727 | ideal I=factorize(f_inf,1); // points at infinity as homogeneous |
---|
1728 | // polynomials |
---|
1729 | int s=size(I); |
---|
1730 | int i; |
---|
1731 | list IP_S=list(); // for singular points at infinity |
---|
1732 | list IP_NS=list(); // for non-singular points at infinity |
---|
1733 | int counter_S; |
---|
1734 | int counter_NS; |
---|
1735 | poly aux; |
---|
1736 | |
---|
1737 | for (i=1;i<=s;i=i+1) |
---|
1738 | { |
---|
1739 | aux=subst(I[i],y,1); |
---|
1740 | if (aux==1) |
---|
1741 | { |
---|
1742 | // the point is (1:0:0) |
---|
1743 | setring r_auxz; |
---|
1744 | poly f_yz=subst(F,x,1); |
---|
1745 | if ( subst(subst(diff(f_yz,y),y,0),z,0)==0 && |
---|
1746 | subst(subst(diff(f_yz,z),y,0),z,0)==0 ) |
---|
1747 | { |
---|
1748 | // the point is singular |
---|
1749 | counter_S=counter_S+1; |
---|
1750 | kill f_yz; |
---|
1751 | setring base_r; |
---|
1752 | IP_S[counter_S]=ideal(I[i],var(3)); |
---|
1753 | } |
---|
1754 | else |
---|
1755 | { |
---|
1756 | // the point is non-singular |
---|
1757 | counter_NS=counter_NS+1; |
---|
1758 | kill f_yz; |
---|
1759 | setring base_r; |
---|
1760 | IP_NS[counter_NS]=ideal(I[i],var(3)); |
---|
1761 | } |
---|
1762 | } |
---|
1763 | else |
---|
1764 | { |
---|
1765 | // the point is (a:1:0) | a is root of aux |
---|
1766 | if (deg(aux)==1) |
---|
1767 | { |
---|
1768 | // the point is rational and no field extension is needed |
---|
1769 | setring r_auxz; |
---|
1770 | poly f_xz=subst(F,y,1); |
---|
1771 | poly aux=imap(base_r,aux); |
---|
1772 | number A=-number(subst(aux,x,0)); |
---|
1773 | map phi=r_auxz,x+A,0,z; |
---|
1774 | poly f_origin=phi(f_xz); |
---|
1775 | |
---|
1776 | if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 && |
---|
1777 | subst(subst(diff(f_origin,z),x,0),z,0)==0 ) |
---|
1778 | { |
---|
1779 | // the point is singular |
---|
1780 | counter_S=counter_S+1; |
---|
1781 | kill f_xz,aux,A,phi,f_origin; |
---|
1782 | setring base_r; |
---|
1783 | |
---|
1784 | IP_S[counter_S]=ideal(I[i],var(3)); |
---|
1785 | } |
---|
1786 | else |
---|
1787 | { |
---|
1788 | // the point is non-singular |
---|
1789 | counter_NS=counter_NS+1; |
---|
1790 | kill f_xz,aux,A,phi,f_origin; |
---|
1791 | setring base_r; |
---|
1792 | IP_NS[counter_NS]=ideal(I[i],var(3)); |
---|
1793 | } |
---|
1794 | } |
---|
1795 | else |
---|
1796 | { |
---|
1797 | // the point is non-rational and a field extension with |
---|
1798 | // minpoly=aux is needed |
---|
1799 | |
---|
1800 | ring r_ext=(char(basering),@a),(x,y,z),lp; |
---|
1801 | poly aux=imap(base_r,aux); |
---|
1802 | minpoly=number(subst(aux,x,@a)); |
---|
1803 | poly F=imap(r_auxz,F); |
---|
1804 | poly f_xz=subst(F,y,1); |
---|
1805 | map phi=r_ext,x+@a,0,z; |
---|
1806 | poly f_origin=phi(f_xz); |
---|
1807 | |
---|
1808 | if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 && |
---|
1809 | subst(subst(diff(f_origin,z),x,0),z,0)==0 ) |
---|
1810 | { |
---|
1811 | |
---|
1812 | // the point is singular |
---|
1813 | counter_S=counter_S+1; |
---|
1814 | setring base_r; |
---|
1815 | |
---|
1816 | kill r_ext; |
---|
1817 | IP_S[counter_S]=ideal(I[i],var(3)); |
---|
1818 | } |
---|
1819 | else |
---|
1820 | { |
---|
1821 | // the point is non-singular |
---|
1822 | counter_NS=counter_NS+1; |
---|
1823 | setring base_r; |
---|
1824 | kill r_ext; |
---|
1825 | IP_NS[counter_NS]=ideal(I[i],var(3)); |
---|
1826 | } |
---|
1827 | } |
---|
1828 | } |
---|
1829 | } |
---|
1830 | kill r_auxz; |
---|
1831 | return(list(IP_S,IP_NS)); |
---|
1832 | } |
---|
1833 | |
---|
1834 | //__________________________________________________________________ |
---|
1835 | static proc transformPlacesAtInfinity(poly f,int k,list InfPlacesF) |
---|
1836 | "USAGE: transformPlacesAtInfinity(f,k,InfPlacesF); f polynomial, |
---|
1837 | k integer, InfPlaces list |
---|
1838 | ASSUME: The base ring must be a ring in three variables, say x,y,z. |
---|
1839 | - f is a homogeneous polynomial and |
---|
1840 | - integer k such that if we subst. in f1 = t^(kn)*f(1/t,y) |
---|
1841 | y*t^k by u, then f1 is a monic polynomial in u with |
---|
1842 | coefficients in k[t]. |
---|
1843 | - InfPlacesF is a list containing the infinite places of F, |
---|
1844 | if the computation was possible, else empty |
---|
1845 | RETURN: A list, say places. |
---|
1846 | - If the computation of the infinite places of F was |
---|
1847 | not possible, then places is empty. |
---|
1848 | - Else, places is list of size m, where m is the number of |
---|
1849 | infinite places of F. |
---|
1850 | places[i] is a list of size 2. places[i][1] an geometric |
---|
1851 | place at infinity, i.e. a point in the chart z = 0, and |
---|
1852 | places[i][2] the corresponding infinite place of F. |
---|
1853 | " |
---|
1854 | { |
---|
1855 | list places; |
---|
1856 | |
---|
1857 | // Kummer's Theorem does not apply |
---|
1858 | // no computation of places at infinity possible |
---|
1859 | if (size(InfPlacesF) == 0) |
---|
1860 | { |
---|
1861 | return(places); |
---|
1862 | } |
---|
1863 | |
---|
1864 | // compute the places of f at z = 0 |
---|
1865 | list InfGeo = infPlaces(f); |
---|
1866 | // no singular places at infinity allowed |
---|
1867 | if ( size(InfGeo[1]) != 0 ) |
---|
1868 | { |
---|
1869 | return(places); |
---|
1870 | } |
---|
1871 | |
---|
1872 | // non-singular places |
---|
1873 | list NSPlaces = InfGeo[2]; |
---|
1874 | if (size(NSPlaces) != size(InfPlacesF)) |
---|
1875 | { |
---|
1876 | ERROR("number of infinite places must be the same"); |
---|
1877 | } |
---|
1878 | |
---|
1879 | int i,j; |
---|
1880 | int d, degtmp; |
---|
1881 | ideal Itemp; |
---|
1882 | list Ltemp; |
---|
1883 | if (size(NSPlaces) == 1) |
---|
1884 | { |
---|
1885 | places[1] = list(NSPlaces[1], InfPlacesF[1]); |
---|
1886 | } |
---|
1887 | else |
---|
1888 | { |
---|
1889 | for (i = 1; i<=size(NSPlaces); i++) |
---|
1890 | { |
---|
1891 | if (reduce(var(1),std(NSPlaces[i])) == 0) |
---|
1892 | // (x,z) is a place at infinity |
---|
1893 | { |
---|
1894 | d = 1; |
---|
1895 | Itemp = NSPlaces[i]; // put it on the end of the list |
---|
1896 | NSPlaces[i] = NSPlaces[size(NSPlaces)]; |
---|
1897 | NSPlaces[size(NSPlaces)] = Itemp; |
---|
1898 | break; |
---|
1899 | } |
---|
1900 | } |
---|
1901 | if (d == 0) |
---|
1902 | // (x,z) is not a place at infinity, simply transform the places |
---|
1903 | { |
---|
1904 | for (i = 1; i<= size(NSPlaces); i++) |
---|
1905 | { |
---|
1906 | Itemp = NSPlaces[i]; |
---|
1907 | degtmp = deg(Itemp[1]); |
---|
1908 | Itemp = Itemp[1], var(1)^(k*degtmp -1); |
---|
1909 | Ltemp = Itemp,var(1)^(k*degtmp); |
---|
1910 | places[i] = list(NSPlaces[i],Ltemp); |
---|
1911 | } |
---|
1912 | } |
---|
1913 | |
---|
1914 | else |
---|
1915 | //(x,z) is a place at infinity, no directly transformation |
---|
1916 | // is possible - transform first the other places, the remaining |
---|
1917 | // places must be (x,z) |
---|
1918 | { |
---|
1919 | for (i = 1; i < size(NSPlaces); i++) |
---|
1920 | { |
---|
1921 | Itemp = NSPlaces[i]; |
---|
1922 | degtmp = deg(Itemp[1]); |
---|
1923 | Itemp = Itemp[1], var(1)^(k*degtmp -1); |
---|
1924 | Ltemp = Itemp,var(1)^(k*degtmp); |
---|
1925 | for (j = 1; j <= size(InfPlacesF); j++) |
---|
1926 | { |
---|
1927 | if (Ltemp[2] == InfPlacesF[j][2]) |
---|
1928 | { |
---|
1929 | if (isEqualId(Ltemp[1], InfPlacesF[j][1])) |
---|
1930 | { |
---|
1931 | places[i] = list(NSPlaces[i], Ltemp); |
---|
1932 | InfPlacesF = delete(InfPlacesF,j); |
---|
1933 | } |
---|
1934 | } |
---|
1935 | } |
---|
1936 | } |
---|
1937 | if (size(InfPlacesF) != 1) |
---|
1938 | { |
---|
1939 | ERROR("more than 1 place left - wrong transformation"); |
---|
1940 | } |
---|
1941 | places[size(NSPlaces)] = list(ideal(x,z),InfPlacesF[1]); |
---|
1942 | } |
---|
1943 | } |
---|
1944 | return(places); |
---|
1945 | } |
---|
1946 | |
---|