Companion program for Joe Horn’s continued fractions article

[dlidstrom]

Hi

This is a companion program for Joe Horn’s article on continued fractions.

Code:


« → a b
  « a →SF →CF b →SF →CF → c d
    « { }
      DO
        IF c HEAD d HEAD DUP2 ≠
        THEN MIN 1 +
        ELSE DROP
        END + c HEAD d HEAD c TAIL 'c' STO d TAIL 'd' STO
      UNTIL ≠
      END CF→
    »
  »
»

You’ll need to type I the programs from the article too.
Give this program two decimal fractions according to the article, and you’ll get back the smallest fraction that evaluates between them.

The article itself can be found here: https://www.hpmuseum.org/forum/thread-7984.html

[John Keith]

For HP 49/50 users, there is a faster System RPL version here. The program converts in both directions, from a fraction to a list or vice versa.

Edit: I just noticed a couple of typos in the linked article. On page 1, near the bottom, the number 777 should be 177.

[dlidstrom]

Hi Keith, thanks for the tip. My program is fast too! It evaluates 3.1415 and 3.1416 to 333/106 in just 1.25s. Being User RPL it is also safe, super easy to inspect, debug, and modify. I kind of prefer User RPL these days myself.

[Joe Horn]

(09-09-2023 08:37 AM)dlidstrom Wrote:  Give this program two decimal fractions according to the article, and you’ll get back the smallest fraction that evaluates between them.

Excellent program, but it seems that some inputs are problematic. E.g. with inputs of 1.4 and 1.5, the program errors out, saying "HEAD Error: Invalid Dimension". It also seems to choke if both inputs are equal.

[dlidstrom]

Hi Joe! Thanks for chiming in and doing some stress tests!
Indeed it fails as it assumes the arrays are same length. I’ll need to think about a fix for those issues.

[Joe Horn]

(09-09-2023 03:10 PM)dlidstrom Wrote:  Hi Joe! Thanks for chiming in and doing some stress tests!
Indeed it fails as it assumes the arrays are same length. I’ll need to think about a fix for those issues.

I hate to make your life more difficult, but another issue arises when the best fraction is one of the inputs. Sometimes your program returns that fraction, but other times not. For example:

Between 1.5 and 1.6, the simplest fraction is 3/2 (same as 1.5), which your program returns.
But between 1.8 and 1.9, the simplest fraction is 9/5, but your program returns 11/6.

In other words, the program should return EITHER the simplest fraction in the open interval (A, B), OR in the closed interval [A, B]. Right now it sometimes does the former and sometimes the latter.

[dlidstrom]

(09-10-2023 12:50 AM)Joe Horn Wrote:  I hate to make your life more difficult, but another issue arises when the best fraction is one of the inputs. Sometimes your program returns that fraction, but other times not. For example:

Between 1.5 and 1.6, the simplest fraction is 3/2 (same as 1.5), which your program returns.
But between 1.8 and 1.9, the simplest fraction is 9/5, but your program returns 11/6.

In other words, the program should return EITHER the simplest fraction in the open interval (A, B), OR in the closed interval [A, B]. Right now it sometimes does the former and sometimes the latter.

This turned out to be much more difficult than I thought. My initial program was small and fast but wrong. Now it has ballooned. But I've tried my best to come up with a more complete program. I say more complete because I'm sure there are issues I haven't thought of or understood completely. Arbitrarily I decided the program can return in the interval (A, B]. That is, it is allowed to return the smallest decimal as the simplest fraction but not the largest decimal.

There are some cases which aren't covered. Such as equal decimals being supplied, and perhaps negative decimals. However, I think I'll leave it here.

Code:


