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



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