P-Adic square roots

From Rosetta Code
Task
P-Adic square roots
You are encouraged to solve this task according to the task description, using any language you may know.


Task.

Convert rational a/b to its approximate p-adic square root. To check the result, square the root and construct rational m/n to compare with radicand a/b.

For rational reconstruction Lagrange's lattice basis reduction algorithm is used.

Recipe: find root x1 modulo p and build a sequence of solutions f(xk) ≡ 0 (mod pk),
using the lifting equation xk+1 = xk + dk * pk  with dk = –(f(xk) / pk) / f ′(x1) (mod p).
The multipliers dk are the successive p-adic digits to find.

If evaluation of f(x) = bx2a overflows, the expansion is cut off and might be too short to retrieve the radicand. Setting a higher precision won't help, using a programming language with built-in large integer support will.


Related task.

p-Adic numbers, basic


Reference.

[1] Solving x2a (mod n)



FreeBASIC

' ***********************************************
'subject: p-adic square roots, Hensel lifting.
'tested : FreeBasic 1.07.0

'The root is squared, approximated by a
'rational, and compared with radicand a/b.


const emx = 48
'exponent maximum

const amx = 700000
'tentative argument maximum

'------------------------------------------------
const Mxd = cdbl(2)^53 - 1
'max. float64 integer

const Pmax = 32749
'max. prime < 2^15


type ratio
   as longint a, b
end type

type padic
declare function sqrt (byref q as ratio, byval sw as integer) as integer
'p-adic square root of q = a/b, set sw to print
declare sub printf (byval sw as integer)
'print expansion, set sw to print rational
declare function crat (byval sw as integer) as ratio
'rational reconstruction

declare sub cmpt (byref a as padic)
'let self:= complement_a
declare sub sqr (byref a as padic)
'let self:= a ^ 2

   as long d(-emx to emx - 1)
   as integer v
end type


'global variables
dim shared as long p1, p = 7
'default prime
dim shared as integer k = 11
'precision

#define min(a, b) iif((a) > (b), b, a)


'------------------------------------------------
'p-adic square root of g = a/b
function padic.sqrt (byref g as ratio, byval sw as integer) as integer
dim as longint a = g.a, b = g.b
dim as longint q, x, pk, pm
dim as long f1, r, s, t
dim i as integer, f as double
sqrt = 0

   if b = 0 then return 1
   if b < 0 then b = -b: a = -a
   if p < 2 or k < 1 then return 1

   'max. short prime
   p = min(p, Pmax)

   if sw then
      'echo numerator, denominator,
      print a;"/";str(b);" + ";
      'prime and precision
      print "O(";str(p);"^";str(k);")"
   end if

   'initialize
   v = 0
   p1 = p - 1
   for i = -emx to emx - 1
      d(i) = 0: next

   if a = 0 then return 0

   'valuation
   do until b mod p
      b \= p: v -= 1
   loop
   do until a mod p
      a \= p: v += 1
   loop

   if (v and 1) = 1 then
      'odd valuation
      print "non-residue mod"; p
      return -1
   end if

   'max. array length
   k = min(k + v, emx - 1)
   k -= v: v shr= 1

   if abs(a) > amx or b > amx then return -1

   if p = 2 then
      '1 / b = b (mod 8)
      'a / b = 1 (mod 8)
      t = a * b
      if (t and 7) - 1 then
         print "non-residue mod 8"
         return -1
      end if

   else
      'find root for small p
      for r = 1 to p1
         q = b * r * r - a
         if q mod p = 0 then exit for
      next r

      if r = p then
         print "non-residue mod"; p
         return -1
      end if

      'f'(r) = 2br
      t = b * r shl 1

      s = 0
      t mod= p
      'modular inverse for small p
      for f1 = 1 to p1
         s += t
         if s > p1 then s -= p
         if s = 1 then exit for
      next f1

      if f1 = p then
         print "impossible inverse mod"
         return -1
      end if
   end if

   'evaluate f(x)
   #macro evalf(x)
      f = b * x * cdbl(x / pk)
      f -= cdbl(a / pk)
      'overflow
      if f > Mxd then exit for
      q = clngint(f)
   #endmacro

   if p = 2 then
      'initialize
      x = 1
      d(v) = 1
      d(v + 1) = 0

      pk = 4
      for i = v + 2 to k - 1 + v
         pk shl= 1
         '2-power overflow
         if pk < 1 then exit for
         evalf(x)
         'next digit
         d(i) = iif(q and 1, 1, 0)
         'lift x
         x += d(i) * (pk shr 1)
      next i

   else
      '-1 / f'(x) mod p
      f1 = p - f1
      x = r
      d(v) = x

      pk = 1
      for i = v + 1 to k - 1 + v
         pm = pk: pk *= p
         if pk \ pm - p then exit for
         evalf(x)
         d(i) = q * f1 mod p
         if d(i) < 0 then d(i) += p
         x += d(i) * pk
      next i
   end if
   k = i - v

   if sw then print "lift:";x;" mod";p;"^";str(k)
end function

'------------------------------------------------
'rational reconstruction
function padic.crat (byval sw as integer) as ratio
dim as integer i, j, t = min(v, 0)
dim as longint s, pk, pm
dim as long q, x, y
dim as double f, h
dim r as ratio

   'weighted digit sum
   s = 0: pk = 1
   for i = t to k - 1 + v
      pm = pk: pk *= p

      if pk \ pm - p then
         'overflow
         pk = pm: exit for
      end if

      s += d(i) * pm '(mod pk)
   next i

   'lattice basis reduction
   dim as longint m(1) = {pk, s}
   dim as longint n(1) = {0, 1}
   'norm(v)^2
   h = cdbl(s) * s + 1
   i = 0: j = 1

   'Lagrange's algorithm
   do
      f = m(i) * (m(j) / h)
      f += n(i) * (n(j) / h)

      'Euclidean step
      q = int(f +.5)
      m(i) -= q * m(j)
      n(i) -= q * n(j)

      f = h
      h = cdbl(m(i)) * m(i)
      h += cdbl(n(i)) * n(i)
      'compare norms
      if h < f then
        'interchange vectors
         swap i, j
      else
         exit do
      end if
   loop

   x = m(j): y = n(j)
   if y < 0 then y = -y: x = -x

   'check determinant
   t = abs(m(i) * y - x * n(i)) = pk

   if t = 0 then
      print "crat: fail"
      x = 0: y = 1

   else
      'negative powers
      for i = v to -1
         y *= p: next

      if sw then
         print x;
         if y > 1 then print "/";str(y);
         print
      end if
   end if

   r.a = x: r.b = y
