I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

P-Adic numbers, basic

From Rosetta Code
Task
P-Adic numbers, basic
You are encouraged to solve this task according to the task description, using any language you may know.
Conversion and addition of p-adic Numbers.


Task.

Convert two rationals to p-adic numbers and add them up. Rational reconstruction is needed to interpret the result.

p-Adic numbers were introduced around 1900 by Hensel. p-Adic expansions (a series of digits 0 ≤ d < p times p-power weights) are finite-tailed and tend to zero in the direction of higher positive powers of p (to the left in the notation used here). For example, the number 4 (100.0) has smaller 2-adic norm than 1/4 (0.01).

If we convert a natural number, the familiar p-ary expansion is obtained: 10 decimal is 1010 both binary and 2-adic. To convert a rational number a/b we perform p-adic long division. If p is actually prime, this is always possible if first the 'p-part' is removed from b (and the p-adic point shifted accordingly). The inverse of b modulo p is then used in the conversion.

Recipe: at each step the most significant digit of the partial remainder (initially a) is zeroed by subtracting a proper multiple of the divisor b. Shift out the zero digit (divide by p) and repeat until the remainder is zero or the precision limit is reached. Because p-adic division starts from the right, the 'proper multiplier' is simply d = partial remainder * 1/b (mod p). The d's are the successive p-adic digits to find.

Addition proceeds as usual, with carry from the right to the leftmost term, where it has least magnitude and just drops off. We can work with approximate rationals and obtain exact results. The routine for rational reconstruction demonstrates this: repeatedly add a p-adic to itself (keeping count to determine the denominator), until an integer is reached (the numerator then equals the weighted digit sum). But even p-adic arithmetic fails if the precision is too low. The examples mostly set the shortest prime-exponent combinations that allow valid reconstruction.


Related task.

p-Adic square roots


Reference.

[1] p-Adic expansions



FreeBASIC[edit]

 
' ***********************************************
'subject: convert two rationals to p-adic numbers,
' add them up and show the result.
'tested : FreeBasic 1.07.0
 
 
'you can change this:
 
const emx = 64
'exponent maximum
 
const dmx = 100000
'approximation loop maximum
 
 
'better not change
'------------------------------------------------
const amx = 1048576
'argument maximum
 
const Pmax = 32749
'max. prime < 2^15
 
 
type ratio
as longint a, b
end type
 
type padic
declare function r2pa (byref q as ratio, byval sw as integer) as integer
'convert q = a/b to p-adic number, set sw to print
declare sub printf (byval sw as integer)
'print expansion, set sw to print rational
declare sub crat ()
'rational reconstruction
 
declare sub add (byref a as padic, byref b as padic)
'let self:= a + b
declare sub cmpt (byref a as padic)
'let self:= complement_a
 
declare function dsum () as long
'weighted digit sum
 
as long d(-emx to emx - 1)
as integer v
end type
 
 
'global variables
dim shared as long p1, p = 7
'default prime
dim shared as integer k = 11
'precision
 
#define min(a, b) iif((a) > (b), b, a)
 
 
'------------------------------------------------
'convert rational a/b to p-adic number
function padic.r2pa (byref q as ratio, byval sw as integer) as integer
dim as longint a = q.a, b = q.b
dim as long r, s, b1
dim i as integer
r2pa = 0
 
if b = 0 then return 1
if b < 0 then b = -b: a = -a
if abs(a) > amx or b > amx then return -1
if p < 2 or k < 1 then return 1
 
'max. short prime
p = min(p, Pmax)
'max. array length
k = min(k, emx - 1)
 
if sw then
'echo numerator, denominator,
print a;"/";str(b);" + ";
'prime and precision
print "O(";str(p);"^";str(k);")"
end if
 
'initialize
v = 0
p1 = p - 1
for i = -emx to emx - 1
d(i) = 0: next
 
if a = 0 then return 0
 
i = 0
'find -exponent of p in b
do until b mod p
b \= p: i -= 1
loop
 
s = 0
r = b mod p
'modular inverse for small p
for b1 = 1 to p1
s += r
if s > p1 then s -= p
if s = 1 then exit for
next b1
 
if b1 = p then
print "r2pa: impossible inverse mod"
return -1
end if
 
