P-Adic numbers, basic: Difference between revisions

Content added Content deleted
mNo edit summary
(Added Wren)
Line 601: Line 601:
577215/6649
577215/6649


</pre>

=={{header|Wren}}==
{{trans|FreeBASIC}}
{{libheader|Wren-dynamic}}
{{libheader|Wren-math}}
<lang ecmascript>import "/dynamic" for Struct
import "/math" for Math

// constants
var EMX = 64 // exponent maximum (if indexing starts at -EMX)
var DMX = 1e5 // approximation loop maximum
var AMX = 1048576 // 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' from Ratio, set 'sw' to print
r2pa(q, sw) {
var a = q.a
var b = q.b
if (b == 0) return 1
if (b < 0) {
b = -b
a = -a
}
if (a.abs > AMX || b > AMX) return -1
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
var i = 0
// find -exponent of P in b
while (b%P == 0) {
b = (b/P).truncate
i = i - 1
}
var s = 0
var r = b % P

// modular inverse for small P
var b1 = 1
while (b1 <= P1) {
s = s + r
if (s > P1) s = s - P
if (s == 1) break
b1 = b1 + 1
}
if (b1 == P) {
System.print("r2pa: impossible inverse mod")
return -1
}
_v = EMX
while (true) {
// find exponent of P in a
while (a%P == 0) {
a = (a/P).truncate
i = i + 1
}

// valuation
if (_v == EMX) _v = i

// upper bound
if (i >= EMX) break

// check precision
if ((i - _v) > K) break

// next digit
_d[i+EMX] = a * b1 % P
if (_d[i+EMX] < 0) _d[i+EMX] = _d[i+EMX] + P

// remainder - digit * divisor
a = a - _d[i+EMX]*b
if (a == 0) break
}
return 0
}

// Horner's rule
dsum() {
var t = Math.min(_v, 0)
var s = 0
for (i in K - 1 + t..t) {
var r = s
s = s * P
if (r != 0 && ((s/r).truncate - P) != 0) {
// overflow
s = -1
break
}
s = s + _d[i+EMX]
}
return s
}

// rational reconstruction
crat() {
var fl = false
var s = this
var j = 0
var i = 1

// denominator count
while (i <= DMX) {
// check for integer
j = K - 1 + _v
while (j >= _v) {
if (s.d[j+EMX] != 0) break
j = j - 1
}
fl = ((j - _v) * 2) < K
if (fl) {
fl = false
break
}

// check negative integer
j = K - 1 + _v
while (j >= _v) {
if (P1 - s.d[j+EMX] != 0) break
j = j - 1
}
fl = ((j - _v) * 2) < K
if (fl) break

// repeatedly add self to s
s = s + this
i = i + 1
}
if (fl) s = s.cmpt

// numerator: weighted digit sum
var x = s.dsum()
var y = i
if (x < 0 || y > DMX) {
System.print("crat: fail")
} else {
// negative powers
i = _v
while (i <= -1) {
y = y * P
i = i + 1
}

// negative rational
if (fl) x = -x
System.write(x)
if (y > 1) System.write("/%(y)")
System.print()
}
}

// 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()
}

// add b to 'this'
+(b) {
var c = 0
var r = Padic.new()
r.v = Math.min(_v, b.v)
for (i in r.v..K + r.v) {
c = c + _d[i+EMX] + b.d[i+EMX]
if (c > P1) {
r.d[i+EMX] = c - P
c = 1
} else {
r.d[i+EMX] = c
c = 0
}
}
return r
}

// complement
cmpt {
var c = 1
var r = Padic.new()
r.v = _v
for (i in _v..K + _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
}
}

var data = [
/* rational reconstruction limits are relative to the precision */
[2, 1, 2, 4, 1, 1],
[4, 1, 2, 4, 3, 1],
[4, 1, 2, 5, 3, 1],
[4, 9, 5, 4, 8, 9],
[-7, 5, 7, 4, 99, 70],
[26, 25, 5, 4, -109, 125],
[49, 2, 7, 6, -4851, 2],
[-9, 5, 3, 8, 27, 7],
[5, 19, 2, 12, -101, 384],
/* three decadic pairs */
[6, 7, 10, 7, -5, 7],
[2, 7, 10, 7, -3, 7],
[34, 21, 10, 9, -39034, 791],
/*familiar digits*/
[11, 4, 2, 43, 679001, 207],
[11, 4, 3, 27, 679001, 207],
[11, 4, 11, 13, 679001, 207],
[-22, 7, 2, 37, 46071, 379],
[-22, 7, 3, 23, 46071, 379],
[-22, 7, 7, 13, 46071, 379],
[-101, 109, 2, 40, 583376, 6649],
[-101, 109, 61, 7, 583376, 6649],
[-101, 109, 32749, 3, 583376, 6649]
]

var sw = 0
var a = Padic.new()
var b = Padic.new()

for (d in data) {
var q = Ratio.new(d[0], d[1])
P = d[2]
K = d[3]
sw = a.r2pa(q, 1)
if (sw == 1) break
a.printf(0)
q.a = d[4]
q.b = d[5]
sw = sw | b.r2pa(q, 1)
if (sw == 1) break
if (sw == 0) {
b.printf(0)
var c = a + b
System.print("+ =")
c.printf(1)
}
System.print()
}</lang>

