(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
You are encouraged to solve this task according to the task description, using any language you may know.

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.

Reference.

## C++

This example displays p-adic numbers in standard mathematical format, consisting of a possibly infinite list of digits extending leftwards from the p-adic point. p-adic numbers are given corrrect to O(prime^40) and rational reconstructions are accurate to O(prime^20).

```#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>

class Rational {
public:
Rational(const int32_t& aNumerator, const int32_t& aDenominator) {
if ( aDenominator < 0 ) {
numerator = -aNumerator;
} else {
numerator = aNumerator;
}

if ( aNumerator == 0 ) {
denominator = 1;
}

const uint32_t divisor = std::gcd(numerator, denominator);
numerator /= divisor;
denominator /= divisor;
}

std::string to_string() const {
return std::to_string(numerator) + " / " + std::to_string(denominator);
}

private:
int32_t numerator;
int32_t denominator;
};

public:
// Create a P_adic number, with p = 'prime', from the given rational 'numerator' / 'denominator'.
P_adic(const uint32_t& prime, int32_t numerator, int32_t denominator) : prime(prime) {
if ( denominator == 0 ) {
throw std::invalid_argument("Denominator cannot be zero");
}

order = 0;

// Process rational zero
if ( numerator == 0 ) {
digits.assign(DIGITS_SIZE, 0);
order = ORDER_MAX;
return;
}

// Remove multiples of 'prime' and adjust the order of the P_adic number accordingly
while ( modulo_prime(numerator) == 0 ) {
numerator /= static_cast<int32_t>(prime);
order += 1;
}

while ( modulo_prime(denominator) == 0 ) {
denominator /= static_cast<int32_t>(prime);
order -= 1;
}

// Standard calculation of P_adic digits
const uint64_t inverse = modulo_inverse(denominator);
while ( digits.size() < DIGITS_SIZE ) {
const uint32_t digit = modulo_prime(numerator * inverse);
digits.emplace_back(digit);

numerator -= digit * denominator;

if ( numerator != 0 ) {
// The denominator is not a power of a prime
uint32_t count = 0;
while ( modulo_prime(numerator) == 0 ) {
numerator /= static_cast<int32_t>(prime);
count += 1;
}

for ( uint32_t i = count; i > 1; --i ) {
digits.emplace_back(0);
}
}
}
}

// Return the sum of this P_adic number with the given P_adic number.
if ( prime != other.prime ) {
}

std::vector<uint32_t> this_digits = digits;
std::vector<uint32_t> other_digits = other.digits;
std::vector<uint32_t> result;

for ( int32_t i = 0; i < -order + other.order; ++i ) {
other_digits.insert(other_digits.begin(), 0);
}

for ( int32_t i = 0; i < -other.order + order; ++i ) {
this_digits.insert(this_digits.begin(), 0);
}

// Standard digit by digit addition
uint32_t carry = 0;
for ( uint32_t i = 0; i < std::min(this_digits.size(), other_digits.size()); ++i ) {
const uint32_t sum = this_digits[i] + other_digits[i] + carry;
const uint32_t remainder = sum % prime;
carry = ( sum >= prime ) ? 1 : 0;
result.emplace_back(remainder);
}

return P_adic(prime, result, all_zero_digits(result) ? ORDER_MAX : std::min(order, other.order));
}

// Return the Rational representation of this P_adic number.
Rational convert_to_rational() {
std::vector<uint32_t> numbers = digits;

// Zero
if ( numbers.empty() || all_zero_digits(numbers) ) {
return Rational(1, 0);
}

// Positive integer
if ( order >= 0 && ends_with(numbers, 0) ) {
for ( int32_t i = 0; i < order; ++i ) {
numbers.emplace(numbers.begin(), 0);
}

return Rational(convert_to_decimal(numbers), 1);
}

// Negative integer
if ( order >= 0 && ends_with(numbers, prime - 1) ) {
negate_digits(numbers);
for ( int32_t i = 0; i < order; ++i ) {
numbers.emplace(numbers.begin(), 0);
}

return Rational(-convert_to_decimal(numbers), 1);
}

// Rational
int32_t denominator = 1;
do {
denominator += 1;
} while ( ! ( ends_with(sum.digits, 0) || ends_with(sum.digits, prime - 1) ) );

const bool negative = ends_with(sum.digits, 6);
if ( negative ) {
negate_digits(sum.digits);
}

int32_t numerator = negative ? -convert_to_decimal(sum.digits) : convert_to_decimal(sum.digits);

if ( order > 0 ) {
numerator *= std::pow(prime, order);
}

if ( order < 0 ) {
denominator *= std::pow(prime, -order);
}

return Rational(numerator, denominator);
}

// Return a string representation of this P_adic number.
std::string to_string() {
std::vector<uint32_t> numbers = digits;

std::string result = "";
for ( int64_t i = numbers.size() - 1; i >= 0; --i ) {
result += std::to_string(digits[i]);
}

if ( order >= 0 ) {
for ( int32_t i = 0; i < order; ++i ) {
result += "0";
}

result += ".0";
} else {
result.insert(result.length() + order, ".");

while ( result[result.length() - 1] == '0' ) {
result = result.substr(0, result.length() - 1);
}
}

return " ..." + result.substr(result.length() - PRECISION - 1);
}

private:
/**
* Create a P_adic, with p = 'prime', directly from a vector of digits.
*
* For example: with 'order' = 0, the vector [1, 2, 3, 4, 5] creates the p-adic ...54321.0,
* 'order' > 0 shifts the vector 'order' places to the left and
* 'order' < 0 shifts the vector 'order' places to the right.
*/
P_adic(const uint32_t& prime, const std::vector<uint32_t>& digits, const int32_t& order)
: prime(prime), digits(digits), order(order) {
}

// Transform the given vector of digits representing a P_adic number
// into a vector which represents the negation of the P_adic number.
void negate_digits(std::vector<uint32_t>& numbers) {
numbers[0] = modulo_prime(prime - numbers[0]);
for ( uint64_t i = 1; i < numbers.size(); ++i ) {
numbers[i] = prime - 1 - numbers[i];
}
}

// Return the multiplicative inverse of the given number modulo 'prime'.
uint32_t modulo_inverse(const uint32_t& number) const {
uint32_t inverse = 1;
while ( modulo_prime(inverse * number) != 1 ) {
inverse += 1;
}
return inverse;
}

// Return the given number modulo 'prime' in the range 0..'prime' - 1.
int32_t modulo_prime(const int64_t& number) const {
const int32_t div = static_cast<int32_t>(number % prime);
return ( div >= 0 ) ? div : div + prime;
}

// The given vector is padded on the right by zeros up to a maximum length of 'DIGITS_SIZE'.
while ( vector.size() < DIGITS_SIZE ) {
vector.emplace_back(0);
}
}

// Return the given vector of base 'prime' integers converted to a decimal integer.
uint32_t convert_to_decimal(const std::vector<uint32_t>& numbers) const {
uint32_t decimal = 0;
uint32_t multiple = 1;
for ( const uint32_t& number : numbers ) {
decimal += number * multiple;
multiple *= prime;
}
return decimal;
}

// Return whether the given vector consists of all zeros.
bool all_zero_digits(const std::vector<uint32_t>& numbers) const {
for ( uint32_t number : numbers ) {
if ( number != 0 ) {
return false;
}
}
return true;
}

// Return whether the given vector ends with multiple instances of the given number.
bool ends_with(const std::vector<uint32_t>& numbers, const uint32_t& number) const {
for ( uint64_t i = numbers.size() - 1; i >= numbers.size() - PRECISION / 2; --i ) {
if ( numbers[i] != number ) {
return false;
}
}
return true;
}

uint32_t prime;
std::vector<uint32_t> digits;
int32_t order;

static const uint32_t PRECISION = 40;
static const uint32_t ORDER_MAX = 1'000;
static const uint32_t DIGITS_SIZE = PRECISION + 5;
};

int main() {
std::cout << "3-adic numbers:" << std::endl;
std::cout << "-2 / 87    => " << padic_one.to_string() << std::endl;
std::cout << "4 / 97     => " << padic_two.to_string() << std::endl;

std::cout << "sum        => " << sum.to_string() << std::endl;
std::cout << "Rational = " << sum.convert_to_rational().to_string() << std::endl;
std::cout << std::endl;

std::cout << "7-adic numbers:" << std::endl;
std::cout << "5 / 8       => " << padic_one.to_string() << std::endl;
std::cout << "353 / 30809 => " << padic_two.to_string() << std::endl;

std::cout << "sum         => " << sum.to_string() << std::endl;
std::cout << "Rational = " << sum.convert_to_rational().to_string() << std::endl;
std::cout << std::endl;
}
```
Output:
```3-adic numbers:
-2 / 87    =>  ...101020111222001212021110002210102011122.2
4 / 97     =>  ...022220111100202001010001200002111122021.0
sum        =>  ...201011000022210220101111202212220210220.2
Rational = 154 / 8439

5 / 8       =>  ...424242424242424242424242424242424242425.0
353 / 30809 =>  ...560462505550343461155520004023663643455.0
sum         =>  ...315035233123101033613062431266421216213.0
Rational = 156869 / 246472
```

## FreeBASIC

```' ***********************************************
'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

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

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

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

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

sw = a.r2pa(q, 1)
if sw = 1 then exit do
a.printf(0)

sw or= b.r2pa(q, 1)
if sw = 1 then exit do
if sw then continue do
b.printf(0)

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

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
}

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
}

c := 0
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
}

c := 1
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
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
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},
{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

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

p-Adic numbers in this implementation are represented in floating point manner, with p-adic unit as mantissa and p-adic norm as an exponent. The base is encoded implicitly at type level, so that combination of p-adics with different bases won't typecheck.

Textual presentation is given in traditional form with infinite part going to the left.

p-Adic arithmetics and conversion between rationals is implemented as instances of `Eq`, `Num`, `Fractional` and `Real` classes, so, they could be treated as usual real numbers (up to existence of some rationals for non-prime bases).

```{-# LANGUAGE KindSignatures, DataKinds  #-}

import Data.Ratio
import Data.List (genericLength)
import GHC.TypeLits

data Padic (n :: Nat) = Null
| Padic { unit :: [Int], order :: Int }

-- valuation of the base
modulo :: (KnownNat p, Integral i) => Padic p -> i
modulo = fromIntegral . natVal

-- Constructor for zero value
pZero :: KnownNat p => Padic p
pZero = Padic (repeat 0) 0

-- Smart constructor, adjusts trailing zeros with the order.
mkPadic :: (KnownNat p, Integral i) => [i] -> Int -> Padic p
mkPadic u k = go 0 (fromIntegral <\$> u)
where
go 17 _ = pZero
go i (0:u) = go (i+1) u
go i u = Padic u (k-i)

mkUnit :: (KnownNat p, Integral i) => [i] -> Padic p
mkUnit u = mkPadic u 0

-- Zero test (up to 1/p^17)
isZero :: KnownNat p => Padic p -> Bool
isZero (Padic u _) = all (== 0) (take 17 u)
isZero _ = False

pNorm :: KnownNat p => Padic p -> Ratio Int
pNorm Null = undefined
pNorm p = fromIntegral (modulo p) ^^ (- order p)

-- test for an integerness up to p^-17
isInteger :: KnownNat p => Padic p -> Bool
isInteger Null = False
isInteger (Padic s k) = case splitAt k s of
([],i) -> length (takeWhile (==0) \$ reverse (take 20 i)) > 3
_ -> False

-- p-adics are shown with 1/p^17 precision
instance KnownNat p => Show (Padic p) where
show Null = "Null"
show (modulo x) ++ "-adic: " ++
(case si of {[] -> "0"; _ -> si})
++ "." ++
(case f of {[] -> "0"; _ -> sf})
where
(f,i) = case compare k 0 of
LT -> ([], replicate (-k) 0 ++ u)
EQ -> ([], u)
GT -> splitAt k (u ++ repeat 0)
sf = foldMap showD \$ reverse \$ take 17 f
si = foldMap showD \$ dropWhile (== 0) \$ reverse \$ take 17 i
el s = if length s > 16 then "…" else ""
showD n = [(['0'..'9']++['a'..'z']) !! n]

instance KnownNat p => Eq (Padic p) where
a == b = isZero (a - b)

instance KnownNat p => Ord (Padic p) where
compare = error "Ordering is undefined fo p-adics."

instance KnownNat p => Num (Padic p) where
fromInteger 0 = pZero
fromInteger n = pAdic (fromInteger n)

where
k = ka `max` kb
(replicate (k-ka) 0 ++ a)
(replicate (k-kb) 0 ++ b)
_ + _ = Null

mkPadic (mulMod (modulo x) a b) (ka + kb)
_ * _ = Null

case map (\y -> modulo x - 1 - y) u of
[] -> pZero
negate _ = Null

abs p = pAdic (pNorm p)

signum = undefined

------------------------------------------------------------
-- conversion from rationals to p-adics

instance KnownNat p => Fractional (Padic p) where

recip Null = Null
| isZero x = Null
| gcd p u /= 1 = Null
| otherwise = mkPadic res (-k)
where
p = modulo x
res = longDivMod p (1:repeat 0) (u:us)

pAdic :: (Show i, Integral i, KnownNat p)
=> Ratio i -> Padic p
where
p = modulo res
(k, q) = getUnit p x
(n, d) = (numerator q, denominator q)
res = maybe Null process \$ recipMod p d

process r = mkPadic (series n) k
where
series n
| n == 0 = repeat 0
| n `mod` p == 0 = 0 : series (n `div` p)
| otherwise =
let m = (n * r) `mod` p
in m : series ((n - m * d) `div` p)

------------------------------------------------------------
-- conversion from p-adics to rationals
-- works for relatively small denominators

instance KnownNat p => Real (Padic p) where
toRational Null = error "no rational representation!"
toRational x@(Padic s k) = res
where
p = modulo x
res = case break isInteger \$ take 10000 \$ iterate (x +) x of
(_,[]) -> - toRational (- x)
(d, i:_) -> (fromBase p (unit i) * (p^(- order i))) % (genericLength d + 1)

fromBase p = foldr (\x r -> r*p + x) 0 .
take 20 . map fromIntegral

--------------------------------------------------------------------------------
-- helper functions

-- extracts p-adic unit from a rational number
getUnit :: Integral i => i -> Ratio i -> (Int, Ratio i)
getUnit p x = (genericLength k1 - genericLength k2, c)
where
(k1,b:_) = span (\n -> denominator n `mod` p == 0) \$
iterate (* fromIntegral p) x
(k2,c:_) = span (\n -> numerator n `mod` p == 0) \$
iterate (/ fromIntegral p) b

-- Reciprocal of a number modulo p (extended Euclidean algorithm).
-- For non-prime p returns Nothing non-invertible element of the ring.
recipMod :: Integral i => i -> i -> Maybe i
recipMod p 1 = Just 1
recipMod p a | gcd p a == 1 = Just \$ go 0 1 p a
| otherwise = Nothing
where
go t _ _ 0 = t `mod` p
go t nt r nr =
let q = r `div` nr
in go nt (t - q*nt) nr (r - q*nr)

-- Addition of two sequences modulo p
where
go 0 [] ys = ys
go 0 xs [] = xs
go s [] ys = go 0 [s] ys
go s xs [] = go 0 xs [s]
go s (x:xs) (y:ys) =
let (q, r) = (x + y + s) `divMod` p
in r : go q xs ys

-- Subtraction of two sequences modulo p
subMod p a (b:bs) = addMod p a \$ (p-b) : ((p - 1 -) <\$> bs)

-- Multiplication of two sequences modulo p
mulMod p as [b] = mulMod p [b] as
mulMod p as bs = case as of
[0] -> repeat 0
[1] -> bs
[a] -> go 0 bs
where
go s [] = [s]
go s (b:bs) =
let (q, r) = (a * b + s) `divMod` p
in r : go q bs
as -> go bs
where
go [] = []
go (b:bs) =
let c:cs = mulMod p [b] as
in c : addMod p (go bs) cs

-- Division of two sequences modulo p
longDivMod p a (b:bs) = case recipMod p b of
Nothing -> error \$
show b ++ " is not invertible modulo " ++ show p
Just r -> go a
where
go [] = []
go (0:xs) = 0 : go xs
go (x:xs) =
let m = (x*r) `mod` p
_:zs = subMod p (x:xs) (mulMod p [m] (b:bs))
in m : go zs
```

Convertation between rationals and p-adic numbers

```:set -XDataKinds
λ> toRational it
2 % 25
λ> λ> 25 :: Padic 10
Null
λ> toRational it
(-12) % 23```

Arithmetic:

```λ> pAdic (12/25) + pAdic (23/56) :: Padic 7
λ> toRational it
1247 % 1400
λ> 12/25 + 23/56 :: Rational
1247 % 1400
λ> let x = 2/7 :: Padic 13
λ> (2*x - x^2) / 3
> toRational it
8 % 49
λ> let x = 2/7 in (2*x - x^2) / 3 :: Rational
8 % 49```

## Java

This example displays p-adic numbers in standard mathematical format, consisting of a possibly infinite list of digits extending leftwards from the p-adic point. p-adic numbers are given correct to O(prime^40) and the rational reconstruction is correct to O(prime^20).

```import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;

public static void main(String[] args) {
System.out.println("-5 / 9    => " + padicOne);
System.out.println("47 / 12   => " + padicTwo);

System.out.println("sum       => " + sum);
System.out.println("Rational = " + sum.convertToRational());
System.out.println();

System.out.println("5 / 8         => " + padicOne);
System.out.println("353 / 30809   => " + padicTwo);

System.out.println("sum           => " + sum);
System.out.println("Rational = " + sum.convertToRational());
}
}

/**
* Create a p-adic, with p = aPrime, from the given rational 'aNumerator' / 'aDenominator'.
*/
if ( aDenominator == 0 ) {
throw new IllegalArgumentException("Denominator cannot be zero");
}

prime = aPrime;
digits = new ArrayList<Integer>(DIGITS_SIZE);
order = 0;

// Process rational zero
if ( aNumerator == 0 ) {
order = MAX_ORDER;
return;
}

// Remove multiples of 'prime' and adjust the order of the p-adic number accordingly
while ( Math.floorMod(aNumerator, prime) == 0 ) {
aNumerator /= prime;
order += 1;
}

while ( Math.floorMod(aDenominator, prime) == 0 ) {
order -= 1;
}

// Standard calculation of p-adic digits
while ( digits.size() < DIGITS_SIZE ) {
final int digit = Math.floorMod(aNumerator * inverse, prime);

if ( aNumerator != 0 ) {
// The denominator is not a power of a prime
int count = 0;
while ( Math.floorMod(aNumerator, prime) == 0 ) {
aNumerator /= prime;
count += 1;
}

for ( int i = count; i > 1; i-- ) {
}
}
}
}

