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)
  • ZpLP stands for Lattice Precision (our package)

Creation of the rings

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

First examples

In [2]:
x = random_element(Z3,5)
x
ZpCR [2]:
...00222
ZpFP [2]:
...00000000000000000222
ZpLP [2]:
...00222

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

In [3]:
3*x
ZpCR [3]:
...002220
ZpFP [3]:
...000000000000000002220
ZpLP [3]:
...002220

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

In [4]:
x + x + x
ZpCR [4]:
...02220
ZpFP [4]:
...000000000000000002220
ZpLP [4]:
...002220
In [5]:
y = 2*x
y
ZpCR [5]:
...01221
ZpFP [5]:
...00000000000000001221
ZpLP [5]:
...01221
In [6]:
x + y
ZpCR [6]:
...02220
ZpFP [6]:
...000000000000000002220
ZpLP [6]:
...002220

The same with multiplication:

In [7]:
x^3
ZpCR [7]:
...002222
ZpFP [7]:
...00000000000220002222
ZpLP [7]:
...002222
In [8]:
x * x * x
ZpCR [8]:
...02222
ZpFP [8]:
...00000000000220002222
ZpLP [8]:
...002222
In [9]:
y = x^2
y
ZpCR [9]:
...21001
ZpFP [9]:
...00000000000000221001
ZpLP [9]:
...21001
In [10]:
x * y
ZpCR [10]:
...02222
ZpFP [10]:
...00000000000220002222
ZpLP [10]:
...002222

Computations with numbers at different precision

In [11]:
x = random_element(Z2,10)
y = random_element(Z2,5)
In [12]:
u = x+y
u
ZpCR [12]:
...01101
ZpFP [12]:
...000000000000000000001011001101
ZpLP [12]:
...01101
In [13]:
v = x-y
v
ZpCR [13]:
...00001
ZpFP [13]:
...000000000000000000001011000001
ZpLP [13]:
...00001

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

In [14]:
u+v
ZpCR [14]:
...01110
ZpFP [14]:
...0000000000000000000010110001110
ZpLP [14]:
...10110001110
In [15]:
2*x
ZpCR [15]:
...10110001110
ZpFP [15]:
...0000000000000000000010110001110
ZpLP [15]:
...10110001110

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 [16]:
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 [17]:
somos(Z2(1,15), Z2(1,15), Z2(1,15), Z2(3,15), 18)
ZpCR [17]:
...11
ZpFP [17]:
...110111001000100100000000000111
ZpLP [17]:
...100000000000111
In [18]:
somos(Z2(1,15), Z2(1,15), Z2(1,15), Z2(3,15), 100)
ZpCR [18]:
PrecisionError: cannot divide by something indistinguishable from zero.
ZpFP [18]:
...110000101000011001001001110001
ZpLP [18]:
...001001001110001

Computations with matrices

Product of many matrices

In [19]:
MS = MatrixSpace(Z2,2)
In [20]:
M = MS(1)
for _ in range(60):
    M *= random_element(MS, prec=8)
In [21]:
M
ZpCR [21]:
[0 0]
[0 0]
ZpFP [21]:
[...10101011001100101101000111001100000000000000000 ...11101110000110101011010000101100000000000000000]
[...00001000100100111111100100100100000000000000000 ...01001111011111011111111001000100000000000000000]
ZpLP [21]:
[...1001100000000000000000 ...0101100000000000000000]
[...0100100000000000000000 ...1000100000000000000000]

Determinants