return r
end function


'print expansion
sub padic.printf (byval sw as integer)
dim as integer i, t = min(v, 0)

   for i = k - 1 + t to t step -1
      print d(i);
      if i = 0 andalso v < 0 then print ".";
   next i
   print

   'rational approximation
   if sw then crat(sw)
end sub

'------------------------------------------------
'let self:= complement_a
sub padic.cmpt (byref a as padic)
dim i as integer, r as padic
dim as long c = 1
with r
  .v = a.v

   for i = .v to k +.v
      c += p1 - a.d(i)
      'carry
      if c > p1 then
        .d(i) = c - p: c = 1
      else
        .d(i) = c: c = 0
      end if
   next i
end with
this = r
end sub

'let self:= a ^ 2
sub padic.sqr (byref a as padic)
dim as long ptr rp, ap = @a.d(a.v)
dim as longint q, c = 0
dim as integer i, j
dim r as padic
with r
  .v = a.v shl 1
   rp = @.d(.v)

   for i = 0 to k
      for j = 0 to i
         c += ap[j] * ap[i - j]
      next j
      'Euclidean step
      q = c \ p
      rp[i] = c - q * p
      c = q
   next i
end with
this = r
end sub


'main
'------------------------------------------------
dim as integer sw
dim as padic a, c
dim as ratio q, r

width 64, 30
cls

'   -7 + O(2^7)
data -7,1, 2,7
data 9,1, 2,8
data 17,1, 2,9
data 497,10496, 2,18
data 10496,497, 2,19

data -577215,664901, 3,23
data 15403,26685, 3,18

data -1,1, 5,8
data 86,25, 5,8
data 2150,1, 5,8

data 2,1, 7,8
data 11696,621467, 7,11
data -27764,11521, 7,11
data -27584,12953, 7,11

data -166420,135131, 11,11
data 14142,135623, 5,15
data -255,256, 257,3

data 0,0, 0,0


print
do
   read q.a,q.b, p,k

   sw = a.sqrt(q, 1)
   if sw = 1 then exit do
   if sw then ? : continue do

   print "sqrt +/-"
   print "...";
   a.printf(0)
   a.cmpt(a)
   print "...";
   a.printf(0)

   c.sqr(a)
   print "sqrt^2"
   print "   ";
   c.printf(0)
   r = c.crat(1)

   '{r = q}
   if q.a * r.b - r.a * q.b then
      print "fail: sqrt^2"
   end if

   print : ?
loop

end
Examples:

-7/1 + O(2^7)
lift: 53 mod 2^7
sqrt +/-
... 0 1 1 0 1 0 1
... 1 0 0 1 0 1 1
sqrt^2
    1 1 1 1 0 0 1
-7


 9/1 + O(2^8)
lift: 253 mod 2^8
sqrt +/-
... 1 1 1 1 1 1 0 1
... 0 0 0 0 0 0 1 1
sqrt^2
    0 0 0 0 1 0 0 1
 9


 17/1 + O(2^9)
lift: 233 mod 2^9
sqrt +/-
... 0 1 1 1 0 1 0 0 1
... 1 0 0 0 1 0 1 1 1
sqrt^2
    0 0 0 0 1 0 0 0 1
 17


 497/10496 + O(2^18)
lift: 92989 mod 2^18
sqrt +/-
... 0 1 0 1 1 0 1 0 1 1 0 0 1 1. 1 1 0 1
... 1 0 1 0 0 1 0 1 0 0 1 1 0 0. 0 0 1 1
sqrt^2
    1 0 0 0 0 0 1 1 0 0. 1 0 0 0 1 0 0 1
 497/10496


 10496/497 + O(2^19)
lift: 148501 mod 2^19
sqrt +/-
... 1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0
... 0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0
sqrt^2
    0 0 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0
 10496/497


-577215/664901 + O(3^23)
lift: 44765673211 mod 3^23
sqrt +/-
... 1 0 2 1 1 1 2 2 1 0 1 1 1 2 1 1 0 1 0 2 0 1 0
... 1 2 0 1 1 1 0 0 1 2 1 1 1 0 1 1 2 1 2 0 2 2 0
sqrt^2
    0 1 2 1 1 0 0 2 1 1 0 0 1 1 0 2 0 1 1 0 1 0 0
-577215/664901


 15403/26685 + O(3^18)
lift: 238961782 mod 3^18
sqrt +/-
... 1 2 1 1 2 2 1 2 2 1 1 1 2 2 1 1 0. 1
... 1 0 1 1 0 0 1 0 0 1 1 1 0 0 1 1 2. 2
sqrt^2
    1 2 1 2 0 1 2 2 0 1 1 0 1 2 2 2. 0 1
 15403/26685


-1/1 + O(5^8)
lift: 280182 mod 5^8
sqrt +/-
... 3 2 4 3 1 2 1 2
... 1 2 0 1 3 2 3 3
sqrt^2
    4 4 4 4 4 4 4 4
-1


 86/25 + O(5^8)
lift: 95531 mod 5^8
sqrt +/-
... 1 1 0 2 4 1 1. 1
... 3 3 4 2 0 3 3. 4
sqrt^2
    0 0 0 0 0 3. 2 1
 86/25


 2150/1 + O(5^8)
lift: 95531 mod 5^8
sqrt +/-
... 1 0 2 4 1 1 1 0
... 3 4 2 0 3 3 4 0
sqrt^2
    0 0 0 3 2 1 0 0
 2150


 2/1 + O(7^8)
lift: 1802916 mod 7^8
sqrt +/-
... 2 1 2 1 6 2 1 3
... 4 5 4 5 0 4 5 4
sqrt^2
    0 0 0 0 0 0 0 2
 2


 11696/621467 + O(7^11)
lift: 967215488 mod 7^11
sqrt +/-
... 3 2 6 5 3 1 2 4 1 4. 1
... 3 4 0 1 3 5 4 2 5 2. 6
sqrt^2
    3 3 1 1 3 3 4 4 5. 1 1
 11696/621467


-27764/11521 + O(7^11)
lift: 453209984 mod 7^11
sqrt +/-
... 1 4 1 4 2 1 3 5 6 2 3
... 5 2 5 2 4 5 3 1 0 4 4
sqrt^2
    5 1 0 2 5 5 5 3 6 6 2
-27764/11521


-27584/12953 + O(7^11)
lift: 1239987302 mod 7^11
sqrt +/-
... 4 2 5 0 4 5 0 1 2 2 1
... 2 4 1 6 2 1 6 5 4 4 6
sqrt^2
    3 2 6 5 3 1 2 4 1 4 1
