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 = 64 'exponent maximum

const amx = 6000 'tentative argument maximum

'------------------------------------------------ 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 f, q, x, pk dim as long f1, r, s, t dim i as integer 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)
  'max. array length
  k = min(k, emx - 1)

  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
  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
        f = b * r * r - a
        if f 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 * x - a
     q = f \ pk
     'overflow
     if f - q * pk then exit for
  #endmacro

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

     pk = 4
     for i = v + 2 to k - 1
        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
        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

  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 r as ratio, f as double dim as long x, y, pk, p1 dim as longint q, s

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

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

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

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

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

     'Euclidean step
     q = int(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
        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,27

data 3141,5926, 3,17 data 2718,281, 3,15

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

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

data 717,8, 11,7 data 1414,213, 41,5 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

system </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: 355133 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^27)
lift: 3818517 mod 2^27
sqrt +/-
... 0 1 1 1 0 1 0 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0
... 1 0 0 0 1 0 1 1 0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0
sqrt^2
    1 1 1 0 0 1 0 1 0 0 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0
 10496/497


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


 2718/281 + O(3^15)
lift: 3882001 mod 3^15
sqrt +/-
... 2 1 0 2 2 0 2 0 0 0 2 2 1 1 0
... 0 1 2 0 0 2 0 2 2 2 0 0 1 2 0
sqrt^2
    0 0 2 2 0 0 1 2 2 0 2 2 1 0 0
 2718/281


-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: 486156 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^10)
lift: 486156 mod 5^10
sqrt +/-
... 1 1 1 0 2 4 1 1 1 0
... 3 3 3 4 2 0 3 3 4 0
sqrt^2
    0 0 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


-2645/28518 + O(7^9)
lift: 281204169 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 + O(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 + O(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 + O(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 + O(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 + O(257^3)
lift: 8867596 mod 257^3
sqrt +/-
... 134 66 68
... 122 190 189
sqrt^2
    255 255 255
-255/256