Setting a maximum demoninator for fractions

[GSL007]

Hello everyone,

I'm sorry if I'm asking a dumb question but I'm totally new to this type of calculator. I purchased an HP35S last year, which I love. With this calculator I can set the maximum denominator for a fraction to be what I wish, 8, 12, 16, whatever up to 4095.

I can't figure out how to set the maximum denominator for fractions in the HP Prime G2. The manual didn't help. Can anyone help or do I misunderstand something?

Thanks, Gary L

[Joe Horn]

(10-02-2023 03:53 PM)GSL007 Wrote:  Hello everyone,

I'm sorry if I'm asking a dumb question but I'm totally new to this type of calculator. I purchased an HP35S last year, which I love. With this calculator I can set the maximum denominator for a fraction to be what I wish, 8, 12, 16, whatever up to 4095.

I can't figure out how to set the maximum denominator for fractions in the HP Prime G2. The manual didn't help. Can anyone help or do I misunderstand something?

Thanks, Gary L

The [a b/c] key in Home view mode doesn't set maximum denominators like the 32SII, 33s, & 35s do. Instead, it uses the current display setting (e.g. FIX 2 or FIX 4 or whatever) to control the accuracy of the fraction approximation. Example: Pi is displayed as 22/7 in FIX 2 mode, but as 201/64 in FIX 3 mode. You can specify the accuracy you want by changing the display mode, but you cannot specify the maximum desired denominator.

The exact(x) function in CAS mode works differently. Instead of using the display mode to control its accuracy, it uses the Epsilon setting in CAS Settings (page 2). For example, if Epsilon = 1E-6, then exact(pi+0.) returns 355/113. Again, you can specify the desired accuracy this way, but not the maximum denominator.

[komame]

(10-02-2023 03:53 PM)GSL007 Wrote:  I can't figure out how to set the maximum denominator for fractions in the HP Prime G2.

HP Prime doesn't have built-in functionality for that, but you can achieve a similar effect with the following program:

Code:

// Closest Fraction
EXPORT CF(x, max_d)
BEGIN
  LOCAL curr_num, best_num, best_denom, best_diff, curr_diff, d;
  best_num := ROUND(x,0);
  best_denom := 1;
  best_diff := ABS(x - best_num);

  FOR d FROM 2 TO max_d DO
    curr_num := ROUND(x * d,0); 
    curr_diff := ABS(x - curr_num / d);

    // If the current fraction is closer to x, update the best fraction
    IF curr_diff < best_diff THEN
      best_diff := curr_diff;
      best_num := curr_num;
      best_denom := d;
    END;
  END;
  exact(best_num/best_denom);
END;

Usage:
CF(fraction, max denominator):

Example:
\(CF(PI, 64) => \cfrac{201}{64}\)
\(CF(1.45, 8) => \cfrac{10}{7}\)
\(CF(0.12345,4) => 0\)

This way, you're not limited to a denominator of 4095; you can enter much larger values (but be cautious, as it may lead to lengthy calculations).

Best regards,
Piotr

[parisse]

Ouch, that's a very slow method. A faster way would be to compute the expansion in continued fraction and try the partial quotients obtained by list truncations.
l:=dfc(pi); dfc2f(l[1..2]); dfc2f(l[1..3]); etc.

[komame]

(10-03-2023 05:54 AM)parisse Wrote:  Ouch, that's a very slow method. A faster way would be to compute the expansion in continued fraction and try the partial quotients obtained by list truncations.
l:=dfc(pi); dfc2f(l[1..2]); dfc2f(l[1..3]); etc.

I know this can be heavily optimized. I wrote it quickly as a PoC, but even for a denominator of 4095 (as Gary mentioned as the limit in HP35S), this program only needs 0.3s to produce the result, so it's not that bad

Nonetheless, I also tried to follow the method you indicated, and I have a question about how to calculate this for PI with a maximum denominator of 5 using your method.
My approach shows
\(\cfrac{16}{5}\) for that,
but how do I do it with your method?