-27584/12953


-166420/135131 + O(11^11)
lift: 34219617965 mod 11^11
sqrt +/-
... 1 3 5 7 0 0 9 10 5 0 5
... 9 7 5 3 10 10 1 0 5 10 6
sqrt^2
    1 4 1 4 2 1 3 5 6 2 3
-166420/135131


 14142/135623 + O(5^15)
lift: 236397452 mod 5^15
sqrt +/-
... 0 0 0 4 4 1 0 0 4 2 0 4 3 0 2
... 4 4 4 0 0 3 4 4 0 2 4 0 1 4 3
sqrt^2
    4 2 1 3 3 4 3 4 3 4 2 3 2 0 4
 14142/135623


-255/256 + O(257^3)
lift: 8867596 mod 257^3
sqrt +/-
... 134 66 68
... 122 190 189
sqrt^2
    255 255 255
-255/256

Haskell

Uses solution for P-Adic_numbers,_basic#Haskell

{-# LANGUAGE KindSignatures, DataKinds  #-}

import Data.Ratio
import Data.List (find)
import GHC.TypeLits
import Padic

pSqrt :: KnownNat p => Rational -> Padic p
pSqrt r = res 
  where
    res = maybe Null mkUnit series 
    (a, b) = (numerator r, denominator r)
    series = case modulo res of
      
      2 | eqMod 4 a 3 -> Nothing
        | not (eqMod 8 a 1) -> Nothing
        | otherwise -> Just $ 1 : 0 : go 8 1
        where
          go pk x =
            let q = ((b*x*x - a) `div` pk) `mod` 2
            in q : go (2*pk) (x + q * (pk `div` 2))

      p -> do
        y <- find (\x -> eqMod p (b*x*x) a) [1..p-1]
        df <- recipMod p (2*b*y)
        let go pk x =
              let f = (b*x*x - a) `div` pk
                  d = (f * (p - df)) `mod` p
              in x `div` (pk `div` p) : go (p*pk) (x + d*pk)
        Just $ go p y

eqMod :: Integral a => a -> a -> a -> Bool
eqMod p a b = a `mod` p == b `mod` p
λ> :set -XDataKinds
λ> pSqrt 2 :: Padic 7
7-adic: 12011266421216213.0
λ> pSqrt 39483088 :: Padic 7
7-adic: 22666526525245311.0
λ> pSqrt 5 :: Padic 7
Null
λ> pSqrt 9 :: Padic 2
2-adic: 11111111111111101.0
λ> it ^ 2
2-adic: 1001.0
λ> toRational it
9 % 1
λ> pSqrt 17 ^ 2 == (17 :: Padic 2)
True
λ> pSqrt 50 ^ 2 == (50 :: Padic 7)
True
λ> pSqrt 52 ^ 2 == (52 :: Padic 7)
False
λ> pSqrt 52 :: Padic 7
Null

Julia

using Nemo, LinearAlgebra

set_printing_mode(FlintPadicField, :terse)

""" convert to Rational (rational reconstruction) """
function toRational(pa::padic)
    rat = lift(QQ, pa)
    r, den = BigInt(numerator(rat)), Int(denominator(rat))
    p, k = Int(prime(parent(pa))), Int(precision(pa))
    N = BigInt(p^k)
    a1, a2 = [N, 0], [r, 1]
    while dot(a1, a1) > dot(a2, a2)
        q = dot(a1, a2) // dot(a2, a2)
        a1, a2 = a2, a1 - round(q) * a2
    end
    if dot(a1, a1) < N
        return (Rational{Int}(a1[1]) // Rational{Int}(a1[2])) // Int(den)
    else
        return Int(r) // den
    end
end

function dstring(pa::padic)
    u, v, n, p, k = pa.u, pa.v, pa.N, pa.parent.p, pa.parent.prec_max
    d = digits(v > 0 ? u * p^v : u, base=p, pad=k)
    return prod([i == k + v && v != 0 ? "$x . " : "$x " for (i, x) in enumerate(reverse(d))])
end

const DATA = [
    [-7, 1, 2, 7],
    [9, 1, 2, 8],
    [17, 1, 2, 9],
    [-1,  1,  5, 8],
    [86, 25,  5, 8],
    [2150, 1,  5, 8],
    [2, 1, 7, 8],
    [3029, 4821, 7, 9],
    [379, 449, 7, 8],
    [717, 8, 11, 7],
    [1414, 213, 41, 5],
    [-255, 256, 257, 3]
]

for (num1, den1, P, K) in DATA
    Qp = PadicField(P, K)
    a = Qp(QQ(num1 // den1))
    c = sqrt(a)
    r = toRational(c * c)
    println(a, "\nsqrt +/-\n", dstring(c), "\n", dstring(-c), "\nCheck sqrt^2:\n", dstring(c * c), "\n", r, "\n")
end
Output:
121 + O(2^7)  
sqrt +/-      
1 1 1 0 1 0 1 
0 0 0 1 0 1 1 
Check sqrt^2: 
1 1 1 1 0 0 1
-7//1

9 + O(2^8)
sqrt +/-
1 1 1 1 1 1 0 1
0 0 0 0 0 0 1 1
Check sqrt^2:
0 0 0 0 1 0 0 1
9//1

17 + O(2^9)
sqrt +/-
0 1 1 1 0 1 0 0 1
1 0 0 0 1 0 1 1 1
Check sqrt^2:
0 0 0 0 1 0 0 0 1
17//1

390624 + O(5^8)
sqrt +/-
3 2 4 3 1 2 1 2
1 2 0 1 3 2 3 3
Check sqrt^2:
4 4 4 4 4 4 4 4
-1//1

86/25 + O(5^6)
sqrt +/-
1 1 0 2 4 1 1 . 1
3 3 4 2 0 3 3 . 4
Check sqrt^2:
0 0 0 0 0 3 . 2 1
86//25

2150 + O(5^10)
sqrt +/-
1 1 0 2 4 1 1 1 0 .
3 3 4 2 0 3 3 4 0 .
Check sqrt^2:
0 0 0 3 2 1 0 0 
2150//1

2 + O(7^8)
sqrt +/-
2 1 2 1 6 2 1 3
4 5 4 5 0 4 5 4
Check sqrt^2:
0 0 0 0 0 0 0 2
2//1

39483088 + O(7^9)
sqrt +/-
5 2 5 2 4 5 3 1 1
1 4 1 4 2 1 3 5 6
Check sqrt^2:
6 5 6 4 1 3 0 2 1
3029//4821

4493721 + O(7^8)
sqrt +/-
0 4 5 0 1 2 2 1
6 2 1 6 5 4 4 6
Check sqrt^2:
5 3 1 2 4 1 4 1
379//449

2435986 + O(11^7)
sqrt +/-
2 10 8 7 0 1 5
8 0 2 3 10 9 6
Check sqrt^2:
1 4 1 4 2 1 3
717//8

36443037 + O(41^5)
sqrt +/-
9 1 38 6 8
31 39 2 34 33
Check sqrt^2:
12 36 31 15 23
1414//213

16908285 + O(257^3)
sqrt +/-
134 66 68
122 190 189
Check sqrt^2:
255 255 255
-255//256

Nim

Translation of: FreeBASIC
Translation of: Wren
import strformat

const
  Emx = 64        # Exponent maximum.
  Amx = 6000      # Argument maximum.
  PMax = 32749    # Prime maximum.

type

  Ratio = tuple[a, b: int]

  Padic = object
    p: int                        # Prime.
    k: int                        # Precision.
    v: int
    d: array[-Emx..(Emx-1), int]

  PadicError = object of ValueError


proc sqrt(pa: var Padic; q: Ratio; sw: bool) =
  ## Return the p-adic square root of q = a/b. Set sw to print.

  var (a, b) = q
  var i, x: int

  if b == 0:
    raise newException(PadicError, &"Wrong rational: {a}/{b}" )
  if b < 0:
    b = -b
    a = -a
  if pa.p  < 2:
    raise newException(PadicError, &"Wrong value for p: {pa.p}")
  if pa.k < 1:
    raise newException(PadicError, &"Wrong value for k: {pa.k}")
  pa.p = min(pa.p, PMax)      # Maximum short prime.

  if sw: echo &"{a}/{b} + 0({pa.p}^{pa.k})"

  # Initialize.
  pa.v = 0
  pa.d.reset()
  if a == 0: return

  # Valuation.
  while b mod pa.p == 0:
    b = b div pa.p
    dec pa.v
  while a mod pa.p == 0:
    a = a div pa.p
    inc pa.v
  if (pa.v and 1) != 0:
    # Odd valuation.
    raise newException(PadicError, &"Non-residue mod {pa.p}.")

  # Maximum array length.
  pa.k = min(pa.k + pa.v, Emx - 1) - pa.v
  pa.v = pa.v shr 1

  if abs(a) > Amx or b > Amx:
    raise newException(PadicError, &"Rational exceeding limits: {a}/{b}.")

  if pa.p == 2:
    # 1 / b = b (mod 8); a / b = 1 (mod 8).
    if (a * b and 7) - 1 != 0:
      raise newException(PadicError, "Non-residue mod 8.")

    # Initialize.
    x = 1
    pa.d[pa.v] = 1
    pa.d[pa.v + 1] = 0
    var pk = 4
    i = pa.v + 2
    while i < pa.k + pa.v:
      pk *= 2
      let f = b * x * x - a
      let q = f div pk
      if f != q * pk: break   # Overflow.
      # Next digit.
      pa.d[i] = if (q and 1) != 0: 1 else: 0
      # Lift "x".
      x += pa.d[i] * (pk shr 1)
      inc i

  else:
    # Find root for small "p".
    var r = 1
    while r < pa.p:
      if (b * r * r - a) mod pa.p == 0: break
      inc r
    if r == pa.p:
      raise newException(PadicError, &"Non-residue mod {pa.p}.")
    let t = (b * r shl 1) mod pa.p
    var s = 0

    # Modular inverse for small "p".
    var f1 = 1
    while f1 < pa.p:
      inc s, t
      if s >= pa.p: dec s, pa.p
      if s == 1: break
      inc f1
    if f1 == pa.p:
      raise newException(PadicError, "Impossible to compute inverse modulo")

    f1 = pa.p - f1
    x = r
    pa.d[pa.v] = x

    var pk = 1
    i = pa.v + 1
    while i < pa.k + pa.v:
      pk *= pa.p
      let f = b * x * x - a
      let q = f div pk
      if f != q * pk: break   # Overflow.
      pa.d[i] = q * f1 mod pa.p
      if pa.d[i] < 0: pa.d[i] += pa.p
      x += pa.d[i] * pk
      inc i

  pa.k = i - pa.v
  if sw: echo &"lift: {x} mod {pa.p}^{pa.k}"


proc crat(pa: Padic; sw: bool): Ratio =
  ## Rational reconstruction.

  # Weighted digit sum.
  var
    s = 0
    pk = 1
  for i in min(pa.v, 0)..<(pa.k + pa.v):
    let pm = pk
    pk *= pa.p
    if pk div pm - pa.p != 0:
      # Overflow.
      pk = pm
      break
    s += pa.d[i] * pm

  # Lattice basis reduction.
  var
    m = [pk, s]
    n = [0, 1]
    i = 0
    j = 1
  s = s * s + 1
  # Lagrange's algorithm.
  while true:
    # Euclidean step.
    var q = ((m[i] * m[j] + n[i] * n[j]) / s).toInt
    m[i] -= q * m[j]
    n[i] -= q * n[j]
    q = s
    s = m[i] * m[i] + n[i] * n[i]
    # Compare norms.
    if s < q: swap i, j   # Interchange vectors.
    else: break

  var x = m[j]
  var y = n[j]
  if y < 0:
    y = -y
    x = -x

  # Check determinant.
  if abs(m[i] * y - x * n[i]) != pk:
    raise newException(PadicError, "Rational reconstruction failed.")

  # Negative powers.
  for i in pa.v..(-1): y *= pa.p

  if sw: echo x, if y > 1: '/' & $y else: ""
  result = (x, y)


func cmpt(pa: Padic): Padic =
  ## Return the complement.
  result = Padic(p: pa.p, k: pa.k, v: pa.v)
  var c = 1
  for i in pa.v..(pa.k + pa.v):
    inc c, pa.p - 1 - pa.d[i]
    if c >= pa.p:
      result.d[i] = c - pa.p
      c = 1
    else:
      result.d[i] = c
      c = 0


func sqr(pa: Padic): Padic =
  ## Return the square of a P-adic number.
  result = Padic(p: pa.p, k: pa.k, v: pa.v * 2)
  var c = 0
  for i in 0..pa.k:
    for j in 0..i:
      c += pa.d[pa.v + j] * pa.d[pa.v + i - j]
    # Euclidean step.
    let q = c div pa.p
    result.d[result.v + i] = c - q * pa.p
    c = q


func `$`(pa: Padic): string =
  ## String representation.
  let t = min(pa.v, 0)
  for i in countdown(pa.k - 1 + t, t):
    result.add $pa.d[i]
    if i == 0 and pa.v < 0: result.add "."
    result.add " "


when isMainModule:

  const Data = [[-7, 1, 2, 7],
                [9, 1, 2, 8],
                [17, 1, 2, 9],
                [497, 10496, 2, 18],
                [10496, 497, 2, 19],
                [3141, 5926, 3, 15],
                [2718,  281, 3, 13],
                [-1,  1,  5, 8],
                [86, 25,  5, 8],
                [2150, 1,  5, 8],
                [2,1, 7, 8],
                [-2645, 28518, 7, 9],
                [3029, 4821, 7, 9],
                [379, 449, 7, 8],
                [717, 8, 11, 7],
                [1414, 213, 41, 5],
                [-255, 256, 257, 3]]

  for d in Data:
    try:
      let q: Ratio = (d[0], d[1])
      var a = Padic(p: d[2], k: d[3])
      a.sqrt(q, true)
      echo "sqrt +/-"
      echo "...", a
      a = a.cmpt()
      echo "...", a
      let c = sqr(a)
      echo "sqrt^2"
      echo "   ", c
      let r = c.crat(true)
      if q.a * r.b - r.a * q.b != 0:
        echo "fail: sqrt^2"
      echo ""
    except PadicError:
      echo getCurrentExceptionMsg()
Output:
-7/1 + 0(2^7)
lift: 53 mod 2^7
sqrt +/-
...0 1 1 0 1 0 1 
...1 0 0 1 0 1 1 
sqrt^2
   1 1 1 1 0 0 1 
-7

9/1 + 0(2^8)
lift: 253 mod 2^8
sqrt +/-
...1 1 1 1 1 1 0 1 
...0 0 0 0 0 0 1 1 
sqrt^2
   0 0 0 0 1 0 0 1 
9

17/1 + 0(2^9)
lift: 233 mod 2^9
sqrt +/-
...0 1 1 1 0 1 0 0 1 
...1 0 0 0 1 0 1 1 1 
sqrt^2
   0 0 0 0 1 0 0 0 1 
17

497/10496 + 0(2^18)
lift: 92989 mod 2^18
sqrt +/-
...0 1 0 1 1 0 1 0 1 1 0 0 1 1. 1 1 0 1 
...1 0 1 0 0 1 0 1 0 0 1 1 0 0. 0 0 1 1 
sqrt^2
   1 0 0 0 0 0 1 1 0 0. 1 0 0 0 1 0 0 1 
497/10496

10496/497 + 0(2^19)
lift: 148501 mod 2^19
sqrt +/-
...1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 
...0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0 
sqrt^2
   0 0 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 
10496/497

3141/5926 + 0(3^15)
lift: 3406966 mod 3^15
sqrt +/-
...2 0 1 0 2 0 0 2 1 1 0 2 2 1 0 
...0 2 1 2 0 2 2 0 1 1 2 0 0 2 0 
sqrt^2
   2 1 1 1 2 2 0 0 0 2 0 1 1 0 0 
3141/5926

2718/281 + 0(3^13)
lift: 693355 mod 3^13
sqrt +/-
...0 2 2 0 2 0 0 0 2 2 1 1 0 
...2 0 0 2 0 2 2 2 0 0 1 2 0 
sqrt^2
   2 2 0 0 1 2 2 0 2 2 1 0 0 
2718/281

-1/1 + 0(5^8)
lift: 280182 mod 5^8
sqrt +/-
...3 2 4 3 1 2 1 2 
...1 2 0 1 3 2 3 3 
sqrt^2
   4 4 4 4 4 4 4 4 
-1

86/25 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 1 0 2 4 1 1. 1 
...3 3 4 2 0 3 3. 4 
sqrt^2
   0 0 0 0 0 3. 2 1 
86/25

2150/1 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 0 2 4 1 1 1 0 
...3 4 2 0 3 3 4 0 
sqrt^2
   0 0 0 3 2 1 0 0 
2150

2/1 + 0(7^8)
lift: 1802916 mod 7^8
sqrt +/-
...2 1 2 1 6 2 1 3 
...4 5 4 5 0 4 5 4 
sqrt^2
   0 0 0 0 0 0 0 2 
2

-2645/28518 + 0(7^9)
lift: 39082527 mod 7^9
sqrt +/-
...6 5 3 1 2 4 1 4. 1 
...0 1 3 5 4 2 5 2. 6 
sqrt^2
   1 1 3 3 4 4 5. 1 1 
-2645/28518

3029/4821 + 0(7^9)
lift: 31104424 mod 7^9
sqrt +/-
...5 2 5 2 4 5 3 1 1 
...1 4 1 4 2 1 3 5 6 
sqrt^2
   6 5 6 4 1 3 0 2 1 
3029/4821

379/449 + 0(7^8)
lift: 555087 mod 7^8
sqrt +/-
...0 4 5 0 1 2 2 1 
...6 2 1 6 5 4 4 6 
sqrt^2
   5 3 1 2 4 1 4 1 
379/449

717/8 + 0(11^7)
lift: 5280093 mod 11^7
sqrt +/-
...2 10 8 7 0 1 5 
...8 0 2 3 10 9 6 
sqrt^2
   1 4 1 4 2 1 3 
717/8

1414/213 + 0(41^5)
lift: 25564902 mod 41^5
sqrt +/-
...9 1 38 6 8 
...31 39 2 34 33 
sqrt^2
   12 36 31 15 23 
1414/213

-255/256 + 0(257^3)
lift: 8867596 mod 257^3
sqrt +/-
...134 66 68 
...122 190 189 
sqrt^2
   255 255 255 
-255/256

Phix

Translation of: FreeBASIC
constant EMX  = 48       // exponent maximum (if indexing starts at -EMX)
constant DMX  = 1e5     // approximation loop maximum
constant AMX  = 700000 // argument maximum
constant PMAX = 32749   // prime maximum
 
// global variables
integer p1 = 0
integer p  = 7    // default prime
integer k  = 11   // precision
 
type Ratio(sequence r)
    return length(r)=2 and integer(r[1]) and integer(r[2])
end type
 
class Padic
    integer v = 0   
    sequence d = repeat(0,EMX*2)
 
    function square_root(Ratio g, integer sw)
        -- p-adic square root of g = a/b
        integer {a,b} = g
        atom f, q, pk, x
        integer f1, r, s, t, i, res = 0
 
        if b = 0 then return 1 end if
        if b < 0 then
            b = -b
            a = -a
        end if
        if p < 2 or k < 1 then return 1 end if
 
        -- max. short prime
        p = min(p, PMAX)
        if sw then
            -- numerator, denominator, prime, precision
            printf(1,"%d/%d + O(%d^%d)\n",{a,b,p,k})
        end if

        -- initialize
        v = 0
        p1 = p - 1
        sequence ntd = repeat(0,2*EMX) -- (new this.d)
        if a = 0 then return 0 end if
 
        -- valuation
        while remainder(b,p)=0 do
            b /= p
            v -= 1
        end while
        while remainder(a,p)=0 do
            a /= p
            v += 1
        end while
 
        if remainder(v,2) then
            -- odd valuation
            printf(1,"(1)non-residue mod %d\n",p)
            return -1
        end if

        -- max. array length
        k = min(k + v, EMX - 1) - v
        v = floor(v/2)
 
        if abs(a) > AMX or b > AMX then return -1 end if
 
        if p = 2 then
            --1 / b = b (mod 8)
            --a / b = 1 (mod 8)
            t = a * b
            if mod(t,8)-1 then
                printf(1,"(2)non-residue mod 8\n")
                return -1
             end if
 
        else
            -- find root for small p
            for r = 1 to p1 do
                f = b * r * r - a
                if mod(f,p) = 0 then exit end if
            end for
 
            if r = p then
                printf(1,"(3)non-residue mod %d\n", p)
                return -1
            end if
 
            -- f'(r) = 2br
            t = b * r * 2
 
            s = 0
            t = mod(t,p)
            -- modular inverse for small p
            for f1 = 1 to p1 do
                s += t
                if s > p1 then s -= p end if
                if s = 1 then exit end if
            end for
 
            if f1 = p then
                printf(1,"impossible inverse mod\n")
                return -1
            end if
        end if
 
        if p = 2 then
            -- initialize
            x = 1
            ntd[v+EMX+1] = 1
            ntd[v+EMX+2] = 0
 
            pk = 4
            for i = v+2 to k-1+v do
                pk *= 2
                f = b * x * x - a
                q = floor(f/pk)
                -- overflow
                if f != q * pk then exit end if
                -- next digit
                ntd[i+EMX+1] = and_bits(q,1)
                -- lift x
                x += ntd[i+EMX+1] * floor(pk/2)
            end for
 
        else
            -- -1 / f'(x) mod p
            f1 = p - f1
            x = r
            ntd[v+EMX+1] = x
 
            pk = 1
            for i = v+1 to k-1 do
                pk *= p
                f = b * x * x - a
                q = floor(f/pk)
                -- overflow
                if f - q * pk then exit end if
                r = mod(q*f1,p)
                if r < 0 then r += p end if
                ntd[i+EMX+1] = r
                x += r * pk
            end for
        end if
        this.d = ntd
        k = i-v
 
        if sw then
            printf(1,"lift: %d mod %d^%d\n",{x,p,k})
        end if  
        return 0
    end function

    function square()
        integer c = 0
        Padic r = new()
        r.v = this.v * 2
        sequence td = this.d,
                 rd = r.d
        for i=0 to k do
            for j=0 to i do
                c += td[v+j+EMX+1] * td[v+i-j+EMX+1]
            end for
            // Euclidean step
            integer q = floor(c/p)
            rd[r.v+i+EMX+1] = c - q*p
            c = q
        end for
        r.d = rd
        return r
    end function

    function complement()
        integer c = 1
        Padic r = new({v})
        sequence rd = r.d
        for i=v to k+v do
            integer dx = i+EMX+1
            c += p1 - this.d[dx]
            if c>p1 then
                rd[dx] = c - p
                c = 1
            else
                rd[dx] = c
                c = 0
            end if
        end for
        r.d = rd
        return r
    end function
 
    function crat(integer sw)
        -- rational reconstruction
        integer i, j, t = min(v, 0)
        Ratio r
        atom f
        integer x, y
        atom p1,pk, q, s
 
        -- weighted digit sum
        s = 0
        pk = 1
        for i = t to k-1+v do
            p1 = pk
            pk *= p
            if floor(pk/p1) - p then
                -- overflow
                pk = p1
                exit
            end if
            s += d[i+EMX+1] * p1 --(mod pk)
        end for
 
        -- lattice basis reduction
        sequence m = {pk, s},
                 n = {0, 1}
        i = 1
        j = 2
        s = s * s + 1 -- norm(v)^2
 
        -- Lagrange's algorithm
        while true do
            f = (m[i] * m[j] + n[i] * n[j]) / s
 
            -- Euclidean step
            q = floor(f +.5)
            m[i] -= q * m[j]
            n[i] -= q * n[j]
 
            q = s
            s = m[i] * m[i] + n[i] * n[i]
            -- compare norms
            if s < q then
                -- interchange vectors
                {i,j} = {j,i}
            else
                exit
            end if
        end while
 
        x = m[j]
        y = n[j]
        if y < 0 then
            y = -y
            x = -x
        end if

        -- check determinant
        t = abs(m[i] * y - x * n[i]) == pk
 
        if not t then
            printf(1,"crat: fail\n")
            x = 0
            y = 1
        else
            -- negative powers
            for i = v to -1 do
                y *= p
            end for
 
            if sw then
--              printf(1,iff(y=1?"%d":"%d/%d"),{x*sgn,y})
                printf(1,iff(y=1?"%d\n":"%d/%d\n"),{x,y})
            end if
        end if
 
        r = {x,y}
        return r
    end function

    procedure prntf(bool sw)
        -- print expansion
        integer t = min(v, 0)
        for i=k-1+t to t by -1 do
            printf(1,"%d",d[i+EMX+1])
            printf(1,iff(i=0 and v<0?". ":" "))
        end for
        printf(1,"\n")
        // rational approximation
        if sw then crat(sw) end if
    end procedure
end class
 
constant tests = {
    {{-7,1},2,7},
--/*
    {{9,1},2,8},
    {{17,1},2,9},
    {{497,10496},2,18},
    {{10496,497},2,19}, 

    {{3141,5926},3,17},
    {{2718,281},3,15},

    {{-1,1},5,8},
    {{86,25},5,8},
    {{2150,1},5,10},

    {{2,1},7,8},
    {{-2645,28518},7,9},
    {{3029,4821},7,9},
    {{379,449},7,8},

    {{717,8},11,7},
    {{1414,213},41,5},
--*/
    {{-255,256},257,3}
}
 
Padic a = new(), c
Ratio q, r

for i=1 to length(tests) do
   {q,p,k} = tests[i]
 
   integer sw = a.square_root(q, 1)
   if sw=1 then exit end if
   if sw=0 then
       printf(1,"square_root +/-\n")
       printf(1,"... ")
       a.prntf(0)
       a = a.complement()
       printf(1,"... ")
       a.prntf(0)
 
       c = a.square()
       printf(1,"square_root^2\n")
       printf(1,"    ")
       c.prntf(0)
       r = c.crat(1)
 
       if q[1] * r[2] - r[1] * q[2] then
          printf(1,"fail: square_root^2\n")
       end if
   end if
   printf(1,"\n")
end for
Output:
-7/1 + O(2^7)
lift: 53 mod 2^7
square_root +/-
... 0 1 1 0 1 0 1
... 1 0 0 1 0 1 1
square_root^2
    1 1 1 1 0 0 1
-7

-255/256 + O(257^3)
lift: 8867596 mod 257^3
square_root +/-
... 134 66 68
... 122 190 189
square_root^2
    255 255 255
-255/256

Wren

Translation of: FreeBASIC
Library: Wren-dynamic
Library: Wren-big
import "./dynamic" for Struct
import "./big" for BigInt

// constants
var EMX  = 64       // exponent maximum (if indexing starts at -EMX)
var AMX  = 6000     // argument maximum
var PMAX = 32749    // prime maximum

// global variables
var P1 = 0
var P  = 7    // default prime
var K  = 11   // precision

var Ratio = Struct.create("Ratio", ["a", "b"])

class Padic {
    // uninitialized
    construct new() {
        _v = 0
        _d = List.filled(2 * EMX, 0) // add EMX to index to be consistent wih FB
    }

    // properties
    v { _v }
    v=(o) { _v = o }
    d { _d }

    // (re)initialize 'this' to the square root of a Ratio, set 'sw' to print
    sqrt(g, sw) {
        var a = g.a
        var b = g.b
        if (b == 0) return 1
        if (b < 0) {
            b = -b
            a = -a
        }
        if (P < 2 || K < 1) return 1
        P = P.min(PMAX)  // maximum short prime
        if (sw != 0) {
            System.write("%(a)/%(b) + ")  // numerator, denominator
            System.print("0(%(P)^%(K))")  // prime, precision
        }

        // (re)initialize
        _v = 0
        P1 = P - 1
        _d = List.filled(2 * EMX, 0)
        if (a == 0) return 0

        //valuation
        while (b%P== 0) {
            b = (b/P).truncate
            _v = _v - 1
        }

        while (a%P == 0) {
            a = (a/P).truncate
            _v = _v + 1
        }

        if ((_v & 1) == 1) {
            // odd valuation
            System.print("non-residue mod %(P)")     
            return -1
        }
        K = (K + _v).min(EMX - 1) - _v  // maximum array length
        _v = (_v/2).truncate

        if (a.abs > AMX || b > AMX) return -1
        var bb = BigInt.new(b) // to avoid overflowing 'f(x) = b * x * x – a'
        var r
        var s
        var t
        var f
        var f1
        if (P == 2) {
            t = a * b
            if ((t & 7) - 1 != 0) {
                System.print("non-residue mod 8")
                return -1
            }
        } else {
            // find root for small P
            r = 1            
            while (r <= P1) {
                f = bb * r * r - a
                if ((f % P) == 0) break
                r = r + 1
            }
            if (r == P) {
                System.print("non-residue mod %(P)")
                return -1
            }
            t = 2 * b * r
            s = 0
            t = t % P

            // modular inverse for small P
            f1 = 1
            while (f1 <= P1) {
                s = s + t
                if (s > P1) s = s - P
                if (s == 1) break
                f1 = f1 + 1
            }
            if (f1 == P) {
                System.print("impossible inverse mod")
                return -1
            }
        }
        var x
        var pk
        var q
        var i
        if (P == 2) {
            // initialize
            x = 1
            _d[_v+EMX] = 1
            _d[_v+1+EMX] = 0
            pk = 4
            i = _v + 2
            while (i <= K - 1 + _v) {
                pk = pk * 2
                f = bb * x * x - a
                q = f / pk
                // overflow
                if (f != q * pk) break
                // next digit
                _d[i+EMX] = ((q & 1) != 0) ? 1 : 0
                // lift x
                x = x + _d[i+EMX]*(pk >> 1)
                i = i + 1
            }

        } else {
            f1 = P - f1
            x = r
            _d[_v+EMX] = x
            pk = 1
            i = _v + 1
            while (i <= K - 1 + _v) {
                pk = pk * P
                f = bb * x * x - a
                q = f / pk
                // overflow
                if (f != q * pk) break
                _d[i+EMX] = q.toSmall * f1 % P
                if (_d[i+EMX] < 0) _d[i+EMX] = _d[i+EMX] + P
                x = x + _d[i+EMX]*pk
                i = i + 1
            }
        }
        K = i - _v
        if (sw != 0) System.print("lift: %(x) mod %(P)^%(K)")
        return 0
    }

    // rational reconstruction
    crat(sw) {
        var t = _v.min(0)
        // weighted digit sum
        var s = 0
        var pk = 1
        for (i in t..K-1+_v) {
            P1 = pk
            pk = pk * P
            if (((pk/P1).truncate - P) != 0) {
                // overflow
                pk = p1
                break
            }
            s = s + _d[i+EMX]*P1
        }

        // lattice basis reduction
        var m = [pk, s]
        var n = [0, 1]
        var i = 0
        var j = 1
        s = s * s + 1
        // Lagrange's algorithm
        while (true) {
            var f = (m[i] * m[j] + n[i] * n[j]) / s
            // Euclidean step
            var q = (f + 0.5).floor
            m[i] = m[i] - q*m[j]
            n[i] = n[i] - q*n[j]
            q = s
            s = m[i] * m[i] + n[i] * n[i]
            // compare norms
            if (s < q) {
                // interchange vectors
                var z = i
                i = j
                j = z
            } else {
                break
            }
        }
        var x = m[j]
        var y = n[j]
        if (y < 0) {
            y = -y
            x = -x
        }

        // check determinant
        t = (m[i]*y - x*n[i]).abs == pk
        if (!t) {
            System.print("crat: fail")
            x = 0
            y = 1
        } else {
            // negative powers
            var i = _v
            while (i <= -1) {
                y = y * P
                i = i + 1
            }
            if (sw != 0) {
                System.write(x)
                if (y > 1) System.write("/%(y)")
                System.print()
            }
        }

        return Ratio.new(x, y)
    }

    // print expansion
    printf(sw) {
        var t = _v.min(0)
        for (i in K - 1 + t..t) {
            System.write(_d[i + EMX])
            if (i == 0 && _v < 0) System.write(".")
            System.write(" ")
        }
        System.print()
        // rational approximation
        if (sw != 0) crat(sw)
    }

    // complement
    cmpt {
        var c = 1
        var r = Padic.new()
        r.v = _v
        for (i in r.v..K + r.v) {
            c = c + P1 - _d[i+EMX]
            if (c > P1) {
                r.d[i+EMX] = c - P
                c = 1
            } else {
                r.d[i+EMX] = c
                c = 0
            }
        }
        return r
    }

    // square
    sqr {
        var c = 0
        var r = Padic.new()
        r.v = _v * 2
        for (i in 0..K) {
            for (j in 0..i) c = c + _d[_v+j+EMX] * _d[_v+i-j+EMX]
            // Euclidean step
            var q = (c/P).truncate
            r.d[r.v+i+EMX] = c - q*P
            c = q
        }
        return r
    }
}

var data = [
    [-7, 1, 2, 7],
    [9, 1, 2, 8],
    [17, 1, 2, 9],
    [497, 10496, 2, 18],
    [10496, 497, 2, 19], 
    [3141, 5926, 3, 15],
    [2718,  281, 3, 13],
    [-1,  1,  5, 8],
    [86, 25,  5, 8],
    [2150, 1,  5, 8],
    [2,1, 7, 8],
    [-2645, 28518, 7, 9],
    [3029, 4821, 7, 9],
    [379, 449, 7, 8],
    [717, 8, 11, 7],
    [1414, 213, 41, 5],
    [-255, 256, 257, 3]
]

var sw = 0
var a = Padic.new()
var c = Padic.new()

for (d in data) {
    var q = Ratio.new(d[0], d[1])
    P = d[2]
    K = d[3]
    sw = a.sqrt(q, 1)
    if (sw == 1) break
    if (sw == 0) {
        System.print("sqrt +/-")
        System.write("...")
        a.printf(0)
        a = a.cmpt
        System.write("...")
        a.printf(0)
        c = a.sqr
        System.print("sqrt^2")
        System.write("   ")
        c.printf(0)
        var r = c.crat(1)
        if (q.a * r.b - r.a * q.b != 0) {
            System.print("fail: sqrt^2")
        }
        System.print()
    }
}
Output:
-7/1 + 0(2^7)
lift: 53 mod 2^7
sqrt +/-
...0 1 1 0 1 0 1 
...1 0 0 1 0 1 1 
sqrt^2
   1 1 1 1 0 0 1 
-7

9/1 + 0(2^8)
lift: 253 mod 2^8
sqrt +/-
...1 1 1 1 1 1 0 1 
...0 0 0 0 0 0 1 1 
sqrt^2
   0 0 0 0 1 0 0 1 
9

17/1 + 0(2^9)
lift: 233 mod 2^9
sqrt +/-
...0 1 1 1 0 1 0 0 1 
...1 0 0 0 1 0 1 1 1 
sqrt^2
   0 0 0 0 1 0 0 0 1 
17

497/10496 + 0(2^18)
lift: 92989 mod 2^18
sqrt +/-
...0 1 0 1 1 0 1 0 1 1 0 0 1 1. 1 1 0 1 
...1 0 1 0 0 1 0 1 0 0 1 1 0 0. 0 0 1 1 
sqrt^2
   1 0 0 0 0 0 1 1 0 0. 1 0 0 0 1 0 0 1 
497/10496

10496/497 + 0(2^19)
lift: 148501 mod 2^19
sqrt +/-
...1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 
...0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0 
sqrt^2
   0 0 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 
10496/497

3141/5926 + 0(3^15)
lift: 3406966 mod 3^15
sqrt +/-
...2 0 1 0 2 0 0 2 1 1 0 2 2 1 0 
...0 2 1 2 0 2 2 0 1 1 2 0 0 2 0 
sqrt^2
   2 1 1 1 2 2 0 0 0 2 0 1 1 0 0 
3141/5926

2718/281 + 0(3^13)
lift: 693355 mod 3^13
sqrt +/-
...0 2 2 0 2 0 0 0 2 2 1 1 0 
...2 0 0 2 0 2 2 2 0 0 1 2 0 
sqrt^2
   2 2 0 0 1 2 2 0 2 2 1 0 0 
2718/281

-1/1 + 0(5^8)
lift: 280182 mod 5^8
sqrt +/-
...3 2 4 3 1 2 1 2 
...1 2 0 1 3 2 3 3 
sqrt^2
   4 4 4 4 4 4 4 4 
-1

86/25 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 1 0 2 4 1 1. 1 
...3 3 4 2 0 3 3. 4 
sqrt^2
   0 0 0 0 0 3. 2 1 
86/25

2150/1 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 0 2 4 1 1 1 0 
...3 4 2 0 3 3 4 0 
sqrt^2
   0 0 0 3 2 1 0 0 
2150

2/1 + 0(7^8)
lift: 1802916 mod 7^8
sqrt +/-
...2 1 2 1 6 2 1 3 
...4 5 4 5 0 4 5 4 
sqrt^2
   0 0 0 0 0 0 0 2 
2

-2645/28518 + 0(7^9)
lift: 39082527 mod 7^9
sqrt +/-
...6 5 3 1 2 4 1 4. 1 
...0 1 3 5 4 2 5 2. 6 
sqrt^2
   1 1 3 3 4 4 5. 1 1 
-2645/28518

3029/4821 + 0(7^9)
lift: 31104424 mod 7^9
sqrt +/-
...5 2 5 2 4 5 3 1 1 
...1 4 1 4 2 1 3 5 6 
sqrt^2
   6 5 6 4 1 3 0 2 1 
3029/4821

379/449 + 0(7^8)
lift: 555087 mod 7^8
sqrt +/-
...0 4 5 0 1 2 2 1 
...6 2 1 6 5 4 4 6 
sqrt^2
   5 3 1 2 4 1 4 1 
379/449

717/8 + 0(11^7)
lift: 5280093 mod 11^7
sqrt +/-
...2 10 8 7 0 1 5 
...8 0 2 3 10 9 6 
sqrt^2
   1 4 1 4 2 1 3 
717/8

1414/213 + 0(41^5)
lift: 25564902 mod 41^5
sqrt +/-
...9 1 38 6 8 
...31 39 2 34 33 
sqrt^2
   12 36 31 15 23 
1414/213

-255/256 + 0(257^3)
lift: 8867596 mod 257^3
sqrt +/-
...134 66 68 
...122 190 189 
sqrt^2
   255 255 255 
-255/256