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