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



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"])
 
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 = 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()
    }
 
    // add b to 'this'
    +(b) {
        var c = 0
        var r = Padic.new()
        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
        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()
}