Also, I have a strange behavior in HOME and in CAS as well (maybe it's intended, but it seems strange to me).
First, I execute dfc:

Code:

v:=dfc(pi);

and I get vector in 'v':
[3, 7, 15, 1, 292, 1, 1, 1, 2]

Then, I calculate the first iteration:

Code:

f:=dfc2f(SUB(v,1,2))

and I get the result
\(\cfrac{22}{7}\)
but at the same time, 'v' has changed to
[3, 7].
Is it supposed to work this way?

PK

[parisse]

SUB does in place modifications if the first argument is an identifier, that's an heritage of the hp48 series.

[Albert Chan]

(10-03-2023 07:58 AM)komame Wrote:  Nonetheless, I also tried to follow the method you indicated, and I have a question about how to calculate this for PI with a maximum denominator of 5 using your method.
My approach shows
\(\cfrac{16}{5}\) for that,
but how do I do it with your method?

pi convergents = 3/1, 22/7, 333/106, 355/113, ...

From best convergent and best semi-convergent, we pick one with least error.
If it is a tie, pick the one with smaller denominator (i.e. convergent)
see Python Fraction.limit_denominator()

Above example, we have 3/1 vs (22-2*3)/(7-2*1) = 16/5
pi - 3/1 ≈ +0.1416
pi - 16/5 ≈ -0.0584      → CF(pi, 5) = 16/5 (smaller error)

Note: code use only HOME numerical calculations, without CAS features.

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, n2,d2, k, t;
  p  := x; q  := 1; // x = p / q
  n1 := 0; n2 := 1; // n convergent
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [n2,d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q;
    t:=d1; d1:=d2; d2:=t+k*d2;
    t:=n1; n1:=n2; n2:=t+k*n2;
    t:=p; p:=q; q:=t;       // SWAP p,q
  END;
  k:=CEILING((d2-dmax)/d1);
  p:=n2-k*n1; q:=d2-k*d1;   // best semi-convergent
  if abs((n1-d1*x)*q) <= abs((p-q*x)*d1) THEN p:=n1; q:=d1; END;
  RETURN [p,q];
END;

We assumed MOD correctly implemented, without rounding errors.

(09-18-2023 11:28 PM)Albert Chan Wrote:  If n≥0, d>0, (euclidean_mod = floor_mod = trunc_mod), and trunc_mod is *exact*
see Eli Bendersky's blog, Computing remainders by doubling, recursive steps are exact.

HOME> CF(3.1416, 9)      → [22, 7], convergent
HOME> CF(3.1416, 99)      → [311, 99], semi-convergent
HOME> CF(3.1416, 999)      → [2862, 911], semi convergent
HOME> CF(3.1416, 9999)      → [3927, 1250], last convergent, return from inside while loop

[Albert Chan]

(10-03-2023 02:44 PM)Albert Chan Wrote:  if abs((n1-d1*x)*q) <= abs((p-q*x)*d1) THEN p:=n1; q:=d1; END;

This test has trouble if numerator is huge, with possibly catastrophic cancellation issue.
We like numerator small, so my next version remove rounded integer part of x, before the test.

I also squared both sides to remove abs, with equivalent test as product of 2 dot products.
HP Prime HOME dot product work with higher precisions, thus more accurate.

Cas> factor(((p-q*x)*d1)^2 - ((n1-d1*x)*q)^2)      → (q*n1-p*d1)*(2*x*q*d1-q*n1-p*d1)

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, d2, k;
  p  := x; q  := 1; // x = p / q
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [ROUND(x*d2,0),d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q;
    k:=d1+k*d2; d1:=d2; d2:=k;
    k:=p; p:=q; q:=k;       // SWAP p,q
  END;
  k:=ROUND(x,0); x:=x-k;    // reduce numerator
  n1:=ROUND(x*d1,0);        // convergent = [n1,d1]
  q:=d2-CEILING((d2-dmax)/d1)*d1;
  p:=ROUND(x*q,0);          // semi-convergent = [p,q]
  IF DOT([p,q],[d1,-n1])*DOT([p,q,2*x],[d1,n1,-d1*q]) >= 0 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;

[komame]

(10-04-2023 01:03 AM)Albert Chan Wrote:  Cas> factor(((p-q*x)*d1)^2 - ((n1-d1*x)*q)^2)      → (q*n1-p*d1)*(2*x*q*d1-q*n1-p*d1)

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, d2, k;
  p  := x; q  := 1; // x = p / q
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [ROUND(x*d2,0),d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q;
    k:=d1+k*d2; d1:=d2; d2:=k;
    k:=p; p:=q; q:=k;       // SWAP p,q
  END;
  k:=ROUND(x,0); x:=x-k;    // reduce numerator
  n1:=ROUND(x*d1,0);        // convergent = [n1,d1]
  q:=d2-CEILING((d2-dmax)/d1)*d1;
  p:=ROUND(x*q,0);          // semi-convergent = [p,q]
  IF DOT([p,q],[d1,-n1])*DOT([p,q,2*x],[d1,n1,-d1*q]) >= 0 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;

Wow, this is incredibly fast regardless of the denominator value. Even for a case like PI with dmax=30,000,000,000; the result is returned instantly.

PK

[Albert Chan]

(10-04-2023 01:03 AM)Albert Chan Wrote:  HP Prime HOME dot product work with higher precisions, thus more accurate.

Note: HOME DOT function extra precision is not coming from CAS side.

HOME> P:=PI                      → 3.14159265359
HOME> Q:=P*P                   → 9.86960440109
HOME> DOT([P,1],[P,-Q])     → 6.5ᴇ−13

P*P (exactly) = 9.8696044010906577398881
This suggested HOME DOT internally use 15 digits, truncated rounding.

CAS DOT is not the same as HOME DOT, even though they share the same name.

CAS> DOT([P,1],[P,-Q])       → 6.25277607469e−13

CAS float is 48 bits (truncated rounding), it is equivalent to IEEE double, stripped 53-48 = 5 bits.

lua> p, q = 3.14159265359, 9.86960440109
lua> p*p - q
6.590283874174929e-013

lua> hpfloat = fn'x: local hi,lo = bit.split(x); bit.join(hi,lo-lo%32)'
lua> p, q = hpfloat(p), hpfloat(q)
lua> hpfloat(p*p - q)
6.430411758628907e-013

lua> hpfloat(hpfloat(p*p) - q)
6.252776074688882e-013

It appeared CAS DOT doesn't use extended precision internally.

[Albert Chan]

(10-04-2023 09:00 AM)komame Wrote:  Even for a case like PI with dmax=30,000,000,000; the result is returned instantly.

For dmax in the millions, code may not be numerically discern which fraction is better.
By design, code will pick the smaller denominator (i.e. convergent)

For above example, it picked wrong, semi-convergent is better, but it is hardly noticeable.

X = PI-3 = 0.14159265359 (exactly)
10225711/72219220      ≈ 0.14159265359000000277      -- error ≈ 2.77E-18
417041790/2945363191 ≈ 0.14159265358999999807      -- error ≈ -1.93E-18

Code:

IF DOT([p,q],[d1,-n1])*DOT([p,q,2*x],[d1,n1,-d1*q]) >= 0 THEN p:=n1; q:=d1; END;

If first DOT has trouble getting ±1, it will return 0.0, due to cancellation.
If first DOT has trouble, so with 2nd DOT (x is between n1/d1 and p/q)

DOT([p,q], [d1,-n1])                = (d1*q) * (p/q - n1/d1)
DOT([p,q,2*x], [d1,n1,-d1*q]) = (d1*q) * ((p/q-x) + (n1/d1-x))

First DOT, we can get correct result, but I choose not to:

DOT([p,q],[d1,-n1]) = (-1)^(while loops)

Difference of 2 CF convergents: n2/d2 - n1/d1 = ±1 / (d1*d2)      → d1*n2 - d2*n1 = ±1
Same result if we replace n2/d2 by semi-convergent p/q

Cas> simplify( DOT([p,q],[d1,-n1]) (p=n2-k*n1, q=d2-k*d1) )      → d1*n2 - d2*n1 = ±1

[Albert Chan]

(10-04-2023 02:33 PM)Albert Chan Wrote:  DOT([p,q], [d1,-n1])                = (d1*q) * (p/q - n1/d1)
DOT([p,q,2*x], [d1,n1,-d1*q]) = (d1*q) * ((p/q-x) + (n1/d1-x))

Perhaps we could use MOD again, to improve accuracy greatly!
Let s = DOT([p,q], [d1,-n1]) = ±1, trivial to get (see previous post)

DOT([p,q,2*x], [d1,n1,-d1*q]) = p*d1 + q*n1 - 2*x*d1*q = s + 2*q*(n1 - d1*x)

We decompose n1 using MOD: n1 = n2 + d2*x:
n1/d1 is a convergent of x    → sign(x) = sign(n1)    → trunc/floor mod, n2 = (n1 MOD x) is exact

x ≈ n1/d1
d1*x ≈ n1 = n2 + d2*x
n2 ≈ (d1-d2)*x

Since |n2| < |x|, (d1-d2) = 0 or 1, no need to use dot product.
Either case, (n2 - (d1-d2)*x) is exact (Sterbenz lemma)

We pick smaller denominator convergent, if product of 2 DOT's ≥ 0:

s * (s + 2*q*(n1 - d1*x)) ≥ 0
1 + 2*s*q*(n1 - d1*x) ≥ 0
s*q*(n1 - d1*x) = s*q*(n2 - (d1-d2)*x) ≥ -1/2

CF, version 3:

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, n2,d2, k, s:=1;
  p  := x; q  := 1; // x = p / q
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [ROUND(x*d2,0),d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q; s:=-s;
    k:=d1+k*d2; d1:=d2; d2:=k;
    k:=p; p:=q; q:=k;       // SWAP p,q
  END;
  k:=ROUND(x,0); x:=x-k;    // reduce numerator
  n1:=ROUND(x*d1,0);        // convergent = [n1,d1]
  q:=d2-CEILING((d2-dmax)/d1)*d1;
  p:=ROUND(x*q,0);          // semi-convergent = [p,q]
  n2:=(n1 MOD x); d2:=(n1-n2)/x;    // n1 = n2 + d2*x
  IF s*q*(n2-(d1-d2)*x) >= -.5 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;

HOME> CF(PI,3000000000)      → [9253131363,2945363191]

Test is now much more accurate: -0.5890726382 ≥ -0.5 is false --> pick semi-convergent

Confirm in Speedcrunch calculator (v 0.12)

x = 0.14159265359
n1 = 10225711
d1 = 72219220
p = 417041790
q = 2945363191

s = p*d1 - q*n1      → -1
s*q*(n1 - d1*x)      → -0.5890726382

[Albert Chan]

(10-04-2023 04:59 PM)Albert Chan Wrote:  We pick smaller denominator convergent, if product of 2 DOT's ≥ 0:

s * (s + 2*q*(n1 - d1*x)) ≥ 0
1 + 2*s*q*(n1 - d1*x) ≥ 0
s*q*(n1 - d1*x) = s*q*(n2 - (d1-d2)*x) ≥ -1/2

But wait! s = sign(x - n1/d1), we don't need to track it!
LHS always negative, we multiply by -1, for equivalent test.

abs(q*(n2 - (d1-d2)*x)) ≤ 1/2

CF, version 4

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, n2,d2, k;
  p  := x; q  := 1; // x = p / q
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [ROUND(x*d2,0),d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q;
    k:=d1+k*d2; d1:=d2; d2:=k;
    k:=p; p:=q; q:=k;       // SWAP p,q
  END;
  k:=ROUND(x,0); x:=x-k;    // reduce numerator
  n1:=ROUND(x*d1,0);        // convergent = [n1,d1]
  q:=d2-CEILING((d2-dmax)/d1)*d1;
  p:=ROUND(x*q,0);          // semi-convergent = [p,q]
  n2:=(n1 MOD x); d2:=(n1-n2)/x;    // n1 = n2 + x*d2
  IF abs(q*(n2-(d1-d2)*x)) <= .5 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;

[Albert Chan]

There is an easier way to build the test!

gap = |n1/d1 - p/q| = 1 / (q*d1)

If best-convergent to x is half gap or less, pick best-convergent.

|n1/d1 - x| ≤ (1/2) / (q*d1)
|q*(n1 - d1*x)| ≤ 1/2

This half-gap test assumed x between n1/d1 and p/q. It needed confirmation.

We assumed "overshooted" convergent n2/d2 exist, q = (d0+k*d1),   d0 ≤ q < d2
(if n1/d1 is x last convergent, CF code would just return it, never touched the test)

Here, we just proof where (n1*d0 - n0*d1) = 1 (similar proof for -1)

n0/d0 < n2/d2 ≤ x < n1/d1

(n1/d1 - x) ≤ (n1/d1 - n2/d2) = 1/(d1*d2)
(n1/d1 - p/q) = (n1*(d0+k*d1)-(n0+k*n1)*d1)/(d1*q) = 1/(d1*q) > 1/(d1*d2)

For (n1*d0 - n0*d1) = 1, we have p/q < x < n1/d1

[GSL007]

Thank you to everyone who replied helping me understand this behaviour. I am still working through some of your answers. Thanks again, Gary.