/* VERSION: 25 February 2007. ======================================================================== This file contains Magma programs for calculating the possible Weil polynomials of curves of a given genus over a given finite field with a given number of rational points. The algorithm we use is a simplified version of the algorithm appearing in Kristin Lauter: Zeta functions of curves over finite fields with many rational points, pp. 167--174 in: Coding Theory, Cryptography, and Related Areas (Johannes Buchmann, Tom Hoholdt, Henning Stichtenoth, Horacio Tapia-Recillas, eds.), Springer-Verlag, Berlin 2000. combined with various results from Everett W. Howe and Kristin Lauter: Improved upper bounds for the number of points on curves over finite fields, arXiv: math.NT/0207101. These programs were designed to be used with Magma version 2.8. If you try to use them with earlier versions of Magma, you may have trouble with... * Coercing a rational number into a p-adic ring. * Having a function return the "throwaway identifier" _. * Using the "Valuation" function. To use this file (which is presumably named CheckQGN.magma), start Magma and load the file by entering load "CheckQGN.magma"; Then, to get a list of possible Weil polynomials of curves over GF(q) of genus g with N points, type CheckQGN(q,g,N); The program will output a list of factored "real Weil polynomials" (see the Howe-Lauter paper for the definition... or see below) of isogeny classes of abelian varieties that might conceivably contain a Jacobian of a curve with N points. Note that "conceivably" here means "not forbidden by theorems that we have incorporated into the program". Here is the notation we will use: Suppose C is a curve of genus g over GF(q). => a[i] will denote the number of degree-i places on C. => N[i] will denote the number of GF(q^i)-rational points on C. Let alpha[i] (for i = 1, ..., 2*g) be the eigenvalues of the Frobenius endomorphism of the Jacobian of C, with alpha[i] and alpha[i+g] being complex conjugates of one another. => p[i] will denote the i-th powersum of the alpha's. => c[i] will denote the coefficient of x^(2*g-i) in the characteristic polynomial of Frobenius of Jac C. If f(x) is the characteristic polynomial of Frobenius, then f(x) = x^g * h(x + q/x) for a polynomial h. Note that h(x)^2 is the characteristic polynomial of (Frobenius + Verschiebung), and that the roots of h are real numbers of absolute value at most 2*sqrt(q). In the paper, we refer to h as the "real Weil polynomial" of the curve. => b[i] will denote the coefficient of x^(g-i) in h. Here are the relations among these quantities: => N[d] is the sum, over the divisors e of d, of e * a[e]. => p[d] = q^d + 1 - N[d]. => The c's and the p's satisfy the Newton relations: n*c[n] + p[n] + sum_{0 < j < n} c[j]*p[n-j] = 0. => The b's and the c's are related by the fact that f(x) = x^g * h(x + q/x). The first part of our program uses these relations to compute the b's in terms of the a's. ======================================================================== VERSION: 13 August 2002. Corrected a sign error that prevented Weil polynomials with multiple roots from being processed correctly. (This sign error had been corrected on the version of the program used to find the results in the associated paper, but it persisted on the version available on the web because of a synchronization error.) VERSION: 18 September 2002 Added code (IsEliminatedByEllipticFactor) to eliminate possible real Weil polynomials if the corresponding isogeny class has an elliptic factor. Modified logical structure of ProcessRealWeilPolynomial (that is, rearranged the order of some if-then-else's ). VERSION: 3 October 2002 Added code to use Smyth's tables of totally positive polynomials with small deficiency to calculate possible Weil polynomials when the defect is small. Required the addition of a routine to calculate the number of places, given a Weil polynomial. VERSION: 23 October 2004 Corrected an error in NumberOfPlaces(q,g,w). If the Weil polynomial w has a linear factor f, then the generator of the number field K = NumberField(f) will *not* be a root of f. VERSION: 25 February 2007 Corrected an error in BuildUpCoefficients. In the part of the procedure where we are looking at polynomials with multiple roots, there is a point where we have to check where a polynomial (called "shift") is actually a constant. We did this by checking whether Degree(shift) was 0. But the zero polynomial has degree -1, so the correct test is to see whether Degree(shift) is less than 1. Fortunately, this error does not change the output of the program in any case considered in the Annales paper. */ // First, load a program that returns a precomputed list of // polynomials of small deficiency. load "DeficientPolynomialList.magma"; // This gives us the variable MAXDEFECT, which tells us the // largest deficiency for which we have computed the totally // positive irreducible polynomials of that deficiency. // It also gives us the function DeficientPolynomialList(R,n) // which returns the totally positive irreducible polynomials // in the ring R with deficiency n, if they have been precomputed. // Now the new routines. forward BuildUpCoefficients; forward ProcessRealWeilPolynomial; forward IsEliminatedByResultant2; forward IsEliminatedByEllipticFactor; forward IsWeilPolynomial; forward AreAllRealRootsInWeilInterval; forward AreAllRealRootsInInterval; forward AreAllRealRootsInIntervalIrreducible; forward NumberOfPlaces; //====================================================================== procedure CheckQGD(q,g,defect: fudgefactor := 0.01, verbose := true) /* fudgefactor tells us how close to Sqrt(q) a certain rational approximation should be. verbose tells us whether to print out possible real Weil polynomials that the program eliminates. */ // Begin by finding a rational approximation to Sqrt(q). // Use brute force instead of continued fractions // because fudgefactor is likely to be pretty large. i := 0; goal := Sqrt(q); repeat i +:=1; sqrtq := Ceiling(goal*i)/i; until (sqrtq^2 ge q) and (sqrtq - goal le fudgefactor); m := Floor(2*Sqrt(q))+2; while (m^2 gt 4*q) do m -:= 1; end while; R:=PolynomialRing(Rationals()); DefList := [DeficientPolynomialList(R,i) : i in [1..defect]]; ShiftedList := [[Evaluate(p,x+m+1): p in dl] : dl in DefList]; IrreducibleList := [[p : p in sl | AreAllRealRootsInInterval(p,2*sqrtq)] : sl in ShiftedList]; for i in [1..defect] do IrreducibleList[i] := [p : p in IrreducibleList[i] | AreAllRealRootsInWeilInterval(p,q)]; end for; ReducibleList := IrreducibleList; for i in [2..defect] do for j in [1.. (i div 2)] do for p in ReducibleList[j], q in ReducibleList[i-j] do pq := p*q; if not pq in ReducibleList[i] then ReducibleList[i] cat:= [pq]; end if; end for; end for; end for; for i in [1..defect] do Sort(~ReducibleList[i]); end for; if defect gt 0 then for p in ReducibleList[defect] do if Degree(p) le g then h := p * (x + m)^(g - Degree(p)); f := Numerator(Evaluate(h,(x^2+q)/x)); if IsWeilPolynomial(f,h,q) then ProcessRealWeilPolynomial([], h, q, verbose); end if; end if; end for; else h := (x + m)^g; ProcessRealWeilPolynomial([], h, q, verbose); end if; end procedure; //====================================================================== procedure CheckQGN(q,g,N: digits := 10, adjustdigits := true, fudgefactor := 0.1, checkdefect := true, verbose := true) /* digits tells us to how many decimal places we should compute real roots of polynomials. adjustdigits tells us to adjust the value of digits upwards (if necessary) so that a certain error term that appears later on will be small. fudgefactor tells us how close to Sqrt(q) a certain rational approximation should be. checkdefect tells us whether we should use the function CheckQGD when N corresponds to a small defect (at most MAXDEFECT). verbose tells us whether to print out possible real Weil polynomials that the program eliminates. */ // First compute the defect, to see whether we should // use CheckQGD. m := Floor(2*Sqrt(q))+2; while (m^2 gt 4*q) do m -:= 1; end while; defect := (q + 1 + g*m) - N; if defect lt 0 then return; end if; if (defect le MAXDEFECT) and checkdefect then CheckQGD(q,g,defect: fudgefactor := fudgefactor, verbose:=verbose); return; end if; // If we don't use CheckQGD then we must use Lauter's algorithm // to list all of the possible real Weil polynomials. // Begin by finding a rational approximation to Sqrt(q). // Use brute force instead of continued fractions // because fudgefactor is likely to be pretty large. i := 0; goal := Sqrt(q); repeat i +:=1; sqrtq := Ceiling(goal*i)/i; until (sqrtq^2 ge q) and (sqrtq - goal le fudgefactor); if adjustdigits then // See the procedure BuildUpCoefficients for details about // why we want digits to meet the following conditions. for number in [1..g-1] do coef := Factorial(g-number-1); errorfactor := (1/2)*(Factorial(g)/Factorial(g-number+1)) * (4*sqrtq + 10.0^-digits)^(number-1); // We want 10^(-2*digits)*errorfactor < coef. newdigits := Ceiling(Log(errorfactor/coef)/Log(10)/2); digits := Max(digits,newdigits); end for; end if; // Now define some rings we will be using, and compute // the b[i] (as mentioned above) in terms of the a[i]. MVR := PolynomialRing(Rationals(),g); // Think of the variables as the a[i] --- that is, a[i] = MVR.i S:=PolynomialRing(MVR); R:=PolynomialRing(Rationals()); // The vector of N's: Nvec := [&+[e * MVR.e : e in Divisors(d)] : d in [1..g]]; // The vector of powersums p: pvec := [q^d + 1 - Nvec[d] : d in [1..g]]; // We calculate the c's recursively using the Newton relations. // First initialize the c's to 0. cvec := [MVR| 0 : i in [1..g]]; // Now use the Newton relations: cvec[1] := -pvec[1]; for n in [2..g] do cvec[n] := (-1/n) * (pvec[n] + &+[pvec[n-j]*cvec[j] : j in [1..(n-1)]]); end for; // Get the characteristic polynomial of Frobenius (the Weil polynomial). weiluniv := (z^(2*g) + q^g) + cvec[g]*z^g + &+[cvec[i] * (z^(2*g-i) + z^i * q^(g-i)) : i in [1..(g-1)]]; // Get the pieces out of which the real Weil polynomial is built: h0 := (z^2 + q)^g; hvec := [(z^2 + q)^(g-i) * z^i : i in [1..(g-1)]]; // Calculate the vector of b's recursively. // First initialize the b's to 0. bvec := [MVR | 0 : i in [1..g]]; // Now use the relation f(x) = x^g * h(x + q/x). weiluniv -:= h0; for i in [1..(g-1)] do bvec[i] := Coefficient(weiluniv,2*g-i); weiluniv -:= bvec[i]*hvec[i]; end for; bvec[g] := Coefficient(weiluniv,g); // OK. Now we know how to compute the b[]'s in terms of the a[]'s. // Now create all possible lists of a[]'s that start with N, // and see whether we can eliminate any of them from consideration. initialavec := [N] cat [0: i in [1..g-1]]; BuildUpCoefficients(q, sqrtq, g, 1, initialavec, bvec, digits, R, verbose); end procedure; //====================================================================== /* Given possible values for a[1], ..., a[n-1], generate a list of possible values of a[n] using Lauter's algorithm, and recurse. If we have values for a[1], ..., a[g] then process the list of counts. */ procedure BuildUpCoefficients(q, sqrtq, g, number, partialavec, bvec, digits, R, verbose) /* q is the field size. sqrtq is a rational number just a bit bigger than Sqrt(q). g is the genus. number is the number of a[i] that have been assigned values. partialavec is the vector of the a[i]'s, with 0's for the unassigned values. (It is *critical* that the unasssigned values be set to 0.) bvec tells us how to compute the b[]'s from the a[]'s. digits tells us to how many decimal places real roots should be computed. R is a polynomial ring over the rationals. */ // First compute the real Weil polynomial (as much as we know of it) // from the vector of a's. top := Min(g,number+1); evalbvec := [Evaluate(bvec[i],partialavec) : i in [1..top]]; x := R.1; partialh := x^g + &+[evalbvec[i] * x^(g-i) : i in [1..top]]; if number lt g then // We have more recursion left to do. derh := Derivative(partialh, g-number-1); derhprime := Derivative(derh); // derhprime is the smallest derivative of h that is completely // determined by the assignment of values to a[1],...,a[number]. // The polynomial derh is determined up to an additive constant; // if a[number+1] is assigned the value k, then derh will become // derh + (g-number-1)! * k. // (That coefficient is worth having handy...) coef := Factorial(g-number-1); // Our goal is to find bounds minval and maxval so that if // the polynomial derh + (g-number-1)! * k has all of its // roots in the Weil interval, then minval <= k <= maxval. // Since k is supposed to be a[number+1], we can also bound // k below by 0 and above by the Weil bound. // If derhprime has a multiple root a, then derh must be // the unique antiderivative of derhprime with a as a root. multipleroots := GCD(derhprime,Derivative(derhprime)); if Degree(multipleroots) gt 0 then // We have multiple roots. // It's likely that there will be a problem, // so start by saying there is a problem... maxval := -1; minval := 1; // ...and fix things up if there isn't one. shift := derh mod multipleroots; // This had better be a number... if Degree(shift) lt 1 then shift := Coefficient(shift,0); // ...and it had better be an integer... if Denominator(shift) eq 1 then // ... and it had better be divisible by coef. shift := Numerator(shift); if shift mod coef eq 0 then value := -shift div coef; // ... and it had better be in the right range... maxbound := (q^(number+1) + 1 + g*Floor(2*sqrtq^(number+1))) / (number+1); if (value ge 0) and (value le maxbound) then // Mirabile dictu! There is no problem! maxval := value; minval := value; end if; end if; end if; end if; else // Now we are in the case where the roots of derhprime are distinct. // We would like an upper bound on the minimum value of derh at its // maxima, and a lower bound on the maximum value of derh at its minima. // Then we could take minval to be the negative of the upper bound // (divided by coef), and maxval to be the negative of the lower bound // (divided by coef). // The strategy: We will ask Magma for the roots of derhprime, // with accuracy 10^-digits. We will interpret these as // rational number r[1], ..., r[k]. We will check to see that // they are ordered correctly (i.e. that the error bars cannot // change the ordering). We will calculate (exactly) the // values derh(r[i]). For each r = r[i] let s be the corresponding // true root of derhprime. We can compute an upper bound for // the difference between derh(r) and derh(s) by using the formula // derh(r) = derh(s) + (r-s)*derhprime(s) // + (1/2)*(r-s)^2*derhprimeprime(t) // for some t lying between r and s. Note that derhprime(s) = 0 // by definition, and that (r-s)^2 is at most 10^(-2*digits) in // absolute value. Now, the absolute value of t is at most // 2*sqrtq + 10^-digits (where sqrtq is the rational approximation // to Sqrt(q) that we are using) so we must bound the value of // derhprimeprime on numbers of this size. But derhprimeprime // is (by assumption) a polynomial all of whose roots are real and // in the interval [-2*sqrtq..2*sqrtq], so derhprimeprime(t) is // bounded by // (leading term of derhprimeprime) * (4*sqrtq + 10^-digits)^deg // where deg = degree(derhprimeprime). // derhprimeprime is the (g-number+1)-th derivative of a monic // polynomial of degree g, so... deg := number - 1; leadingterm := Factorial(g) div Factorial(g-number+1); errorbound := (1/2) * 10^(-2*digits) * leadingterm * (4*sqrtq + 10^-digits)^deg; extrema := Roots(derhprime,RealField(): Digits:=digits); // Verify that the extrema are real, simple, and at most // sqrtq in absolute value (within limits of accuracy). flag := true; // true means everything is OK. for ex in extrema do if not IsReal(ex[1]) then flag := false; break; end if; if ex[2] gt 1 then flag eq false; break; end if; if AbsoluteValue(ex[1]) gt 2*sqrtq + 10^-digits then flag := false; break; end if; end for; if flag then // We're OK. Continue. orderedextrema := Sort([-Round(10^digits*ex[1])/10^digits : ex in extrema]); orderedextrema := [-ex : ex in orderedextrema]; // So now the extrema are listed in decreasing order. // Did we compute to enough accuracy? if &or [orderedextrema[i] - orderedextrema[i+1] lt 2*10^-digits : i in [1..#orderedextrema - 1]] then print "WARNING: INSUFFICIENT PRECISION IN ROOT FINDING."; printf "Increase value of digits, currently %0.\n\n", digits; end if; oddextrema := [orderedextrema[i] : i in [1..#orderedextrema] | IsOdd(i) ]; evenextrema := [orderedextrema[i] : i in [1..#orderedextrema] | IsEven(i)]; if #evenextrema gt 0 then minval := Max([-Evaluate(derh,ex) : ex in evenextrema]) - errorbound; else minval := 0; end if; minval := Ceiling(Max(0,minval/coef)); maxval := Min([-Evaluate(derh,ex) : ex in oddextrema]) + errorbound; maxval := maxval/coef; // maxval cannot be greater than what is allowed by the // the Weil-Serre estimate on the number of points. maxbound := (q^(number+1) + 1 + g*Floor(2*sqrtq^(number+1))) / (number+1); maxval := Floor(Min(maxbound, maxval)); // (We could improve things by replacing the Weil-Serre bound here // with the "explicit formulae" bound.) else // We're not OK. Give bad values for maxval and minval. maxval := -1; minval := 1; end if; end if; // Now recurse. for nexta in [minval..maxval] do partialavec[number+1] := nexta; BuildUpCoefficients(q, sqrtq, g, number+1, partialavec, bvec, digits, R, verbose); end for; else // If we are at this point, then the procedure has been called with // values assigned to all of the a[i]. // Process the completely determined real Weil polynomial. ProcessRealWeilPolynomial(partialavec, partialh, q, verbose); end if; end procedure; //====================================================================== /* Given a real Weil polynomial, check to see whether any of the arguments from our paper can eliminate the polynomial from consideration. */ procedure ProcessRealWeilPolynomial(avec, h, q, verbose) /* If avec is nonempty, it is the list of numbers of places of various degrees. If avec is empty, then we will compute it from h when it is needed. h is the real Weil polynomial. */ x := Parent(h).1; g := Degree(h); weil := Numerator(Evaluate(h, (x^2 + q)/x)); isweil, reason := IsWeilPolynomial(weil,h,q); if isweil then factorsh := Factorization(h); if #factorsh eq 1 then if #avec eq 0 then avec := NumberOfPlaces(q,g,weil); end if; if &and [a ge 0 : a in avec] then print avec; print factorsh; print "Power of a simple abelian variety.\n"; end if; else // Make a matrix of the values of the resultants of the // irreducible factors of h. resmat := [ [AbsoluteValue(Resultant(h0[1],h1[1])) : h1 in factorsh]: h0 in factorsh]; // min will be the minimal value of the resultant. // minlist will indicate how this minimal value is attained. min := 0; minlist := []; // Loop over nonempty subsets of the factors not containing the // last factor. for i in [1..2^(#factorsh-1)-1] do // list0 will be the factors in the given subset; // list1 will be the other factors. // r will be the resultant of the product of the elements in list0 // with the product of the elements in list1. r := 1; list0 := []; list1 := []; for j in [1..#factorsh] do if IsOdd(i div 2^(j-1)) then list0 cat:= [j]; r *:= &*[Rationals()| resmat[j][k] : k in list1]; else list1 cat:= [j]; r *:= &*[Rationals()| resmat[j][k] : k in list0]; end if; end for; if (min eq 0) or (r lt min) then min := r; minlist := list0; end if; end for; if min eq 1 then if verbose then print avec, factorsh, "ELIMINATED: resultant=1 method."; printf "Splitting = %o\n\n", minlist; end if; else // Minimal resultant is greater than 1. bool := false; if min eq 2 then // Compute avec if it is empty: if #avec eq 0 then avec := NumberOfPlaces(q,g,weil); end if; bool, splitting, reasons := IsEliminatedByResultant2(q,avec,h); end if; if bool then if verbose then print avec, factorsh, "ELIMINATED: resultant=2 method."; printf "Splitting = %o\nReasons: %o\n\n", splitting, reasons; end if; else // Not eliminated by resultant 1 or 2 methods. // What about elliptic factor method? // Compute the number of points on C if we don't already // know it. if #avec gt 0 then N := avec[1]; else N := q + 1 + Coefficient(h,g-1); end if; bool, splitting := IsEliminatedByEllipticFactor(q,N,h); if bool then if verbose then print avec, factorsh, "ELIMINATED: elliptic factor method."; printf "Splitting = %o\n\n", splitting; end if; else // We've been unable to eliminate this one. // Compute avec if we don't know it. if #avec eq 0 then avec := NumberOfPlaces(q,g,weil); end if; if &and [a ge 0: a in avec] then print avec; print factorsh; printf "Smallest resultant is %o.\n", min; printf "Resultant matrix is %o\n\n", resmat; end if; end if; end if; end if; end if; else if verbose and reason eq "p-adic" then print avec, Factorization(h), "ELIMINATED: Not Weil polynomial.\n"; end if; end if; end procedure; //====================================================================== /* Determine whether a Weil polynomial is eliminated by the resultant = 2 argument. */ function IsEliminatedByResultant2(q,avec,h); /* avec is the list of numbers of places of various degrees. h is the real Weil polynomial. Return the boolean value. If is true, return the splitting and the reasons for elimination as well. */ g := Degree(h); R := Parent(h); factorsh := Factorization(h); // Loop over nonempty subsets of the factors not containing the // last factor. for i in [1..2^(#factorsh - 1) - 1] do // h0 and h1 will be the minimal polynomials of F+V on the two factors // of the abelian variety. // hh0 and hh1 will be the real Weil polynomials of the two factors. // split0 will keep track of which factors go into h0. h0 := R!1; h1 := R!1; hh0 := R!1; hh1 := R!1; split0 := []; for j in [1..#factorsh] do if IsOdd(i div 2^(j-1)) then split0 cat:= [j]; h0 *:= factorsh[j][1]; hh0 *:= factorsh[j][1]^factorsh[j][2]; else h1 *:= factorsh[j][1]; hh1 *:= factorsh[j][1]^factorsh[j][2]; end if; end for; r := AbsoluteValue(Resultant(h0,h1)); if r eq 2 then // Make sure there is a problem with both factors. // The boolean variables will be "true" if there *is* a problem. bool0 := false; bool1 := false; // ram will be the minimal ramification contribution in // the Riemann-Hurwitz formula, if C is a double cover. ram := &+[Integers()|d : d in [1..#avec] | IsOdd(d) and IsOdd(avec[d])]; // Check to see whether Riemann-Hurwitz allows a genus-g curve // to cover a curve with genus equal to degree(hh0) or degree(hh1), // given that there are at least ram ramification points. if g lt 2*Degree(hh0)-1+ram/2 then bool0 := true; reason0 := "Riemann-Hurwitz"; end if; if g lt 2*Degree(hh1)-1+ram/2 then bool1 := true; reason1 := "Riemann-Hurwitz"; end if; // Now check to see whether the number of points on the top curve // is at most twice the number of the bottom curve. if (q + 1 + Coefficient(h, g-1)) gt 2*(q + 1 + Coefficient(hh0, Degree(hh0)-1)) then bool0 := true; reason0 := "point counts"; end if; if (q + 1 + Coefficient(h, g-1)) gt 2*(q + 1 + Coefficient(hh1, Degree(hh1)-1)) then bool1 := true; reason1 := "point counts"; end if; // If there is a problem with both factors, return true. if bool0 and bool1 then reasons := reason0 cat ", " cat reason1; return true, split0, reasons; end if; end if; end for; return false, _, _; end function; //====================================================================== /* Determine whether a Weil polynomial is eliminated by the elliptic factor argument. */ function IsEliminatedByEllipticFactor(q,N,h); /* N is the number of rational points on the covering curve. h is the real Weil polynomial. Return the boolean value. If is true, return the splitting. */ g := Degree(h); R := Parent(h); factorsh := Factorization(h); // Loop over the factors. for i in [1..#factorsh] do h0 := factorsh[i][1]; e0 := factorsh[i][2]; if (e0 eq 1) and (Degree(h0) eq 1) then h1 := R!1; for j in [1..#factorsh] do if j ne i then h1 *:= factorsh[j][1]; end if; end for; // h0 and h1 are the minimal polynomials of F+V on the two factors // of the abelian variety. r := AbsoluteValue(Resultant(h0,h1)); sizeofE := q + 1 + Coefficient(h0,0); if N gt r * sizeofE then // We have a problem --- too many points on top curve. return true, i; end if; end if; end for; return false, _; end function; //====================================================================== /* Given a polynomial f, a prime power q, and a polynomial h such that f(x) = h(x + q/x) * x^degree(h), determine whether or not f is a Weil q-polynomial. If not, return the reason: either "real" if there roots of f are not on circle of radius sqrt(q), or "p-adic" if the Honda-Tate exponents don't work out. */ function IsWeilPolynomial(f,h,q) if not AreAllRealRootsInWeilInterval(h,q) then return false, "real"; end if; fq := Factorization(q)[1]; p := fq[1]; exp := fq[2]; facf := Factorization(f); // I don't know what p-adic precision Magma requires in order // to factor a polynomial over a p-adic ring. Certainly the // discriminant of the polynomial should be nonzero up to // the precision chosen, but apparently more is required. // The following seems to work... but is not guaranteed! // Fortunately, we don't need to call IsWeilPolynomial very // often, so the inefficiency of extra precision is not // such a large price to pay. precisions := [exp*Degree(fe[1]) + Valuation(Discriminant(fe[1]),p) : fe in facf]; padicprecision := 10 + Max(precisions); T := PolynomialRing(pAdicRing(p,padicprecision)); for fe in facf do localf := T!fe[1]; HondaTateExponent := 1; // 1 is the initial value... it will change. faclocalf := Factorization(localf); for fl in faclocalf do if fl[2] gt 1 then print "WARNING: BAD PRECISION in IsWeilPolynomial"; // There should never be multiple roots. end if; val := Valuation(Coefficient(fl[1],0)); val0 := Valuation(Coefficient(fl[1],0)*0); if val eq val0 then print "WARNING: POSSIBLE BAD PRECISION in IsWeilPolynomial"; end if; localHondaTateExponent := Denominator(val/exp); HondaTateExponent := LCM(HondaTateExponent, localHondaTateExponent); end for; if 0 ne (fe[2] mod HondaTateExponent) then return false, "p-adic"; end if; end for; return true, _; end function; //====================================================================== /* Given an irreducible monic polynomial h with rational coefficients and a positive rational number r, determine whether or not all of the complex roots of h are real and lie in the interval [-r, r]. We use Sturm's theorem (see Theorem 4.1.10 of Henri Cohen's _A Course in Computational Algebraic Number Theory_). Our implementation is naive. */ function AreAllRealRootsInIntervalIrreducible(h,r) if Degree(h) eq 1 then root := -Coefficient(h,0); return AbsoluteValue(root) le r; end if; remainderseq := [h, Derivative(h)]; i := 1; repeat i +:= 1; newterm := -(remainderseq[i-1] mod remainderseq[i]); remainderseq cat:= [newterm]; until Degree(newterm) eq 0; valseq1 := [Evaluate(rs,-r) : rs in remainderseq]; valseq2 := [Evaluate(rs, r) : rs in remainderseq]; signseq1 := [Sign(vs) : vs in valseq1 | vs ne 0]; if #signseq1 ne Degree(h) + 1 then return false; end if; signseq2 := [Sign(vs) : vs in valseq2 | vs ne 0]; allchange := &and[signseq1[i] ne signseq1[i+1] : i in [1..#signseq1-1]]; nochange := &and[signseq2[i] eq signseq2[i+1] : i in [1..#signseq2-1]]; return allchange and nochange; end function; //====================================================================== /* Given a polynomial h with rational coefficients and a positive rational number r, determine whether or not all of the complex roots of h are real and lie in the interval [-r, r]. To check this, simply check this for all of the monic irreducible factors of h. */ function AreAllRealRootsInInterval(h,r) fach := Factorization(h); for hh in fach do if not AreAllRealRootsInIntervalIrreducible(hh[1],r) then return false; end if; end for; return true; end function; //====================================================================== /* Given a polynomial h with rational coefficients and a positive rational number q, determine whether or not all of the complex roots of h are real and lie in the interval [-2*sqrt(q), 2*sqrt(q)]. */ function AreAllRealRootsInWeilInterval(h,q) bool, qq := IsSquare(q); if bool then // q is a square, and qq its square root. return AreAllRealRootsInInterval(h,2*qq); end if; // That was the easy case. Now for the hard case. // Compute the polynomial whose roots are the squares of the roots of h. x := Parent(h).1; h1 := &+[x^i*Coefficient(h,i) : i in [0..Degree(h)] | IsEven(i)]; h2 := &+[x^i*Coefficient(h,i) : i in [0..Degree(h)] | IsOdd(i)]; hh := h1^2 - h2^2; newh := &+[x^i*Coefficient(hh,2*i) : i in [0..Degree(hh) div 2]]; // newh is the polynomial whose roots are the squares of the roots of h. // We want to check whether newh's roots are all in [0..4*q]. newnewh := Evaluate(newh, x + 2*q); return AreAllRealRootsInInterval(newnewh, 2*q); end function; //====================================================================== /* Given a prime power q, a genus g, and a Weil polynomial w, compute the number of degree-i places, for i = 1, ..., g. */ function NumberOfPlaces(q,g,w) Nvec := [q^i + 1 : i in [1..g]]; factorsofw := Factorization(w); for fe in factorsofw do f := fe[1]; e := fe[2]; K:=NumberField(f); if Degree(f) gt 1 then pi := K.1; else pi:=Roots(f)[1][1]; end if; for i in [1..g] do Nvec[i] -:= e * Trace(pi^i); end for; end for; avec := [Numerator( (&+[Nvec[d]*MoebiusMu(n div d) : d in Divisors(n)]) / n) : n in [1..g]]; return avec; end function;