I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

You are encouraged to solve this task according to the task description, using any language you may know.

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.

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, bend type type padicdeclare function sqrt (byref q as ratio, byval sw as integer) as integer'p-adic square root of q = a/b, set sw to printdeclare sub printf (byval sw as integer)'print expansion, set sw to print rationaldeclare function crat (byval sw as integer) as ratio'rational reconstruction declare sub cmpt (byref a as padic)'let self:= complement_adeclare sub sqr (byref a as padic)'let self:= a ^ 2    as long d(-emx to emx - 1)   as integer vend type  'global variablesdim shared as long p1, p = 7'default primedim shared as integer k = 11'precision #define min(a, b) iif((a) > (b), b, a)  '------------------------------------------------'p-adic square root of g = a/bfunction padic.sqrt (byref g as ratio, byval sw as integer) as integerdim as longint a = g.a, b = g.bdim as longint q, x, pk, pmdim as long f1, r, s, tdim i as integer, f as doublesqrt = 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 reconstructionfunction padic.crat (byval sw as integer) as ratiodim as integer i, j, t = min(v, 0)dim as longint s, pk, pmdim as long q, x, ydim as double f, hdim 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 = yreturn rend function  'print expansionsub 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_asub padic.cmpt (byref a as padic)dim i as integer, r as padicdim as long c = 1with 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 iend withthis = rend sub 'let self:= a ^ 2sub padic.sqr (byref a as padic)dim as long ptr rp, ap = @a.d(a.v)dim as longint q, c = 0dim as integer i, jdim r as padicwith 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 iend withthis = rend sub  'main'------------------------------------------------dim as integer swdim as padic a, cdim as ratio q, r width 64, 30cls '   -7 + O(2^7)data -7,1, 2,7data 9,1, 2,8data 17,1, 2,9data 497,10496, 2,18data 10496,497, 2,19 data -577215,664901, 3,23data 15403,26685, 3,18 data -1,1, 5,8data 86,25, 5,8data 2150,1, 5,8 data 2,1, 7,8data 11696,621467, 7,11data -27764,11521, 7,11data -27584,12953, 7,11 data -166420,135131, 11,11data 14142,135623, 5,15data -255,256, 257,3 data 0,0, 0,0  printdo   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

```

## 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    endend 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 maximumconstant AMX  = 700000 // argument maximumconstant PMAX = 32749   // prime maximum // global variablesinteger p1 = 0integer p  = 7    // default primeinteger 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 procedureend 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(), cRatio 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-math
Library: Wren-big
`import "/dynamic" for Structimport "/math" for Mathimport "/big" for BigInt // constantsvar EMX  = 64       // exponent maximum (if indexing starts at -EMX)var AMX  = 6000     // argument maximumvar PMAX = 32749    // prime maximum // global variablesvar P1 = 0var P  = 7    // default primevar 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 = 0var 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
```