/**
* Return the sum of this p-adic number and the given p-adic number.
*/
if ( prime != aOther.prime ) {
}

List<Integer> result = new ArrayList<Integer>();

for ( int i = 0; i < -order + aOther.order; i++ ) {
}

for ( int i = 0; i < -aOther.order + order; i++ ) {
}

// Standard digit by digit addition
int carry = 0;
for ( int i = 0; i < Math.min(digits.size(), aOther.digits.size()); i++ ) {
final int sum = digits.get(i) + aOther.digits.get(i) + carry;
final int remainder = Math.floorMod(sum, prime);
carry = ( sum >= prime ) ? 1 : 0;
}

// Reverse the changes made to the digits
for ( int i = 0; i < -order + aOther.order; i++ ) {
aOther.digits.removeFirst();
}

for ( int i = 0; i < -aOther.order + order; i++ ) {
digits.removeFirst();
}

return new Padic(prime, result, allZeroDigits(result) ? MAX_ORDER : Math.min(order, aOther.order));
}

/**
* Return the Rational representation of this p-adic number.
*/
public Rational convertToRational() {
List<Integer> numbers = new ArrayList<Integer>(digits);

// Zero
if ( numbers.isEmpty() || allZeroDigits(numbers) ) {
return new Rational(0, 1);
}

// Positive integer
if ( order >= 0 && endsWith(numbers, 0) ) {
for ( int i = 0; i < order; i++ ) {
}

return new Rational(convertToDecimal(numbers), 1);
}

// Negative integer
if ( order >= 0 && endsWith(numbers, prime - 1) ) {
negateList(numbers);
for ( int i = 0; i < order; i++ ) {
}

return new Rational(-convertToDecimal(numbers), 1);
}

// Rational
int denominator = 1;
do {
denominator += 1;
} while ( ! ( endsWith(sum.digits, 0) || endsWith(sum.digits, prime - 1) ) );

final boolean negative = endsWith(sum.digits, prime - 1);
if ( negative ) {
negateList(sum.digits);
}

int numerator = negative ? -convertToDecimal(sum.digits) : convertToDecimal(sum.digits);

if ( order > 0 ) {
numerator *= Math.pow(prime, order);
}

if ( order < 0 ) {
denominator *= Math.pow(prime, -order);
}

return new Rational(numerator, denominator);
}

