ZpL: a p-adic precision package

Legend

  • ZpCR stands for Capped Relative precision: this is roughly the analogue of interval arithmetic (a $p$-adic number is represented by a subset of the shape $x + O(p^n)$ and precision is tracked by operating of these subsets)
  • ZpFP stands for Floating Point arithmetic: this is the analogue of usual floating point arithmetic (just a certain subset of $p$-adic numbers are available and results are rounded so that they always lie in this subset)
  • ZpLC stands for Lattice Precision with cap (provided by our package)
  • ZpLF stands for Lattice Precision combined with Floating Point arithmetics (provided by our package)

Creation of the rings

In [1]:
Z2 = Zp(2, prec=30, print_mode="digits")
Q2 = Z2.fraction_field()
In [2]:
Z3 = Zp(3, prec=20, print_mode="digits")

First examples

In [3]:
x = random_element(Z3,5)
x
ZpCR [3]:
0.000 s
...02112
ZpFP [3]:
0.000 s
...00000000000000002112
ZpLC [3]:
0.000 s
...02112
ZpLF [3]:
0.000 s
...02112

When we multiply by $p$, we gain one digit of absolute precision:

In [4]:
3*x
ZpCR [4]:
0.000 s
...021120
ZpFP [4]:
0.000 s
...000000000000000021120
ZpLC [4]:
0.004 s
...021120
ZpLF [4]:
0.004 s
...021120

The ZpL package sees this gain of precision even if the computation is split into several steps:

In [5]:
x + x + x
ZpCR [5]:
0.004 s
...21120
ZpFP [5]:
0.000 s
...000000000000000021120
ZpLC [5]:
0.000 s
...021120
ZpLF [5]:
0.000 s
...021120
In [6]:
y = 2*x
y
ZpCR [6]:
0.000 s
...12001
ZpFP [6]:
0.000 s
...00000000000000012001
ZpLC [6]:
0.000 s
...12001
ZpLF [6]:
0.000 s
...12001
In [7]:
x + y
ZpCR [7]:
0.000 s
...21120
ZpFP [7]:
0.000 s
...000000000000000021120
ZpLC [7]:
0.000 s
...021120
ZpLF [7]:
0.000 s
...021120

The same with multiplication:

In [8]:
x^3
ZpCR [8]:
0.000 s
...022122
ZpFP [8]:
0.000 s
...00000000120222022122
ZpLC [8]:
0.000 s
...022122
ZpLF [8]:
0.000 s
...022122
In [9]:
x * x * x
ZpCR [9]:
0.000 s
...22122
ZpFP [9]:
0.000 s
...00000000120222022122
ZpLC [9]:
0.000 s
...022122
ZpLF [9]:
0.000 s
...022122
In [10]:
y = x^2
y
ZpCR [10]:
0.000 s
...00021
ZpFP [10]:
0.000 s
...00000000000020100021
ZpLC [10]:
0.000 s
...00021
ZpLF [10]:
0.000 s
...00021
In [11]:
x * y
ZpCR [11]:
0.000 s
...22122
ZpFP [11]:
0.000 s
...00000000120222022122
ZpLC [11]:
0.000 s
...022122
ZpLF [11]:
0.000 s
...022122

Computations with numbers at different precision

In [12]:
x = random_element(Z2,10)
y = random_element(Z2,5)
In [13]:
u = x+y
u
ZpCR [13]:
0.000 s
...10101
ZpFP [13]:
0.000 s
...000000000000000000000111010101
ZpLC [13]:
0.000 s
...10101
ZpLF [13]:
0.000 s
...10101
In [14]:
v = x-y
v
ZpCR [14]:
0.000 s
...00111
ZpFP [14]:
0.000 s
...000000000000000000000111000111
ZpLC [14]:
0.000 s
...00111
ZpLF [14]:
0.000 s
...00111

Now we compute $u + v$, which is also $2x$. Here is what we get:

In [15]:
u+v
ZpCR [15]:
0.000 s
...11100
ZpFP [15]:
0.000 s
...00000000000000000000001110011100
ZpLC [15]:
0.000 s
...01110011100
ZpLF [15]:
0.004 s
...01110011100
In [16]:
2*x
ZpCR [16]:
0.000 s
...01110011100
ZpFP [16]:
0.000 s
...00000000000000000000001110011100
ZpLC [16]:
0.000 s
...01110011100
ZpLF [16]:
0.000 s
...01110011100

