Algebraic approach to p-adic conversion of rational numbers

Information Processing North-Holland

Letters 18 (1984) 167-171

ALGEBRAIC Alfonso

APPROACH

TO pADIC

30 March 1984

CONVERSION

OF RATIONAL

NUMBERS

MIOLA

Istltufo di Analisi dei Ststemi ed Informatica del C. N. R., 00185 Roma, Italy Communicated by L. Bolliet Received March 1983 Revised July 1983

Keywords:

Exact arithmetic, p-adic arithmetic, approximate rational extended Euclid’s algorithm, continued fraction expansion

arithmetic.

p-adic

representation

of rational

numbers,

1. Introduction The problem of error free computations has been investigated by many authors in the recent past. A natural well-known approach to this problem is based on the use of infinite precision integer and rational arithmetics [ll]. However, this kind of arithmetics is always very space- and time-consuming. A significant alternative that lies between this infinite precision approach and the floating point arithmetic, is the approximate rational arithmetic that assures low complexity operations together with powerful capabilities. In the papers by Horn [lo] and Matula and Kornerup [14] in fact the method used to approximate a real number is based on the algorithm for finding successive convergents (i.e., rational numbers) of the continued fraction expansion of the given number. This method is very well known the antiquity [7], and, since an approximate value of the rational p/q is accurate to within l/q’, it gives similar accuracy as from a double precision arithmetic. Another very interesting approximate rational arithmetic is based on the use of the so-called p-adic arithmetic, first introduced by Hensel [9]. This approach has been followed by Krishnamurthy, Rao and Subramanian [13] and by Hehner and Horspool [S]. The first three authors defined a particular kind of p-adic representation with fixed 0020-0190/84/$3.00

0 1984, Elsevier Science Publishers

length numbers, that they call Hensel-codes. Subsequently Gregory [3,4,5] studied this finite segment p-adic arithmetic and proposed two algorithms for direct and inverse transformations of rational numbers into their corresponding Hensel-codes. However, these algorithms are not very general (i.e., the value of the number p, in the p-adic system, is fixed once and for all in every computation) and rather inefficient (i.e., a table look-up method is used). The main goal of this paper is to present new algorithms for the direct and inverse mappings. We shall recognize the relationship between this p-adic approach and the approximation method based on continued fraction expansion of rational numbers. The study of this relationship will then provide us with a natural and efficient method for computing the homomorphic images of Henselcodes. Then the two new algorithms for direct and inverse transformations will strengthen the p-adic approach for a computationally significant approximate rational arithmetic. In fact, it will now be possible to define efficient algorithms for all the p-adic arithmetic operations, also allowing to change easily the value of p from one computation to another. The new algorithms solve some of the open problems in p-adic approximate arithmetics and therefore they will encourage further research to also succeed the overflow problem. We now can say that the overflow appears at

B.V. (North-Holland)

167

Volume

18, Number

3

INFORMATION

PROCESSING

the end of the computation when the inverse mapping can generate an answer which is not the true one. However, the obtained value and the true one belong to the same residue class. At the same time the basic algebraic background, common to all the approaches that we have mentioned here, will be evident together with the natural similarity between numbers and polynomials. Moreover, this similarity suggests also to investigate the relationship of this approximation method for numbers with the polynomial approximation of functions (described in [2,6,15.16]).