/**
* Return a string representation of this p-adic.
*/
public String toString() {
List<Integer> numbers = new ArrayList<Integer>(digits);
Collections.reverse(numbers);
String numberString = numbers.stream().map(String::valueOf).collect(Collectors.joining());
StringBuilder builder = new StringBuilder(numberString);

if ( order >= 0 ) {
for ( int i = 0; i < order; i++ ) {
builder.append("0");
}

builder.append(".0");
} else {
builder.insert(builder.length() + order, ".");

while ( builder.toString().endsWith("0") ) {
builder.deleteCharAt(builder.length() - 1);
}
}

return " ..." + builder.toString().substring(builder.length() - PRECISION - 1);
}

// PRIVATE //

/**
* Create a p-adic, with p = 'aPrime', directly from a list of digits.
*
* With 'aOrder' = 0, the list [1, 2, 3, 4, 5] creates the p-adic ...54321.0
* 'aOrder' > 0 shifts the list 'aOrder' places to the left and
* 'aOrder' < 0 shifts the list 'aOrder' places to the right.
*/
prime = aPrime;
order = aOrder;
}

/**
* Return the multiplicative inverse of the given decimal number modulo 'prime'.
*/
private int moduloInverse(int aNumber) {
int inverse = 1;
while ( Math.floorMod(inverse * aNumber, prime) != 1 ) {
inverse += 1;
}

return inverse;
}

