Rational Arithmetic

Roger Hui




Introduction

Computing with rational numbers opens new vistas. For example:

      +∘÷ \ 10 ⍴ 1
1 2 1.5 1.66667 1.6 1.625 1.61538 1.61905 1.61765 1.61818

      ⍕Q (+Q)∘(÷Q) \ 10 ⍴ qi 1
1 2 3r2 5r3 8r5 13r8 21r13 34r21 55r34 89r55

We presents a facility for extended-precision rational arithmetic: a representation for rational numbers, a function qi (“rational input”) to specify rational numbers conveniently, and a monadic operator Q for working with rational numbers,

This text, the slides, and the code are available for download from the Dyalog ’20 webpage.

An earlier version of this text appeared in the Dyalog APL Chat Forum post Rational Arithmetic [Hui 2020]. A more detailed description of the implementation can be found in Rational Arithmetic Implementation Notes [Hui 2020b] and in Rational Arithmetic — RNG [Hui 2020c]. The original design was changed, informed by the analysis in A Design Exercise [Hui 2020d].

The operator “rats”in the dfns workspace [Scholes 2014] also implements rational arithmetic; the current facility is more general and more capable. There was a Conference Edition of Dyalog APL distributed at the 2011 Dyalog User Conference [Hui 2011], in which rational arithmetic was implemented as part of the interpreter. The facility was never released as part of the official interpreter due to difficulties with making exact and inexact numbers play nice together.
 

Representation

A rational number is a scalar of an enclosed 3-element vector, whose items are:

• a vector of the decimal digits of the numerator
a vector of the decimal digits of the denominator
the sign, ¯1 , 0 , or 1

The numerator and the denominator are relatively prime. The vectors of digits are required to have no leading 0s (except for rational 0 , represented as ⊂(,¨0 1),0). Therefore, the representation is unique: rational numbers p and q are equal if and only if p≡q .

Because a rational number so represented is a scalar, “structural” ( , , etc.) and “selection” ( , [;] , , etc.) functions just work on arrays of rational numbers.

For example:

   Hilbert←{÷1+∘.+⍨⍳⍵}

   Hilbert 5          ⍝ an array of floating point numbers
1        0.5      0.333333 0.25     0.2     
0.5      0.333333 0.25     0.2      0.166667
0.333333 0.25     0.2      0.166667 0.142857
0.25     0.2      0.166667 0.142857 0.125   
0.2      0.166667 0.142857 0.125    0.111111

   qi Hilbert 5       ⍝ an array of rational numbers
┌───────┬───────┬───────┬───────┬───────┐
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│1│1│││1│2│1│││1│3│1│││1│4│1│││1│5│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│2│1│││1│3│1│││1│4│1│││1│5│1│││1│6│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│3│1│││1│4│1│││1│5│1│││1│6│1│││1│7│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│4│1│││1│5│1│││1│6│1│││1│7│1│││1│8│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
├───────┼───────┼───────┼───────┼───────┤
│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│┌─┬─┬─┐│
││1│5│1│││1│6│1│││1│7│1│││1│8│1│││1│9│1││
│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│└─┴─┴─┘│
└───────┴───────┴───────┴───────┴───────┘
      ⍕Q qi Hilbert 5    ⍝ rational formatted output
  1 1r2 1r3 1r4 1r5
1r2 1r3 1r4 1r5 1r6
1r3 1r4 1r5 1r6 1r7
1r4 1r5 1r6 1r7 1r8
1r5 1r6 1r7 1r8 1r9

The Function qi

The function qi (“rational input”) provides a convenient way to specify rational numbers. It applies independently to items of the argument; each item can be one of several types:

•   an ordinary number per APL, e.g. 1234, ¯2.75, 1.234e¯10
the decimal digits of the numerator, e.g. 1 2 3 4 5 6 7 8 (the denominator is assumed to be 1)
a character vector of the number per APL, but without the usual limits, e.g. '1.2345e6789' . This is a way to specify exactly a large or tiny rational number.

qi has an optional left argument of the comparison tolerance to use (if different from the current ⎕ct value) in rationalizing a floating point number.

Examples:

   qi ¯2.75
┌──────────┐
│┌───┬─┬──┐│
││1 1│4│¯1││
│└───┴─┴──┘│
└──────────┘
   ⍕Q qi ¯2.75
