(12C) Decimal to Fraction

[Albert Chan]

(08-08-2018 06:05 PM)Thomas Klemm Wrote:  

Code:

def dec2frac(f, c):
    if c < 2 and int(f) == f: return
    u, v, F = 0, 1, f
    while v < c and f:
        g = 1 / f
        f, _ = modf(g)
        u, v = v, u + v * int(g)
    if v > c: v -= u * ceil((v - c) / u)
    q = u if diff(F, u) < diff(F, v) else v
    p = round(q * F, 0)
    return p, q

Hi, Thomas

I noticed some issues with above Python code.

dec2frac(2.0, 10) crashes with ZeroDivisionError. "if c < 2 and" should be "if c < 2 or"
-> since it is translated from calculator code, you might want to check this.

Although not a bug in Python, when code running in a calculator, v can be huge.
-> v might be much bigger than c, corrected v below c might be off.

Below is the patched code:

Code:

def dec2frac(f, c):
    if c < 2 or int(f) == f: return
    u, v, F = 0, 1, f
    while f:
        f = 1/f
        w = u + v * int(f)  # check before update
        if w > c: break     # guaranteed v <= c
        u, v, f = v, w, f - int(f)
    q = u + (c - u)//v * v  # biggest quotient
    if diff(F, q) >= diff(F, v): q = v
    return int(q * F + 0.5), q

[Thomas Klemm]

(08-08-2018 10:26 PM)Albert Chan Wrote:  dec2frac(2.0, 10) crashes with ZeroDivisionError. "if c < 2 and" should be "if c < 2 or"
-> since it is translated from calculator code, you might want to check this.

Nah, Joe is of course correct:

Code:

\<< DUP2 @ Must be two arguments.  Exit now if max denominator < 2,
  IF 1 > SWAP FP AND @ or if decimal fraction is an integer.

It's only me that should study De Morgan's laws.

Quote:Although not a bug in Python, when code running in a calculator, v can be huge.
-> v might be much bigger than c, corrected v below c might be off.

Below is the patched code:

Code:

def dec2frac(f, c):
    if c < 2 or int(f) == f: return
    u, v, F = 0, 1, f
    while f:
        f = 1/f
        w = u + v * int(f)  # check before update
        if w > c: break     # guaranteed v <= c
        u, v, f = v, w, f - int(f)
    q = u + (c - u)//v * v  # biggest quotient
    if diff(F, q) >= diff(F, v): q = v
    return int(q * F + 0.5), q

Not sure if I understand that correctly but that's something Joe might want to change in his RPL program as well. So this is rather a "jump forwards"?

Thanks for the feedback
Thomas

[Thomas Klemm]

(08-08-2018 09:44 PM)Dieter Wrote:  Yes - as already mentioned:

Oopsie Woopsie. I should really read your posts more faithfully.

Cheers
Thomas

[Gamo]

This routine is to check for the FIX setting.

If 0.000 the result will show 3

This might be helpful to use in the decimal to fraction program.

Code:


01  3
02 1/x
03 RND
04 LSTx
05  -
06 ENTER
07  x
08 √x
09 1/x
10 LN
11  1
12  0
13 LN
14  ÷
15 INTG

Gamo

[Thomas Klemm]

This program for the HP-11C includes Dieter's and Albert's suggestions:

Code:

  *LBL A        ; c f
   STO 1        ; c f
   CLx          ; 0 f
   STO 2        ; u=0 f
   R↓           ; f
   STO 0        ; f
   1            ; 1 f
   STO 3        ; v=1 f
  *LBL 0        ; while
   R↓           ; f
   x=0          ; f = 0 ?
   GTO 1        ; done
   1/x          ; 1/f
   FRAC         ; f'=frac(1/f)
   RCL 2        ; u f'
   LSTx         ; 1/f u f'
   INT          ; ⌊1/f⌋ u f'
   RCL 3        ; v ⌊1/f⌋ u f'
   STO 2        ; u'=v
   ×            ; v*⌊1/f⌋ u f'
   x<>y         ; u v*⌊1/f⌋ f'
   +            ; v'=u+v*⌊1/f⌋ f'
   RCL 1        ; c v' f'
   x≤y          ; c ≤ v' ?
   GTO 1        ; done
   R↓           ; v' f'
   STO 3        ; v' f'
   GTO 0        ; while
  *LBL 1        ; done
   LSTx         ; u c
   STO 2        ; u c
   -            ; c-u
   RCL 3        ; v c-v
   ÷            ; (c-u)/v
   INT          ; ⌊(c-u)/v⌋
   RCL 3        ; v ⌊(c-u)/v⌋
   ×            ; v*⌊(c-u)/v⌋
   STO + 2      ; u += v*⌊(c-u)/v⌋
   RCL 2        ; u
   GSB 4        ; diff
   RCL 3        ; v
   GSB 4        ; diff
   x≤y
   GTO 2
   RCL 2        ; u
   GTO 3
  *LBL 2
   RCL 3        ; v
  *LBL 3        ; round
   ENTER        ; d d
   ENTER        ; d d d
   RCL 0        ; f d d
   ×            ; f*d d
   .            ;
   5            ; 0.5 f*d d
   +            ; f*d+0.5 d
   INT          ; round(f*d) d
   RTN          ; n d
  *LBL 4        ; diff
   GSB 3        ; n d
   x<>y         ; d n
   ÷            ; n/d
   RCL 0        ; f n/d
   -            ; n/d-f
   ABS          ; |n/d-f|
   RTN          ; diff

[Albert Chan]

Hi Thomas,

I was in the middle of typing this, when I see you fixed the code.
Anyway, below illustrated why I request the fix.

Try 0.10322, approx. with fraction, denominator <= 1e8

Using Mathematica, Continued Fraction coefficient = [0,9,1,2,4,1,6,2,1,4,3]

With HP-12C, due to rounding error, I get [0,9,1,2,4,1,6,2,1,4,2,1,64,1,4,3,1,74139, ...]

-> original setup, last 2 quotients = 68773261, 5098833198335
-> the big quotient rounded down in HP-12C to 5098833198000

-> interpolated quotient = 5098833198000 - 74139 * 68773261 = 52401000 (HP-12C rounded)
-> fraction = 5408831 / 52401000 = 0.103219996

The fraction is wrong.
Even with bad continued fraction coefficients, value should be very close to 0.10322

This is what the patched version does to above:
-> fraction = 5408837 / 52401056 = 0.1032200000

[Dieter]

(08-09-2018 12:51 AM)Gamo Wrote:  This routine is to check for the FIX setting.

If 0.000 the result will show 3

This might be helpful to use in the decimal to fraction program.

I cannot see how this may be related to the topic of this thread, but anyway it's a nice little utility that you may post in a separate thread.

The program calculates round(1/3) – 1/3. This is always negative, so there is no need for an ABS via ENTER × √x: a CHS will do. Or simply reverse the order: 1/3 – round(1/3) is always positive. I assume that after this usually the base-10-logarithm is calculated, which the 12C does not offer. But instead of dividing the natural log by ln 10 (which is slow) here any value between 2,193 and 2,436 may be used. For instance √5 = 2,236. This is faster than a log and saves one step.

So here is an updated version. I added an ENTER at the start, just to be sure...

Code:

01 ENTER
02 3
03 1/x
04 ENTER
05 RND
06 -
07 1/x
08 LN
09 5
10 √x
11 ÷
12 INTG
13 GTO 00

FIX 4 [R/S] => 4,0000

Re. line 07/08: LN CHS may by a tiiiny bit faster. ;-)

Dieter

[Dieter]

(08-09-2018 02:17 AM)Albert Chan Wrote:  Try 0.10322, approx. with fraction, denominator <= 1e8
...
This is what the patched version does to above:
-> fraction = 5408837 / 52401056 = 0.1032200000

That's why I think the algorithm should stop if the fraction matches the input.
My 35s program returns 5161/50000.

Dieter

[Gamo]

0.10322 on HP Prime also returns 5161/50000

Gamo

[Dieter]

(08-08-2018 10:26 PM)Albert Chan Wrote:  Below is the patched code:

Code:

def dec2frac(f, c):
    if c < 2 or int(f) == f: return
    u, v, F = 0, 1, f
    while f:
        f = 1/f
        w = u + v * int(f)  # check before update
        if w > c: break     # guaranteed v <= c
        u, v, f = v, w, f - int(f)
    q = u + (c - u)//v * v  # biggest quotient
    if diff(F, q) >= diff(F, v): q = v
    return int(q * F + 0.5), q

I have transcoded this for the 35s. Let me add a few remarks:

First I think that q can also be calculated via q = c  –  (c–u) mod v. Can you confirm this?
Also I do not quite understand why earlier code (and Joe's as well) used the CEIL function while here it's the FLOOR (or INT) function. Thomas even calculated 1+INT() which in most cases is the same as CEIL... but not for integer arguments.

Then I wonder how exactly this Python diff() command works. Its arguments here are just the decimal input and a denominator – but not the corresponding numerator. So how does this work in detail? For instance, what is the value for diff(pi, 106) and diff(pi, 113)?

In a calculator implementation first the numerators for the two possible results (denominators v and q) are calculated and then the errors are compared. These errors are evaluated by abs(input – numerator/denominator). This may lead to inaccuracies due to digit cancellation. Think of the √2 example (577/408 vs. 816/577). Is there a better way to compare the two possible fractions?

Finally the code currently rejects integer input (which would lead to a division by zero). But can't this be trapped? For instance by leaving the while-loop if f=0? Or simply return v=1 if F is an integer?

Dieter

[Thomas Klemm]

Maybe I missed something but is there a reason for not using Farey sequences?

Python

Code:

def farey(q, n):
    a = 0
    b = c = d = 1
    while True:
        e, f = a+c, b+d
        k = f * q
        if e == k:
            return e, f
        if e < k:
            a, b = e, f
        else:
            c, d = e, f
        if b+d > n:
            return e, f

HP-11C

Code:

LBL A
STO 1       ; n q
R↓          ; q
STO 0       ; q
CLx         ; 0
STO 2       ; a = 0
1           ; 1
STO 3       ; b = 1
STO 4       ; c = 1
STO 5       ; d = 1
RCL 3       ; 1 1
+           ; 2
LBL 0       ; while
STO 7       ; f = b+d
RCL 0       ; q f
*           ; k
RCL 2       ; a k
RCL 4       ; c a k
+           ; a+c k
STO 6       ; e=a+c k
x=y         ; e = k ?
GTO 3       ; done
x≤y?        ; e < k ?
GTO 1
RCL 6       ; e
STO 4       ; c = e
RCL 7       ; f
STO 5       ; d = f
GTO 2
LBL 1
RCL 6       ; e
STO 2       ; a = e
RCL 7       ; f
STO 3       ; b = f
LBL 2
RCL 1       ; n
RCL 3       ; b n
RCL 5       ; d b n
+           ; b+d n
x≤y         ; b+d ≤ n ?
GTO 0       ; while
LBL 3       ; done
RCL 6       ; e
RCL 7       ; f
RTN

This should also be easier to convert for the HP-12C since no subroutines are used.

However the decimal number \(q\) to transform must be \(0 < q < 1\).

Example

0.10322 ENTER EEX 8 [A]
Y: 5,161
X: 50,000


Cheers
Thomas

[Thomas Klemm]

(08-09-2018 09:27 AM)Dieter Wrote:  Then I wonder how exactly this Python diff() command works. Its arguments here are just the decimal input and a denominator – but not the corresponding numerator. So how does this work in detail?

I assume that this is just the same function that I used:

Code:

def diff(f, d):
    n = round(d * f, 0)
    return abs(f - n / d)

Quote:For instance, what is the value for diff(pi, 106) and diff(pi, 113)?

>>> from math import pi
>>> diff(pi, 106)
8.32196275291075e-05
>>> diff(pi, 113)
2.667641894049666e-07


HTH
Thomas

[Dieter]

(08-09-2018 09:30 AM)Thomas Klemm Wrote:  Maybe I missed something but is there a reason for not using Farey sequences?
...
However the decimal number \(q\) to transform must be \(0 < q < 1\).

This should be easy to handle (use 1/input if input>1 and swap numerator and denominator afterwards).
But this also means that you define the largest possible numerator. #-)

I tried this algorithm on the 12C. Here is what I got for pi, using 1/pi as the input:

3,141592653 [1/x] 500 => 113/355 (as expected)
3,141592653 [1/x] 600 => 113/355 (as expected)
3,141592653 [1/x] 700 => 219/688 (?!?)

What's this?

Dieter

[Thomas Klemm]

(08-09-2018 10:30 AM)Dieter Wrote:  What's this?

Oh, I see. It's just that: \(\frac{113}{355} < \frac{1}{\pi} < \frac{219}{688}\)

Thus we can't simply return \(e, f\):

Code:

        if b+d > n:
            return e, f

Instead we should check whether \(\frac{a}{b}\) or \(\frac{c}{d}\) is closer to \(\frac{1}{\pi}\).

Here's the output of: a, b, c, d

Code:

>>> farey(1/pi, 1000)
0 1 1 2
0 1 1 3
1 4 1 3
2 7 1 3
3 10 1 3
4 13 1 3
5 16 1 3
6 19 1 3
7 22 1 3
7 22 8 25
7 22 15 47
7 22 22 69
7 22 29 91
7 22 36 113
7 22 43 135
7 22 50 157
7 22 57 179
7 22 64 201
7 22 71 223
7 22 78 245
7 22 85 267
7 22 92 289
7 22 99 311
7 22 106 333
113 355 106 333
113 355 219 688
(219, 688)

As we know \(\frac{113}{355}\) is closer for a long time. We have to wait until \(\frac{33215}{104348}\) to get closer:

Code:

>>> farey(1/pi, 200000)
0 1 1 2
0 1 1 3
1 4 1 3
2 7 1 3
3 10 1 3
4 13 1 3
5 16 1 3
6 19 1 3
7 22 1 3
7 22 8 25
7 22 15 47
7 22 22 69
7 22 29 91
7 22 36 113
7 22 43 135
7 22 50 157
7 22 57 179
7 22 64 201
7 22 71 223
7 22 78 245
7 22 85 267
7 22 92 289
7 22 99 311
7 22 106 333
113 355 106 333
113 355 219 688
113 355 332 1043
113 355 445 1398
113 355 558 1753
113 355 671 2108
113 355 784 2463
113 355 897 2818
113 355 1010 3173
113 355 1123 3528
113 355 1236 3883
113 355 1349 4238
113 355 1462 4593
113 355 1575 4948
113 355 1688 5303
113 355 1801 5658
113 355 1914 6013
113 355 2027 6368
113 355 2140 6723
113 355 2253 7078
113 355 2366 7433
113 355 2479 7788
113 355 2592 8143
113 355 2705 8498
113 355 2818 8853
113 355 2931 9208
113 355 3044 9563
113 355 3157 9918
113 355 3270 10273
113 355 3383 10628
113 355 3496 10983
113 355 3609 11338
113 355 3722 11693
113 355 3835 12048
113 355 3948 12403
113 355 4061 12758
113 355 4174 13113
113 355 4287 13468
113 355 4400 13823
113 355 4513 14178
113 355 4626 14533
113 355 4739 14888
113 355 4852 15243
113 355 4965 15598
113 355 5078 15953
113 355 5191 16308
113 355 5304 16663
113 355 5417 17018
113 355 5530 17373
113 355 5643 17728
113 355 5756 18083
113 355 5869 18438
113 355 5982 18793
113 355 6095 19148
113 355 6208 19503
113 355 6321 19858
113 355 6434 20213
113 355 6547 20568
113 355 6660 20923
113 355 6773 21278
113 355 6886 21633
113 355 6999 21988
113 355 7112 22343
113 355 7225 22698
113 355 7338 23053
113 355 7451 23408
113 355 7564 23763
113 355 7677 24118
113 355 7790 24473
113 355 7903 24828
113 355 8016 25183
113 355 8129 25538
113 355 8242 25893
113 355 8355 26248
113 355 8468 26603
113 355 8581 26958
113 355 8694 27313
113 355 8807 27668
113 355 8920 28023
113 355 9033 28378
113 355 9146 28733
113 355 9259 29088
113 355 9372 29443
113 355 9485 29798
113 355 9598 30153
113 355 9711 30508
113 355 9824 30863
113 355 9937 31218
113 355 10050 31573
113 355 10163 31928
113 355 10276 32283
113 355 10389 32638
113 355 10502 32993
113 355 10615 33348
113 355 10728 33703
113 355 10841 34058
113 355 10954 34413
113 355 11067 34768
113 355 11180 35123
113 355 11293 35478
113 355 11406 35833
113 355 11519 36188
113 355 11632 36543
113 355 11745 36898
113 355 11858 37253
113 355 11971 37608
113 355 12084 37963
113 355 12197 38318
113 355 12310 38673
113 355 12423 39028
113 355 12536 39383
113 355 12649 39738
113 355 12762 40093
113 355 12875 40448
113 355 12988 40803
113 355 13101 41158
113 355 13214 41513
113 355 13327 41868
113 355 13440 42223
113 355 13553 42578
113 355 13666 42933
113 355 13779 43288
113 355 13892 43643
113 355 14005 43998
113 355 14118 44353
113 355 14231 44708
113 355 14344 45063
113 355 14457 45418
113 355 14570 45773
113 355 14683 46128
113 355 14796 46483
113 355 14909 46838
113 355 15022 47193
113 355 15135 47548
113 355 15248 47903
113 355 15361 48258
113 355 15474 48613
113 355 15587 48968
113 355 15700 49323
113 355 15813 49678
113 355 15926 50033
113 355 16039 50388
113 355 16152 50743
113 355 16265 51098
113 355 16378 51453
113 355 16491 51808
113 355 16604 52163
113 355 16717 52518
113 355 16830 52873
113 355 16943 53228
113 355 17056 53583
113 355 17169 53938
113 355 17282 54293
113 355 17395 54648
113 355 17508 55003
113 355 17621 55358
113 355 17734 55713
113 355 17847 56068
113 355 17960 56423
113 355 18073 56778
113 355 18186 57133
113 355 18299 57488
113 355 18412 57843
113 355 18525 58198
113 355 18638 58553
113 355 18751 58908
113 355 18864 59263
113 355 18977 59618
113 355 19090 59973
113 355 19203 60328
113 355 19316 60683
113 355 19429 61038
113 355 19542 61393
113 355 19655 61748
113 355 19768 62103
113 355 19881 62458
113 355 19994 62813
113 355 20107 63168
113 355 20220 63523
113 355 20333 63878
113 355 20446 64233
113 355 20559 64588
113 355 20672 64943
113 355 20785 65298
113 355 20898 65653
113 355 21011 66008
113 355 21124 66363
113 355 21237 66718
113 355 21350 67073
113 355 21463 67428
113 355 21576 67783
113 355 21689 68138
113 355 21802 68493
113 355 21915 68848
113 355 22028 69203
113 355 22141 69558
113 355 22254 69913
113 355 22367 70268
113 355 22480 70623
113 355 22593 70978
113 355 22706 71333
113 355 22819 71688
113 355 22932 72043
113 355 23045 72398
113 355 23158 72753
113 355 23271 73108
113 355 23384 73463
113 355 23497 73818
113 355 23610 74173
113 355 23723 74528
113 355 23836 74883
113 355 23949 75238
113 355 24062 75593
113 355 24175 75948
113 355 24288 76303
113 355 24401 76658
113 355 24514 77013
113 355 24627 77368
113 355 24740 77723
113 355 24853 78078
113 355 24966 78433
113 355 25079 78788
113 355 25192 79143
113 355 25305 79498
113 355 25418 79853
113 355 25531 80208
113 355 25644 80563
113 355 25757 80918
113 355 25870 81273
113 355 25983 81628
113 355 26096 81983
113 355 26209 82338
113 355 26322 82693
113 355 26435 83048
113 355 26548 83403
113 355 26661 83758
113 355 26774 84113
113 355 26887 84468
113 355 27000 84823
113 355 27113 85178
113 355 27226 85533
113 355 27339 85888
113 355 27452 86243
113 355 27565 86598
113 355 27678 86953
113 355 27791 87308
113 355 27904 87663
113 355 28017 88018
113 355 28130 88373
113 355 28243 88728
113 355 28356 89083
113 355 28469 89438
113 355 28582 89793
113 355 28695 90148
113 355 28808 90503
113 355 28921 90858
113 355 29034 91213
113 355 29147 91568
113 355 29260 91923
113 355 29373 92278
113 355 29486 92633
113 355 29599 92988
113 355 29712 93343
113 355 29825 93698
113 355 29938 94053
113 355 30051 94408
113 355 30164 94763
113 355 30277 95118
113 355 30390 95473
113 355 30503 95828
113 355 30616 96183
113 355 30729 96538
113 355 30842 96893
113 355 30955 97248
113 355 31068 97603
113 355 31181 97958
113 355 31294 98313
113 355 31407 98668
113 355 31520 99023
113 355 31633 99378
113 355 31746 99733
113 355 31859 100088
113 355 31972 100443
113 355 32085 100798
113 355 32198 101153
113 355 32311 101508
113 355 32424 101863
113 355 32537 102218
113 355 32650 102573
113 355 32763 102928
113 355 32876 103283
113 355 32989 103638
113 355 33102 103993
33215 104348 33102 103993
(33215, 104348)

Sorry for the slip.
Thomas

[Albert Chan]

(08-09-2018 09:27 AM)Dieter Wrote:  I think that q can also be calculated via q = c  –  (c-u) mod v. Can you confirm this?

q = u + (c-u)//v * v
= c - ((c-u) - (c-u)//v * v)
= c - (c-u) mod v -- confirmed

Quote:I do not quite understand why earlier code (and Joe's as well) used the CEIL function why here it's the FLOOR (or INT) function.

Thomas code use ceiling because v were over-shoot passed c:
q = v - k u <= c
k u >= v - c
k >= (v - c) / u
k = ceiling((v-c) / u) -- to maximize q

My patch, both u, v are below c:
q = u + k v <= c
k <= (c - u) / v
k = (c - u) // v -- to maximize q

Quote:Thomas even used 1+INT() which in most cases is the same as CEIL... but not for integer arguments.

1+INT() for ceiling is wrong
A simple example is 0.00005, u, v = 1, 20000
Since u = 1, 1+INT() always over-shoot ceiling by 1

[Albert Chan]

(08-09-2018 09:30 AM)Thomas Klemm Wrote:  Maybe I missed something but is there a reason for not using Farey sequences?
...
However the decimal number \(q\) to transform must be \(0 < q < 1\).

Perhaps Continued Fraction route is faster (and less restriction) ?

(08-09-2018 10:30 AM)Dieter Wrote:  This should be easy to handle (use 1/input if input>1 and swap numerator and denominator afterwards).
But this also means that you define the largest possible numerator. #-)

Is it better do let Farey sequence handle fraction of input ?
Numerator can be easily fixed later.

Probably not even need fixing, mixed fraction is a "feature" :-)

[Albert Chan]

I removed \(0 < x < 1\) restriction of farey fraction.

Code:

>>> farey(0.10322, 1e8)   # fractions match, stop
(5161, 50000)
>>> farey(-0.10322, 1e8)  # negatives ok
(-5161, 50000)
>>> farey(1/pi, 700)      # pick the best estimate
(113, 355)
>>> farey(pi, 100)        # abs(x) > 1 ok
(311, 99)
>>> farey(3, 100)         # integer ok
(3, 1)

Code:

def farey(x, limit=100):
    'return Farey fraction (n,d) closest to x with denominator <= limit'
    n1, d1 = int(x), 1      
    if n1 == x: return n1, d1
    if x < 0: n1 -= 1       # lo = floor(x)
    n2, d2 = n1+1, 1        # hi = lo + 1
    while 1:
        n, d = n1 + n2, d1 + d2
        if d > limit: break
        if   x * d > n: n1, d1 = n, d
        elif x * d < n: n2, d2 = n, d
        else: return n, d   # exact match
    pick_hi = (n2/d2 - x) < (x - n1/d1)
    return (n2, d2) if pick_hi else (n1, d1)

[Thomas Klemm]

(08-09-2018 03:39 PM)Albert Chan Wrote:  I removed \(0 < x < 1\) restriction of farey fraction.

Very nice! This is my attempt to translate your program for the HP-11C:

Code:

*LBL A
 STO 1       ; limit x
 R↓          ; x
 STO 0       ; x
 INT         ; ⌊x⌋
 STO 2       ; lo_n
 1           ; 1 lo_n
 STO 3       ; lo_d=1
 STO 4       ; hi_n=1
 R↓          ; lo_n
 RCL 0       ; x lo_n
 x=y         ; x = lo_n ?
 GTO 4       ; return lo_n lo_d
 Clx         ; 0
 STO 5       ; hi_d=0
*LBL 0       ; while
 RCL 2       ; lo_n
 RCL 4       ; hi_n lo_n
 +           ; lo_n+hi_n
 STO 6       ; n=lo_n+hi_n
 RCL 1       ; limit
 RCL 3       ; lo_d limit
 RCL 5       ; hi_d lo_d limit
 +           ; lo_d+hi_d limit
 STO 7       ; d=lo_d+hi_d limit
 x>y         ; d > limit ?
 GTO 3       ; break
 RCL 6       ; n
 RCL 0       ; x n
 RCL 7       ; d x n
 *           ; x*d n
 x≤y         ; x*d ≤ n ?
 GTO 1
 RCL 6       ; n
 STO 2       ; lo_n
 RCL 7       ; d
 STO 3       ; lo_d
 GTO 0       ; while
*LBL 1       ; x*d n
 x=y         ; x*d = n ?
 GTO 2       ; return n d
 RCL 6       ; n
 STO 4       ; hi_n
 RCL 7       ; d
 STO 5       ; hi_d
 GTO 0       ; while
*LBL 2       ; return n d
 RCL 6       ; n
 RCL 7       ; d
 RTN
*LBL 3       ; break
 RCL 0       ; x
 RCL 4       ; hi_n x
 RCL 5       ; hi_d hi_n x
 /           ; hi_n/hi_d x
 -           ; x-hi_n/hi_d
 ABS         ; ∆_hi=|x-hi_n/hi_d|
 RCL 0       ; x ∆_hi
 RCL 2       ; lo_n x ∆_hi
 RCL 3       ; lo_d lo_n x ∆_hi
 /           ; lo_n/lo_d x ∆_hi
 -           ; x-lo_n/lo_d ∆_hi
 ABS         ; ∆_lo=|x-lo_n/lo_d| ∆_hi
 x≤y         ; ∆_lo ≤ ∆_hi ?
 GTO 4       ; return lo_n lo_d
 RCL 4       ; hi_n
 RCL 5       ; hi_d hi_n
 RTN
*LBL 4
 RCL 2       ; lo_n
 RCL 3       ; lo_d lo_n
 RTN

Quote:

Code:

>>> farey(-0.10322, 1e8)  # negatives ok
(-5161, 50000)

That's not implemented yet. So I'm leaving this as an exercise for the dear reader.

Kind regards
Thomas

[Albert Chan]

(08-09-2018 06:38 PM)Thomas Klemm Wrote:  

(08-09-2018 03:39 PM)Albert Chan Wrote:  I removed \(0 < x < 1\) restriction of farey fraction.

Very nice! This is my attempt to translate your program for the HP-11C

I had revised my Python code after your HP-11C translation. Now:
1. hi = lo + 1 instead of 1/0 (infinity), thus saved 1 iteration
2. Both abs() calls are removed

For reference, this was how I removed the restriction (see Farey Pair): Continued Fractions without Tears

[Dieter]

(08-09-2018 06:38 PM)Thomas Klemm Wrote:  Very nice! This is my attempt to translate your program for the HP-11C:

Service: here is a translation for the 12C.

Not available tests have been substituted, n_max is stored in the n-register, so that R1 is free and R2...R7 can move down to R1...R6. The output order has been reversed so that the numerator is returned in X and the denominator in Y.

Code:

01 STO 0
02 INTG
03 STO 1
04 1
05 STO 2
06 STO 3
07 RCL 0
08 FRAC
09 x=0?
10 GTO 66
11 Clx
12 STO 4
13 RCL 1
14 RCL 3
15 +
16 STO 5
17 RCL n
18 RCL 2
19 RCL 4
20 +
21 STO 6
22 x≤y?
23 GTO 25
24 GTO 47
25 RCL 5
26 RCL 0
27 RCL 6
28 *
29 x≤y?
30 GTO 36
31 RCL 5
32 STO 1
33 RCL 6
34 STO 2
35 GTO 13
36 -
37 x=0?
38 GTO 44
39 RCL 5
40 STO 3
41 RCL 6
42 STO 4
43 GTO 13
44 RCL 6
45 RCL 5
46 GTO 00
47 RCL 0
48 RCL 3
49 RCL 4
50 /
51 -
52 ENTER
53 x
54 RCL 0
55 RCL 1
56 RCL 2
57 /
58 -
59 ENTER
60 x
61 x≤y?
62 GTO 66
63 RCL 4
64 RCL 3
65 GTO 00
66 RCL 2
67 RCL 1
68 GTO 00

The max. denominator is set with the [n] key.
Output: numerator [X↔Y] denominator.

99 [n]
3,141592654 [R/S] => 311 [X↔Y] 99

2 [√x] [R/S] => 140 [X↔Y] 99

1 [e^x] [R/S] => 193 [X↔Y] 71

20 [n]
1 [e^x] [R/S] => 49 [X↔Y] 18

700 [n]
2 [√x] [R/S] => 816 [X↔Y] 577

999 [n]
3,141592654 [R/S] => 355 [X↔Y] 113

Now I wonder how long all this takes on an original hardware 12C... ;-)

Note:
You may replace line 23 with an x=0? test. This yields an x≤y? x=0? combination, and since x is always > 0 at this point the result logically is the desired x>y? test. ;-)

Dieter