P-Adic square roots
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.
- Task.
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.
- 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 = 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
- 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)
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
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 + 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 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 + v
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: 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^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: 46453687 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: 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^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: 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 + 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
Wren
<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
K = Math.min(K, EMX-1) // maximum array length
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
}
_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) {
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) {
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
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+t) {
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, 27], [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]
]
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: 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 + 0(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 + 0(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 + 0(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 + 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: 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 + 0(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 + 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: 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 + 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