¯11r4

   qi ⊂1 2 3 4 5 6 7 8
┌─────────────────────┐
│┌───────────────┬─┬─┐│
││1 2 3 4 5 6 7 8│1│1││
│└───────────────┴─┴─┘│
└─────────────────────┘
   ⍕Q qi ⊂1 2 3 4 5 6 7 8
12345678

   qi ○1      ⍝ ⎕ct=1e¯14
┌───────────────────────────────┐
│┌─────────────┬─────────────┬─┐│
││5 4 1 9 3 5 1│1 7 2 5 0 3 3│1││
│└─────────────┴─────────────┴─┘│
└───────────────────────────────┘
   18 ⍕ 5419351÷1725033
 3.141592653589815___

   0 qi ○1    ⍝ ⎕ct=0
┌─────────────────────────────────────┐
│┌─────────────────┬───────────────┬─┐│
││2 4 5 8 5 0 9 2 2│7 8 2 5 6 7 7 9│1││
│└─────────────────┴───────────────┴─┘│
└─────────────────────────────────────┘
   18 ⍕ 245850922÷78256779
 3.141592653589793___

   18 ⍕ ○1
 3.141592653589793___

   t←qi ⊂'1.2345e6789'
   ⍴¨¨t        ⍝ # decimal digits in numerator & denominator
┌─────────┐
│┌────┬─┬┐│
││6790│1│││
│└────┴─┴┘│
└─────────┘
   ¯10 ⍕Q t    ⍝ exponential format
1.234500000e6789

The Operator Q

f Q derives a function which works on rational numbers and produces a rational result except as noted otherwise.

Scalar Functions

• monadic and dyadic: + - × ÷ ⌊ ⌈ | !
dyadic only: * ∨ ∧ < ≤ = ≠ ≥ >
monadic only:?

If f and g are dyadic functions in the above list, then f Q/ , f Q⌿ , f Q\ , f Q⍀ , (f Q).(g Q) , and ∘.(f Q) also work on rational arguments.

The boolean functions < ≤ = ≠ ≥ > return a boolean result, an ordinary 0 or 1 instead of a rational 0 or 1 . Likewise ×⍵ returns an ordinary ¯1 , 0 , or 1 .
 

Non-Scalar Functions