The SOMOS 4 sequence

The SOMOS 4 sequence is the sequence defined by the recurrence $$u_{n+4} = \frac{u_{n+1}u_{n+3} + u_{n+2}^2}{u_n}.$$ It turns out that its computation is highly unstable.
For this reason, using the ZpL package can here be decisive.

In [17]:
def somos(u0, u1, u2, u3, n):
    a, b, c, d = u0, u1, u2, u3
    for _ in range(4, n+1):
        a, b, c, d = b, c, d, (b*d + c*c) / a
    return d
In [18]:
somos(Z2(1,15), Z2(1,15), Z2(1,15), Z2(3,15), 18)
ZpCR [18]:
0.004 s
...11
ZpFP [18]:
0.000 s
...110111001000100100000000000111
ZpLC [18]:
0.052 s
...100000000000111
ZpLF [18]:
0.032 s
...100000000000111
In [19]:
somos(Z2(1,15), Z2(1,15), Z2(1,15), Z2(3,15), 100)
ZpCR [19]:
0.000 s
PrecisionError: cannot divide by something indistinguishable from zero.
ZpFP [19]:
0.000 s
...110000101000011001001001110001
ZpLC [19]:
0.364 s
...001001001110001
ZpLF [19]:
0.228 s
...001001001110001

Computations with matrices

Product of many matrices

In [20]:
MS = MatrixSpace(Z2,2)
In [21]:
M = MS(1)
for _ in range(60):
    M *= random_element(MS, prec=8)
In [22]:
M
ZpCR [22]:
0.000 s
[ ...010100000000000000 ...1111000000000000000]
[   ...0011000000000000   ...10010000000000000]
ZpFP [22]:
0.000 s
[ ...00101100101000111110111100010100000000000000 ...011010011000000010010100111111000000000000000]
[   ...100110110001111110010001100011000000000000   ...0111001111111001101010101110010000000000000]
ZpLC [22]:
0.000 s
[ ...00010100000000000000 ...111111000000000000000]
[   ...100011000000000000   ...1110010000000000000]
ZpLF [22]:
0.000 s
[ ...00010100000000000000 ...111111000000000000000]
[   ...100011000000000000   ...1110010000000000000]

Determinants

In [23]:
D = diagonal_matrix([ Z3(1,5), Z3(9,5), Z3(27,5), Z3(81,5) ])
D
ZpCR [23]:
0.040 s
[...00001        0        0        0]
[       0 ...00100        0        0]
[       0        0 ...01000        0]
[       0        0        0 ...10000]
ZpFP [23]:
0.000 s
[    ...00000000000000000001                           0                           0                           0]
[                          0   ...0000000000000000000100                           0                           0]
[                          0                           0  ...00000000000000000001000                           0]
[                          0                           0                           0 ...000000000000000000010000]
ZpLC [23]:
0.004 s
[...00001        0        0        0]
[       0 ...00100        0        0]
[       0        0 ...01000        0]
[       0        0        0 ...10000]
ZpLF [23]:
0.000 s
[...00001        0        0        0]
[       0 ...00100        0        0]
[       0        0 ...01000        0]
[       0        0        0 ...10000]

The determinant of $D$ is known at higher precision. This is due to the particular shape of $D$.

In [24]:
D.determinant()
ZpCR [24]:
0.012 s
...1000000000
ZpFP [24]:
0.004 s
...00000000000000000001000000000
ZpLC [24]:
0.152 s
...1000000000
ZpLF [24]:
0.056 s
...1000000000

This phenomenon is preserved when $D$ is modified by (left and right) multiplication by some matrix. However the software may or may not see it, depending on the method it uses for tracking precision

