Anonymous user
P-Adic numbers, basic: Difference between revisions
→{{header|Nim}}
(→{{header|Raku}}: update) |
|||
Line 1,304:
1961835256 + O(7^11)
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>
| |||