⊤Q ⍵ applies to rational scalar positive integer and produces its binary representation. ( does not have a monadic definition. One possibility is ⊤⍵ ←→ 2⊤⍵ , that is, bits or binary representation. (Similarly, ⊥⍵ would be 2⊥⍵ , binary value.) This is in fact the definition used in A Formal Description of SYSTEM/360 [Falkoff et al., 1964, Table 1]. It was actually implemented in APL\360 for a time but was later removed due to space limitations [Hui 1992]. For example, ⊤Q qi 123 ←→ 1 1 1 1 0 1 1 . ⍺ ⊤Q ⍵ applies to rational scalar positive integer and produces its base-representation, for an ordinary integer . For example, 5 ⊤Q qi 1234567 ←→ 3 0 4 0 0 1 2 3 2 .

⍕Q ⍵ formats rational array , with an r separating the numerator and the denominator (if not 1). ⍺ ⍕Q ⍵ formats rational arrayin fixed-point (if ⍺≥0) or exponential (if ⍺<0) format. The left argumentmust be a single integer: ⍺≥0 specifies the number of digits after the decimal point; ⍺<0 specifies the number of decimal digits in exponential notation.

⌹Q is analoguous tofor rational arguments. As usual, a DOMAIN ERROR is signalled if is non-singular (but not as usual, if the error is signalledreally is singular).
 

Applications

Computing with rational numbers opens new vistas.
 

Continued Fractions

   +∘÷\ 16⍴1
1 2 1.5 1.66667 1.6 1.625 1.61538 1.61905 1.61765 1.61818 
      1.61798 1.61806 1.61803 1.61804 1.61803 1.61803 

   ⍕Q (+Q)∘(÷Q)\ 16 ⍴ qi 1
1 2 3r2 5r3 8r5 13r8 21r13 34r21 55r34 89r55 144r89 233r144 
      377r233 610r377 987r610 1597r987 

How about that: our friends the Fibonacci numbers.

   A001203← 3 7 15 1 292 1 1 1 2 1 3 1 14  ⍝ OEIS A001203
   30 ⍕Q ⍪ (+Q)∘(÷Q)\ qi A001203
3.000000000000000000000000000000
3.142857142857142857142857142857
3.141509433962264150943396226415
3.141592920353982300884955752212
3.141592653011902604072261494774
3.141592653921421044708715941593
3.141592653467436705520454785349
3.141592653618936623397500301411
3.141592653581077771204419306582
3.141592653591403978482542414219
3.141592653589389171543687321707
3.141592653589815383241943777307
3.141592653589792659375626945712

Square Root via Newton’s Method

If W is a number whose square root is sought andis an approximation of the root, Newton’s method improves the approximation by computing the average of , and W divided by . By using rational arithmetic, a sequence of rational approximations obtains.

sqrt←{  ⍝ apply Newton’s method ⍺ times to find sqrt ⍵
 q2←qi 2
 W←⍵
 f←{q2÷Q⍨⍵+Q W÷Q ⍵}
 ⍺ {1≥⍺:⍵ ⋄ (⍺-1)∇ ⍵,f ¯1↑⍵} qi 1
}

   ⍕Q t←8 sqrt qi 2
1 3r2 17r12 577r408 665857r470832 886731088897r627013566048 
      1572584048032918633353217r1111984844349868137938112 ...

The following code displays the rational approximations to 40 decimal places, and the differences between 2 and the squares of the approximations. The differences illustrate that Newton’s method doubles the number of correct digits with each iteration.

   (40 ⍕Q ⍪t), ' ', ¯6 ⍕Q ⍪ (qi 2) -Q t ×Q t
1.0000000000000000000000000000000000000000  1.00000e0  
1.5000000000000000000000000000000000000000 ¯2.50000e¯1 
1.4166666666666666666666666666666666666667 ¯6.94444e¯3 
1.4142156862745098039215686274509803921569 ¯6.00730e¯6 
1.4142135623746899106262955788901349101166 ¯4.51095e¯12
1.4142135623730950488016896235025302436150 ¯2.54358e¯24
1.4142135623730950488016887242096980785697 ¯8.08728e¯49
1.4142135623730950488016887242096980785697 ¯8.17550e¯98

Matrix Inverse

Hilbert’s matrix is known to be invertible, but is notorious for being nearly singular.

      A← qi Hilbert 7
      ⍕Q A
  1 1r2 1r3  1r4  1r5  1r6  1r7
1r2 1r3 1r4  1r5  1r6  1r7  1r8
1r3 1r4 1r5  1r6  1r7  1r8  1r9
1r4 1r5 1r6  1r7  1r8  1r9 1r10
1r5 1r6 1r7  1r8  1r9 1r10 1r11
1r6 1r7 1r8  1r9 1r10 1r11 1r12
1r7 1r8 1r9 1r10 1r11 1r12 1r13

      B← ⌹Q A
      ⍕Q A (+Q).(×Q) B
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1

Inverting the 7-by-7 Hilbert matrix as 64-bit floating point numbers:

   X←Hilbert 7
   Y←⌹X

   ¯3 ⍕ X +.× Y
 1.00E0    1.46E¯11  1.16E¯10 4.66E¯10  9.31E¯10 ¯9.31E¯10  0.00E0  
¯9.09E¯13  1.00E0   ¯3.49E¯10 2.33E¯9   9.31E¯10 ¯4.66E¯10  0.00E0  
 0.00E0   ¯7.28E¯12  1.00E0   9.31E¯10  9.31E¯10 ¯1.40E¯9   4.66E¯10
 2.27E¯13 ¯2.18E¯11 ¯1.16E¯10 1.00E0    1.86E¯9  ¯2.33E¯9   2.33E¯10
 4.55E¯13 ¯2.91E¯11 ¯5.82E¯11 4.66E¯10  1.00E0   ¯2.79E¯9   5.82E¯10
 0.00E0    2.91E¯11 ¯3.49E¯10 9.31E¯10 ¯1.40E¯9   1.00E0   ¯4.66E¯10
¯4.55E¯13 ¯1.46E¯11 ¯1.16E¯10 9.31E¯10 ¯1.86E¯9   0.00E0    1.00E0

How does it look for a matrix of random numbers?

   ⎕rl←6*5    ⍝ for reproducible random numbers
   X←¯10+?6 6⍴25

   ¯3 ⍕ X +.× ⌹X
  1.00E0    3.55E¯15 ¯2.84E¯14 ¯1.78E¯14  2.49E¯14 ¯2.22E¯15
  5.55E¯17  1.00E0    2.84E¯14  2.13E¯14 ¯7.11E¯15  0.00E0  
 ¯6.94E¯17  3.55E¯15  1.00E0   ¯2.13E¯14  2.13E¯14  2.11E¯15
  2.78E¯17  3.55E¯15  7.11E¯15  1.00E0    0.00E0   ¯8.88E¯16
 ¯1.39E¯16 ¯1.07E¯14 ¯1.42E¯14  1.42E¯14  1.00E0    8.88E¯16
  3.47E¯16 ¯9.77E¯15 ¯1.78E¯14  2.13E¯14 ¯2.31E¯14  1.00E0  

   A← qi X
   ⍕Q A (+Q).(×Q) ⌹Q A
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1


Addition examples of using rational numbers can be found in the APL Chat Forum posts modpower [Hui 2020e] and Miller-Rabin Primality Testing [Hui 2020f].
 

In Closing …

A more detailed description of the rational arithmetic implementation can be found in Rational Arithmetic Implementation Notes [Hui 2020b] and in Rational Arithmetic — RNG [Hui 2020c]. The original design was changed, informed by the analysis in A Design Exercise [Hui 2020d].

There was a “Conference Edition” of Dyalog APL distributed at the 2011 Dyalog User Conference [Hui 2011], in which rational arithmetic was implemented as part of the interpreter. The facility was never released as part of the official interpreter due to difficulties with making exact and inexact numbers play nice together.

The operator “rats”in the dfns workspace [Scholes 2014] also implements rational arithmetic; the current facility is more general and more capable.
 

Postscript

We have shown that it is possible to compute rational approximations to non-rational functions. Another approach is to have floating point arithmetic to a specified number of digits (“big floats”). The following examples are from work in progress.

MAXDIGITS specifies the number of significant floating point decimal digits. Like the rational arithmetic facility, there is a specification for the representation, a function fi (“floating point input”) to specify big floats conveniently, and a monadic operator F for working with big floats. The examples illustrate the computation of Euler’s number (*1) using a Taylor series. (The number of terms in the series, 60, is estimated by !⍣¯1 .)

   MAXDIGITS ← 80

   fi 3.14159e¯20
┌───────────────────┐
│┌───────────┬───┬─┐│
││3 1 4 1 5 9│¯25│1││
│└───────────┴───┴─┘│
└───────────────────┘

   ¯20 ⍕ +/ ÷ ! ⍳60
 2.718281828459046____E0

   ¯80 ⍕F +F/ ÷F !F fi ⍳60
2.71828182845904523536028747135266249775724709369995957 … 476e0
   ' ',⎕d[A001113]   ⍝ OEIS A001113
 271828182845904523536028747135266249775724709369995957 … 4759 …

   !⍣¯1 ⊢ 1e80
58.92

References

•  Falkoff, A.D., K.E. Iverson, E.H. Sussenguth, A Formal Description of SYSTEM/360, IBM Systems Journal, volume 3, number 3, 1964.
Hui, R.K.W., Josephus Problem, J Wiki Essay, 2005-10-17.
Hui, R.K.W., Rational Numbers, Dyalog User Conference 2011, 2011-10-03.
Hui, R.K.W., Rational Arithmetic, Dyalog APL Chat Forum, 2020-07-07. (a)
Hui, R.K.W., Rational Arithmetic Implementation Notes, Dyalog APL Chat Forum, 2020-07-13. (b)
Hui, R.K.W., Rational Arithmetic — RNG, Dyalog APL Chat Forum, 2020-07-17. (c)
Hui, R.K.W., A Design Exercise, Dyalog APL Chat Forum, 2020-07-30. (d)
Hui, R.K.W., modpower, Dyalog APL Chat Forum, 2020-08-06. (e)
Hui, R.K.W., Miller-Rabin Primality Testing, Dyalog APL Chat Forum, 2020-08-20. (f)
Scholes, J.M., Rational Number, Dfns, 2014-08-06.


created:  2020-07-07 21:20
updated: