summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorXavier Caruso <xavier.caruso@univ-rennes1.fr>2017-07-22 00:09:18 -0400
committerXavier Caruso <xavier.caruso@univ-rennes1.fr>2017-07-22 00:09:18 -0400
commit4b2af06c55734697998f25742eea3659334144ae (patch)
tree6a0d6bdce67aa5c3787ac4b1bc829731d36e2f57
parentPseudo-code for lattice precision (diff)
First more-or-less working implementation
-rw-r--r--src/sage/rings/padics/all.py2
-rw-r--r--src/sage/rings/padics/lattice_precision.py367
-rw-r--r--src/sage/rings/padics/padic_base_generic.py5
3 files changed, 296 insertions, 78 deletions
diff --git a/src/sage/rings/padics/all.py b/src/sage/rings/padics/all.py
index 809b4aa..9665553 100644
--- a/src/sage/rings/padics/all.py
+++ b/src/sage/rings/padics/all.py
@@ -7,3 +7,5 @@ from .padic_generic import local_print_mode
from .pow_computer import PowComputer
from .pow_computer_ext import PowComputer_ext_maker
from .discrete_value_group import DiscreteValueGroup
+
+from .lattice_precision import ZpLP
diff --git a/src/sage/rings/padics/lattice_precision.py b/src/sage/rings/padics/lattice_precision.py
index 7f1e98d..640ca37 100644
--- a/src/sage/rings/padics/lattice_precision.py
+++ b/src/sage/rings/padics/lattice_precision.py
@@ -1,110 +1,323 @@
-# The parent
-############
+"""
+Toy implementation of lattice precision for p-adic numbers.
-class pAdicFieldLattice(pAdicRingBaseGeneric):
+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)
+"""
+