/**
* Transform the given list of digits representing a p-adic number
* into a list which represents the negation of the p-adic number.
*/
for ( int i = 1; i < aDigits.size(); i++ ) {
}
}

/**
* Return the given list of base 'prime' integers converted to a decimal integer.
*/
private int convertToDecimal(List<Integer> aNumbers) {
int decimal = 0;
int multiple = 1;
for ( int number : aNumbers ) {
decimal += number * multiple;
multiple *= prime;
}

return decimal;
}

/**
* Return whether the given list consists of all zeros.
*/
private static boolean allZeroDigits(List<Integer> aList) {
return aList.stream().allMatch( i -> i == 0 );
}

/**
* The given list is padded on the right by zeros up to a maximum length of 'PRECISION'.
*/
private static void padWithZeros(List<Integer> aList) {
while ( aList.size() < DIGITS_SIZE ) {
}
}

/**
* Return whether the given list ends with multiple instances of the given number.
*/
for ( int i = aDigits.size() - 1; i >= aDigits.size() - PRECISION / 2; i-- ) {
return false;
}
}

return true;
}

private static class Rational {

public Rational(int aNumerator, int aDenominator) {
if ( aDenominator < 0 ) {
numerator = -aNumerator;
} else {
numerator = aNumerator;
}

if ( aNumerator == 0 ) {
denominator = 1;
}

final int gcd = gcd(numerator, denominator);
numerator /= gcd;
denominator /= gcd;
}

public String toString() {
return numerator + " / " + denominator;
}

private int gcd(int aOne, int aTwo) {
if ( aTwo == 0 ) {
return Math.abs(aOne);
}
return gcd(aTwo, Math.floorMod(aOne, aTwo));
}

private int numerator;
private int denominator;

}

