rats ← {tolerance←⎕ct} ##.rational nums     ⍝ Rational approximation to real ⍵.

Argument [nums] is a numeric,  possibly nested, array of shape S.  The result is
an array of shape: (2,S) with pairs of whole numbers, along the first axis, with
a quotient which is  ⍺-tolerably  close  to  the  corresponding  argument  item.
[rational] guarantees to return the smallest pair of such numbers.  For example,
if we choose pi (an irrational number) as argument, then with the default ⎕CT of
1e¯14, we see:

        rational ○1             ⍝ tolerable rational approximation to pi
    5419351 1725033

        (○1) = ÷/ rational ○1   ⍝ quotient of pair is ⎕CT-equal to argument.
    1
        ⎕ct←1e¯10               ⍝ coarser tolerance.

        rational ○1             ⍝ coarser tolerance yields smaller pair.
    208341 66317

        ⊢ numbs ← +/¨1⊂ 2 2 2⍴ +\10*-⍳8         ⍝ nested numeric array.
    ┌─────────────────┬───────────────────┐
    │0.1     0.111    │0.11     0.1111    │
    │0.11111 0.1111111│0.111111 0.11111111│
    └─────────────────┴───────────────────┘

        disp rational numbs                     ⍝ rational pairs.
    ┌───────────────┬─────────────────┐
    │    1     111  │    11     1111  │
    │11111 1111111  │111111 11111111  │
    ├───────────────┼─────────────────┤
    │    10     1000│    100     10000│
    │100000 10000000│1000000 100000000│
    └───────────────┴─────────────────┘

        numbs ≡ ÷⌿ rational numbs               ⍝ round-trip
    1

The  integer  pair  could  be reconstituted directly from the continued fraction
returned  by  function →cfract←. The inner reduction in the coding that follows,
represents the expression transformation:

  ···+1     ···+1        ···+1        ···+1         ┐ ···+ ┌ yz+1
      ─         ─            ─            ─         │      │ ─────────  => ···
      x+1       x+1    ┐   ┌ x+z    ┐   ┌ x(yz+1)+z ├───→──┤ x(yz+1)+z
        ─         ─    ├─→─┤   ──── ├─→─┤ ───────── │      └
      ┌ y+1 ┐   ┌ yz+1 │   └   yz+1 ┘   └   yz+1    ┘
      │   ─ ├─→─┤ ──── │
      └   z ┘   └   z  ┘


    rational←{              ⍝ Rational approximation from continued fraction.
        ⌽↑{                 ⍝ reversed accumulation of,
            top bot←⍵       ⍝ tailmost,
            ⌽top 0+⍺ 1×bot  ⍝ fractions.
        }/(cfract ⍵),1      ⍝ continued fraction representation.
    }

or the equivalent one-liner:

    rational←{⌽↑{↑⍺{⌽⍺ 0+⍺⍺ 1×⍵}/⍵}/(cfract ⍵),⊂1 1}

For  amusement, [rational]  could be recast as a derived function, although it's
hard to see a practical benefit of doing so:

    rational ← ↑∘(÷∘⊂∘(1∘∨)⍨∘(1∘(,⍨∘⊂))⍨)       ⍝ derived fn

Historical note
---------------
Prior to the implementation of primitive function GCD (∨) this more complex cod-
ing, which was slower by a factor of around 2, was used:

    rational←{                  ⍝ Rational number near real ⍵.
        ⍺←⎕CT ⋄ ⎕CT←⍺           ⍝ default comparison tolerance.
        real←⍵                  ⍝ "real" number.
        real{                   ⍝ starting with real number.
            An←⌊⍺               ⍝ nth term of continued fraction.
            Cn←1 An+.×⍵         ⍝ nth convergent pair.
            real=÷/Cn:Cn        ⍝ tolerably close rational: done.
            (÷⍺-An)∇ 1 0↓⍵⍪Cn   ⍝ otherwise: next term.
        }2 2⍴0 1 1 0            ⍝ convergents: C¯2, C¯1.
    }

A Perfectly Accurate (⎕ct-intolerant) Version
---------------------------------------------
If we ignore comparison tolerance, all IEEE floating point numbers are rational,
with a denominator of a power of 2 (yeah?).

