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.
- 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
- 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
Haskell
Uses solutin for P-Adic_numbers,_basic#Haskell
<lang haskell>{-# LANGUAGE KindSignatures, DataKinds #-}
import Data.Ratio import Data.List (find) import GHC.TypeLits import Padic
pSqrt :: KnownNat p => Rational -> Padic p pSqrt r = res
where res = maybe Null mkUnit series (a, b) = (numerator r, denominator r) series = case modulo res of 2 | eqMod 4 a 3 -> Nothing | not (eqMod 8 a 1) -> Nothing | otherwise -> Just $ 1 : 0 : go 8 1 where go pk x = let q = ((b*x*x - a) `div` pk) `mod` 2 in q : go (2*pk) (x + q * (pk `div` 2))
p -> do y <- find (\x -> eqMod p (b*x*x) a) [1..p-1] df <- recipMod p (2*b*y) let go pk x = let f = (b*x*x - a) `div` pk d = (f * (p - df)) `mod` p in x `div` (pk `div` p) : go (p*pk) (x + d*pk) Just $ go p y
eqMod :: Integral a => a -> a -> a -> Bool eqMod p a b = a `mod` p == b `mod` p</lang>
λ> :set -XDataKinds λ> pSqrt 2 :: Padic 7 7-adic: 12011266421216213.0 λ> pSqrt 39483088 :: Padic 7 7-adic: 22666526525245311.0 λ> pSqrt 5 :: Padic 7 Null λ> pSqrt 9 :: Padic 2 2-adic: 11111111111111101.0 λ> it ^ 2 2-adic: 1001.0 λ> toRational it 9 % 1 λ> pSqrt 17 ^ 2 == (17 :: Padic 2) True λ> pSqrt 50 ^ 2 == (50 :: Padic 7) True λ> pSqrt 52 ^ 2 == (52 :: Padic 7) False λ> pSqrt 52 :: Padic 7 Null
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
Nim
<lang Nim>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()</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
Phix
<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
<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