P-Adic square roots

P-Adic square roots is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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 x₁ modulo p and build a sequence of solutions f(x_k) ≡ 0 (mod p^k),
using the lifting equation x_k+1 = x_k + d_k * p^k  with d_k = –(f(x_k) / p^k) / f ′(x₁) (mod p).
The multipliers d_k are the successive p-adic digits to find.

If evaluation of f(x) = bx² – a 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 x² ≡ a (mod n)

FreeBASIC

<lang 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

1. 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 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
        pk *= p
        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 </lang>

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

Wren

Translation of: FreeBASIC

Library: Wren-dynamic

Library: Wren-math

Library: Wren-big

<lang ecmascript>import "/dynamic" for Struct import "/math" for Math 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 = Math.min(P, 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 = Math.min(K + _v, 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 = Math.min(_v, 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 = Math.min(_v, 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()
   }

}</lang>

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