billw@hpcvra.cv.hp.com. (William C Wickes) (04/03/91)
From the HP 48 Design Team:
->Q OR NOT ->Q, THAT IS THE QUESTION
Thanks to Joseph Horn, the HP28/48 community now has additional fast
rationalization functionality for their machines (see "HP 48 Improved
->Q".) We congratulate him for this elegant solution! There are a
number of interesting points regarding DEC2FRACs relation to ->Q
which deserve some mention. Perhaps these will stimulate discussion
and further contributions.
DEC2FRAC is an addition and not a replacement for ->Q because its aim
is rather different. DEC2FRAC's aim is to produce the fraction with
the _smallest error_ and a given denominator size. ->Q's aim is to
produce the fraction with the _smallest denominator_ and a given error
size. This is best illustrated with the golden mean, '(1+\sqr(5))/2'.
Given this and 1E11, DEC2FRAC returns '139583862445/86267571272' while
->Q returns (in STD mode) '514229/317811'. When evaluated, these
return _exactly_ the same floating point number. The idea of ->Q is
to "round to simpler" rather than "round to closer." This is useful
where you have what you believe to be a "simple" fraction contaminated
with roundoff or other noise.
The HP48 manuals did little to clarify this issue, especially since we
wanted to make no specific claim with only a conjecture of correctness. In
this regard, posting a proof that the algorithm fills in all possible
fractions would be a great boon to the community.
While for many "interesting" numbers the continued fraction sequence as
generated by DUP IP SWAP FP INV ... is exact, even though the floating
point representation is not, in other cases, you get "noise" unexpectedly
early. Does this have any significance regarding the "safety" of generating
the numerator of the fraction from the denominator by division?
While ->Q keeps around the entire continued fraction sequence entailing
some speed penalty, there may be other statistics about the sequence
(such as periodicity) which can be of interest in "deciphering" a
floating point number. Suggestions along these lines are most welcome.
Thanks again to J. Horn for a stimulating and useful post.
Included below is a program, NEW2Q which follows Horn's algorithm, but
uses exit conditions like those of ->Q.
@ NEW2Q: Version of ->Q based on J.K.Horn's Algorithm
@
@ Input:
@
@ 2: Decimal Number to be converted to a fraction
@ 1: Number of decimal places of zeros required in the error.
@
@ Output:
@
@ 1: 'Numerator/Denominator'
@
@ Example:
@
@ What's the simplest fraction which agrees with sqrt(3) to 4 decimal places?
@ 3 \sqr 4 NEQ2Q returns '97/56'
@ '97/56-\sqr(3)' EVAL returns .00009294957
@ ^^^^ note 4 zeros.
@
@ Checksum and bytes:
@ #3992d
@ 620.5
%%HP: T(3)A(R)F(.);
\<< DUP2
IF 1 > SWAP FP AND
THEN OVER XPON 1 - @ calculate the input exponent.
\<< \-> X 'IFTE(X==0,-500,XPON(X))' \>> @ define a 0-tolerant XPON.
\-> f c x expo
\<< 0 1 1 f DUP IP SWAP FP @ set recursion initial cond.s.
WHILE
OVER 5 PICK / f - ABS expo EVAL @ calculate expon. of error
x SWAP - c < @ and compare with input.
OVER AND @ if not close enough and
@ the remainder's non-zero
REPEAT
INV DUP FP SWAP IP @ then calculate next iterate.
\-> B0 B1 A0 A1 R B
\<< B1 'B*B1+B0' EVAL
A1 'B*A1+A0' EVAL
R
\>>
END
DROP SWAP DROP SWAP
DUP 4 ROLL - DUP f * 0 RND @ calc. "missing" den. and num.
\-> N D D0 N0
\<<
IF 'x-expo(ABS(f-N0/D0))<c' @ if "missing" frac. is not
THEN N D @ good enough, use original.
ELSE N0 D0
IF 'N0\=/N' @ If it is really new,
THEN 200 .2 BEEP @ then beep.
END
END
\>>
\>> @ We're done, now clean up.
IF DUP ABS 1 >
THEN # 352318d SYSEVAL
ELSE DROP
END
ELSE DROP
END
\>>cloos@acsu.buffalo.edu (James H. Cloos) (04/11/91)
In article <25590130@hpcvra.cv.hp.com.> billw@hpcvra.cv.hp.com. (William C Wickes) writes: >From the HP 48 Design Team: > [ETC] > >Included below is a program, NEW2Q which follows Horn's algorithm, but >uses exit conditions like those of ->Q. > >@ NEW2Q: Version of ->Q based on J.K.Horn's Algorithm >@ >@ Input: >@ >@ 2: Decimal Number to be converted to a fraction >@ 1: Number of decimal places of zeros required in the error. >@ >@ Output: >@ >@ 1: 'Numerator/Denominator' >@ >@ Example: >@ >@ What's the simplest fraction which agrees with sqrt(3) to 4 decimal places? >@ 3 \sqr 4 NEQ2Q returns '97/56' >@ '97/56-\sqr(3)' EVAL returns .00009294957 >@ ^^^^ note 4 zeros. >@ >@ Checksum and bytes: >@ #3992d >@ 620.5 > >%%HP: T(3)A(R)F(.); >\<< DUP2 > IF 1 > SWAP FP AND > THEN OVER XPON 1 - @ calculate the input exponent. > \<< \-> X 'IFTE(X==0,-500,XPON(X))' \>> @ define a 0-tolerant XPON. > \-> f c x expo > \<< 0 1 1 f DUP IP SWAP FP @ set recursion initial cond.s. > WHILE > OVER 5 PICK / f - ABS expo EVAL @ calculate expon. of error > x SWAP - c < @ and compare with input. > OVER AND @ if not close enough and > @ the remainder's non-zero > REPEAT > INV DUP FP SWAP IP @ then calculate next iterate. > \-> B0 B1 A0 A1 R B > \<< B1 'B*B1+B0' EVAL > A1 'B*A1+A0' EVAL > R > \>> > END > DROP SWAP DROP SWAP > DUP 4 ROLL - DUP f * 0 RND @ calc. "missing" den. and num. > \-> N D D0 N0 > \<< > IF 'x-expo(ABS(f-N0/D0))<c' @ if "missing" frac. is not > THEN N D @ good enough, use original. > ELSE N0 D0 > IF 'N0\=/N' @ If it is really new, > THEN 200 .2 BEEP @ then beep. > END > END > \>> > \>> @ We're done, now clean up. > IF DUP ABS 1 > > THEN # 352318d SYSEVAL > ELSE DROP > END > ELSE DROP > END >\>> When I saw this article, I immediately typed the program in (the set-up here makes that easier for ascii programs than using Kermit ;^( ), assigned it to key 33.2 (where \->Q normally is), and have been happily using it ever since. I did run into a couple of problems today, though. They are listed below with suggested work arounds. After playing around with the recently posted LSQ, I had left the calculator in Numeric Evaluation mode (ie, flag -3 was set). I discovered that this makes for a bit of a problem in using NEW2Q, a problem that does not show up in \->Q. For the curious, the routine at 5603Eh (351218d) is just the symbolic divide function; the same one that is called by the user language / function when it sees that its arguments are both either Global Names or Algebraics. If numeric evaluation is in effect, you might as well call \->NUM right after the SYSEVAL above; the effect is the same. The end result is that you get another real to the stack rather than the symbolic you were expecting. If the rational is exact, it'll look just like the level 1 argument is simply DROPed. The solution is to save the flags in a local, clear -3, and restore the flags at the end. The other problem is more of a bug. If the value you are trying to rationalize, to the number of decimal places specified in level 1, comes out as an integer, then you will get a 0-divided-by-0 error. Specifically, in this case the locals 'N0' and 'D0' are both 0. To solve this, I took the following section: \-> N D N0 D0 \<< IF 'x-expo(ABS(f-N0/D0))<c' THEN N D ELSE N0 D0 IF 'N0\=/N' THEN 200 .2 BEEP END END \>> and changed it to: \-> N D N0 D0 \<< IF D0 THEN IF 'x-expo(ABS(f-N0/D0))<c' THEN N D ELSE N0 D0 IF 'N0\=/N' THEN 200 .2 BEEP END END ELSE N D END \>> While looking at the code in my 48 just now to make sure I got the above changes correct, I discovered another small change I made. As written, if the real is negative, you can get results like '5/-2', which I don't exacly like. The change I made to ensure that the result would look like '-(5/2)' instead, is from: IF DUP ABS 1 > THEN # 352318d SYSEVAL ELSE DROP END to: IF DUP ABS 1 . THEN DUP2 * SIGN ROT ABS ROT ABS # 5603Eh SYSEVAL * ELSE DROP END I used * after the DUP2 rather than / because multiplication is usually faster than division, but as the speed difference in this case is so minimal, change it to / if you think it looks more understandable. On a slightly different topic, if you would like a \->ZQ function using this algorithm (ZQ meaning integer+rational), change the section immediately above to look like: IF DUP ABS 1 > THEN DUP2 / DUP SIGN SWAP ABS IP 4 ROLL ABS 4 PICK ABS MOD 4 ROLL ABS # 352318d SYSEVAL + * ELSE DROP END I'm presently calling this version of the function NW2ZQ, as it fits better in the menu. (As a side note, I would tend to expect that if these were written in the internal -JimC -- James H. Cloos, Jr. Phone: +1 716 673-1250 cloos@ACSU.Buffalo.EDU Snail: PersonalZipCode: 14048-0772, USA cloos@ub.UUCP Quote: <>