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) = bx2 – 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 x2 ≡ 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, 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 </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

Julia

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

</lang>

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

Phix

Translation of: FreeBASIC

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

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

<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
Retrieved from "https://rosettacode.org/wiki/P-Adic_square_roots?oldid=307909"