Coverage for local/lib/python2.7/site-packages/sage/rings/number_field/totallyreal_data.pyx : 75%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" Enumeration of Totally Real Fields
AUTHORS:
- Craig Citro and John Voight (2007-11-04): Type checking and other polishing. - John Voight (2007-10-09): Improvements: Smyth bound, Lagrange multipliers for b. - John Voight (2007-09-19): Various optimization tweaks. - John Voight (2007-09-01): Initial version. """
#***************************************************************************** # Copyright (C) 2007 William Stein and John Voight # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #*****************************************************************************
from __future__ import absolute_import, print_function
from libc.math cimport sqrt from cysignals.memory cimport sig_malloc, sig_free
from sage.libs.gmp.mpz cimport * from sage.rings.integer cimport Integer
# Other global variables
from libc.math cimport lrint, floor, ceil, fabs, round
#********************************************************************* # Auxiliary routine # Hermite constant, naive Newton-Raphson, and a specialized Lagrange # multiplier solver. #*********************************************************************
r""" This function returns the nth Hermite constant
The nth Hermite constant (typically denoted `\gamma_n`), is defined to be
.. MATH::
\max_L \min_{0 \neq x \in L} ||x||^2
where `L` runs over all lattices of dimension `n` and determinant `1`.
For `n \leq 8` it returns the exact value of `\gamma_n`, and for `n > 9` it returns an upper bound on `\gamma_n`.
INPUT:
- n -- integer
OUTPUT:
- (an upper bound for) the Hermite constant gamma_n
EXAMPLES::
sage: hermite_constant(1) # trivial one-dimensional lattice 1.0 sage: hermite_constant(2) # Eisenstein lattice 1.1547005383792515 sage: 2/sqrt(3.) 1.15470053837925 sage: hermite_constant(8) # E_8 2.0
.. NOTE::
The upper bounds used can be found in [CS]_ and [CE]_.
REFERENCES:
.. [CE] Henry Cohn and Noam Elkies, New upper bounds on sphere packings I, Ann. Math. 157 (2003), 689--714.
.. [CS] \J.H. Conway and N.J.A. Sloane, Sphere packings, lattices and groups, 3rd. ed., Grundlehren der Mathematischen Wissenschaften, vol. 290, Springer-Verlag, New York, 1999.
AUTHORS:
- John Voight (2007-09-03) """
# Exact values are known for gamma_n. elif n <= 36: gamma = [2.13263235569928, 2.26363016185702, 2.39334691240146, 2.52178702088414, 2.64929462619823, 2.77580405570023, 2.90147761892077, 3.02639364467182, 3.15067928476872, 3.27433066745617, 3.39744386110070, 3.52006195697466, 3.64224310140724, 3.76403701226104, 3.88547626036618, 4.00659977840648, 4.12744375027069, 4.24804458298177, 4.36843113799634, 4.48863097933934, 4.60866759008263, 4.72856660611662, 4.84834821242630, 4.96803435811402, 5.08764086822471, 5.20718687262715, 5.32668836123079, 5.44615801810606][n-9] else: # Mordell's inequality. gamma = 5.44615801810606**((n-1.)/35)
cdef double eval_seq_as_poly(int *f, int n, double x): r""" Evaluates the sequence a, thought of as a polynomial with
.. MATH::
f[n]*x^n + f[n-1]*x^(n-1) + ... + f[0]. """ cdef double s, xp
# Horner's method: With polynomials of small degree, we shouldn't # expect asymptotic methods to be any faster.
cdef double newton(int *f, int *df, int n, double x0, double eps): r""" Find the real root x of f (with derivative df) near x0 with provable precision eps, i.e. |x-z| < eps where z is the actual root. The sequence a corresponds to the polynomial f with
.. MATH::
f(x) = x^n + a[n-1]*x^(n-1) + ... + a[0].
This function (for speed reasons) has no error checking and no guarantees are made as to the convergence; a naive Newton-Raphson method is used. """ cdef double x, rdx, dx, fx
# In truly optimized code, one could tune by automatically # iterating a certain number of times based on the size of dx to # save on a few comparisons. # This savings should be almost negligible...?
# Small hack for improved performance elsewhere: if it is close to an # integer, give it full precision as an integer.
# Now ensure that either f(x-eps) or f(x+eps) has opposite sign # as f(x), which implies that |x-z| < eps. dx = eval_seq_as_poly(f,n,x)/eval_seq_as_poly(df,n-1,x) x -= dx fx = eval_seq_as_poly(f,n,x)
cdef void newton_in_intervals(int *f, int *df, int n, double *beta, double eps, double *rts): r""" Find the real roots of f in the intervals specified by beta:
(beta[0],beta[1]), (beta[1],beta[2]), ..., (beta[n-1], beta[n])
Calls newton_kernel, so same provisos apply---in particular, each interval should contain a unique simple root. Note the derivative *df is passed but is recomputed--this is just a way to save a malloc and free for each call. """ cdef int i
cpdef lagrange_degree_3(int n, int an1, int an2, int an3): r""" Private function. Solves the equations which arise in the Lagrange multiplier for degree 3: for each 1 <= r <= n-2, we solve
r*x^i + (n-1-r)*y^i + z^i = s_i (i = 1,2,3)
where the s_i are the power sums determined by the coefficients a. We output the largest value of z which occurs. We use a precomputed elimination ideal.
EXAMPLES::
sage: ls = sage.rings.number_field.totallyreal_data.lagrange_degree_3(3,0,1,2) sage: [RealField(10)(x) for x in ls] [-1.0, -1.0] sage: sage.rings.number_field.totallyreal_data.lagrange_degree_3(3,6,1,2) # random [-5.8878, -5.8878]
TESTS:
Check that :trac:`13101` is solved::
sage: sage.rings.number_field.totallyreal_data.lagrange_degree_3(4,12,19,42) [0.0, -1.0] """ cdef double zmin, zmax, val cdef double *roots_data cdef long coeffs[7] cdef int r, rsq, rcu cdef int nr, nrsq, nrcu cdef int s1, s1sq, s1cu, s1fo, s2, s2sq, s2cu, s3, s3sq
# Newton's relations.
# common subexpressions
## x^6 r*nrcu + 5*r*nrsq + 7*r*nr + 3*r + nrcu + \ 3*nrsq + 3*nr + 1
## x^5 12*r*s1 - 6*nrsq*s1 - 12*nr*s1 - 6*s1
## x^4 3*r*nrsq*s2 + 15*r*nr*s1sq - 6*r*nr*s2 + 18*r*s1sq - \ 3*r*s2 - 3*nrcu*s2 + 3*nrsq*s1sq - 6*nrsq*s2 + \ 18*nr*s1sq - 3*nr*s2 + 15*s1sq
## x^3 6*rsq*nr*s3 + 12*rsq*s1*s2 - 2*r*nrcu*s3 + \ 6*r*nrsq*s1*s2 - 6*r*nrsq*s3 - 4*r*nr*s1cu + \ 12*r*nr*s1*s2 - 4*r*nr*s3 - 12*r*s1cu + 12*r*s1*s2 + \ 12*nrsq*s1*s2 - 12*nr*s1cu + 12*nr*s1*s2 - \ 20*s1cu
## x^2 6*rsq*s1sq*s2 + 3*rsq*s2sq + 6*r*nrsq*s1*s3 - \ 3*r*nrsq*s2sq - 6*r*nr*s1sq*s2 + 12*r*nr*s1*s3 + \ 3*r*nr*s2sq + 3*r*s1fo - 18*r*s1sq*s2 + \ 3*nrcu*s2sq - 6*nrsq*s1sq*s2 + 3*nrsq*s2sq + \ 3*nr*s1fo - 18*nr*s1sq*s2 + 15*s1fo
## x^1 12*r*nr*s1sq*s3 - 6*r*nr*s1*s2sq + 12*r*s1cu*s2 - \ 6*nrsq*s1*s2sq + 12*nr*s1cu*s2 - 6*s1*s1fo
## x^0 6*rsq*nr*s1*s2*s3 + rsq*nr*s2cu + 3*rsq*s1sq*s2sq + \ r*nrcu*s3sq - 6*r*nrsq*s1*s2*s3 + r*nrsq*s2cu + \ 4*r*nr*s1cu*s3 + 3*r*nr*s1sq*s2sq - \ 3*r*s1fo*s2 - nrcu*s2cu + \ 3*nrsq*s1sq*s2sq - 3*nr*s1fo*s2 + \ s1sq*s1fo
cdef long primessq[46]
r""" Returns the largest a such that a^2 divides d and a has prime divisors < 200.
EXAMPLES::
sage: from sage.rings.number_field.totallyreal_data import int_has_small_square_divisor sage: int_has_small_square_divisor(500) 100 sage: is_prime(691) True sage: int_has_small_square_divisor(691) 1 sage: int_has_small_square_divisor(691^2) 1 """
cdef int i cdef Integer asq
cdef int eval_seq_as_poly_int(int *f, int n, int x): r""" Evaluates the sequence a, thought of as a polynomial with
.. MATH::
f[n]*x^n + f[n-1]*x^(n-1) + ... + f[0]. """ cdef int s, xp
cdef double eps_abs, phi, sqrt2
cdef int easy_is_irreducible(int *a, int n): r""" Very often, polynomials have roots in {+/-1, +/-2, +/-phi, sqrt2}, so we rule these out quickly. Returns 0 if reducible, 1 if inconclusive. """ cdef int s, t, st, sgn, i
# Check if a has a root in {1,-1,2,-2}.
# Check if f has factors x^2-x-1, x^2+x-1, x^2-2, respectively. # Note we only call the ZZx constructor if we're almost certain to reject.
""" Used solely for testing easy_is_irreducible.
EXAMPLES::
sage: sage.rings.number_field.totallyreal_data.easy_is_irreducible_py(pari('x^2+1')) 1 sage: sage.rings.number_field.totallyreal_data.easy_is_irreducible_py(pari('x^2-1')) 0 """ cdef int a[10]
#**************************************************************************** # Main class and routine #****************************************************************************
# Global precision to find roots; this should probably depend on the # architecture in some way. Algorithm gives provably correct results # for any eps, but an optimal value of eps will be neither too large # (which gives trivial bounds on coefficients) nor too small (which # spends needless time computing higher precision on the roots). cdef double eps_global
cdef class tr_data: r""" This class encodes the data used in the enumeration of totally real fields.
We do not give a complete description here. For more information, see the attached functions; all of these are used internally by the functions in totallyreal.py, so see that file for examples and further documentation. """
r""" Initialization routine (constructor).
INPUT:
n -- integer, the degree B -- integer, the discriminant bound a -- list (default: []), the coefficient list to begin with, where a[len(a)]*x^n + ... + a[0]x^(n-len(a))
OUTPUT:
the data initialized to begin enumeration of totally real fields with degree n, discriminant bounded by B, and starting with coefficients a.
EXAMPLES::
sage: T = sage.rings.number_field.totallyreal_data.tr_data(2,100) sage: T.printa() k = 0 a = [0, -1, 1] amax = [0, 0, 1] beta = [...] gnk = [...] """
cdef int i
# Initialize constants.
# Declare the coefficients of the polynomials (and max such). # df is memory set aside for the derivative, as # used in Newton iteration above.
# beta is an array of arrays (of length n) which list the # roots of the derivatives. # gnk is the collection of (normalized) derivatives.
# Initialize variables. # No starting input, all polynomials will be found; initialize to zero. elif len(a) <= n+1: # First few coefficients have been specified. # The value of k is the largest index of the coefficients of a which is # currently unknown; e.g., if k == -1, then we can iterate # over polynomials, and if k == n-1, then we have finished iterating. if a[len(a)-1] != 1: raise ValueError("a[len(a)-1](=%s) must be 1 so polynomial is monic" % a[len(a)-1])
k = n-len(a) self.k = k a = [0]*(k+1) + a for i from 0 <= i < n+1: self.a[i] = a[i] self.amax[i] = a[i]
# Bounds come from an application of Lagrange multipliers in degrees 2,3. self.b_lower = -1./n*(self.a[n-1] + (n-1.)*sqrt((1.*self.a[n-1])**2 - 2.*(1+1./(n-1))*self.a[n-2])) self.b_upper = -1./n*(self.a[n-1] - (n-1.)*sqrt((1.*self.a[n-1])**2 - 2.*(1+1./(n-1))*self.a[n-2])) if k < n-3: bminmax = lagrange_degree_3(n,a[n-1],a[n-2],a[n-3]) if bminmax: self.b_lower = bminmax[0] self.b_upper = bminmax[1]
# Annoying, but must reverse coefficients for NumPy. gnk = [int(binomial(j,k+2))*a[j] for j in range(k+2,n+1)] gnk.reverse() import numpy rts = numpy.roots(gnk).tolist() rts.sort() self.beta[(k+1)*(n+1)+0] = self.b_lower for i from 0 <= i < n-k-2: self.beta[(k+1)*(n+1)+(i+1)] = rts[i] self.beta[(k+1)*(n+1)+(n-k-1)] = self.b_upper
# Now to really initialize gnk. gnk = [0] + [binomial(j,k+1)*a[j] for j in range (k+2,n+1)] for i from 0 <= i < n-k: self.gnk[(k+1)*n+i] = gnk[i] else: # Bad input! raise ValueError("a has length %s > n+1" % len(a))
def __dealloc__(self): r""" Destructor. """
def increment(self, verbose=False, haltk=0, phc=False): r""" This function 'increments' the totally real data to the next value which satisfies the bounds essentially given by Rolle's theorem, and returns the next polynomial as a sequence of integers.
The default or usual case just increments the constant coefficient; then inductively, if this is outside of the bounds we increment the next higher coefficient, and so on.
If there are no more coefficients to be had, returns the zero polynomial.
INPUT:
- verbose -- boolean to print verbosely computational details - haltk -- integer, the level at which to halt the inductive coefficient bounds - phc -- boolean, if PHCPACK is available, use it when k == n-5 to compute an improved Lagrange multiplier bound
OUTPUT:
The next polynomial, as a sequence of integers
EXAMPLES::
sage: T = sage.rings.number_field.totallyreal_data.tr_data(2,100) sage: T.increment() [-24, -1, 1] sage: for i in range(19): _ = T.increment() sage: T.increment() [-3, -1, 1] sage: T.increment() [-25, 0, 1] """ cdef int *f_out cdef int i
raise MemoryError("unable to allocate coefficient list")
cdef void incr(self, int *f_out, int verbose, int haltk, int phc): r""" This function 'increments' the totally real data to the next value which satisfies the bounds essentially given by Rolle's theorem, and returns the next polynomial in the sequence f_out.
The default or usual case just increments the constant coefficient; then inductively, if this is outside of the bounds we increment the next higher coefficient, and so on.
If there are no more coefficients to be had, returns the zero polynomial.
INPUT:
- f_out -- an integer sequence, to be written with the coefficients of the next polynomial - verbose -- boolean to print verbosely computational details - haltk -- integer, the level at which to halt the inductive coefficient bounds - phc -- boolean, if PHCPACK is available, use it when k == n-5 to compute an improved Lagrange multiplier bound
OUTPUT:
None. The return value is stored in the variable f_out. """
cdef int n, np1, k, i, j, nk, kz cdef int *gnkm cdef int *gnkm1 cdef double *betak cdef double bl, br, akmin, akmax, tmp_dbl cdef bint maxoutflag
# If k == -1, we have a full polynomial, so we add 1 to the constant coefficient. # Can't have constant coefficient zero! else: print(" ", end="") for i from 0 <= i < np1: print(self.a[i], end="") print(">", end="") for i from 0 <= i < np1: print(self.amax[i], end="") print("")
# Already reached maximum, so "carry the 1" to find the next value of k. # If we are working through an initialization routine, treat that. self.a[k] += 1 if self.a[k] > self.amax[k]: k += 1 while k <= n-1 and self.a[k] >= self.amax[k]: k += 1 self.a[k] += 1 self.gnk[n*k] = 0 k -= 1
# If in the previous step we finished all possible values of # the lastmost coefficient, so we must compute bounds on the next coefficient. # Recall k == n-1 implies iteration is complete. # maxoutflag flags a required abort along the way
# Recall k == -1 means all coefficients are good to go. print(k, ":", end="") for i from 0 <= i < np1: print(self.a[i], end="") print("")
# We only know the value of a[n-1], the trace. Need to apply # basic bounds from Hunter's theorem and Siegel's theorem, with # improvements due to Smyth to get bounds on a[n-2]. bl = 1./2*(1-1./n)*(1.*self.a[n-1])**2 \
# If maximum is already greater than the minimum, break! print(" ", end="") for i from 0 <= i < np1: print(self.a[i], end="") print(">", end="") for i from 0 <= i < np1: print(self.amax[i], end="") print("")
# Knowing a[n-1] and a[n-2] means we can apply bounds from # the Lagrange multiplier in degree 2, which can be solved # immediately.
# Initialize the second derivative.
print(" ", '%.2f' % self.beta[k * np1 + 1]) else: # Compute the roots of the derivative. &self.beta[(k+1)*np1], eps_global, &self.beta[k*np1+1]) print(" ", end="") for i from 0 <= i < n-k-1: print('%.2f' % self.beta[k * np1 + 1 + i], end="") print("")
# This happens reasonably infrequently, so calling # the Python routine should be sufficiently fast... # Could just take self.gnk(k+2)*n+i, but this may not be initialized... print(" gnk has multiple factor!")
# Bounds come from an application of Lagrange multipliers in degrees 2,3. # New bounds from Lagrange multiplier in degree 3. # New bounds using phc/Lagrange multiplier in degree 4. bminmax = __lagrange_bounds_phc(n, 4, [self.a[i] for i from 0 <= i <= n]) if len(bminmax) > 0: self.b_lower = bminmax[0] self.b_upper = bminmax[1] else: maxoutflag = 1 break
print(" [LM bounds:", '%.2f' % self.b_lower, '%.2f' % self.b_upper, end="") tb = sqrt((1.*self.a[n-1])**2 - 2.*self.a[n-2]) print("vs. +/-", '%.2f' % tb, ']')
# Compute next g_(n-(k+1)), k times the formal integral of g_(n-k).
# Compute upper and lower bounds which guarantee one retains # a polynomial with all real roots. akmin = -eval_seq_as_poly(gnkm, n-k, betak[nk+1]) \ # Use the fact that f(z) <= f(x)+|f'(x)|eps if |x-z| < eps # for sufficiently small eps, f(z) = 0, and f''(z) < 0. tmp_dbl = -eval_seq_as_poly(gnkm, n-k, betak[nk+1-2*i]) \
akmax = -eval_seq_as_poly(gnkm, n-k, betak[nk]) \ # Similar calculus statement here. tmp_dbl = -eval_seq_as_poly(gnkm, n-k, betak[nk-2*i]) \
# Can replace alpha by -alpha, so if all # "odd" coefficients are zero, may assume next # "odd" coefficient is positive. # Can't have constant coefficient zero!
print(" ", end="") for i from 0 <= i < np1: print(self.a[i], end="") print(">", end="") for i from 0 <= i < np1: print(self.amax[i], end="") print("")
else:
# k == n-1, so iteration is complete; return the zero polynomial (of degree n+1).
def printa(self): """ Print relevant data for self.
EXAMPLES::
sage: T = sage.rings.number_field.totallyreal_data.tr_data(3,2^10) sage: T.printa() k = 1 a = [0, 0, -1, 1] amax = [0, 0, 0, 1] beta = [...] gnk = [...]
""" |