In [25]:
MS = D.parent()
P = random_element(MS, prec=5)
Q = random_element(MS, prec=5)
M = P*D*Q
M
ZpCR [25]:
0.000 s
[ ...00120  ...20020  ...20110 ...200100]
[ ...12000  ...21200  ...20200  ...10100]
[ ...21000  ...11100  ...22100  ...22200]
[ ...20021  ...01101  ...21002  ...20220]
ZpFP [25]:
0.004 s
[  ...000000011020110000120   ...000000002100101020020   ...000000011002210020110  ...0000000010121011200100]
[...00000000020011211112000  ...0000000011121221221200  ...0000000021100012120200  ...0000000012221210210100]
[...00000000100010222221000  ...0000000012111001011100  ...0000000100012210022100  ...0000000021210022222200]
[   ...00000011220000220021    ...00000010010221001101    ...00000012020102121002   ...000000011100122120220]
ZpLC [25]:
0.408 s
[ ...00120  ...20020  ...20110 ...200100]
[ ...12000  ...21200  ...20200  ...10100]
[ ...21000  ...11100  ...22100  ...22200]
[ ...20021  ...01101  ...21002  ...20220]
ZpLF [25]:
0.428 s
[ ...00120  ...20020  ...20110 ...200100]
[ ...12000  ...21200  ...20200  ...10100]
[ ...21000  ...11100  ...22100  ...22200]
[ ...20021  ...01101  ...21002  ...20220]
In [26]:
M.determinant()
ZpCR [26]:
0.004 s
0
ZpFP [26]:
0.000 s
...000000010121201222210000000000
ZpLC [26]:
0.560 s
...10000000000
ZpLF [26]:
0.376 s
...10000000000
In [27]:
P.determinant() * D.determinant() * Q.determinant()
ZpCR [27]:
0.004 s
...10000000000
ZpFP [27]:
0.000 s
...001101010121201222210000000000
ZpLC [27]:
0.676 s
...10000000000
ZpLF [27]:
0.568 s
...10000000000

Notice that ZpCR cannot decide whether the determinant $M$ vanishes. This may have important consequences if, at some point, we will test if $M$ is invertible or not and jump to one branch or another one depending on the result.

In [28]:
M.charpoly()
ZpCR [28]:
0.004 s
(...00000000000000000001)*x^4 + (...10120)*x^3 + (...02100)*x^2 + (0)*x + (0)
ZpFP [28]:
0.000 s
...00000000000000000001*x^4 + ...222222011120011010120*x^3 + ...0022121020011212002100*x^2 + ...00012200222121121011000000*x + ...000000010121201222210000000000
ZpLC [28]:
0.356 s
...00000000000000000001*x^4 + ...10120*x^3 + ...02100*x^2 + ...1000000*x + ...10000000000
ZpLF [28]:
0.360 s
...00000000000000000001*x^4 + ...10120*x^3 + ...02100*x^2 + ...1000000*x + ...10000000000

Dodgson condensation algorithm

Dodgson condensation is a recursive method for computing the determinant of a matrix. A description of the method can be found here.

In [29]:
def dodgson(M):   # M must be a square matrix
    n = M.nrows() - 1
    A = M
    B = matrix(n, n, 
               [ [ M[i,j]*M[i+1,j+1] - M[i+1,j]*M[i,j+1] 
                   for j in range(n) ] for i in range(n) ])
    for m in range(n-1, 0, -1):
        A, B = B, matrix(m, m, 
                         [ [ (B[i,j]*B[i+1,j+1] - B[i+1,j]*B[i,j+1]) / A[i+1,j+1] 
                             for j in range(m) ] for i in range(m) ])
    return B[0,0]
In [30]:
MS = MatrixSpace(Z2, 6)
M = random_element(MS, prec=10)
In [31]:
dodgson(M)
ZpCR [31]:
0.012 s
...1110
ZpFP [31]:
0.008 s
...0101110010011100100110100001110
ZpLC [31]:
1.004 s
...0100001110
ZpLF [31]:
0.528 s
...0100001110
In [32]:
M.determinant()
ZpCR [32]:
0.004 s
...0100001110
ZpFP [32]:
0.000 s
...0001110010011100100110100001110
ZpLC [32]:
1.980 s
...0100001110
ZpLF [32]:
1.668 s
...0100001110

Computations with polynomials

