summaryrefslogtreecommitdiffstats
path: root/src/sage/rings/padics/lattice_precision.py
diff options
context:
space:
mode:
Diffstat (limited to 'src/sage/rings/padics/lattice_precision.py')
-rw-r--r--src/sage/rings/padics/lattice_precision.py832
1 files changed, 538 insertions, 294 deletions
diff --git a/src/sage/rings/padics/lattice_precision.py b/src/sage/rings/padics/lattice_precision.py
index 640ca37..1808afe 100644
--- a/src/sage/rings/padics/lattice_precision.py
+++ b/src/sage/rings/padics/lattice_precision.py
@@ -1,323 +1,567 @@
-"""
-Toy implementation of lattice precision for p-adic numbers.
-
-Here is a small demo::
-
- sage: R = ZpLP(3)
- sage: x = R(1,10)
-
-Of course, when we multiply by 3, we gain one digit of absolute
-precision::
-
- sage: 3*x
- 3 + O(3^11)
-
-The lattice precision machinery sees this even if we decompose
-the computation into several steps::
-
- sage: y = x+x
- sage: y
- 2 + O(3^10)
- sage: x+y
- 3 + O(3^11)
-
-The same works for the multiplication::
-
- sage: z = x^2
- sage: z
- 1 + O(3^10)
- sage: x*z
- 1 + O(3^11)
-
-This comes more funny when we are working with elements given
-at different precisions::
-
- sage: R = ZpLP(2)
- sage: x = R(1,10)
- sage: y = R(1,5)
- sage: z = x+y; z
- 2 + O(2^5)
- sage: t = x-y; t
- 0 + O(2^5)
- sage: z+t # observe that z+t = 2*x
- 2 + O(2^11)
- sage: z-t # observe that z-t = 2*y
- 2 + O(2^6)
-
- sage: x = R(28888,15)
- sage: y = R(204,10)
- sage: z = x/y; z
- 242 + O(2^9)
- sage: z*y # which is x
- 28888 + O(2^15)
-
-The SOMOS sequence is the sequence defined by the recurrence::
-
-..MATH::
-
- u_n = \frac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}}
-
-It is known for its numerical instability.
-On the one hand, one can show that if the initial values are
-invertible in `\mathbb{Z}_p` and known at precision `O(p^N)`
-then all the next terms of the SOMOS sequence will be known
-at the same precision as well.
-On the other hand, because of the division, when we unroll
-the recurrence, we loose a lot of precision. Observe::
-
- sage: R = Zp(2, print_mode='terse')
- sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 4 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 13 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 55 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 21975 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 6639 + O(2^13)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 7186 + O(2^13)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 569 + O(2^13)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 253 + O(2^13)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 4149 + O(2^13)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 2899 + O(2^12)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 3072 + O(2^12)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 349 + O(2^12)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 619 + O(2^12)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 243 + O(2^12)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 3 + O(2^2)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 2 + O(2^2)
-
-If instead, we use the lattice precision, everything goes well::
-
- sage: R = ZpLP(2)
- sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 4 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 13 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 55 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 21975 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 23023 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 31762 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 16953 + O(2^15)
- sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
- 16637 + O(2^15)
-
- sage: for _ in range(100):
- ....: a,b,c,d = b,c,d,(b*d+c*c)/a
- sage: a
- 15519 + O(2^15)
- sage: b
- 32042 + O(2^15)
- sage: c
- 17769 + O(2^15)
- sage: d
- 20949 + O(2^15)
-"""
-
import _weakref as weakref
-from sage.misc.prandom import randint
+from sage.misc.misc import walltime
+
+from sage.structure.sage_object import SageObject
+from sage.structure.unique_representation import UniqueRepresentation
+
from sage.rings.integer_ring import ZZ
-from sage.rings.padics.generic_nodes import pAdicRingBaseGeneric
-from sage.rings.padics.padic_generic_element import pAdicGenericElement
+from sage.rings.rational_field import QQ
+from sage.rings.infinity import Infinity
+
+
+# Class pRational
+#################
+
+# It is a temporary class implementing rational numbers
+# as approximations of p-adic numbers
+
+class pRational:
+ def __init__(self, p, x, exponent=0, valuation=None):
+ self.p = p
+ if x in ZZ:
+ self.x = ZZ(x)
+ else:
+ self.x = x
+ self.exponent = exponent
+ self._valuation = valuation
+
+ def __repr__(self):
+ if self.exponent == 0:
+ return str(self.x)
+ else:
+ return "%s^%s * %s" % (self.p, self.exponent, self.x)
+
+ def reduce(self, prec):
+ if prec is Infinity:
+ return self
+ x = self.x
+ exp = self.exponent
+ if x.parent() is ZZ:
+ if prec > exp:
+ x = x % (self.p ** (prec-exp))
+ else:
+ x = 0
+ elif x.parent() is QQ:
+ num = x.numerator()
+ denom = x.denominator()
+ valdenom = denom.valuation(self.p)
+ denom /= self.p ** valdenom
+ exp -= valdenom
+ modulo = self.p ** (prec - exp)
+ # probably we should use Newton iteration instead
+ # (but it is actually slower for now - Python implementation)
+ _, inv, _ = denom.xgcd(modulo)
+ x = (num*inv) % modulo
+ if self.x == 0:
+ val = Infinity
+ else:
+ val = self._valuation
+ return self.__class__(self.p, x, exp, valuation=val)
+
+ def normalize(self):
+ if self.x == 0:
+ self.exponent = 0
+ else:
+ val = self.valuation()
+ exp = self.exponent
+ self.x /= self.p ** (val-exp)
+ if self.x in ZZ:
+ self.x = ZZ(self.x)
+ self.exponent = val
+
+ def valuation(self):
+ if self._valuation is None:
+ valx = self.x.valuation(self.p)
+ self._valuation = self.exponent + valx
+ return self._valuation
+
+ def is_p_power(self):
+ self.normalize()
+ return self.x == 1
+
+ def is_zero(self):
+ return self.x == 0
+
+ def __add__(self, other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if self._valuation is None or other._valuation is None:
+ val = None
+ elif self._valuation < other._valuation:
+ val = self._valuation
+ elif self._valuation > other._valuation:
+ val = other._valuation
+ else:
+ val = None
+ if sexp < oexp:
+ return self.__class__(p, self.x + other.x * p**(oexp-sexp), sexp, valuation=val)
+ else:
+ return self.__class__(p, self.x * p**(sexp-oexp) + other.x, oexp, valuation=val)
+
+ def __sub__(self, other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if self._valuation is None or other._valuation is None:
+ val = None
+ elif self._valuation < other._valuation:
+ val = self._valuation
+ elif self._valuation > other._valuation:
+ val = other._valuation
+ else:
+ val = None
+ if sexp < oexp:
+ return self.__class__(p, self.x - other.x * p**(oexp-sexp), sexp, valuation=val)
+ else:
+ return self.__class__(p, self.x * p**(sexp-oexp) - other.x, oexp, valuation=val)
+
+ def __neg__(self):
+ return self.__class__(self.p, -self.x, self.exponent, valuation=self._valuation)
+
+ def __mul__(self, other):
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = self._valuation + other._valuation
+ return self.__class__(self.p, self.x * other.x, self.exponent + other.exponent, valuation=val)
+
+ def __div__(self, other):
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = self._valuation - other._valuation
+ return self.__class__(self.p, self.x / other.x, self.exponent - other.exponent, valuation=val)
-class pAdicRingLattice(pAdicRingBaseGeneric):
+ def __lshift__(self, n):
+ if self._valuation is None:
+ val = None
+ else:
+ val = self._valuation + n
+ return self.__class__(self.p, self.x, self.exponent + n, valuation=val)
+
+ def __rshift__(self, n):
+ return self << (-n)
+
+ def unit_part(self):
+ if self.is_zero():
+ raise ValueError("the unit part of zero is not defined")
+ p = self.p
+ val = self.valuation()
+ x = self.x / (p ** (val-self.exponent))
+ return self.__class__(p, x, 0, valuation=0)
+
+ def xgcd(self,other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if sexp < oexp:
+ a = ZZ(self.x)
+ b = ZZ(other.x * (p ** (oexp-sexp)))
+ exp = sexp
+ else:
+ a = ZZ(self.x * (p ** (sexp-oexp)))
+ b = ZZ(other.x)
+ exp = oexp
+ d, u, v = a.xgcd(b)
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = min(self._valuation, other._valuation)
+ d = self.__class__(p, d, exp, valuation=val)
+ u = self.__class__(p, u)
+ v = self.__class__(p, v)
+ return d, u, v
+
+ def value(self):
+ return (self.p ** self.exponent) * self.x
+
+ def list(self, prec):
+ if self.x not in ZZ:
+ raise NotImplementedError
+ val = self.valuation()
+ p = self.p
+ x = ZZ(self.x * p**(self.exponent - val))
+ l = [ ]; i = val
+ while i < prec:
+ x, digit = x.quo_rem(p)
+ l.append(digit)
+ i += 1
+ return l
+
+
+# Class PrecisionLattice
+########################
+
+def list_of_padics(elements):
+ from sage.rings.padics.padic_lattice_element import pAdicLatticeElement
+ if isinstance(elements, pAdicLatticeElement):
+ return [ weakref.ref(elements) ]
+ if not isinstance(elements, list):
+ elements = list(elements)
+ ans = [ ]
+ for x in elements:
+ ans += list_of_padics(x)
+ return ans
+
+def format_history(tme, status, timings):
+ status = ''.join(status)
+ if timings:
+ if tme < 0.000001:
+ s = " --- "
+ elif tme >= 10:
+ s = " >= 10s "
+ else:
+ s = "%.6fs" % tme
+ return s + " " + status
+ else:
+ return status
+
+class PrecisionLattice(UniqueRepresentation, SageObject):
# Internal variables:
- # . self._prec_cap
+ # . self._cap
# a cap for the (working) precision
- # meaning that the precision lattice always contains p^(self._prec_cap)
+ # meaning that the precision lattice always contains p^(self._cap)
# . self._elements
# list of weak references of elements in this parent
- # . self._precision
+ # . self._lattice
# an upper triangular matrix over ZZ representing the
# lattice of precision
# (its columns are indexed by self._elements)
-
- def __init__(self, p, prec=50, print_mode={'mode':'terse'}, label=None):
- if label is None:
- self._label = None
- else:
- self._label = str(label)
+ # . self._absolute_precisions
+ def __init__(self, p, label):
+ self._p = p
+ self._label = label
self._elements = [ ]
- self._precision = { }
- self._prec_cap = prec
- self._modulus = p**prec
- pAdicRingBaseGeneric.__init__(self, p, prec, print_mode, '', pAdicLatticeElement)
-
- def _prec_type(self):
- return 'lattice'
+ self._lattice = { }
+ self._absolute_precisions = { }
+ self._marked_for_deletion = [ ]
+ self._approx_zero = pRational(p, ZZ(0))
+ # History
+ self._history_init = None
+ self._history = None
def label(self):
return self._label
- def precision_cap(self):
- return self._prec_cap
-
def _repr_(self):
- s = "%s-adic Field with lattice precision" % self.prime()
- if self._label is not None:
- s += ", called '%s'" % self._label
- return s
+ n = len(self._elements)
+ if self._label is None:
+ if n > 1:
+ return "Precision Lattice on %s objects" % len(self._elements)
+ else:
+ return "Precision Lattice on %s object" % len(self._elements)
+ else:
+ if n > 1:
+ return "Precision Lattice on %s objects (label: %s)" % (len(self._elements), self._label)
+ else:
+ return "Precision Lattice on %s object (label: %s)" % (len(self._elements), self._label)
+
+ def prime(self):
+ return self._p
- def _new_element(self, x, prec, dx):
+ def reduce(self, index=0, partial=False):
+ n = len(self._elements)
+ if index >= n-1:
+ return
+ if partial:
+ # Partial reduction
+ # Cost: O(m^2) with m = n-index
+ tme = walltime()
+ diffval = (n-index) * [0]
+ for j in range(n-1, index, -1):
+ col = self._lattice[self._elements[j]]
+ prec = col[j].valuation() - diffval[j-index]
+ for i in range(index,j):
+ col[i] = col[i].reduce(prec)
+ col[i].normalize() # seems to be faster then
+ dval = col[i].valuation() - prec
+ if dval < diffval[i-index]:
+ diffval[i-index] = dval
+ # We update history
+ if self._history is not None:
+ self._history.append(('partial reduce', index, walltime(tme)))
+ else:
+ # Full Hermite reduction
+ # Cost: O(m^3) with m = n-index
+ tme = walltime()
+ for j in range(index+1, n):
+ # In what follows, we assume that col[j] is a power of p
+ col = self._lattice[self._elements[j]]
+ valpivot = col[j].valuation()
+ for i in range(index, j):
+ reduced = col[i].reduce(valpivot)
+ scalar = (col[i] - reduced) >> valpivot
+ if scalar.is_zero(): continue
+ col[i] = reduced
+ col[i].normalize()
+ for j2 in range(j+1, n):
+ col2 = self._lattice[self._elements[j2]]
+ col2[i] -= scalar*col2[i]
+ col2[i].normalize()
+ # We update history
+ if self._history is not None:
+ self._history.append(('full reduce', index, walltime(tme)))
+
+ def new_element(self, x, dx, bigoh, dx_mode='linear_combinaison'):
+ # First we delete some elements marked for deletion
+ if self._marked_for_deletion:
+ self.del_elements(thresold=50)
+
+ # Then we add the new element
+ tme = walltime()
+ p = self._p
n = len(self._elements)
- x_ref = weakref.ref(x, self._del_element)
+ x_ref = weakref.ref(x, self.mark_for_deletion)
self._elements.append(x_ref)
- if prec is None or prec > self._prec_cap:
- prec = self._prec_cap
- modulus = self.prime()**prec
- col = n * [ZZ(0)]
- for ref,scalar in dx:
- c = self._precision[ref]
- for i in range(len(c)):
- col[i] += scalar * c[i]
+ col = n * [self._approx_zero]
+ if dx_mode == 'linear_combinaison':
+ for elt,scalar in dx:
+ ref = weakref.ref(elt)
+ if not isinstance(scalar, pRational):
+ scalar = pRational(p, scalar)
+ c = self._lattice[ref]
+ for i in range(len(c)):
+ col[i] += scalar * c[i]
+ elif dx_mode == 'values':
+ for elt,scalar in dx:
+ ref = weakref.ref(elt)
+ if not isinstance(scalar, pRational):
+ scalar = pRational(p, scalar)
+ i = len(self._lattice[ref]) - 1
+ col[i] = scalar
+ else:
+ raise ValueError("dx_mode must be either 'linear_combinaison' or 'values'")
for i in range(n):
- col[i] = ZZ(col[i]) % modulus
- # We should handle the case where col[i] is not an integer but I'm lazy...
- col.append(modulus)
- self._precision[x_ref] = col
-
- def _del_element(self, ref):
+ col[i] = col[i].reduce(bigoh)
+ col.append(pRational(p, ZZ(1), bigoh))
+ self._lattice[x_ref] = col
+
+ # We compute the absolute precision of the new element and cache it
+ self._absolute_precisions[x_ref] = min([ c.valuation() for c in col ])
+
+ # We update history
+ if self._history is not None:
+ self._history.append(('add', None, walltime(tme)))
+
+ def _set_precision(self, x, column={}):
+ p = self._p
+ x_ref = weakref.ref(x)
+ index = len(self._lattice[x_ref]) - 1
+ n = len(self._elements)
+ col = n * [self._approx_zero]
+ self._lattice[x_ref] = col
+ self._absolute_precisions[x_ref] = min([ c.valuation() for c in col ])
+
+ def mark_for_deletion(self, ref):
+ tme = walltime()
try:
- index = len(self._precision[ref]) - 1
+ index = len(self._lattice[ref]) - 1
except IndexError:
return
- del self._elements[index]
- del self._precision[ref]
+ self._marked_for_deletion.append(index)
+ if self._history is not None:
+ self._history.append(('mark', index, walltime(tme)))
- # Now, we echelonize...
- p = self.prime()
+ def del_elements(self, thresold=None):
+ p = self._p
n = len(self._elements)
- modulus = self._modulus
- for i in range(index,n):
- col = self._precision[self._elements[i]]
- d, u, v = col[i].xgcd(col[i+1])
- up, vp = col[i+1]/d, -col[i]/d
- col[i] = d
- del col[i+1]
- for j in range(i+1,n):
- col = self._precision[self._elements[j]]
- col[i], col[i+1] = (u*col[i] + v*col[i+1]) % modulus, (up*col[i] + vp*col[i+1]) % modulus
- # Actually, the matrix we obtained this way is not in Hermite normal form
- # because its entries above the diagonal are not necessarily completely reduced
- # (though they are at least modulo p^precision_cap)
- # Should we reduce them further (this affects seriously the complexity)?
- def precision_absolute(self, x):
+ self._marked_for_deletion.sort(reverse=True)
+ count = 0
+ for index in self._marked_for_deletion:
+ if thresold is not None and index < n - thresold: break
+ n -= 1; count += 1
+
+ tme = walltime()
+ del self._lattice[self._elements[index]]
+ del self._elements[index]
+
+ # Now, we echelonize
+ for i in range(index,n):
+ col = self._lattice[self._elements[i]]
+ vali = col[i].valuation()
+ valj = col[i+1].valuation()
+ d, u, v = col[i].xgcd(col[i+1])
+ up, vp = col[i+1]/d, col[i]/d
+ col[i] = d
+ del col[i+1]
+ for j in range(i+1,n):
+ col = self._lattice[self._elements[j]]
+ col[i], col[i+1] = u*col[i] + v*col[i+1], up*col[i] - vp*col[i+1]
+
+ # We update history
+ if self._history is not None:
+ self._history.append(('del', index, walltime(tme)))
+
+ # And we reduce a bit
+ # (we do not perform a complete reduction because it is costly)
+ self.reduce(index, partial=True)
+
+ del self._marked_for_deletion[:count]
+
+ def lift_to_precision(self, x, prec):
ref = weakref.ref(x)
- p = self.prime()
- return min([ c.valuation(p) for c in self._precision[ref] ])
-
- def working_precision(self, x):
- # Do something better here
- return self._prec_cap
-
-ZpLP = pAdicRingLattice
-
-
-# The elements
-##############
-
-class pAdicLatticeElement(pAdicGenericElement):
- # Internal variable:
- # . self._approximation
- # an approximation of this p-adic number
- # it is defined modulo p^(working_precision)
-
- def __init__(self, parent, x, prec=None, dx={}):
- self._parent = parent
-