2. p-adic arithmetic The p-adic number arithmetic was first introduced by Hensel [9] and it has been studied by many authors from several points of view (see. for instance, [12]). Let (Y be a rational number, and let p be an integer, that usually is assumed to be prime. The number 01 can be represented by the sequence of digits

with 73

C a,p’

(2)

i=n

and and

O
for i=n, negative

[13,3,81. The a,‘s are computed a, = ;mod

19X4

is represented by the periodic sequence (2 3 1 3 1 . .)5 with n = 1. In fact the computation of the values of a, for a given 1y can be done by successive mod p operations. This p-adic representation leads to a powerful arithmetic that has very useful and significant properties for error free computation [4,&l 21. 2. I. Fixed

length p-u&

numher.s

The computer representation of a p-adic number remains a crucial question, even if these numbers have not infinite sequence of digits. Therefore, we would like to deal only with fixed length p-adic numbers, namely p-adic sequence as in (1). only with, say, r p-adic digits. Since a rational number (Y can be uniquely expressed as Q

=

!t

pll,

d with p and n as in (2) and c, d and p pairwise relatively prime, it follows from (2) that

p,

a, =

Then the number (Ycan be expressed (as proposed in [13,3]) by the value e,, = n that represents its exponent, and by the finite sequence of the first r p-adic digits of representation (4) corresponding to ncr-1

n+l, n+2 . . . . . a,,#0 or zero, as proposed in

m, =

as follows:

that is a sort of mantissa. For a given p we can then represent cx by the so-called Hensel-code.

a-

a,_,

C I=”

alp’-“.

(5)

mod p

P

fori=n+l, n+2,.... Expression (2) is the p-adic representation of cy and the values a, for i = n, n + 1,. . , are called p-adic digits. It must be quoted here that sequence (1) is either finite (namely a, = 0 for all i greater than a certain value) or periodic. For instance, we can easily see that for p = 5 the number 44 is represented by the finite sequence (4 3 l), with n = 0. and the number l/15 168

30 March

(1)

a,a,+,a,+,...

a=

LETTERS

H(P. r. a> = (m,?, a).

(6)

In order now to deal with this kind of representation we have to devise algorithms: (i) to compute efficiently the digits a,. (ii) to reconstruct a rational number cy from its Hensel-code (6). Gregory [5] proposed two methods for (i) and (ii). The first simply noticed that the value of m, as in (5), for given (Y, p and r. can be obtained as m,, = (c/d)

mod pr.

(7)

Volume

18, Number

3

INFORMATION

PROCESSING

The method proposed for the second problem recognizes that a unique answer is only possible if some constraint is posed on the expected rational number. Gregory suggested to solve problem (ii) for a so-called order-N Farey fraction: namely a rational number a/b such that 0 < ]a] < N and 0 < b < N for given N. Under this assumption he devised a table look-up method to find the desired rational number among previously computed values. Unfortunately, it happens that both these methods suffer for a lack of a more appropriate algebraic approach, and more dramatically the proposed method for the inverse mapping looses a crucial property of the p-adic approach; namely, the possibility to repeat any computation with different values of p, or more generally to let the value of p varying from one problem to the next.

3. Algebraic background: The extended Euclid’s algorithm, continued fractions, diophantine equations

The computational properties of infinite precision integer arithmetic of modular and p-adic arithmetics have been deeply studied in the literature, together with all their algebraic bases [l,ll]. It is also very important to underline the fundamental contributions of computer algebra research in the understanding of the integers and polynomials similarity, especially in the case of p-adic computation. In this sense the work by Gustavson and Yun [6,16] and by Wang [15] draws the exact shape of how to handle the problem of approximation by the p-adic approach. From this algebraic algorithmic point of view a basic observation can be done on all the arithmetics proposed as mentioned in Section 1. The computing methods used by the authors we have cited, make a very intensive use of the extended Euclid’s algorithm, as is easy to see. That allows us to stress the unified view of all these methods and also to recognize the natural relation of these methods with the fundamental number theory of the continued fraction approximation [7,11].

LETTERS

30 March 1984

We will now give a short overview of these basic algebraic and number-theoretic elements in order to devise a more appropriate solution of the two problems stated in the previous paragraph. The extended Euclid’s algorithm (EEA) [l,ll], integers, computes x for given a,, a, nonnegative and y, as well as the greatest common divisor (gcd) of a, and a,, such that a,x + a,y = d = gcd{a,, The computation

a,}.

is carried

(8) out by using

three

vectors (a,, a, ,..., a,), (x0, x,,...,x,) and (yO, y,, . . , y,) such that at each step the following relations hold: a ,+, = a,-, - qa, xi+1 = xi_1 - qxi Yi+ I

=YI-1

for given a,, a,, for x0 = 1, x, = 0,

where q = [ai_,/ai], a, = 0 with a,-, a, }, and also such that the property a,xi

(9)

fory,=O,y,=l,

-9Yi

= gcd{a,,

+ a,y, = a,

(10)

holds for all i=O, l,..., n- 1. The values y,, _ 1 are then the solution of (8). Eq. (8) has the general form a,x + a,y = c,

x,-,,

(11)

that is known as a diophantine equation in the integer unknowns x and y, for given integers a,,, a,, c. If gcd{a,, a, } divides c, then the solution X, 4 of (11) is related to the solution x,_ ,, y,_, of (8) (see, for instance, [9]) in the following way: X = mx,_,

+ ka,/gcd{a,,

a,},

Y = my,-l

+ ka,/gcd{a,,,

a,},

for any integer

(12)

k and m such that c = m gcd{a,,

a, 1. On the other hand the EEA and solution (12) of (11) have also interesting relation with the continued fraction expansion of numbers. In fact the problem of solving (8) can also be formulated as the calculation of approximation to a,/a, since, in (8),

~C’=$<’ X

1

X2

(13) 169

Volume

18. Number

3

INFORMATION

if we assume 0 < x -C a, [7]. Therefore -y/x is an approximation of a,/a, that is quite close to so/a, and has a very small denominator. From this observation we can also verify that the successive values of q in the EEA are the successive terms of the continued fraction expansion of a,Ja, : 1 a,/a,

= ql +

(14)

1 qz+

q,+

The successive convergents of a,/a, are then the values -y/x, for i 2 0. with x,, y, computed throughout the EEA as in (9) (see [7,11]). These unified number-theoretic elements can now be used to solve the two problems posed in Section 2.1. 4. Conversion

of Hensel-codes

Theorem

1. Given (Y = c/d and pr such that c, d and pr are pairwise relatively prime, the value of m, as in (5) is computed by the EEA applied to pr and d and it is

m, = cy mod pr,

(15)

where y is the second output of the EEA.

unknown

(16)

for k integer. This equation is a diophantine equation solvable by the EEA with solutions as in (12). In this case we have

c and

30 March

d, for given

mcr and

lYX4

p’, with

m,

< P’. The EEA with pr and m, as input will compute the sequence of x, an y, and therefore the continued fraction convergents of m,,/p’, according to the results recalled in the previous section. If we now assume that 2cd < pr and we look for a rational number that is an order-N Farey fraction with N = dp’/2, then the approximation leads to a unique value c/d. In fact the following claims can be proved. Lemma.

Jf

1 ‘
P --x Iq

then p/q

is a convergent

Proof. For the proof.

2. If c/d

m
to x.

see [7, Section

10.71.

•I

is such that

pr,

(18)

with 0 < ICI< N, 0 < d < N, and N = Jp’/2. then, .since EEA applied to p’ and rn_ computes the convergents of m ,/p’, there exists an i such that the i th convergent determines exact(y c/d = a ,/y,. Proof. Expression

(18) can be rewritten

as

!T+!+’ P

Proof. The value of m, as in (5) can be obtained as a solution of (7) and also be written as

m,d = c

LETTERS

Theorem

Let cy be a rational number. The computation of m, as in (5) can be done according to the following theorem.

kp’+

PROCESSING

(19)

dp’

for some k, and -k/d m,/p’ with an error of

is an approximation

to

(20)

Y = cy,, , + kp’,

(17)

and the assertion follows by recalling is less than pr. q

that here m,

We can now apply the lemma to prove that -k/d is a convergent of m,/pr and it can be computed by EEA (which computes all the convergents). Then the value c/d can be obtained by EEA as a/y, for the i such that ]y,] < N. q

The inverse mapping, namely to recover cx from its Hensel-code, consists of solving (15) for the

Using this theorem we can easily define an algorithm (a slight variation of the EEA) to solve our inverse mapping problem, and to make availa-

X = cx,_,

170

+ kd,

Volume 18. Number

INFORMATION

3

PROCESSING

ble a significant approximate arithmetic that maintains the power of the p-adic approach as specified in [8] and the strong algebraic and number-theoretical bases.

111A.V. Aho, J.E. Hopcroft

I31

Let us compute now, for example, the expression H(p, r, (w), as in (6), for the rational number CY= 2/7, with p = 5 and r = 4, according to Theorem 1 of Section 4. In this case we have to solve the equation. k. 625 + m, .7 = 2.

[41 ISI WI

(21) [71

By applying the EEA we do have the following sequences as in (9): (625, 7, 2, 1, 0) for the ai’s and (0, 1, - 89, 268, - 625) for the yi’s. Then the solution of (21) with respect to m, is

PI

m,, = (2.268)

[91

mod 625 = 536,

(22)

which corresponds to 4121 base 5. The Hensel-code we are looking for will then be given by H(5,4,

2/7)

= (0.1214,

0).

DOI 1111

(23)

Let us now examplify the inverse mapping procedure by computing the appropriate rational number c/d, with Hensel-code as in (23), via Theorem 2. Now we have that the value of N is 17, and if we apply the EEA to p’ = 625 and 4121, = 536 we obtain the following sequences as in (9): (625, 536, 89, 2, 1, 0) for the a,‘s and (0, 1, - 1, 7, -309, 625) for the yi’s. Then the expected ratio: a/y, = 2/7, is obtained at the second step of the EEA.

30 March 1984

References

PI 5. Example of direct computation and inverse mapping

LETTERS

1121 u31

[I41

VW

[I61

and J.D. Ullman, The Design and Analysis of Computer Algorithms (Addison-Wesley, Reading, MA, 1974). C. Brezinski, Pad&type approximation and general orthogonal polynomials, ISNM Vol. 50 (Birkhaiiser, Basel, 1980). R.T. Gregory, The use of finite-segment p-adic arithmetic for exact computations, BIT 18 (1978) 282-300. R.T. Gregory, Error-Free Computation (Krieger, Huntington, WV, 1980). R.T. Gregory. Error-free computation with rational numbers, BIT 21 (1981) 194-202. F.G. Gustavson and D.Y.Y. Yun, Fast computation of Pade approximations and Toeplitz systems of equations via the extended Euclidean algorithm, IBM Res. Rept. RC7551, 1979. G.H. Hardy and E.M. Wright, An Introduction to the Theory of Numbers (Clarendon Press, Oxford, 4th ed., 1960). E.C.R. Hehner and R.N.S. Horspool, A new representation of the rational numbers for fast easy arithmetic, SIAM J. Comput. 8 (1979) 124-134. K. Hensel, Theorie der Algebraischen Zahlen (Teubner, Leipzig, 1908). B.K.P. Horn, Rational arithmetic for minicomputers, Software Prac. Exper. 8 (1978) 171-176. D.E. Knuth, The Art of Computer Programming: Vol. II, Seminumerical Algorithms (Addison-Wesley, Reading, MA, 1978). N. Koblitz, p-Adic Numbers, p-Adic Analysis, and Zeta Functions (Springer. New York, 1977). E.V. Krishnamurthy, T.M. Rao and K. Subramanian, Finite segment p-adic systems with applications to exact computation, Proc. Indian Acad. Sci. 81A (1975) 58-79. D.W. Matula and P. Kornerup, Approximate rational arithmetic systems: Analysis of recovery of simple fractions during expression evalution, Lecture Notes in Comput. Sci. 72 (Springer, Berlin, 1979). P.S. Wang, A p-adic algorithm for univariate partial fractions, Proc. 1981 ACM Symp. on Symbolic and Algebraic Computation (ACM Inc., New York, 1981). D.Y.Y. Yun, Algebraic algorithms using p-adic constructions, Proc. ACM SYMSAC ‘76. 1976.

171