pair ← {tolerance} ##.rational numb ⍝ Rational number near real ⍵.
The result is the "best" pair of integers whose quotient is (⎕CT-tolerably) equ-
al to original argument. In this context, "best" means smallest. For example, if
we choose pi (an irrational number) as argument, then with default ⎕CT of 1e¯14,
we see:
⎕ct←1e¯14 ⍝ default value of ⎕CT.
rational ○1 ⍝ tolerable rational approximation to pi
5419351 1725033
(○1) = ÷/ rational ○1 ⍝ quotient of pair is ⎕CT-equal to original.
1
A coarser tolerance results in a smaller pair of numbers.
⎕ct←1e¯10 ⍝ coarser tolerance
rational ○1 ⍝ → smaller pair.
208341 66317
(○1)=÷/ rational ○1 ⍝ quotient of pair is always ⎕CT-equal to original.
1
(
NB: In Dyalog V11.0 and above we can do better by using primitive function ∨
(GCD). The function is pervasive in that for a (possibly nested) argument of
shape (S), it returns a result of shape (2,S), with rational pairs along the
first axis:
rational ← {⎕ml←0 ⋄ ↑⍵ 1÷⊂1∨⍵} ⍝ rational (V11 coding) (V)
⎕← 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↓
└~─────────────→┴~───────────────→┘
)
Technical note:
The 2-item result vector is generated from a continued fraction representation
of the scalar numeric argument. See also →cfract←.
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.
}
See: http://en.wikipedia.org/wiki/Continued_fraction for details.
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}
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 ⎕ML ⎕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.
bsign '',¨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
disp ratf 1234 ⍝ whole number.
┌→───┬─┐
│1234│1│
└───→┴→┘
disp ratf 0.1 ⍝ ÷10 is not a finite binary number.
┌→───────────────┬─────────────────┐
│3602879701896397│36028797018963968│
└───────────────→┴────────────────→┘
disp ratf 2*¯16 ⍝ reciprocal of power of 2 is spot-on.
┌→┬─────┐
│1│65536│
└→┴────→┘
disp 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.
Examples:
rational 0.75 ⍝ three quarters.
3 4
rational ÷3 ⍝ one third.
1 3
rational ¯0.1234 ⍝ negative real => negative numerator.
¯617 5000
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
See also: factors sieve gcd cfract efract
See also: hexf nats
Back to: contents
Back to: Dyalog APL
Trouble seeing APL font?