v = emx
do
'find exponent of p in a
do until a mod p
a \= p: i += 1
loop
 
'valuation
if v = emx then v = i
 
'upper bound
if i >= emx then exit do
'check precision
if (i - v) > k then exit do
 
'next digit
d(i) = a * b1 mod p
if d(i) < 0 then d(i) += p
 
'remainder - digit * divisor
a -= d(i) * b
loop while a
end function
 
'------------------------------------------------
'Horner's rule
function padic.dsum () as long
dim as integer i, t = min(v, 0)
dim as long r, s = 0
 
for i = k - 1 + t to t step -1
r = s: s *= p
if r andalso s \ r - p then
'overflow
s = -1: exit for
end if
s += d(i)
next i
 
return s
end function
 
#macro pint(cp)
for j = k - 1 + v to v step -1
if cp then exit for
next j
fl = ((j - v) shl 1) < k
#endmacro
 
'rational reconstruction
sub padic.crat ()
dim as integer i, j, fl
dim as padic s = this
dim as long x, y
 
'denominator count
for i = 1 to dmx
'check for integer
pint(s.d(j))
if fl then fl = 0: exit for
 
'check negative integer
pint(p1 - s.d(j))
if fl then exit for
 
'repeatedly add self to s
s.add(s, this)
next i
 
if fl then s.cmpt(s)
 
'numerator: weighted digit sum
x = s.dsum: y = i
 
if x < 0 or y > dmx then
print "crat: fail"
 
else
'negative powers
for i = v to -1
y *= p: next
 
'negative rational
if fl then x = -x
 
print x;
if y > 1 then print "/";str(y);
print
end if
end sub
 
 
'print expansion
sub padic.printf (byval sw as integer)
dim as integer i, t = min(v, 0)
 
for i = k - 1 + t to t step -1
print d(i);
if i = 0 andalso v < 0 then print ".";
next i
print
 
'rational approximation
if sw then crat
end sub
 
'------------------------------------------------
'carry
#macro cstep(dt)
if c > p1 then
dt = c - p: c = 1
else
dt = c: c = 0
end if
#endmacro
 
'let self:= a + b
sub padic.add (byref a as padic, byref b as padic)
dim i as integer, r as padic
dim as long c = 0
with r
.v = min(a.v, b.v)
 
for i = .v to k +.v
c += a.d(i) + b.d(i)
cstep(.d(i))
next i
end with
this = r
end sub
 
'let self:= complement_a
sub padic.cmpt (byref a as padic)
dim i as integer, r as padic
dim as long c = 1
with r
.v = a.v
 
for i = .v to k +.v
c += p1 - a.d(i)
cstep(.d(i))
next i
end with
this = r
end sub
 
 
'main
'------------------------------------------------
dim as integer sw
dim as padic a, b, c
dim q as ratio
 
width 64, 30
cls
 
'rational reconstruction
'depends on the precision -
'until the dsum-loop overflows.
data 2,1, 2,4
data 1,1
 
data 4,1, 2,4
data 3,1
 
data 4,1, 2,5
data 3,1
 
' 4/9 + O(5^4)
data 4,9, 5,4
data 8,9
 
data 26,25, 5,4
data -109,125
 
data 49,2, 7,6
data -4851,2
 
data -9,5, 3,8
data 27,7
 
data 5,19, 2,12
data -101,384
 
'two 'decadic' pairs
data 2,7, 10,7
data -1,7
 
data 34,21, 10,9
data -39034,791
 
'familiar digits
data 11,4, 2,43
data 679001,207
 
data -8,9, 23,9
data 302113,92
 
data -22,7, 3,23
data 46071,379
 
data -22,7, 32749,3
data 46071,379
 
data 35,61, 5,20
data 9400,109
 
data -101,109, 61,7
data 583376,6649
 
data -25,26, 7,13
data 5571,137
 
data 1,4, 7,11
data 9263,2837
 
data 122,407, 7,11
data -517,1477
 
'more subtle
data 5,8, 7,11
data 353,30809
 
data 0,0, 0,0
 
 
print
do
read q.a,q.b, p,k
 
