P-Adic square roots: Difference between revisions

=={{header|Nim}}==

{{trans|FreeBASIC}}

{{trans|Wren}}

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

{{out}}

<pre>-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</pre>