P-Adic numbers, basic: Difference between revisions
Content added Content deleted
(→{{header|Raku}}: update) |
|||
Line 1,304: | Line 1,304: | ||
1961835256 + O(7^11) |
1961835256 + O(7^11) |
||
6 6 4 2 1 2 1 6 2 1 3 -25145//36122 |
6 6 4 2 1 2 1 6 2 1 3 -25145//36122 |
||
</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Go}} |
|||
Translation of Go with some modifications, especially using exceptions when an error occurs. |
|||
<lang Nim>import math, strformat |
|||
const |
|||
Emx = 64 # Exponent maximum (if indexing starts at -Emx). |
|||
Dmx = 100000 # Approximation loop maximum. |
|||
Amx = 1048576 # 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 r2pa(pa: var Padic; q: Ratio; sw: int) = |
|||
## Convert "q" to p-adic number, set "sw" to print. |
|||
var (a, b) = q |
|||
if b == 0: |
|||
raise newException(PadicError, &"Wrong rational: {a}/{b}" ) |
|||
if b < 0: |
|||
b = -b |
|||
a = -a |
|||
if abs(a) > Amx or b > Amx: |
|||
raise newException(PadicError, &"Rational exceeding limits: {a}/{b}") |
|||
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. |
|||
pa.k = min(pa.k, Emx - 1) # Maximum array length. |
|||
if sw != 0: |
|||
echo &"{a}/{b} + 0({pa.p}^{pa.k})" |
|||
# Initialize. |
|||
pa.v = 0 |
|||
pa.d.reset() |
|||
if a == 0: return |
|||
var i = 0 |
|||
# Find -exponent of "p" in "b". |
|||
while b mod pa.p == 0: |
|||
b = b div pa.p |
|||
dec i |
|||
var s = 0 |
|||
var r = b mod pa.p |
|||
# Modular inverse for small "p". |
|||
var b1 = 1 |
|||
while b1 < pa.p: |
|||
inc s, r |
|||
if s >= pa.p: dec s, pa.p |
|||
if s == 1: break |
|||
inc b1 |
|||
if b1 == pa.p: |
|||
raise newException(PadicError, "Impossible to compute reverse modulo") |
|||
pa.v = Emx |
|||
while true: |
|||
# Find exponent of "p" in "a". |
|||
while a mod pa.p == 0: |
|||
a = a div pa.p |
|||
inc i |
|||
# Valuation. |
|||
if pa.v == Emx: pa.v = i |
|||
# Upper bound. |
|||
if i >= Emx: break |
|||
# Check precision. |
|||
if i - pa.v > pa.k: break |
|||
# Next digit. |
|||
pa.d[i] = floorMod(a * b1, pa.p) |
|||
# Remainder - digit * divisor. |
|||
dec a, pa.d[i] * b |
|||
if a == 0: break |
|||
func dsum(pa: Padic): int = |
|||
## Horner's rule. |
|||
let t = min(pa.v, 0) |
|||
for i in countdown(pa.k - 1 + t, t): |
|||
var r = result |
|||
result *= pa.p |
|||
if r != 0 and (result div r - pa.p) != 0: |
|||
return -1 # Overflow. |
|||
inc result, pa.d[i] |
|||
func `+`(pa, pb: Padic): Padic = |
|||
## Add two p-adic numbers. |
|||
assert pa.p == pb.p and pa.k == pb.k |
|||
result.p = pa.p |
|||
result.k = pa.k |
|||
var c = 0 |
|||
result.v = min(pa.v, pb.v) |
|||
for i in result.v..(pa.k + result.v): |
|||
inc c, pa.d[i] + pb.d[i] |
|||
if c >= pa.p: |
|||
result.d[i] = c - pa.p |
|||
c = 1 |
|||
else: |
|||
result.d[i] = c |
|||
c = 0 |
|||
func cmpt(pa: Padic): Padic = |
|||
## Return the complement. |
|||
var c = 1 |
|||
result.p = pa.p |
|||
result.k = pa.k |
|||
result.v = pa.v |
|||
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 crat(pa: Padic): string = |
|||
## Rational reconstruction. |
|||
var s = pa |
|||
# Denominator count. |
|||
var i = 1 |
|||
var fl = false |
|||
while i <= Dmx: |
|||
# Check for integer. |
|||
var j = pa.k - 1 + pa.v |
|||
while j >= pa.v: |
|||
if s.d[j] != 0: break |
|||
dec j |
|||
fl = (j - pa.v) * 2 < pa.k |
|||
if fl: |
|||
fl = false |
|||
break |
|||
# Check negative integer. |
|||
j = pa.k - 1 + pa.v |
|||
while j >= pa.v: |
|||
if pa.p - 1 - s.d[j] != 0: break |
|||
dec j |
|||
fl = (j - pa.v) * 2 < pa.k |
|||
if fl: break |
|||
# Repeatedly add "pa" to "s". |
|||
s = s + pa |
|||
inc i |
|||
if fl: s = s.cmpt() |
|||
# Numerator: weighted digit sum. |
|||
var x = s.dsum() |
|||
var y = i |
|||
if x < 0 or y > Dmx: |
|||
raise newException(PadicError, &"Error during rational reconstruction: {x}, {y}") |
|||
# Negative powers. |
|||
for i in pa.v..(-1): y *= pa.p |
|||
# Negative rational. |
|||
if fl: x = -x |
|||
result = $x |
|||
if y > 1: result.add &"/{y}" |
|||
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 " " |
|||
proc print(pa: Padic; sw: int) = |
|||
echo pa |
|||
# Rational approximation. |
|||
if sw != 0: echo pa.crat() |
|||
when isMainModule: |
|||
# Rational reconstruction depends on the precision |
|||
# until the dsum-loop overflows. |
|||
const Data = [[2, 1, 2, 4, 1, 1], |
|||
[4, 1, 2, 4, 3, 1], |
|||
[4, 1, 2, 5, 3, 1], |
|||
[4, 9, 5, 4, 8, 9], |
|||
[26, 25, 5, 4, -109, 125], |
|||
[49, 2, 7, 6, -4851, 2], |
|||
[-9, 5, 3, 8, 27, 7], |
|||
[5, 19, 2, 12, -101, 384], |
|||
# Two decadic pairs. |
|||
[2, 7, 10, 7, -1, 7], |
|||
[34, 21, 10, 9, -39034, 791], |
|||
# Familiar digits. |
|||
[11, 4, 2, 43, 679001, 207], |
|||
[-8, 9, 23, 9, 302113, 92], |
|||
[-22, 7, 3, 23, 46071, 379], |
|||
[-22, 7, 32749, 3, 46071, 379], |
|||
[35, 61, 5, 20, 9400, 109], |
|||
[-101, 109, 61, 7, 583376, 6649], |
|||
[-25, 26, 7, 13, 5571, 137], |
|||
[1, 4, 7, 11, 9263, 2837], |
|||
[122, 407, 7, 11, -517, 1477], |
|||
# More subtle. |
|||
[5, 8, 7, 11, 353, 30809]] |
|||
for d in Data: |
|||
try: |
|||
var a, b = Padic(p: d[2], k: d[3]) |
|||
r2pa(a, (d[0], d[1]), 1) |
|||
print(a, 0) |
|||
r2pa(b, (d[4], d[5]), 1) |
|||
print(b, 0) |
|||
echo "+ =" |
|||
print(a + b, 1) |
|||
echo "" |
|||
except PadicError: |
|||
echo getCurrentExceptionMsg()</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 |
|||
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 |
|||
2/7 + 0(10^7) |
|||
5 7 1 4 2 8 6 |
|||
-1/7 + 0(10^7) |
|||
7 1 4 2 8 5 7 |
|||
+ = |
|||
2 8 5 7 1 4 3 |
|||
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 |
|||
-8/9 + 0(23^9) |
|||
2 12 17 20 10 5 2 12 17 |
|||
302113/92 + 0(23^9) |
|||
5 17 5 17 6 0 10 12. 2 |
|||
+ = |
|||
18 12 3 4 11 3 0 6. 2 |
|||
2718281/828 |
|||
-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(32749^3) |
|||
28070 18713 23389 |
|||
46071/379 + 0(32749^3) |
|||
4493 8727 10145 |
|||
+ = |
|||
32563 27441 785 |
|||
314159/2653 |
|||
35/61 + 0(5^20) |
|||
2 3 2 3 0 2 4 1 3 3 0 0 4 0 2 2 1 2 2 0 |
|||
9400/109 + 0(5^20) |
|||
3 1 4 4 1 2 3 4 4 3 4 1 1 3 1 1 2 4 0 0 |
|||
+ = |
|||
1 0 2 2 2 0 3 1 3 1 4 2 0 3 3 3 4 1 2 0 |
|||
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 |
|||
-25/26 + 0(7^13) |
|||
2 6 5 0 5 4 4 0 1 6 1 2 2 |
|||
5571/137 + 0(7^13) |
|||
3 2 4 1 4 5 4 2 2 5 5 3 5 |
|||
+ = |
|||
6 2 2 2 3 3 1 2 4 4 6 6 0 |
|||
141421/3562 |
|||
1/4 + 0(7^11) |
|||
1 5 1 5 1 5 1 5 1 5 2 |
|||
9263/2837 + 0(7^11) |
|||
6 5 6 6 0 3 2 0 4 4 1 |
|||
+ = |
|||
1 4 1 4 2 1 3 5 6 2 3 |
|||
39889/11348 |
|||
122/407 + 0(7^11) |
|||
6 2 0 3 0 6 2 4 4 4 3 |
|||
-517/1477 + 0(7^11) |
|||
1 2 3 4 3 5 4 6 4 1. 1 |
|||
+ = |
|||
3 2 6 5 3 1 2 4 1 4. 1 |
|||
-27584/90671 |
|||
5/8 + 0(7^11) |
|||
4 2 4 2 4 2 4 2 4 2 5 |
|||
353/30809 + 0(7^11) |
|||
2 3 6 6 3 6 4 3 4 5 5 |
|||
+ = |
|||
6 6 4 2 1 2 1 6 2 1 3 |
|||
47099/10977 |
|||
</pre> |
</pre> |
||