« → ←a ←b
  « ←a →SF →CF DUP SIZE ←b →SF →CF DUP SIZE ROT SWAP DUP2 -
    « → ah at al bh bt bl
      «
        CASE 'ah<0'
          THEN bh "0" →TAG +
          END 'bh<0'
          THEN ah "1" →TAG +
          END '←a<←b AND al≤0.'
          THEN ah "2" →TAG +
          END '←a<←b AND bl≤0.'
          THEN ah "3" →TAG + at +
          END '←b<←a AND al≤0.'
          THEN bh "4" →TAG + bt +
          END '←b<←a AND bl≤0.'
          THEN bh "5" →TAG +
          END 'ah<bh'
          THEN ah 1 + "6" →TAG +
          END 'bh<ah'
          THEN bh 1 + "7" →TAG +
          END ah + at HEAD at TAIL al 1 -
            bt HEAD bt TAIL bl 1 - ←proc EVAL
        END
      »
    » → a b al bl c ←proc
    «
      CASE 'c<0'
        THEN 'a' a -1 c NEG
        END 'b' b -1 c
      END NDUPN →LIST + SWAP STO { } a HEAD a TAIL al 1 -
      b HEAD b TAIL bl 1 - ←proc EVAL DUP CF→ DUP →NUM
    »
  »
»

# A222h
887.5

[Albert Chan]

Perhaps a simple linear search?
Denominator to seek normally small anyway.
Because we never skip, denominator guaranteed smallest!

This has the advantage of only needing fast float. (It is hard to get exact CF coefs. with float)

Code:

function frac_between(a,b)
    if a>b then a,b = b,a end
    local d = 1
    while (a*d) > floor(b*d) do d=d+1 end
    return ceil(a*d), d
end

lua> frac_between(3.1415, 3.1416)
333      106
lua> frac_between(1.8, 1.9)
9      5
lua> frac_between(1.8, 1.8)
9      5
lua> frac_between(1.8, 1.7)
7      4
lua> frac_between(1.8, 1.6)
5      3
lua> frac_between(1.8, 1.5)
3      2
lua> frac_between(-pi, -pi*.99)
-25      8

Update: ceil(a*d) <= floor(b*d) test simplified to (a*d) <= floor(b*d)

[Joe Horn]

(09-13-2023 10:50 PM)Albert Chan Wrote:  Perhaps a simple linear search?
Denominator to seek normally small anyway.
Because we never skip, denominator guaranteed smallest!

Yes, that works ok when the denominator is small, but for a general-purpose math engine it fails terribly when the denominator is large. Dlidstrom's program finds the simplest fraction between 0.999999999977 and 0.999999999978 very fast (namely, 43478260869/43478260870), but linear searching would take way too long... at least on an HP 50g.

[Albert Chan]

(09-14-2023 03:40 AM)Joe Horn Wrote:  Dlidstrom's program finds the simplest fraction between 0.999999999977 and 0.999999999978 very fast (namely, 43478260869/43478260870), but linear searching would take way too long... at least on an HP 50g.

Although math is correct, I doubt user really want fraction with denominator of billions.
We might as well just use 0.999999999977 or 0.999999999978

Anyway, here is the sem-convergent search way.
We need to be less agressive with CF, and use RHS bolded result.

CF(1.8) = [1; 1, 4] = [1; 1, 4, Inf] = [1; 1, 3, 0, 0, 0 ..., 0, 1]

In case both numbers have same CF coef. list length, we keep them sorted.
(and maintain sort order going into recursion)

Code:

from fractions import Fraction

def d_convergent(a, b, d1=1, d2=0):
    ia, ib = a//1, b//1 # assumed a <= b
    if ia == a:  return d2, d1+(ia-1)*d2
    if ia != ib: return d2, d1+ia*d2
    return d_convergent(1/(b-ib), 1/(a-ia), d2, d1+ia*d2)

def frac_between(a,b):
    if a>b: a,b = b,a
    d, d2 = d_convergent(a,b)
    while 1:
        d = d + d2      # semi-convergents
        if (a*d) <= (b*d//1): return Fraction(-(-a*d//1), d)

>>> frac_between(Fraction('0.999999999977'), Fraction('0.999999999978'))
Fraction(43478260869, 43478260870)

CF(10/7) = [1;2,3]
CF(13/9) = [1;2,4]

>>> frac_between(Fraction(10, 7), Fraction(13, 9))
Fraction(10, 7)

Update: ceil(a*d) ≤ floor(b*d) test simplified to (a*d) ≤ floor(b*d)

(a*d) ≤ ceil(a*d) ≤ floor(b*d) ≤ (b*d)

--> a ≤ ceil(a*d)/d ≤ floor(b*d)/d ≤ b

[dlidstrom]

Hi Albert! Are you running this Python program on an HP Prime with beta firmware?