P-Adic numbers, basic: Difference between revisions
→{{header|Phix}}
(→{{header|Wren}}: Aligned examples with revised FB ones.) |
|||
Line 591:
47099/10977
</pre>
=={{header|Phix}}==
{{libheader|Phix/Class}}
<lang Phix>// constants
constant EMX = 64 // exponent maximum (if indexing starts at -EMX)
constant DMX = 1e5 // approximation loop maximum
constant AMX = 1048576 // argument maximum
constant PMAX = 32749 // prime maximum
// global variables
integer p1 = 0
integer p = 7 // default prime
integer k = 11 // precision
type Ratio(sequence r)
return length(r)=2 and integer(r[1]) and integer(r[2])
end type
procedure pad_to(string fmt, sequence data, integer len)
fmt = sprintf(fmt,data)
puts(1,fmt&repeat(' ',len-length(fmt)))
end procedure
class Padic
integer v = 0
sequence d = repeat(0,EMX*2)
// (re)initialize 'this' from Ratio, set 'sw' to print
function r2pa(Ratio q, integer sw)
integer {a,b} = q
if b=0 then return 1 end if
if b<0 then
b = -b
a = -a
end if
if abs(a)>AMX or b>AMX then return -1 end if
if p<2 or k<1 then return 1 end if
p = min(p, PMAX) // maximum short prime
k = min(k, EMX-1) // maximum array length
if sw!=0 then
-- numerator, denominator, prime, precision
pad_to("%d/%d + O(%d^%d)",{a,b,p,k},30)
end if
// (re)initialize
v = 0
p1 = p - 1
sequence ntd = repeat(0,2*EMX) -- (new this.d)
if a=0 then return 0 end if
// find -exponent of p in b
integer i = 0
while remainder(b,p)=0 do
b /= p
i -= 1
end while
integer s = 0,
r = remainder(b,p)
// modular inverse for small P
integer b1 = 1
while b1<=p1 do
s += r
if s>p1 then s -= p end if
if s=1 then exit end if
b1 += 1
end while
if b1=p then
printf(1,"r2pa: impossible inverse mod")
return -1
end if
v = EMX
while true do
// find exponent of P in a
while remainder(a,p)=0 do
a /= p
i += 1
end while
// valuation
if v=EMX then v = i end if
// upper bound
if i>=EMX then exit end if
// check precision
if i-v>k then exit end if
// next digit
integer rdx = remainder(a*b1,p)
if rdx<0 then rdx += p end if
if rdx<0 or rdx>=p then ?9/0 end if -- sanity chk
ntd[i+EMX+1] = rdx
// remainder - digit * divisor
a -= rdx*b
if a=0 then exit end if
end while
this.d = ntd
return 0
end function
// Horner's rule
function dsum()
integer t = min(v, 0),
s = 0
for i=k-1+t to t by -1 do
integer r = s
s *= p
if r!=0 and floor(s/r)-p!=0 then
// overflow
s = -1
exit
end if
s += d[i+EMX+1]
end for
return s
end function
// add b to 'this'
function add(Padic b)
integer c = 0
Padic r = new({min(v,b.v)})
sequence rd = r.d
for i=r.v to k+r.v do
integer dx = i+EMX+1
c += d[dx] + b.d[dx]
if c>p1 then
rd[dx] = c - p
c = 1
else
rd[dx] = c
c = 0
end if
end for
r.d = rd
return r
end function
// complement
function complement()
integer c = 1
Padic r = new({v})
sequence rd = r.d
for i=v to k+v do
integer dx = i+EMX+1
c += p1 - this.d[dx]
if c>p1 then
rd[dx] = c - p
c = 1
else
rd[dx] = c
c = 0
end if
end for
r.d = rd
return r
end function
// rational reconstruction
procedure crat()
integer sgn = 1
Padic s = this
integer j = 0,
i = 1
// denominator count
while i<=DMX do
// check for integer
j = k-1+v
while j>=v and s.d[j+EMX+1]=0 do
j -= 1
end while
if ((j-v)*2)<k then exit end if
// check for negative integer
j = k-1+v
while j>=v and p1-s.d[j+EMX+1]=0 do
j -= 1
end while
if ((j-v)*2)<k then
s = s.complement()
sgn = -1
exit
end if
// repeatedly add self to s
s = s.add(this)
i += 1
end while
// numerator: weighted digit sum
integer x = s.dsum(),
y = i
if x<0 or y>DMX then
printf(1,"crat: fail")
else
// negative powers
for i=v to -1 do
y *= p
end for
pad_to(iff(y=1?"%d":"%d/%d"),{x*sgn,y},26)
printf(1,"+ = ")
end if
end procedure
// print expansion
procedure prntf(bool sw)
integer t = min(v, 0)
// rational approximation
if sw!=0 then crat() end if
for i=k-1+t to t by -1 do
printf(1,"%d",d[i+EMX+1])
printf(1,iff(i=0 and v<0?". ":" "))
end for
printf(1,"\n")
end procedure
end class
sequence 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}},
-- all tested, but let's keep the output reasonable:
-- {{-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}},
-- /* four decadic pairs */
-- {{6, 7}, 10, 7, {-5, 7}},
-- {{2, 7}, 10, 7, {-3, 7}},
-- {{2, 7}, 10, 7, {-1, 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}},
-- {{-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}}
}
integer sw = 0,qa,qb
Padic a = new()
Padic b = new()
for i=1 to length(data) do
{Ratio q, p, k, Ratio q2} = data[i]
sw = a.r2pa(q, 1)
if sw=1 then exit end if
a.prntf(0)
sw = sw or b.r2pa(q2, 1)
if sw=1 then exit end if
if sw=0 then
b.prntf(0)
Padic c = a.add(b)
c.prntf(1)
end if
printf(1,"\n")
end for</lang>
{{out}}
<pre>
2/1 + O(2^4) 0 0 1 0
1/1 + O(2^4) 0 0 0 1
3 + = 0 0 1 1
4/1 + O(2^4) 0 1 0 0
3/1 + O(2^4) 0 0 1 1
-2/2 + = 0 1 1 1
4/1 + O(2^5) 0 0 1 0 0
3/1 + O(2^5) 0 0 0 1 1
7 + = 0 0 1 1 1
4/9 + O(5^4) 4 2 1 1
8/9 + O(5^4) 3 4 2 2
4/3 + = 3 1 3 3
5/8 + O(7^11) 4 2 4 2 4 2 4 2 4 2 5
353/30809 + O(7^11) 2 3 6 6 3 6 4 3 4 5 5
47099/10977 + = 6 6 4 2 1 2 1 6 2 1 3
</pre>
| |||