In [22]:
D = diagonal_matrix([ Z3(1,5), Z3(9,5), Z3(27,5), Z3(81,5) ])
D
ZpCR [22]:
[...00001        0        0        0]
[       0 ...00100        0        0]
[       0        0 ...01000        0]
[       0        0        0 ...10000]
ZpFP [22]:
[    ...00000000000000000001                           0                           0                           0]
[                          0   ...0000000000000000000100                           0                           0]
[                          0                           0  ...00000000000000000001000                           0]
[                          0                           0                           0 ...000000000000000000010000]
ZpLP [22]:
[...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 [23]:
D.determinant()
ZpCR [23]:
...1000000000
ZpFP [23]:
...00000000000000000001000000000
ZpLP [23]:
...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 [24]:
MS = D.parent()
P = random_element(MS, prec=5)
Q = random_element(MS, prec=5)
M = P*D*Q
M
ZpCR [24]:
[ ...21020  ...22200  ...02011  ...20121]
[ ...11220  ...10200  ...01021  ...00001]
[ ...20000  ...12100  ...01200  ...20100]
[...221200  ...12200  ...01210  ...02110]
ZpFP [24]:
[   ...000000000222212121020   ...0000000001021020022200     ...00000000221022202011     ...00000000101212020121]
[   ...000000010200001011220   ...0000000011120021010200     ...00000010200110101021     ...00000002222201100001]
[...000000000002200020020000   ...0000000010100111112100   ...0000000002210202101200   ...0000000002020212220100]
[  ...0000000011202110221200   ...0000000020001021212200    ...000000011201221201210    ...000000010021211202110]
ZpLP [24]:
[ ...21020  ...22200  ...02011  ...20121]
[ ...11220  ...10200  ...01021  ...00001]
[ ...20000  ...12100  ...01200  ...20100]
[...221200  ...12200  ...01210  ...02110]
In [25]:
M.determinant()
ZpCR [25]:
0
ZpFP [25]:
...00000021000002111222000000000
ZpLP [25]:
...2000000000
In [26]:
P.determinant() * D.determinant() * Q.determinant()
ZpCR [26]:
...2000000000
ZpFP [26]:
...01012021000002111222000000000
ZpLP [26]:
...2000000000

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 [27]:
M.charpoly()
ZpCR [27]:
(...00000000000000000001)*x^4 + (...10000)*x^3 + (0)*x^2 + (...200000)*x + (0)
ZpFP [27]:
...00000000000000000001*x^4 + ...000222222120100021010000*x^3 + ...0001011101102010221200000*x^2 + ...2220101021012200020200000*x + ...00000021000002111222000000000
ZpLP [27]:
...00000000000000000001*x^4 + ...10000*x^3 + ...0200000*x + ...2000000000

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 [28]:
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 [29]:
MS = MatrixSpace(Z2, 6)
M = random_element(MS, prec=10)
In [30]:
dodgson(M)
ZpCR [30]:
...01111
ZpFP [30]:
...110011000001111011111000101111
ZpLP [30]:
...1000101111
In [31]:
M.determinant()
ZpCR [31]:
...1000101111
ZpFP [31]:
...110011000001111011111000101111
ZpLP [31]:
...1000101111

Computations with polynomials

In [32]:
S.<x> = PolynomialRing(Z2)
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 [32]:
(...000000000000000000000000000001)*x^5 + (...10111100)*x^4 + (...11100111)*x^3 + (...01110101)*x^2 + (...00001011)*x + (...10100110)
ZpFP [32]:
...000000000000000000000000000001*x^5 + ...00000000000000000000000010111100*x^4 + ...000000000000000000000011100111*x^3 + ...000000000000000000000001110101*x^2 + ...000000000000000000000000001011*x + ...0000000000000000000000010100110
ZpLP [32]:
...000000000000000000000000000001*x^5 + ...10111100*x^4 + ...11100111*x^3 + ...01110101*x^2 + ...00001011*x + ...10100110

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 [33]:
def euclidean(A,B):
    while B != 0:
        A, B = B, A % B
    return A.monic()
In [34]:
euclidean(D*P, D*Q)
ZpCR [34]:
(...1)*x^6 + (...10)*x^5 + (...1)*x^4 + (...1)*x^3 + (...1)*x^2 + (...100)*x + (...100)
ZpFP [34]:
...000000000000000000000000000001
ZpLP [34]:
...000000000000000000000000000001*x^5 + ...10111100*x^4 + ...11100111*x^3 + ...01110101*x^2 + ...00001011*x + ...10100110

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$.