In [33]:
S.<x> = PolynomialRing(Q2)
P = random_element(S, degree=10, prec=5)
Q = random_element(S, degree=10, prec=5)
D = x^5 + random_element(S, degree=4, prec=8)
D
ZpCR [33]:
0.004 s
(...000000000000000000000000000001)*x^5 + (...00110010)*x^4 + (...00011111)*x^3 + (...00010010)*x^2 + (...00011111)*x + (...10101000)
ZpFP [33]:
0.000 s
...000000000000000000000000000001*x^5 + ...0000000000000000000000000110010*x^4 + ...000000000000000000000000011111*x^3 + ...0000000000000000000000000010010*x^2 + ...000000000000000000000000011111*x + ...000000000000000000000000010101000
ZpLC [33]:
0.076 s
...000000000000000000000000000001*x^5 + ...00110010*x^4 + ...00011111*x^3 + ...00010010*x^2 + ...00011111*x + ...10101000
ZpLF [33]:
0.080 s
...000000000000000000000000000001*x^5 + ...00110010*x^4 + ...00011111*x^3 + ...00010010*x^2 + ...00011111*x + ...10101000

With high probability, $P$ and $Q$ are coprime, implying that the gcd of $DP$ is $DQ$ is $D$.
However, we observe that we do not always get this expected result:

In [34]:
def euclidean(A,B):
    while B != 0:
        A, B = B, A % B
    return A.monic()
In [35]:
euclidean(D*P, D*Q)
ZpCR [35]:
0.004 s
(0)*x^13 + (0)*x^12 + (0)*x^11 + (...01)*x^10 + (0)*x^9 + (0)*x^8 + (0)*x^7 + (0)*x^6 + (...10)*x^5 + (0)*x^4 + (0)*x^3 + (0)*x^2 + (0)*x + (0)
ZpFP [35]:
0.004 s
...000000000000000000000000000001
ZpLC [35]:
4.024 s
...000000000000000000000000000001*x^5 + ...0010*x^4 + ...0011111*x^3 + ...0010*x^2 + ...1111*x + ...0101000
ZpLF [35]:
2.812 s
...000000000000000000000000000001*x^5 + ...0010*x^4 + ...0011111*x^3 + ...0010*x^2 + ...1111*x + ...0101000

More precisely, ZpCR often found a gcd of higher degree. Indeed the Euclidean algorithm stops prematurely because the software believes that some early remainder vanishes due to the lack of precision.
The situation with ZpFP is the exact opposite: since the precision is too high (always the maximal cap), the software does not notice that the remainder vanishes. Therefore the Euclidean algorithm does not stop until it obtains a remainder of degree $0$.

Gröbner basis

In [36]:
R.<x,y,z> = PolynomialRing(Q2, order = 'invlex')
F = [Q2(2,10)*x + Q2(1,10)*z, 
     Q2(1,10)*x^2 + Q2(1,10)*y^2 - Q2(2,10)*z^2,
     Q2(4,10)*y^2 + Q2(1,10)*y*z + Q2(8,10)*z^2]
In [37]:
from sage.rings.polynomial.toy_buchberger import buchberger_improved
g = buchberger_improved(ideal(F)); g.sort()
g
ZpCR [37]:
0.016 s
[x^3, x*y + ...1100010*x^2, y^2 + ...11001*x^2, z + ...0000000010*x]
ZpFP [37]:
0.016 s
[x^3, x*y + ...1111111111111111111111111100010*x^2, y^2 + ...011111111111111111111111111001*x^2, z + ...0000000000000000000000000000010*x]
ZpLC [37]:
1.656 s
[x^3, x*y + ...111100010*x^2, y^2 + ...1111111001*x^2, z + ...0000000010*x]
ZpLF [37]:
1.868 s
[x^3, x*y + ...111100010*x^2, y^2 + ...1111111001*x^2, z + ...0000000010*x]

$p$-adic differential equations

In [38]:
def edp_simple(g, h, N):
    S = g.parent()
    u = S.gen()
    for i in range(N):
        u -= h(u) * (u.derivative()/h(u) - g).integral()
        u = u[:2^(i+1)]
    return u
In [39]:
N = 4
S.<t> = PowerSeriesRing(Q2, 2^N)
In [40]:
h = 1 + t + t^3
y = t + t^2 * random_element(S,10)
g = y.derivative() / h(y)
In [41]:
u = edp_simple(g, h, N)
u[15]
ZpCR [41]:
0.056 s
...10011
ZpFP [41]:
0.008 s
...011000000000000000000011110011
ZpLC [41]:
41.623 s
...0011110011
ZpLF [41]:
25.746 s
...0011110011