diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index ad4218e..86ba7a7 100644 --- a/realalg/base_algebraic.py +++ b/realalg/base_algebraic.py @@ -39,7 +39,7 @@ def __init__(self, coefficients, index=-1): # List of integers and integer inde raise ValueError(f'Polynomial {self.sp_polynomial} has no real roots') self.sp_place = real_roots[index] self._accuracy = 0 - self._intervals = None + self._intervals = dict() # pylint:disable=use-dict-literal self.lmbda = None # To be created by instances. def __str__(self): @@ -52,19 +52,24 @@ def __reduce__(self): return (self.__class__, (self.coefficients,)) def __eq__(self, other): return self.coefficients == other.coefficients and self.index == other.index + def find_root_as_interval(self, precision): + ''' Return an Interval around self.lmbda with at least the requested precision. ''' + s = str(sp.N(self.sp_place, precision)) + return Interval.from_string(s, precision) def intervals(self, accuracy): ''' Return intervals around self.lmbda**i with at least the requested accuracy. ''' assert isinstance(accuracy, Integral) assert accuracy > 0 - if accuracy > self._accuracy: + if accuracy not in self._intervals: precision = int(accuracy + self.degree*self.log_bound + 1) + 1 # Cheap ceil. - s = str(sp.N(self.sp_place, precision)) - interval = Interval.from_string(s, precision) - self._intervals = [interval**i for i in range(self.degree)] - assert all(I.accuracy >= accuracy for I in self._intervals) - self._accuracy = accuracy - return [I.simplify(accuracy+1) for I in self._intervals] + interval = self.find_root_as_interval(precision) + intervals = [interval**i for i in range(self.degree)] + assert all(I.accuracy >= accuracy for I in intervals) + intervals = [I.simplify(accuracy+1) for I in intervals] + self._intervals[accuracy] = intervals + self._accuracy = max(accuracy, self._accuracy) + return self._intervals[accuracy] @total_ordering class BaseRealAlgebraic(ABC): @@ -81,6 +86,7 @@ def __init__(self, field, rep): if not self.coefficients: self.coefficients = [Fraction(0, 1)] self.length = sum(LOG_2 + log_plus(coefficient.numerator) + log_plus(coefficient.denominator) + index * self.field.length for index, coefficient in enumerate(self.coefficients)) + self._intervals = dict() # pylint:disable=use-dict-literal def __str__(self): return str(self.N()) def __repr__(self): @@ -172,10 +178,13 @@ def degree(self): def interval(self, accuracy=8): ''' Return an interval around self with at least the requested accuracy. ''' - intermediate_accuracy = int(accuracy + max(log_plus(coefficient) for coefficient in self.coefficients) + len(self.coefficients)) + 1 - interval = sum(coeff * interval for coeff, interval in zip(self.coefficients, self.field.intervals(intermediate_accuracy))) - assert interval.accuracy >= accuracy - return interval.simplify(accuracy+1) + if accuracy not in self._intervals: + intermediate_accuracy = int(accuracy + max(log_plus(coefficient) for coefficient in self.coefficients) + len(self.coefficients)) + 1 + interval = sum(coeff * interval for coeff, interval in zip(self.coefficients, self.field.intervals(intermediate_accuracy))) + assert interval.accuracy >= accuracy + self._intervals[accuracy] = interval.simplify(accuracy + 1) + return self._intervals[accuracy] + def N(self, accuracy=8): ''' Return a string approximating self to at least ``accuracy`` digits. ''' return self.interval(accuracy).midpoint() diff --git a/realalg/cypari2_algebraic.py b/realalg/cypari2_algebraic.py index a57a987..fa405e9 100644 --- a/realalg/cypari2_algebraic.py +++ b/realalg/cypari2_algebraic.py @@ -2,9 +2,13 @@ ''' A module for representing and manipulating real algebraic numbers and the fields that they live in using cypari2. ''' from fractions import Fraction +from math import log2 import numpy as np import cypari2 # pylint: disable=import-error from .base_algebraic import BaseRealNumberField, BaseRealAlgebraic +from .interval import Interval + +LOG_10 = log2(10) cp = cypari2.Pari() cp_x = cp('x') @@ -21,6 +25,15 @@ def __init__(self, coefficients, index=-1): # List of integers and / or Fractio super().__init__(coefficients, index) self.cp_polynomial = cp_polynomial(self.coefficients) self.lmbda = self([0, 1]) + + def find_root_as_interval(self, precision): + bit_prec = int(LOG_10 * precision) + 1 + old_bit_prec = cp.get_real_precision_bits() + cp.set_real_precision_bits(bit_prec) + roots = list(self.cp_polynomial.polrootsreal(precision=bit_prec)) + s = str(roots[self.index]) + cp.set_real_precision_bits(old_bit_prec) + return Interval.from_string(s, precision) def __call__(self, coefficients): return RealAlgebraic(self, cp_polynomial(coefficients).Mod(self.cp_polynomial)) diff --git a/realalg/cypari_algebraic.py b/realalg/cypari_algebraic.py index 7c2b154..6913059 100644 --- a/realalg/cypari_algebraic.py +++ b/realalg/cypari_algebraic.py @@ -2,9 +2,13 @@ ''' A module for representing and manipulating real algebraic numbers and the fields that they live in using cypari. ''' from fractions import Fraction +from math import log2 import numpy as np import cypari # pylint: disable=import-error from .base_algebraic import BaseRealNumberField, BaseRealAlgebraic +from .interval import Interval + +LOG_10 = log2(10) cp = cypari.pari cp_x = cp('x') @@ -21,6 +25,15 @@ def __init__(self, coefficients, index=-1): # List of integers and / or Fractio super().__init__(coefficients, index) self.cp_polynomial = cp_polynomial(self.coefficients) self.lmbda = self([0, 1]) + + def find_root_as_interval(self, precision): + bit_prec = int(LOG_10 * precision) + 1 + old_bit_prec = cp.get_real_precision_bits() + cp.set_real_precision_bits(bit_prec) + roots = list(self.cp_polynomial.polrootsreal(precision=bit_prec)) + s = str(roots[self.index]) + cp.set_real_precision_bits(old_bit_prec) + return Interval.from_string(s, precision) def __call__(self, coefficients): return RealAlgebraic(self, cp_polynomial(coefficients).Mod(self.cp_polynomial))