sw = a.r2pa(q, 1)
if sw = 1 then exit do
a.printf(0)
 
read q.a,q.b
 
sw or= b.r2pa(q, 1)
if sw = 1 then exit do
if sw then continue do
b.printf(0)
 
c.add(a, b)
print "+ ="
c.printf(1)
 
print : ?
loop
 
system
 
Examples:
 2/1 + O(2^4)
 0 0 1 0
 1/1 + O(2^4)
 0 0 0 1
+ =
 0 0 1 1
 3


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


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


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


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


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


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


 5/19 + O(2^12)
 0 0 1 0 1 0 0 0 0 1 1 1
-101/384 + O(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 + O(10^7)
 5 7 1 4 2 8 6
-1/7 + O(10^7)
 7 1 4 2 8 5 7
+ =
 2 8 5 7 1 4 3
 1/7


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


 11/4 + O(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 + O(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 + O(23^9)
 2 12 17 20 10 5 2 12 17
 302113/92 + O(23^9)
 5 17 5 17 6 0 10 12. 2
+ =
 18 12 3 4 11 3 0 6. 2
 2718281/828


-22/7 + O(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 + O(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 + O(32749^3)
 28070 18713 23389
 46071/379 + O(32749^3)
 4493 8727 10145
+ =
 32563 27441 785
 314159/2653


 35/61 + O(5^20)
 2 3 2 3 0 2 4 1 3 3 0 0 4 0 2 2 1 2 2 0
 9400/109 + O(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 + O(61^7)
 33 1 7 16 48 7 50
 583376/6649 + O(61^7)
 33 1 7 16 49 34. 35
+ =
 34 8 24 3 57 23. 35
 577215/6649


-25/26 + O(7^13)
 2 6 5 0 5 4 4 0 1 6 1 2 2
 5571/137 + O(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 + O(7^11)
 1 5 1 5 1 5 1 5 1 5 2
 9263/2837 + O(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 + O(7^11)
 6 2 0 3 0 6 2 4 4 4 3
-517/1477 + O(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 + 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
+ =
 6 6 4 2 1 2 1 6 2 1 3
 47099/10977

Go[edit]

Translation of: FreeBASIC
package main
 
import "fmt"
 
// constants
const EMX = 64 // exponent maximum (if indexing starts at -EMX)
const DMX = 100000 // approximation loop maximum
const AMX = 1048576 // argument maximum
const PMAX = 32749 // prime maximum
 
// global variables
var p1 = 0
var p = 7 // default prime
var k = 11 // precision
 
func abs(a int) int {
if a >= 0 {
return a
}
return -a
}
 
func min(a, b int) int {
if a < b {
return a
}
return b
}
 
type Ratio struct {
a, b int
}
 
type Padic struct {
v int
d [2 * EMX]int // add EMX to index to be consistent wih FB
}
 
// (re)initialize receiver from Ratio, set 'sw' to print
func (pa *Padic) r2pa(q Ratio, sw int) int {
a := q.a
b := q.b
if b == 0 {
return 1
}
if b < 0 {
b = -b
a = -a
}
if abs(a) > AMX || b > AMX {
return -1
}
if p < 2 || k < 1 {
return 1
}
p = min(p, PMAX) // maximum short prime
k = min(k, EMX-1) // maxumum array length
if sw != 0 {
fmt.Printf("%d/%d + ", a, b) // numerator, denominator
fmt.Printf("0(%d^%d)\n", p, k) // prime, precision
}
 
// (re)initialize
pa.v = 0
p1 = p - 1
pa.d = [2 * EMX]int{}
if a == 0 {
return 0
}
i := 0
 
// find -exponent of p in b
for b%p == 0 {
b = b / p
i--
}
s := 0
r := b % p
 
// modular inverse for small p
b1 := 1
for b1 <= p1 {
s += r
if s > p1 {
s -= p
}
if s == 1 {
break
}
b1++
}
if b1 == p {
fmt.Println("r2pa: impossible inverse mod")
return -1
}
pa.v = EMX
for {
// find exponent of P in a
for a%p == 0 {
a = a / p
i++
}
 
// valuation
if pa.v == EMX {
pa.v = i
}
 
// upper bound
if i >= EMX {
break
}
 
// check precision
if (i - pa.v) > k {
break
}
 
// next digit
pa.d[i+EMX] = a * b1 % p
if pa.d[i+EMX] < 0 {
pa.d[i+EMX] += p
}
 
// remainder - digit * divisor
a -= pa.d[i+EMX] * b
if a == 0 {
break
}
}
return 0
}
 
// Horner's rule
func (pa *Padic) dsum() int {
t := min(pa.v, 0)
s := 0
for i := k - 1 + t; i >= t; i-- {
r := s
s *= p
if r != 0 && (s/r-p != 0) {
// overflow
s = -1
break
}
s += pa.d[i+EMX]
}
return s
}
 
// add b to receiver
func (pa *Padic) add(b Padic) *Padic {
c := 0
r := Padic{}
r.v = min(pa.v, b.v)
for i := r.v; i <= k+r.v; i++ {
c += pa.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 of receiver
func (pa *Padic) cmpt() *Padic {
c := 1
r := Padic{}
r.v = pa.v
for i := pa.v; i <= k+pa.v; i++ {
c += p1 - pa.d[i+EMX]
if c > p1 {
r.d[i+EMX] = c - p
c = 1
} else {
r.d[i+EMX] = c
c = 0
}
}
return &r
}
 
// rational reconstruction
func (pa *Padic) crat() {
fl := false
s := pa
j := 0
i := 1
 
// denominator count
for i <= DMX {
// check for integer
j = k - 1 + pa.v
for j >= pa.v {
if s.d[j+EMX] != 0 {
break
}
j--
}
fl = ((j - pa.v) * 2) < k
if fl {
fl = false
break
}
 
// check negative integer
j = k - 1 + pa.v
for j >= pa.v {
if p1-s.d[j+EMX] != 0 {
break
}
j--
}
fl = ((j - pa.v) * 2) < k
if fl {
break
}
 
// repeatedly add self to s
s = s.add(*pa)
i++
}
if fl {
s = s.cmpt()
}
 
// numerator: weighted digit sum
x := s.dsum()
y := i
if x < 0 || y > DMX {
fmt.Println(x, y)
fmt.Println("crat: fail")
} else {
// negative powers
i = pa.v
for i <= -1 {
y *= p
i++
}
 
// negative rational
if fl {
x = -x
}
fmt.Print(x)
if y > 1 {
fmt.Printf("/%d", y)
}
fmt.Println()
}
}
 
// print expansion
func (pa *Padic) printf(sw int) {
t := min(pa.v, 0)
for i := k - 1 + t; i >= t; i-- {
fmt.Print(pa.d[i+EMX])
if i == 0 && pa.v < 0 {
fmt.Print(".")
}
fmt.Print(" ")
}
fmt.Println()
// rational approximation
if sw != 0 {
pa.crat()
}
}
 
func main() {
data := [][]int{
/* rational reconstruction depends on the precision
until the dsum-loop overflows */

{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},
}
 
sw := 0
a := Padic{}
b := Padic{}
 
for _, d := range data {
q := Ratio{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)
c := a.add(b)
fmt.Println("+ =")
c.printf(1)
}
fmt.Println()
}
}
Output:
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

Julia[edit]

Uses the Nemo abstract algebra library. The Nemo library's rational reconstruction function gives up quite easily, so another alternative to FreeBasic's crat() using vector products is below.

using Nemo, LinearAlgebra
 
set_printing_mode(FlintPadicField, :terse)
 
""" convert to Rational (rational reconstruction) """
function toRational(pa::padic)
rat = lift(QQ, pa)
r, den = BigInt(numerator(rat)), Int(denominator(rat))
p, k = Int(prime(parent(pa))), Int(precision(pa))
N = BigInt(p^k)
a1, a2 = [N, 0], [r, 1]
while dot(a1, a1) > dot(a2, a2)
q = dot(a1, a2) // dot(a2, a2)
a1, a2 = a2, a1 - BigInt(round(q)) * a2
end
if dot(a1, a1) < N
return (Rational{Int}(a1[1]) // Rational{Int}(a1[2])) // Int(den)
else
return Int(r) // den
end
end
 
function dstring(pa::padic)
u, v, n, p, k = pa.u, pa.v, pa.N, pa.parent.p, pa.parent.prec_max
d = digits(v > 0 ? u * p^v : u, base=pa.parent.p, pad=k)
return prod([i == k + v && v != 0 ? "$x . " : "$x " for (i, x) in enumerate(reverse(d))])
end
 
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],
 
# Base 10 10-adic p-adics are not allowed by Nemo library -- p must be a prime
 
# 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 (num1, den1, P, K, num2, den2) in DATA
Qp = PadicField(P, K)
a = Qp(QQ(num1 // den1))
b = Qp(QQ(num2 // den2))
c = a + b
r = toRational(c)
println(a, "\n", dstring(a), "\n", b, "\n", dstring(b), "\n+ =\n", c, "\n", dstring(c), " $r\n")
end
 
Output:
2 + O(2^5)
0 0 1 0   
1 + O(2^4)
0 0 0 1   
+ =       
3 + O(2^4)
0 0 1 1    3//1

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

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

556 + O(5^4)
4 2 1 1
487 + O(5^4)
3 4 2 2
+ =
418 + O(5^4)
3 1 3 3    4//3

26/25 + O(5^2)
0 1 . 0 1
516/125 + O(5^1)
4 . 0 3 1
+ =
21/125 + O(5^1)
0 . 0 4 1    21//125

58849 + O(7^6)
3 3 3 4 0 0
56399 + O(7^6)
3 2 3 3 0 0
+ =
115248 + O(7^6)
6 6 0 0 0 0    0//1

5247 + O(3^8)
2 1 0 1 2 1 0 0
3753 + O(3^8)
1 2 0 1 1 0 0 0
+ =
2439 + O(3^8)
1 0 1 0 0 1 0 0    72//35

647 + O(2^12)
0 0 1 0 1 0 0 0 0 1 1 1
2697/128 + O(2^5)
1 0 1 0 1 . 0 0 0 1 0 0 1
+ =
3593/128 + O(2^5)
1 1 1 0 0 . 0 0 0 1 0 0 1    3593//128

11/4 + O(2^41)
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
2379619371607 + O(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
+ =
722384464231/4 + O(2^41)
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    330311//1618052

200128073495 + O(23^9)
2 12 17 20 10 5 2 12 17
450288240894/23 + O(23^8)
5 17 5 17 6 0 10 12 . 2
+ =
1450928608353/23 + O(23^8)
18 12 3 4 11 3 0 6 . 2    1450928608353//23

40347076637 + O(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
69303290076 + O(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
+ =
15507187886 + O(3^23)
0 1 1 1 1 0 0 0 2 0 1 1 1 1 2 0 2 2 0 0 0 0 2    15507187886//1

30105603673496 + O(32749^3)
28070 18713 23389
4819014836161 + O(32749^3)
4493 8727 10145
+ =
34924618509657 + O(32749^3)
32563 27441 785    314159//2653

51592217117060 + O(5^20)
2 3 2 3 0 2 4 1 3 3 0 0 4 0 2 2 1 2 2 0
64744861847850 + O(5^20)
3 1 4 4 1 2 3 4 4 3 4 1 1 3 1 1 2 4 0 0
+ =
20969647324285 + O(5^20)
1 0 2 2 2 0 3 1 3 1 4 2 0 3 3 3 4 1 2 0    577215//6649

1701117681882 + O(61^7)
33 1 7 16 48 7 50
1701117687235/61 + O(61^6)
33 1 7 16 49 34 . 35
+ =
1758782693344/61 + O(61^6)
34 8 24 3 57 23 . 35    1758782693344//61

40991504402 + O(7^13)
2 6 5 0 5 4 4 0 1 6 1 2 2
46676457609 + O(7^13)
3 2 4 1 4 5 4 2 2 5 5 3 5
+ =
87667962011 + O(7^13)
6 2 2 2 3 3 1 2 4 4 6 6 0    141421//3562

494331686 + O(7^11)
1 5 1 5 1 5 1 5 1 5 2
1936205041 + O(7^11)
6 5 6 6 0 3 2 0 4 4 1
+ =
453209984 + O(7^11)
1 4 1 4 2 1 3 5 6 2 3    39889//11348

1778136580 + O(7^11)
6 2 0 3 0 6 2 4 4 4 3
384219886/7 + O(7^10)
1 2 3 4 3 5 4 6 4 1 . 1
+ =
967215488/7 + O(7^10)
3 2 6 5 3 1 2 4 1 4 . 1    967215488//7

1235829215 + O(7^11)
4 2 4 2 4 2 4 2 4 2 5
726006041 + O(7^11)
2 3 6 6 3 6 4 3 4 5 5
+ =
1961835256 + O(7^11)
6 6 4 2 1 2 1 6 2 1 3    -25145//36122

Phix[edit]

Library: Phix/Class
// 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
Output:
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

Raku[edit]

# 20210225 Raku programming solution
 
class Padic { has ($.p is default(2), @.v is default([])) is rw ;
 
method r2pa (Rat $x is copy, \p, \d) { # Reference: math.stackexchange.com/a/1187037
self.p = p;
while +self.v < d {
my %d = ^p Z=> (( $x «-« ^p ) »/» p )».&{ .denominator % p }; # .kv
for %d.keys { self.v.unshift: $_ and last if %d{$_} != 0 }
$x = ($x - self.v.first) / p ;
}
self
}
 
method add (Padic \x, \p) {
my $div = 0;
return Padic.new:
p => p,
v => reverse gather for ( self.v.reverse Z x.v.reverse )».sum
{ take ( $_ + $div ) % p and $div = $_ div p } # .polymod
}
 
method gist { # en.wikipedia.org/wiki/P-adic_number#Notation
my %H = (0..9) Z=> ('₀'..'₉'); # (0x2080 .. 0x2089);
'⋯ ' ~ self.v ~ ' ' ~ [~] self.p.comb».&{ %H{$_} }
}
}
 
my @T;
for my \D = (
#`[[[[[ these are not working
< -7/5 99/70 7 4> ,
< 26/25 -109/125 5 5 >,
< 49/2 -4851/2 7 6 >,
< -9/5 27/7 3 8>,
< 5/19 -101/384 2 12>,
< 6/7 -5/7 10 7 >,
< 2/7 -3/7 10 7 >,
< 2/7 -1/7 10 7 >,
<34/21 -39034/791 10 9 >,
< 11/4 679001/207 2 43>,
< 11/4 679001/207 3 27 >,
< -22/7 46071/379 2 37 >,
< -22/7 46071/379 3 23 >,
< -22/7 46071/379 7 13 >,
< -101/109 583376/6649 2 40>,
< -101/109 583376/6649 61 7>,
< -101/109 583376/6649 32749 3>,
< -25/26 5571/137 7 13>,
< 122/407 -517/1477 7 11>,
# ]]]]]
#`[[[[[ working cases
< 2/1 1/1 2 4>,
< 4/1 3/1 2 4>,
< 4/1 3/1 2 5>,
< 4/9 8/9 5 4>,
< 11/4 679001/207 11 13 >,
< 1/4 9263/2837 7 11 >,
# ]]]]]
 
< 5/8 353/30809 7 11 >,
) -> \D {
given @T[0] = Padic.new { say D[0]~' = ', .r2pa: D[0],D[2],D[3] }
given @T[1] = Padic.new { say D[1]~' = ', .r2pa: D[1],D[2],D[3] }
given @T[0].add: @T[1], D[2] {
given @T[2] = Padic.new { .r2pa: D[0]+D[1], D[2], D[3] }
say "Addition result = ", $_.gist, ( $_.gist eqv @T[2].gist ) ?? ''
!! ' but ' ~ (D[0]+D[1]).nude.join('/') ~ ' = ' ~ @T[2].gist
}
}
 
Output:
5/8 = ⋯ 4 2 4 2 4 2 4 2 4 2 5 ₇
353/30809 = ⋯ 2 3 6 6 3 6 4 3 4 5 5 ₇
Addition result = ⋯ 6 6 4 2 1 2 1 6 2 1 3 ₇

Wren[edit]

Translation of: FreeBASIC
Library: Wren-dynamic
Library: Wren-math
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 depends on the precision
until the dsum-loop overflows */

[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]
]
 
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()
}
Output:
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