private List<Integer> digits;
private int order;

private final int prime;

private static final int MAX_ORDER = 1_000;
private static final int PRECISION = 40;
private static final int DIGITS_SIZE = PRECISION + 5;

}
```
Output:
```3-adic numbers:
-5 / 9    =>  ...22222222222222222222222222222222222222.11
47 / 12   =>  ...020202020202020202020202020202020202101.2
sum       =>  ...20202020202020202020202020202020202101.01
Rational = 121 / 36

5 / 8         =>  ...424242424242424242424242424242424242425.0
353 / 30809   =>  ...560462505550343461155520004023663643455.0
sum           =>  ...315035233123101033613062431266421216213.0
Rational = 156869 / 246472
```

## Julia

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

""" convert to Rational (rational reconstruction) """
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

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

## Nim

Translation of: Go

Translation of Go with some modifications, especially using exceptions when an error is encountered.

```import math, strformat

const
Emx = 64        # Exponent maximum.
Dmx = 100000    # Approximation loop maximum.
Amx = 1048576   # Argument maximum.
PMax = 32749    # Prime maximum.

type

Ratio = tuple[a, b: int]

p: int                        # Prime.
k: int                        # Precision.
v: int
d: array[-Emx..(Emx-1), int]

proc r2pa(pa: var Padic; q: Ratio; sw: bool) =
## 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: 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 inverse 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

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

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

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

