From 91d5949af93c12497af5078c9846376cd2e2c7a9 Mon Sep 17 00:00:00 2001 From: Nathan Dunfield Date: Mon, 1 Jun 2026 21:10:00 -0500 Subject: [PATCH 01/11] Better caching of BaseRealNumberField.intervals: 18% speedup --- realalg/__init__.py | 1 + realalg/base_algebraic.py | 20 ++++++++++++-------- realalg/cypari2_algebraic.py | 1 + 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/realalg/__init__.py b/realalg/__init__.py index f0daf54..a2e376a 100644 --- a/realalg/__init__.py +++ b/realalg/__init__.py @@ -14,6 +14,7 @@ RealNumberField = module.RealNumberField RealAlgebraic = module.RealAlgebraic eigenvectors = module.eigenvectors + print(f'Using {interface}') break # We could add some code here to find out which interface was loaded. except ImportError: diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index ad4218e..fe58b8a 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() self.lmbda = None # To be created by instances. def __str__(self): @@ -52,19 +52,23 @@ 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): + 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) + invervals = [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): diff --git a/realalg/cypari2_algebraic.py b/realalg/cypari2_algebraic.py index a57a987..9c3ede6 100644 --- a/realalg/cypari2_algebraic.py +++ b/realalg/cypari2_algebraic.py @@ -18,6 +18,7 @@ class RealNumberField(BaseRealNumberField): __engine = 'cypari2' def __init__(self, coefficients, index=-1): # List of integers and / or Fractions, integer index + print(coefficients, index) super().__init__(coefficients, index) self.cp_polynomial = cp_polynomial(self.coefficients) self.lmbda = self([0, 1]) From b82010965a15790bfe093327fea9ca41265b2557 Mon Sep 17 00:00:00 2001 From: Nathan Dunfield Date: Mon, 1 Jun 2026 22:14:48 -0500 Subject: [PATCH 02/11] cypari2 backend: use PARI to find roots, about a 10-times speedup over sympy --- realalg/__init__.py | 1 - realalg/cypari2_algebraic.py | 11 ++++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/realalg/__init__.py b/realalg/__init__.py index a2e376a..f0daf54 100644 --- a/realalg/__init__.py +++ b/realalg/__init__.py @@ -14,7 +14,6 @@ RealNumberField = module.RealNumberField RealAlgebraic = module.RealAlgebraic eigenvectors = module.eigenvectors - print(f'Using {interface}') break # We could add some code here to find out which interface was loaded. except ImportError: diff --git a/realalg/cypari2_algebraic.py b/realalg/cypari2_algebraic.py index 9c3ede6..b7294ce 100644 --- a/realalg/cypari2_algebraic.py +++ b/realalg/cypari2_algebraic.py @@ -5,6 +5,7 @@ import numpy as np import cypari2 # pylint: disable=import-error from .base_algebraic import BaseRealNumberField, BaseRealAlgebraic +from .interval import Interval cp = cypari2.Pari() cp_x = cp('x') @@ -18,10 +19,18 @@ class RealNumberField(BaseRealNumberField): __engine = 'cypari2' def __init__(self, coefficients, index=-1): # List of integers and / or Fractions, integer index - print(coefficients, index) 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(3.32192809488737 * 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)) From 0f41e728de76c71de874a5a6bc0fdce21d007135 Mon Sep 17 00:00:00 2001 From: Nathan Dunfield Date: Tue, 2 Jun 2026 13:34:41 -0500 Subject: [PATCH 03/11] cypari backend: same improvement as for cypari2 --- realalg/cypari_algebraic.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/realalg/cypari_algebraic.py b/realalg/cypari_algebraic.py index 7c2b154..6d34d9a 100644 --- a/realalg/cypari_algebraic.py +++ b/realalg/cypari_algebraic.py @@ -5,6 +5,8 @@ import numpy as np import cypari # pylint: disable=import-error from .base_algebraic import BaseRealNumberField, BaseRealAlgebraic +from .interval import Interval + cp = cypari.pari cp_x = cp('x') @@ -21,6 +23,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(3.32192809488737 * 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)) From cca994cf4722a35ee725ec776c058a66cae75228 Mon Sep 17 00:00:00 2001 From: Nathan Dunfield Date: Tue, 2 Jun 2026 14:29:48 -0500 Subject: [PATCH 04/11] Added benchmark script --- tests/flipper_speed.py | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 tests/flipper_speed.py diff --git a/tests/flipper_speed.py b/tests/flipper_speed.py new file mode 100644 index 0000000..bd0ccd2 --- /dev/null +++ b/tests/flipper_speed.py @@ -0,0 +1,4 @@ +import flipper +F = flipper.create_triangulation([(~20, ~16, ~17), (~19, ~18, ~15), (~14, ~12, 20), (~13, 15, ~8), (~11, 19, 18), (~10, 16, 17), (~9, 10, ~4), (~7, ~3, 9), (~6, 14, 13), (~5, 12, 11), (~2, 7, 8), (~1, 6, 5), (~0, 3, 4), (0, 2, 1)]) +h = F.encode([{0: 3}, 9, 1, 5, 7, 17, 20, 19, 16, 10, 18, 9, 11, 7, 5, 17, 20, 14, 5, 13, 8, 2, 5, 1, 6, 13]) +h.canonical().stratum() From 784228414ba29e2932e9322a2f15df1df87a5cf2 Mon Sep 17 00:00:00 2001 From: Nathan Dunfield Date: Mon, 15 Jun 2026 10:52:44 -0500 Subject: [PATCH 05/11] Big speedup by caching BaseRealAlgebraic.inverval --- realalg/base_algebraic.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index fe58b8a..ef56d9c 100644 --- a/realalg/base_algebraic.py +++ b/realalg/base_algebraic.py @@ -85,6 +85,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() def __str__(self): return str(self.N()) def __repr__(self): @@ -176,10 +177,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() From 352d1f8d0a1919f1d5eba7987e8ba58ece8ea042 Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:05:44 +0100 Subject: [PATCH 06/11] Derive constant --- realalg/cypari2_algebraic.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/realalg/cypari2_algebraic.py b/realalg/cypari2_algebraic.py index b7294ce..fa405e9 100644 --- a/realalg/cypari2_algebraic.py +++ b/realalg/cypari2_algebraic.py @@ -2,11 +2,14 @@ ''' 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') @@ -24,7 +27,7 @@ def __init__(self, coefficients, index=-1): # List of integers and / or Fractio self.lmbda = self([0, 1]) def find_root_as_interval(self, precision): - bit_prec = int(3.32192809488737 * precision) + 1 + 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)) From 79d9da90caaf8fab9e92335e378aeeb528c3817f Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:07:23 +0100 Subject: [PATCH 07/11] Derive constant --- realalg/cypari_algebraic.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/realalg/cypari_algebraic.py b/realalg/cypari_algebraic.py index 6d34d9a..6913059 100644 --- a/realalg/cypari_algebraic.py +++ b/realalg/cypari_algebraic.py @@ -2,11 +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') @@ -25,7 +27,7 @@ def __init__(self, coefficients, index=-1): # List of integers and / or Fractio self.lmbda = self([0, 1]) def find_root_as_interval(self, precision): - bit_prec = int(3.32192809488737 * precision) + 1 + 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)) From 6c38e27af9e464cf36194f199445fa95e6d68abd Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:09:30 +0100 Subject: [PATCH 08/11] Typo --- realalg/base_algebraic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index ef56d9c..0620c01 100644 --- a/realalg/base_algebraic.py +++ b/realalg/base_algebraic.py @@ -65,7 +65,7 @@ def intervals(self, accuracy): 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) - invervals = [I.simplify(accuracy+1) 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] From abf0ece28aa48b34b535ff314e902738425fdc13 Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:18:22 +0100 Subject: [PATCH 09/11] Remove trailing whitespace --- realalg/base_algebraic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index 0620c01..88754f1 100644 --- a/realalg/base_algebraic.py +++ b/realalg/base_algebraic.py @@ -182,7 +182,7 @@ def interval(self, accuracy=8): 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] + return self._intervals[accuracy] def N(self, accuracy=8): ''' Return a string approximating self to at least ``accuracy`` digits. ''' From 1e01876056f451d6bce7f5a8ec5148784001f824 Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:20:45 +0100 Subject: [PATCH 10/11] Delete tests/flipper_speed.py Flipper_speed was not a unit test --- tests/flipper_speed.py | 4 ---- 1 file changed, 4 deletions(-) delete mode 100644 tests/flipper_speed.py diff --git a/tests/flipper_speed.py b/tests/flipper_speed.py deleted file mode 100644 index bd0ccd2..0000000 --- a/tests/flipper_speed.py +++ /dev/null @@ -1,4 +0,0 @@ -import flipper -F = flipper.create_triangulation([(~20, ~16, ~17), (~19, ~18, ~15), (~14, ~12, 20), (~13, 15, ~8), (~11, 19, 18), (~10, 16, 17), (~9, 10, ~4), (~7, ~3, 9), (~6, 14, 13), (~5, 12, 11), (~2, 7, 8), (~1, 6, 5), (~0, 3, 4), (0, 2, 1)]) -h = F.encode([{0: 3}, 9, 1, 5, 7, 17, 20, 19, 16, 10, 18, 9, 11, 7, 5, 17, 20, 14, 5, 13, 8, 2, 5, 1, 6, 13]) -h.canonical().stratum() From 0c5f50c85cbc491eae4a85a0df922acf0970ddfb Mon Sep 17 00:00:00 2001 From: Mark Bell Date: Sat, 20 Jun 2026 00:23:45 +0100 Subject: [PATCH 11/11] Lint --- realalg/base_algebraic.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/realalg/base_algebraic.py b/realalg/base_algebraic.py index 88754f1..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 = dict() + self._intervals = dict() # pylint:disable=use-dict-literal self.lmbda = None # To be created by instances. def __str__(self): @@ -53,6 +53,7 @@ def __reduce__(self): 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) @@ -85,7 +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() + self._intervals = dict() # pylint:disable=use-dict-literal def __str__(self): return str(self.N()) def __repr__(self):