Coverage for local/lib/python2.7/site-packages/sage/rings/number_field/bdd_height.py : 97%

r""" 

Elements of bounded height in number fields 

Sage functions to list all elements of a given number field with height less 

than a specified bound. 

AUTHORS: 

- John Doyle (2013): initial version 

- David Krumm (2013): initial version 

REFERENCES: 

.. [Doyle-Krumm] John R. Doyle and David Krumm, Computing algebraic numbers 

of bounded height, :arxiv:`1111.4963` (2013). 

""" 

#***************************************************************************** 

# Copyright (C) 2013 John Doyle and David Krumm 

# 

# Distributed under the terms of the GNU General Public License (GPL) 

# 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 print_function 

from six.moves import range 

from copy import copy 

from itertools import product 

from sage.rings.real_mpfr import RealField 

from sage.rings.number_field.unit_group import UnitGroup 

from sage.modules.free_module_element import vector 

from sage.matrix.constructor import column_matrix 

from sage.rings.rational_field import QQ 

from sage.functions.other import ceil 

from sage.geometry.polyhedron.constructor import Polyhedron 

from sage.structure.proof.all import number_field 

from sage.libs.pari.all import pari 

def bdd_norm_pr_gens_iq(K, norm_list): 

r""" 

Compute generators for all principal ideals in an imaginary quadratic field 

`K` whose norms are in ``norm_list``. 

The only keys for the output dictionary are integers n appearing in 

``norm_list``. 

The function will only be called with `K` an imaginary quadratic field. 

The function will return a dictionary for other number fields, but it may be 

incorrect. 

INPUT: 

- `K` - an imaginary quadratic number field 

- ``norm_list`` - a list of positive integers 

OUTPUT: 

- a dictionary of number field elements, keyed by norm 

EXAMPLES: 

In `QQ(i)`, there is one principal ideal of norm 4, two principal ideals of 

norm 5, but no principal ideals of norm 7:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_gens_iq 

sage: K.<g> = NumberField(x^2 + 1) 

sage: L = range(10) 

sage: bdd_pr_ideals = bdd_norm_pr_gens_iq(K, L) 

sage: bdd_pr_ideals[4] 

[2] 

sage: bdd_pr_ideals[5] 

[-g - 2, -g + 2] 

sage: bdd_pr_ideals[7] 

[] 

There are no ideals in the ring of integers with negative norm:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_gens_iq 

sage: K.<g> = NumberField(x^2 + 10) 

sage: L = range(-5,-1) 

sage: bdd_pr_ideals = bdd_norm_pr_gens_iq(K,L) 

sage: bdd_pr_ideals 

{-5: [], -4: [], -3: [], -2: []} 

Calling a key that is not in the input ``norm_list`` raises a KeyError:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_gens_iq 

sage: K.<g> = NumberField(x^2 + 20) 

sage: L = range(100) 

sage: bdd_pr_ideals = bdd_norm_pr_gens_iq(K, L) 

sage: bdd_pr_ideals[100] 

Traceback (most recent call last): 

... 

KeyError: 100 

""" 

gens = dict() 

for n in norm_list: 

gens[n] = K.elements_of_norm(n) 

return gens 

def bdd_height_iq(K, height_bound): 