{{out}}
<pre>
2/1 + 0(2^4)
0 0 1 0
1/1 + 0(2^4)
0 0 0 1
+ =
0 0 1 1
3

4/1 + 0(2^4)
0 1 0 0
3/1 + 0(2^4)
0 0 1 1
+ =
0 1 1 1
-2/2

4/1 + 0(2^5)
0 0 1 0 0
3/1 + 0(2^5)
0 0 0 1 1
+ =
0 0 1 1 1
7

4/9 + 0(5^4)
4 2 1 1
8/9 + 0(5^4)
3 4 2 2
+ =
3 1 3 3
4/3

-7/5 + 0(7^4)
2 5 4 0
99/70 + 0(7^4)
0 5 0. 5
+ =
6 2 0. 5
1/70

26/25 + 0(5^4)
0 1. 0 1
-109/125 + 0(5^4)
4. 0 3 1
+ =
0. 0 4 1
21/125

49/2 + 0(7^6)
3 3 3 4 0 0
-4851/2 + 0(7^6)
3 2 3 3 0 0
+ =
6 6 0 0 0 0
-2401

-9/5 + 0(3^8)
2 1 0 1 2 1 0 0
27/7 + 0(3^8)
1 2 0 1 1 0 0 0
+ =
1 0 1 0 0 1 0 0
72/35

5/19 + 0(2^12)
0 0 1 0 1 0 0 0 0 1 1 1
-101/384 + 0(2^12)
1 0 1 0 1. 0 0 0 1 0 0 1
+ =
1 1 1 0 0. 0 0 0 1 0 0 1
1/7296

6/7 + 0(10^7)
7 1 4 2 8 5 8
-5/7 + 0(10^7)
5 7 1 4 2 8 5
+ =
2 8 5 7 1 4 3
1/7

2/7 + 0(10^7)
5 7 1 4 2 8 6
-3/7 + 0(10^7)
1 4 2 8 5 7 1
+ =
7 1 4 2 8 5 7
-1/7

34/21 + 0(10^9)
9 5 2 3 8 0 9 5 4
-39034/791 + 0(10^9)
1 3 9 0 6 4 4 2 6
+ =
0 9 1 4 4 5 3 8 0
-16180/339

11/4 + 0(2^43)
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0. 1 1
679001/207 + 0(2^43)
0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 1 1 1 0 0 0 0 0 1 0 1 0 0 1 0 1 0 1 1 1
+ =
0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 1 1 1 0 0 0 0 0 1 0 1 0 0 1 0 1 1 0 0 1. 1 1
2718281/828

11/4 + 0(3^27)
2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 1 2
679001/207 + 0(3^27)
1 1 0 2 2 0 1 2 2 1 2 1 1 0 2 2 1 0 1 1 0 0 2 2 2. 0 1
+ =
0 2 0 0 1 1 1 0 1 2 1 2 0 1 2 0 0 1 0 1 2 1 2 1 1. 0 1
2718281/828

11/4 + 0(11^13)
8 2 8 2 8 2 8 2 8 2 8 3 0
679001/207 + 0(11^13)
8 7 9 5 6 10 6 3 6 4 2 10 9
+ =
5 10 6 8 4 2 3 6 3 7 0 2 9
2718281/828

-22/7 + 0(2^37)
1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 1 0
46071/379 + 0(2^37)
1 1 1 1 1 1 0 1 0 1 0 0 1 1 0 0 0 1 0 1 0 0 1 1 1 1 0 0 0 1 0 1 1 0 1 0 1
+ =
1 0 0 0 1 1 1 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 0 0 1 0 1 0 1 1 1 1 1 0 1 1
314159/2653

-22/7 + 0(3^23)
1 0 2 1 2 0 1 0 2 1 2 0 1 0 2 1 2 0 1 0 2 0 2
46071/379 + 0(3^23)
2 0 1 2 1 2 1 2 2 1 2 1 0 0 2 2 0 1 1 2 1 0 0
+ =
0 1 1 1 1 0 0 0 2 0 1 1 1 1 2 0 2 2 0 0 0 0 2
314159/2653

-22/7 + 0(7^13)
6 6 6 6 6 6 6 6 6 6 6 3. 6
46071/379 + 0(7^13)
6 4 1 6 6 5 1 2 2 1 3 2 4
+ =
4 1 6 6 5 1 2 2 1 3 2 0. 6
314159/2653

-101/109 + 0(2^40)
0 1 1 1 1 1 1 0 1 1 0 1 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 0 1 0 1 1 0 0 1 0 0 1 1 1
583376/6649 + 0(2^40)
1 0 1 0 0 1 0 1 1 1 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 0 0 0
+ =
0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 0 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 0 1 1 1
577215/6649

-101/109 + 0(61^7)
33 1 7 16 48 7 50
583376/6649 + 0(61^7)
33 1 7 16 49 34. 35
+ =
34 8 24 3 57 23. 35
577215/6649

-101/109 + 0(32749^3)
5107 21031 15322
583376/6649 + 0(32749^3)
5452 13766 16445
+ =
10560 2048 31767
577215/6649
</pre>
</pre>
Retrieved from "https://rosettacode.org/wiki/P-Adic_numbers,_basic"