P-Adic square roots: Difference between revisions
Added Wren
mNo edit summary |
(Added Wren) |
||
Line 596:
-255/256
</pre>
=={{header|Wren}}==
{{trans|FreeBASIC}}
{{libheader|Wren-dynamic}}
{{libheader|Wren-math}}
{{libheader|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
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>
{{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: 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
</pre>
| |||