HP Prime program: rounding to a fraction

[Patrice]

Theses 2 programs are rounding a decimal value to the value of a fraction using Farey series. Use 'a b/c' key to display the fraction.

EXPORT FareyMax(Vl, DMax)

// Round a Vl to the best fraction with denominator < DMax

BEGIN

LOCAL VlE, Tmp;

LOCAL DbN,DbD, FnN,FnD, RsN,RsD,RsE;

VlE:= ABS(Vl);

DbN:= INT(VlE); DbD:=1; FnN:=DbN+1; FnD:=1;

RsN:= ROUND(VlE,0); RsD:= 1; RsE:= ABS(VlE-(RsN/RsD));

WHILE DbD+FnD <= DMax DO

  Tmp:= (DbN+FnN)/(DbD+FnD);

  IF RsE > ABS(VlE-Tmp) THEN RsN:= (DbN+FnN); RsD:= (DbD+FnD); RsE:= ABS(VlE-(RsN/RsD)); END;

  IF Tmp < VlE THEN DbN:= (DbN+FnN); DbD:= (DbD+FnD); ELSE FnN:= (DbN+FnN); FnD:= (DbD+FnD); END;

END;

RETURN SIGN(Vl)*RsN/RsD;

END;

EXPORT FareyDelta(Vl, Error)

// round Vl to the smallest fraction with |Vl - fraction| < Error

BEGIN

LOCAL VlE, Tmp;

LOCAL DbN,DbD, FnN,FnD, RsN,RsD,RsE;

VlE:= ABS(Vl);

DbN:= INT(VlE); DbD:=1; FnN:=DbN+1; FnD:=1;

RsN:= ROUND(VlE,0); RsD:= 1; RsE:= ABS(VlE-(RsN/RsD));

WHILE RsE > Error DO

  Tmp:= (DbN+FnN)/(DbD+FnD);

  IF RsE > ABS(VlE-Tmp) THEN RsN:= (DbN+FnN); RsD:= (DbD+FnD); RsE:= ABS(VlE-(RsN/RsD)); END;

  IF Tmp < VlE THEN DbN:= (DbN+FnN); DbD:= (DbD+FnD); ELSE FnN:= (DbN+FnN); FnD:= (DbD+FnD); END;

END;

RETURN SIGN(Vl)*RsN/RsD;

END;

[Joe Horn]

Machine roundoff error causes FareyDelta followed by the [a b/c] button to fail when Error is too small.

Example:
FareyDelta(pi,1E-11) [a b/c] --> 312689/99532

However, the correct answer is 833719/265381, as can be calculated in CAS by setting epsilon:=1E-11 and executing exact(approx(pi)).

This problem can be avoided by having FareyDelta return RsN and RsD explicitly (rather than dividing them, and attempting to reconstruct them with the fraction button), in which case it returns the correct result for the above example.

-Joe-

Edited: 30 Oct 2013, 1:15 a.m.

[Patrice]

Hi Joe,

Do you prefer this way ? no play with 'a b/c'

EXPORT FareyMax(Vl, DMax)

// Round a Vl to the best fraction with denominator < DMax

BEGIN

LOCAL VlE, Tmp;

LOCAL DbN,DbD, FnN,FnD, RsN,RsD,RsE;

VlE:= ABS(Vl);

DbN:= INT(VlE); DbD:=1; FnN:=DbN+1; FnD:=1;

RsN:= ROUND(VlE,0); RsD:= 1; RsE:= ABS(VlE-(RsN/RsD));

WHILE DbD+FnD <= DMax DO

  Tmp:= (DbN+FnN)/(DbD+FnD);

  IF RsE > ABS(VlE-Tmp) THEN RsN:= (DbN+FnN); RsD:= (DbD+FnD); RsE:= ABS(VlE-(RsN/RsD)); END;

  IF Tmp < VlE THEN DbN:= (DbN+FnN); DbD:= (DbD+FnD); ELSE FnN:= (DbN+FnN); FnD:= (DbD+FnD); END;

END;

RETURN EXPR("QUOTE("+SIGN(Vl)*RsN+"/"+RsD+")");

END;

EXPORT FareyDelta(Vl, Error)

// round Vl to the smallest fraction with |Vl - fraction| < Error

BEGIN

LOCAL VlE, Tmp;

LOCAL DbN,DbD, FnN,FnD, RsN,RsD,RsE;

VlE:= ABS(Vl);

DbN:= INT(VlE); DbD:=1; FnN:=DbN+1; FnD:=1;

RsN:= ROUND(VlE,0); RsD:= 1; RsE:= ABS(VlE-(RsN/RsD));

WHILE RsE > Error DO

  Tmp:= (DbN+FnN)/(DbD+FnD);

  IF RsE > ABS(VlE-Tmp) THEN RsN:= (DbN+FnN); RsD:= (DbD+FnD); RsE:= ABS(VlE-(RsN/RsD)); END;

  IF Tmp < VlE THEN DbN:= (DbN+FnN); DbD:= (DbD+FnD); ELSE FnN:= (DbN+FnN); FnD:= (DbD+FnD); END;

END;

RETURN EXPR("QUOTE("+SIGN(Vl)*RsN+"/"+RsD+")");

END;

[Joe Horn]

Quote:

---

RETURN EXPR("QUOTE("+SIGN(Vl)*RsN+"/"+RsD+")");

---

Ah, a thing of beauty! Very nice indeed.