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



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