Using →hexf←, this version decodes the bits from the internal IEEE double float-
ing-point representation and uses →nats← to provide a perfectly (⎕CT=0) accurate
rational  pair.  The  two-item  result is a pair of character vectors of decimal
digits.

    ratf←{⎕IO ⎕CT←0                     ⍝ Exact rational from IEEE double.
        digs←16↑⎕D,⎕A                   ⍝ hex digits.
        bits←,⍉2 2 2 2⊤digs⍳hexf ⍵      ⍝ binary floating number.
        split←(⍳⍴bits)∊0 1 12           ⍝ split fields: sign exponent mantissa.
        bsign bexp bmant←split⊂bits     ⍝ bit-vector fields.
        sign←bsign/'¯'                  ⍝ negative sign.
        exp←¯1022+2⊥bexp                ⍝ signed numeric binary exponent.
        deco←{↑⍺{⍺+nats ⍺⍺×nats ⍵}/⌽⍵}  ⍝ accurate ⍺-decode.
        top←2 deco 1,bmant              ⍝ numerator.
        bot←2*nats 53-exp               ⍝ denominator.
        sign '',¨top{                   ⍝ rational pair, char vectors à la nats.
            1∊'13579'∊,↑¯1↑¨⍺ ⍵:⍺ ⍵     ⍝ either number is odd: done.
            (⍺÷nats 2)∇ ⍵÷nats 2        ⍝ both even: cancel 2s.
        }bot
    }

Notice that some friendly decimal numbers such as 0.1 are not representable with
perfect  accuracy in binary.  This means that the closest number to 0.1 that may
be represented in 53 bits, is the rational number:

    3602879701896397 ÷ 36028797018963968
0.1

    ratf 1234                           ⍝ whole number.
┌────┬─┐
│1234│1│
└────┴─┘

    ratf 0.1                            ⍝ ÷10 is not a finite binary number.
┌────────────────┬─────────────────┐
│3602879701896397│36028797018963968│
└────────────────┴─────────────────┘

    ratf 2*¯16                          ⍝ reciprocal of power of 2 is spot-on.
┌─┬─────┐
│1│65536│
└─┴─────┘

    ratf ¯0.75                          ⍝ negative sign fixed to numerator.
┌──┬─┐
│¯3│4│
└──┴─┘

Notice  that  the denominator of the rational pair returned by ratf is always an
integer power of 2.

Ref: http://en.wikipedia.org/wiki/Continued_fraction

Examples:

    rational 0.75             ⍝ three quarters.
3 4
    rational ÷3               ⍝ one third.
1 3

    rational ¯0.1234          ⍝ negative real => negative numerator.
¯617 5000

    rational 0                ⍝ 0 => 0÷1 (this is handy for →cfract←).
0 1

    1e¯10 rational 2*÷2       ⍝ coarse rational approximation to sqrt(2).
114243 80782

    ⎕←pi←rational ○1          ⍝ tolerable approxmiation to pi.
5419351 1725033

    (○1)=÷/pi                 ⍝ ... compares within ⎕ct.
1
    gcd/pi                    ⍝ pair is relatively prime.
1
    gcd/⎕←rational 9÷16       ⍝ pair of relatively prime numbers:
9 16
1
    rational 11÷29            ⍝ pair matches rational of quotient ...
11 29
                              ⍝ ... for all relatively prime pairs:

    {⍵≡rational÷/⍵}¨ ∘.,⍨ sieve 2 to 100
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0

    ⍉ rational +∘÷\ cfract ○1   ⍝ successive rational approximations to Pi.
      3       1
     22       7
    333     106
    355     113
 103993   33102
 104348   33215
 208341   66317
 312689   99532
 833719  265381
1146408  364913
5419351 1725033

⍝ Looking at the third row above, we see that if we divide the number-of-the-
⍝ beast (666) by the NYC area code (212) we get a reasonable appoximation to Pi.

See also: factors sieve gcd cfract efract pco
See also: hexf nats
See also: numbers

Back to: contents

Back to: Workspaces