r""" 

Compute all elements in the imaginary quadratic field `K` which have 

relative multiplicative height at most ``height_bound``. 

The function will only be called with `K` an imaginary quadratic field. 

If called with `K` not an imaginary quadratic, the function will likely 

yield incorrect output. 

ALGORITHM: 

This is an implementation of Algorithm 5 in [Doyle-Krumm]. 

INPUT: 

- `K` - an imaginary quadratic number field 

- ``height_bound`` - a real number 

OUTPUT: 

- an iterator of number field elements 

EXAMPLES:: 

sage: from sage.rings.number_field.bdd_height import bdd_height_iq 

sage: K.<a> = NumberField(x^2 + 191) 

sage: for t in bdd_height_iq(K,8): 

....: print(exp(2*t.global_height())) 

1.00000000000000 

1.00000000000000 

1.00000000000000 

4.00000000000000 

4.00000000000000 

4.00000000000000 

4.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

There are 175 elements of height at most 10 in `QQ(\sqrt(-3))`:: 

sage: from sage.rings.number_field.bdd_height import bdd_height_iq 

sage: K.<a> = NumberField(x^2 + 3) 

sage: len(list(bdd_height_iq(K,10))) 

175 

The only elements of multiplicative height 1 in a number field are 0 and 

the roots of unity:: 

sage: from sage.rings.number_field.bdd_height import bdd_height_iq 

sage: K.<a> = NumberField(x^2 + x + 1) 

sage: list(bdd_height_iq(K,1)) 

[0, a + 1, a, -1, -a - 1, -a, 1] 

A number field has no elements of multiplicative height less than 1:: 

sage: from sage.rings.number_field.bdd_height import bdd_height_iq 

sage: K.<a> = NumberField(x^2 + 5) 

sage: list(bdd_height_iq(K,0.9)) 

[] 

""" 

if height_bound < 1: 

return 

yield K(0) 

roots_of_unity = K.roots_of_unity() 

for zeta in roots_of_unity: 

yield zeta 

# Get a complete set of ideal class representatives 

class_group_reps = [] 

class_group_rep_norms = [] 

for c in K.class_group(): 

a = c.ideal() 

class_group_reps.append(a) 

class_group_rep_norms.append(a.norm()) 

class_number = len(class_group_reps) 

# Find principal ideals of bounded norm 

possible_norm_set = set([]) 

for n in range(class_number): 

for m in range(1, height_bound + 1): 

possible_norm_set.add(m*class_group_rep_norms[n]) 

bdd_ideals = bdd_norm_pr_gens_iq(K, possible_norm_set) 

# Distribute the principal ideals 

generator_lists = [] 

for n in range(class_number): 

this_ideal = class_group_reps[n] 

this_ideal_norm = class_group_rep_norms[n] 

gens = [] 

for i in range(1, height_bound + 1): 

for g in bdd_ideals[i*this_ideal_norm]: 

if g in this_ideal: 

gens.append(g) 

generator_lists.append(gens) 

# Build all the output numbers 

for n in range(class_number): 

gens = generator_lists[n] 

s = len(gens) 

for i in range(s): 

for j in range(i + 1, s): 

if K.ideal(gens[i], gens[j]) == class_group_reps[n]: 

new_number = gens[i]/gens[j] 

for zeta in roots_of_unity: 

yield zeta*new_number 

yield zeta/new_number 

def bdd_norm_pr_ideal_gens(K, norm_list): 

r""" 

Compute generators for all principal ideals in a number field `K` whose 

norms are in ``norm_list``. 

INPUT: 

- `K` - a number field 

- ``norm_list`` - a list of positive integers 

OUTPUT: 

- a dictionary of number field elements, keyed by norm 

EXAMPLES: 

There is only one principal ideal of norm 1, and it is generated by the 

element 1:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_ideal_gens 

sage: K.<g> = QuadraticField(101) 

sage: bdd_norm_pr_ideal_gens(K, [1]) 

{1: [1]} 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_ideal_gens 

sage: K.<g> = QuadraticField(123) 

sage: bdd_norm_pr_ideal_gens(K, range(5)) 

{0: [0], 1: [1], 2: [-g - 11], 3: [], 4: [2]} 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_norm_pr_ideal_gens 

sage: K.<g> = NumberField(x^5 - x + 19) 

sage: b = bdd_norm_pr_ideal_gens(K, range(30)) 

sage: key = ZZ(28) 

sage: b[key] 

[157*g^4 - 139*g^3 - 369*g^2 + 848*g + 158, g^4 + g^3 - g - 7] 

""" 

negative_norm_units = K.elements_of_norm(-1) 

gens = dict() 

if len(negative_norm_units) == 0: 

for n in norm_list: 

if not n: 

gens[n] = [K.zero()] 

else: 

gens[n] = K.elements_of_norm(n) + K.elements_of_norm(-n) 

else: 

for n in norm_list: 

gens[n] = K.elements_of_norm(n) 

return gens 

def integer_points_in_polytope(matrix, interval_radius): 

r""" 

Return the set of integer points in the polytope obtained by acting on a 

cube by a linear transformation. 

Given an r-by-r matrix ``matrix`` and a real number ``interval_radius``, 

this function finds all integer lattice points in the polytope obtained by 

transforming the cube [-interval_radius,interval_radius]^r via the linear 

map induced by ``matrix``. 

INPUT: 

- ``matrix`` - a square matrix of real numbers 

- ``interval_radius`` - a real number 

OUTPUT: 

- a list of tuples of integers 

EXAMPLES: 

Stretch the interval [-1,1] by a factor of 2 and find the integers in the 

resulting interval:: 

sage: from sage.rings.number_field.bdd_height import integer_points_in_polytope 

sage: m = matrix([2]) 

sage: r = 1 

sage: integer_points_in_polytope(m,r) 

[(-2), (-1), (0), (1), (2)] 

Integer points inside a parallelogram:: 

sage: from sage.rings.number_field.bdd_height import integer_points_in_polytope 

sage: m = matrix([[1, 2],[3, 4]]) 

sage: r = RealField()(1.3) 

sage: integer_points_in_polytope(m,r) 

[(-3, -7), (-2, -5), (-2, -4), (-1, -3), (-1, -2), (-1, -1), (0, -1), (0, 0), (0, 1), (1, 1), (1, 2), (1, 3), (2, 4), (2, 5), (3, 7)] 

Integer points inside a parallelepiped:: 

sage: from sage.rings.number_field.bdd_height import integer_points_in_polytope 

sage: m = matrix([[1.2,3.7,0.2],[-5.3,-.43,3],[1.2,4.7,-2.1]]) 

sage: r = 2.2 

sage: L = integer_points_in_polytope(m,r) 

sage: len(L) 

4143 

If ``interval_radius`` is 0, the output should include only the zero tuple:: 

sage: from sage.rings.number_field.bdd_height import integer_points_in_polytope 

sage: m = matrix([[1,2,3,7],[4,5,6,2],[7,8,9,3],[0,3,4,5]]) 

sage: integer_points_in_polytope(m,0) 

[(0, 0, 0, 0)] 

""" 

T = matrix; d = interval_radius; r = T.nrows() 

# Find the vertices of the given box 

box_vertices = [vector(x) for x in product([-d, d], repeat=r)] 

# Transform the vertices 

T_trans = T.transpose() 

transformed_vertices = [v*T_trans for v in box_vertices] 

# Create polyhedron from transformed vertices and find integer points inside 

return list(Polyhedron(transformed_vertices, base_ring=QQ).integral_points()) 

def bdd_height(K, height_bound, precision=53, LLL=False): 

r""" 

Computes all elements in the number field `K` which have relative 

multiplicative height at most ``height_bound``. 

The algorithm requires arithmetic with floating point numbers; 

``precision`` gives the user the option to set the precision for such 

computations. 

It might be helpful to work with an LLL-reduced system of fundamental 

units, so the user has the option to perform an LLL reduction for the 

fundamental units by setting ``LLL`` to True. 

Certain computations may be faster assuming GRH, which may be done 

globally by using the number_field(True/False) switch. 

The function will only be called for number fields `K` with positive unit 

rank. An error will occur if `K` is `QQ` or an imaginary quadratic field. 

ALGORITHM: 

This is an implementation of the main algorithm (Algorithm 3) in 

[Doyle-Krumm]. 

INPUT: 

- ``height_bound`` - real number 

- ``precision`` - (default: 53) positive integer 

- ``LLL`` - (default: False) boolean value 

OUTPUT: 

- an iterator of number field elements 

.. WARNING:: 

In the current implementation, the output of the algorithm cannot be 

guaranteed to be correct due to the necessity of floating point 

computations. In some cases, the default 53-bit precision is 

considerably lower than would be required for the algorithm to 

generate correct output. 

.. TODO:: 

Should implement a version of the algorithm that guarantees correct 

output. See Algorithm 4 in [Doyle-Krumm] for details of an 

implementation that takes precision issues into account. 

EXAMPLES: 

There are no elements of negative height:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = NumberField(x^5 - x + 7) 

sage: list(bdd_height(K,-3)) 

[] 

The only nonzero elements of height 1 are the roots of unity:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = QuadraticField(3) 

sage: list(bdd_height(K,1)) 

[0, -1, 1] 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = QuadraticField(36865) 

sage: len(list(bdd_height(K,101))) # long time (4 s) 

131 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = NumberField(x^3 - 197*x + 39) 

sage: len(list(bdd_height(K, 200))) # long time (5 s) 

451 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = NumberField(x^6 + 2) 

sage: len(list(bdd_height(K,60,precision=100))) # long time (5 s) 

1899 

:: 

sage: from sage.rings.number_field.bdd_height import bdd_height 

sage: K.<g> = NumberField(x^4 - x^3 - 3*x^2 + x + 1) 

sage: len(list(bdd_height(K,10,LLL=true))) 

99 

""" 

B = height_bound 

r1, r2 = K.signature(); r = r1 + r2 -1 

if B < 1: 

return 

yield K(0) 

roots_of_unity = K.roots_of_unity() 

if B == 1: 

for zeta in roots_of_unity: 

yield zeta 

return 

RF = RealField(precision) 

embeddings = K.places(prec=precision) 

logB = RF(B).log() 

def log_map(number): 

r""" 

Computes the image of an element of `K` under the logarithmic map. 

""" 

x = number 

x_logs = [] 

for i in range(r1): 

sigma = embeddings[i] 

x_logs.append(abs(sigma(x)).log()) 

for i in range(r1, r + 1): 

tau = embeddings[i] 

x_logs.append(2*abs(tau(x)).log()) 

return vector(x_logs) 

def log_height_for_generators(n, i, j): 

r""" 

Computes the logarithmic height of elements of the form `g_i/g_j`. 

""" 

gen_logs = generator_logs[n] 

Log_gi = gen_logs[i]; Log_gj = gen_logs[j] 

arch_sum = sum([max(Log_gi[k], Log_gj[k]) for k in range(r + 1)]) 

return (arch_sum - class_group_rep_norm_logs[n]) 

def packet_height(n, pair, u): 

r""" 

Computes the height of the element of `K` encoded by a given packet. 

""" 

gen_logs = generator_logs[n] 

i = pair[0] ; j = pair[1] 

Log_gi = gen_logs[i]; Log_gj = gen_logs[j] 

Log_u_gi = Log_gi + unit_log_dictionary[u] 

arch_sum = sum([max(Log_u_gi[k], Log_gj[k]) for k in range(r + 1)]) 

return (arch_sum - class_group_rep_norm_logs[n]) 

class_group_reps = [] 

class_group_rep_norms = [] 

class_group_rep_norm_logs = [] 

for c in K.class_group(): 

a = c.ideal() 

a_norm = a.norm() 

class_group_reps.append(a) 

class_group_rep_norms.append(a_norm) 

class_group_rep_norm_logs.append(RF(a_norm).log()) 

class_number = len(class_group_reps) 

# Get fundamental units and their images under the log map 

fund_units = UnitGroup(K).fundamental_units() 

fund_unit_logs = [log_map(fund_units[i]) for i in range(r)] 

unit_prec_test = fund_unit_logs 

try: 

[l.change_ring(QQ) for l in unit_prec_test] 

except ValueError: 

raise ValueError('Precision too low.') # QQ(log(0)) may occur if precision too low 

# If LLL is set to True, find an LLL-reduced system of fundamental units 

if LLL: 

cut_fund_unit_logs = column_matrix(fund_unit_logs).delete_rows([r]) 

lll_fund_units = [] 

for c in pari(cut_fund_unit_logs).qflll().sage().columns(): 

new_unit = 1 

for i in range(r): 

new_unit *= fund_units[i]**c[i] 

lll_fund_units.append(new_unit) 

fund_units = lll_fund_units 

fund_unit_logs = [log_map(_) for _ in fund_units] 

unit_prec_test = fund_unit_logs 

try: 

[l.change_ring(QQ) for l in unit_prec_test] 

except ValueError: 

raise ValueError('Precision too low.') # QQ(log(0)) may occur if precision too low 

# Find generators for principal ideals of bounded norm 

possible_norm_set = set([]) 

for n in range(class_number): 

for m in range(1, B + 1): 

possible_norm_set.add(m*class_group_rep_norms[n]) 

bdd_ideals = bdd_norm_pr_ideal_gens(K, possible_norm_set) 

# Distribute the principal ideal generators 

generator_lists = [] 

generator_logs = [] 

for n in range(class_number): 

this_ideal = class_group_reps[n] 

this_ideal_norm = class_group_rep_norms[n] 

gens = [] 

gen_logs = [] 

for i in range(1, B + 1): 

for g in bdd_ideals[i*this_ideal_norm]: 

if g in this_ideal: 

gens.append(g) 

gen_logs.append(log_map(g)) 

generator_lists.append(gens) 

generator_logs.append(gen_logs) 

# Compute the lists of relevant pairs and corresponding heights 

gen_height_dictionary = dict() 

relevant_pair_lists = [] 

for n in range(class_number): 

relevant_pairs = [] 

gens = generator_lists[n] 

s = len(gens) 

for i in range(s): 

for j in range(i + 1, s): 

if K.ideal(gens[i], gens[j]) == class_group_reps[n]: 

relevant_pairs.append([i, j]) 

gen_height_dictionary[(n, i, j)] = log_height_for_generators(n, i, j) 

relevant_pair_lists.append(relevant_pairs) 

# Find the bound for units to be computed 

gen_height_list = [gen_height_dictionary[x] for x in gen_height_dictionary.keys()] 

if len(gen_height_list) == 0: 

d = logB 

else: 

d = logB + max(gen_height_list) 

# Create the matrix whose columns are the logs of the fundamental units 

S = column_matrix(fund_unit_logs).delete_rows([r]) 

try: 

T = S.inverse() 

except ZeroDivisionError: 

raise ValueError('Precision too low.') 

# Find all integer lattice points in the unit polytope 

U = integer_points_in_polytope(T, ceil(d)) 

U0 = []; L0 = [] 

# Compute unit heights 

unit_height_dictionary = dict() 

unit_log_dictionary = dict() 

Ucopy = copy(U) 

for u in U: 

u_log = sum([u[j]*fund_unit_logs[j] for j in range(r)]) 

unit_log_dictionary[u] = u_log 

u_height = sum([max(u_log[k], 0) for k in range(r + 1)]) 

unit_height_dictionary[u] = u_height 

if u_height <= logB: 

U0.append(u) 

if u_height > d: 

Ucopy.remove(u) 

U = Ucopy 

# Sort U by height 

U = sorted(U, key=lambda u: unit_height_dictionary[u]) 

U_length = len(U) 

all_unit_tuples = set(copy(U0)) 

# Check candidate heights 

for n in range(class_number): 

relevant_pairs = relevant_pair_lists[n] 

for pair in relevant_pairs: 

i = pair[0] ; j = pair[1] 

gen_height = gen_height_dictionary[(n, i, j)] 

u_height_bound = logB + gen_height 

for k in range(U_length): 

u = U[k] 

u_height = unit_height_dictionary[u] 

if u_height <= u_height_bound: 

candidate_height = packet_height(n, pair, u) 

if candidate_height <= logB: 

L0.append([n, pair, u]) 

all_unit_tuples.add(u) 

else: 

break 

# Use previous data to build all necessary units 

units_dictionary = dict() 

for u in all_unit_tuples: 

unit = K(1) 

for j in range(r): 

unit *= (fund_units[j])**(u[j]) 

units_dictionary[u] = unit 

# Build all the output numbers 

for u in U0: 

unit = units_dictionary[u] 

for zeta in roots_of_unity: 

yield zeta*unit 

for packet in L0: 

n = packet[0] ; pair = packet[1] ; u = packet[2] 

i = pair[0] ; j = pair[1] 

relevant_pairs = relevant_pair_lists[n] 

gens = generator_lists[n] 

unit = units_dictionary[u] 

c = unit*gens[i]/gens[j] 

for zeta in roots_of_unity: 

yield zeta*c 

yield zeta/c