## 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}"

## String representation.
let t = min(pa.v, 0)
for i in countdown(pa.k - 1 + t, t):
if i == 0 and pa.v < 0: 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],
[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]), true)
print(a, 0)
r2pa(b, (d[4], d[5]), true)
print(b, 0)
echo "+ ="
print(a + b, 1)
echo ""
echo getCurrentExceptionMsg()
```
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
```

## Phix

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

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

integer c = 0
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
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
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
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
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

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

```# 20210225 Raku programming solution

#!/usr/bin/env raku

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 ;
\$x += p**d if \$x < 0 ;  # complement

my \$lowerest = 0;
my (\$num,\$den) = \$x.nude;
while (\$den % p) == 0 { \$den /= p and \$lowerest-- }
\$x = \$num / \$den;

while +self.v < d {
my %d = ^p Z=> (( \$x «-« ^p ) »/» p )».&{ .denominator % p }; # .kv
for %d.keys { self.v.{\$lowerest++} = \$_ and last if %d{\$_} != 0 }
\$x = (\$x - self.v.{\$lowerest-1}) / p ;
}
self
}

my \$div = 0;
my \$lowerest = (self.v.keys.sort({.Int}).first,
x.v.keys.sort({.Int}).first  ).min ;
p => self.p,
v => gather for ^d {
my \$power = \$lowerest + \$_;
given ((self.v.{\$power}//0)+(x.v.{\$power}//0)+\$div).polymod(x.p)
{ take (\$power, .[0]).Slip and \$div = .[1] }
}
}

method gist {
# my %H = (0..9) Z=> ('₀'..'₉'); # (0x2080 .. 0x2089);
# '⋯ ' ~ self.v ~ ' ' ~ [~] self.p.comb».&{ %H{\$_} }

# express as a series
my %H = ( 0…9 ,'-') Z=> ( '⁰','¹','²','³','⁴'…'⁹','⁻');
[~] self.v.keys.sort({.Int}).map: {
' + ' ~ self.v.{\$_} ~ '*' ~ self.p ~ [~] \$_.comb».&{ %H{\$_}} }
}
}

my @T;
for my \D = (
#`[[ these are not working
< 26/25 -109/125 5 4 >,
< 6/7 -5/7 10 7 >,
< 2/7 -3/7 10 7 >,
< 2/7 -1/7 10 7 >,
< 34/21 -39034/791 10 9 >,
#]]
#`[[[[[ Works
< 11/4 679001/207 2 43>,
< 11/4 679001/207 3 27 >,
< 5/19 -101/384 2 12>,
< -22/7 46071/379 7 13 >,
< -7/5 99/70 7 4> ,
< -101/109 583376/6649 61 7>,
< 122/407 -517/1477 7 11>,

< 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 >,
< 49/2 -4851/2 7 6 >,
< -9/5 27/7 3 8>,
< -22/7 46071/379 2 37 >,
< -22/7 46071/379 3 23 >,

< -101/109 583376/6649 2 40>,
< -101/109 583376/6649 32749 3>,
< -25/26 5571/137 7 13>,
#]]]]]

< 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[2] = Padic.new { .r2pa: D[0]+D[1], D[2], D[3] }
say "Addition result = ", \$_.gist; #
unless ( \$_.v.Str eq @T[2].v.Str ) {
say 'but ' ~ (D[0]+D[1]).nude.join('/') ~ ' = ' ~ @T[2].gist
}
}
}
```
Output:
```5/8 =  + 5*7⁰ + 2*7¹ + 4*7² + 2*7³ + 4*7⁴ + 2*7⁵ + 4*7⁶ + 2*7⁷ + 4*7⁸ + 2*7⁹ + 4*7¹⁰
353/30809 =  + 5*7⁰ + 5*7¹ + 4*7² + 3*7³ + 4*7⁴ + 6*7⁵ + 3*7⁶ + 6*7⁷ + 6*7⁸ + 3*7⁹ + 2*7¹⁰
Addition result =  + 3*7⁰ + 1*7¹ + 2*7² + 6*7³ + 1*7⁴ + 2*7⁵ + 1*7⁶ + 2*7⁷ + 4*7⁸ + 6*7⁹ + 6*7¹⁰
```

## Wren

Translation of: FreeBASIC
Library: Wren-dynamic
```import "./dynamic" for Struct

// 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"])

// 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 = P.min(PMAX)  // maximum short prime
K = K.min(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 = _v.min(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 = _v.min(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()
}

+(b) {
var c = 0
r.v = _v.min(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
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],
[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

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