(50G) PDQ Algorithm in SRPL and URPL

[Joe Horn]

PDQ Algorithm in HP 50g System RPL and User RPL, by Joe Horn

Here are two HP 50g versions of the "PDQ Algorithm", fully described HERE. The first version is in System RPL (more or less), and the second is 100% User RPL. They are identical in operation except for speed: the System RPL version is roughly twice as fast as the User RPL version. Example:

5. √ 2 PDQ --> 29/13
SRPL version: 0.37 seconds
URPL version: 0.78 seconds

Differences from the Prime version linked above:

(1) Inputs:
Stack level 2: Number to be approximated (as a real decimal, integer ratio, or string).
Stack level 1: Tolerance (as a real decimal, or integer ratio).
Allowing the inputs to be in the form of integer ratios or strings is what lets the user tap the "infinite precision" of the algorithm (see examples below).

(2) Outputs:
Stack level 2: Exact best fraction
Stack level 1: Exact error (tagged "N" or "X")
"N" means "Normal" (principal) convergent
"X" means "eXtra" (intermediate) convergent

Examples:

3.14159265359 7 PDQ
--> 75948/24175, X:-9602153/96700000000000
What this means: The input means, "What is the simplest fraction within ±1/10^7 of 3.14159265359?" The answer is 75948/24175, which is an intermediate convergent of the input (that's what the "X" indicates). The answer differs from the input by exactly -9602153/96700000000000.

'31415926535897932384626/10000000000000000000000' '1/800' PDQ
--> 179/57, X:exact error (ratio of two large integers)

"3.1415926535897932384626" .00125 PDQ
--> same as previous example

(3) 'ic' and 'err' global variables are not created.

(4) Since the 50g uses 12-digit BCD, and Prime's CAS uses 48-bit binary floating point, real inputs on both machines will usually be slightly different, which can cause PDQ's results to differ on both machines. When their inputs are actually identical (e.g. when they are exact ratios of integers), PDQ's results will also be identical on both machines.

PDQ in System RPL
BYTES (when assembled): 852. #A4ABh

Code:

%%HP: T(3)A(R)F(.);
"::
  CK2NOLASTWD
  STRIPTAGS
  STRIPTAGSl2
  '
  ::
    DUP
    xTYPE
    %9
    %-
    xTHEN
    %0
    xELSE
    ::
      DUP
      x\->LST
      DUP
      %1
      '
      xTYPE
      xDOSUBS
      xR>I
      {
        ZINT 28
        ZINT 28
        ZINT 18
      }
      xSAME
      SWAPDUP
      xSIZE
      xGET
      '
      x/
      xSAME
      xAND
    ;
    xTHEN
    xFXND
    xELSE
    ::
      DUP
      xFXND
      xMOD
      xTHEN
      ::
        % 11.
        OVER
        %EXPONENT
        %-
        %ALOG
        SWAPOVER
        %*
        xR>I
        SWAP
        xR>I
        xSIMP2
      ;
      xELSE
      ZINT 1
    ;
  ;
  xRPN->
  LAM tof
  x<<
  OVER
  xTYPE
  %2
  xSAME
  xTHEN
  ::
    SWAPDUP
    \".\"
    xPOS
    2DUP
    %1
    %-
    %1
    SWAP
    xSUB
    3PICK
    3PICK
    ZINT 1
    x+
    OVER
    xSIZE
    xSUB
    x+
    xSTR>
    ROT
    xSIZE
    ROT
    x-
    xR>I
    xALOG
    x/
    SWAP
  ;
  xIF
  OVER
  xFXND
  xMOD
  ZINT 0
  x#?
  xTHEN
  ::
    DUP
    ZINT 1
    x>=?
    OVER
    xTYPE
    % 28.
    xSAME
    xAND
    xTHEN
    ::
      xNEG
      xALOG
    ;
    xIFEND
    LAM tof
    EVAL
    UNROT
    xABS
    UNROT
    LAM tof
    EVAL
    2DUP
    ZINT 0
    ZINT 1
    xRPN->
    LAM a
    LAM b
    LAM n0
    LAM d0
    LAM n
    LAM d
    LAM cd
    LAM pd
    x<<
    ZINT 1
    LAM d0
    LAM a
    FPTR2 ^QMul
    BEGIN
    LAM pd
    LAM cd
    DUP
    '
    LAM pd
    STOLAM
    LAM n
    LAM d
    DUP
    '
    LAM n
    STOLAM
    FPTR2 ^IDIV2
    '
    LAM d
    STOLAM
    FPTR2 ^QMul
    FPTR2 ^QAdd
    '
    LAM cd
    STOLAM
    SWAP
    FPTR2 ^QNeg
    SWAP
    LAM b
    LAM d
    FPTR2 ^QMul
    OVER
    LAM cd
    FPTR2 ^QMul
    Z<=
    UNTIL
    DUP
    LAM cd
    FPTR2 ^QMul
    LAM b
    LAM d
    FPTR2 ^QMul
    FPTR2 ^QSub
    SWAP
    LAM pd
    FPTR2 ^QMul
    LAM b
    LAM n
    FPTR2 ^QMul
    FPTR2 ^QAdd
    FPTR2 ^IDIV2
    DROP
    LAM n
    OVER
    FPTR2 ^QMul
    LAM d
    FPTR2 ^QAdd
    ROT
    FPTR2 ^QMul
    OVER
    FPTR2 ^QNeg
    LAM pd
    FPTR2 ^QMul
    LAM cd
    FPTR2 ^QAdd
    2DUP
    LAM n0
    FPTR2 ^QMul
    FPTR2 ^QAdd
    LAM d0
    ROT
    FPTR2 ^QMul
    SWAPOVER
    FPTR2 ^QDiv
    UNROT
    FPTR2 ^QDiv
    ROT
    ZINT 0
    Z<>
    ITE
    \"X\"
    \"N\"
    >TAG
    ABND
  ;
  xELSE
  ::
    DROP
    TAG N
    ZINT 0
  ;
  ABND
;
@"

PDQ in User RPL
BYTES: 818.5 #7596h

Code:

%%HP: T(3)A(R)F(.);
\<< DTAG SWAP DTAG SWAP
  \<<
    IF
      IF DUP TYPE 9. -
      THEN 0.
      ELSE DUP \->LST DUP 1.
        \<< TYPE
        \>> DOSUBS R\->I { 28 28 18 } SAME SWAP DUP SIZE GET { / } 1. GET SAME AND
      END
    THEN FXND
    ELSE
      IF DUP FXND MOD
      THEN 11. OVER XPON - ALOG SWAP OVER * R\->I SWAP R\->I SIMP2
      ELSE 1
      END
    END
  \>> \-> tof
  \<<
    IF OVER TYPE 2. SAME
    THEN SWAP DUP "." POS DUP2 1. - 1. SWAP SUB PICK3 PICK3 1 + OVER SIZE SUB + STR\-> ROT SIZE ROT - R\->I ALOG / SWAP
    END
    IF OVER FXND MOD 0 \=/
    THEN
      IF DUP 1 \>= OVER TYPE 28. SAME AND
      THEN NEG ALOG
      END tof EVAL ROT tof EVAL DUP2 0 1 \-> a b n0 d0 n d cd pd
      \<< 1 d0 a *
        DO pd cd DUP 'pd' STO n d DUP 'n' STO IDIV2 'd' STO * + 'cd' STO SWAP NEG SWAP
        UNTIL b d * OVER cd * \<=
        END DUP cd * b d * - SWAP pd * b n * + IQUOT n OVER * d + ROT * OVER NEG pd * cd + DUP2 n0 * + d0 ROT * SWAP OVER / UNROT /
        IF ROT 0 \=/
        THEN "X"
        ELSE "N"
        END \->TAG
      \>>
    ELSE DROP :N: 0
    END
  \>>
\>>

[brickviking]

(05-04-2015 04:48 PM)Joe Horn Wrote:  PDQ Algorithm in HP 50g System RPL and User RPL, by Joe Horn

Here are two HP 50g versions of the "PDQ Algorithm", fully described HERE. The first version is in System RPL (more or less), and the second is 100% User RPL.
...

I'm completely new, and I've been trying to get this to work, so far I'm not having any luck.
First, I tried cut-and-pasting this into an editor (vim), set the line-endings to unix, then saved that as a file on a SD card. I stuck the card back into the HP, in the filer, I navigated to the right place, put my cursor on

Code:

"..." PDQ  STRNG  2414

and chose COPY, to put it into a directory in main memory. Then, I went to that location, and hit EVAL, which put it into my stack as 1:
I then typed out ASM, and got back the message asm Error: Invalid File.
So, what did I actually need to do instead?

(Post 4)

Regards, BrickViking

[rprosperi]

(08-26-2015 09:19 PM)brickviking Wrote:  So, what did I actually need to do instead?

The way the listings are posted, with the "%%HP: T(3)A(R)F(.);" at the top, they are intended to be transferred to the 50g using the HP ConnKit. The %%HP.... directives are instructions to tell the ConnKit how to translate the file when downloading.

If you copy as a string via the SD card, remove this top line as it has no meaning to the SysRPL compiler.

[brickviking]

Okay. What step do I do after I've transferred it to the calculator? It's in my Home directory at the moment, as a STRNG object, 2382 bytes long. Do I put the cursor on it in the Filer, and hit Eval?

This would put the object into my stack, after which, I hit ON and then I should type ASM ? Doing this shows me (after MASD finishes) with the following:

Code:


1: External External
    External  * 1d
    External * 10868d
...... continues on with more code

* means some funny circle with spikes coming off it, I can't type it here.

Now, do I store that into a varname, like .. ooo, I dunno. ->PDQ ?

Ah, found it in UsrRPLTut. PDQ Enter, STO>

Cheers for your patience.

(Post 6)

Regards, BrickViking

[Gerald H]

What is the command "x\->LST"?

[Joe Horn]

(08-27-2015 05:19 AM)Gerald H Wrote:  What is the command "x\->LST"?

Also known as "ROMPTR 100 A", it's one of the built-in Development Library commands. It's similar to →LIST, but it also converts programs and algebraic objects to lists. If you type 256.02 MENU, you'll see it on the F5 key. If you keep library 256 attached (a good idea), you'll also be able to see it by pressing [APPS] [up-arrow] [ENTER] [NXT]. If you set flag -86, your 50g will automatically attach library 256 at every warmstart.

In the AUR, →LST is described on page 6-6 (which is page 504 in the PDF version).

[Gerald H]

Thanks, Joe, don't know how I could have missed the command previously, seems extremely useful.

pdq programme works nicely on 50G, nicely optimised for User RPL.

[Han]

Hi Joe,

Thank you for contributing your PDQ code for both the HP Prime and 50G. Just as a matter of curiosity, though, I was wondering why your SysRPL code was not "purge" SysRPL and instead is a combination of of SysRPL and regular UserRPL.

Han

[Joe Horn]

(09-11-2015 06:08 PM)Han Wrote:  Thank you for contributing your PDQ code for both the HP Prime and 50G. Just as a matter of curiosity, though, I was wondering why your SysRPL code was not "purge" SysRPL and instead is a combination of of SysRPL and regular UserRPL.

Assuming that you meant "pure" SysRPL, the answer is that my only goal was to improve the speed. The inner loop, where most execution time is spent, is pure SysRPL, so my goal was achieved. Converting the rest of the program into pure SysRPL is left as an exercise for people with more free time and patience than I have.

[Gerald H]

(09-12-2015 03:58 AM)Joe Horn Wrote:  

(09-11-2015 06:08 PM)Han Wrote:  Thank you for contributing your PDQ code for both the HP Prime and 50G. Just as a matter of curiosity, though, I was wondering why your SysRPL code was not "purge" SysRPL and instead is a combination of of SysRPL and regular UserRPL.

Assuming that you meant "pure" SysRPL, the answer is that my only goal was to improve the speed. The inner loop, where most execution time is spent, is pure SysRPL, so my goal was achieved. Converting the rest of the program into pure SysRPL is left as an exercise for people with more free time and patience than I have.

While I doubt that I have more patience than Joe I have increased the extent of sys commands in the programme, however some of the user commands are flexible in just the way needed in the programme & concise.

I will publish in due course & hope someone else beats me to it.

[cdmackay]

(05-04-2015 04:48 PM)Joe Horn Wrote:  PDQ Algorithm in HP 50g System RPL and User RPL, by Joe Horn

belated thank you Joe.

Is there a SysRPL version for HP48 anywhere? I can't seem to find one…

or should I just use the UserRPL version above?

thanks very much.

[Joe Horn]

(02-23-2019 09:09 PM)cdmackay Wrote:  Is there a SysRPL version for HP48 anywhere? I can't seem to find one…

or should I just use the UserRPL version above?

thanks very much.

Unfortunately, the User RPL program above cannot be used on the HP 48, because it uses long integers, which is an object type which was introduced into RPL with the HP 49G (after the HP 48). However, you might enjoy a less powerful version (limited to the HP 48's range of 12-digit mantissas) which was posted by Bill Wickes to the comp.sys.handhelds newsgroup on 2 April 1991, and included on Goodies Disk #3. Its commented source code follows. Don't worry about the SYSEVAL; it should work in any HP 48.

Code:

%%HP:T(3)F(.);
@ by the HP 48 Design Team
@
@ NEW2Q: Version of ->Q based on J.K.Horn's Algorithm,
@ but uses exit conditions like those of ->Q.
@
@ 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 û 4 NEQ2Q returns '97/56'
@   '97/56-û3' EVAL returns .00009294957
@                                 ^^^^  note 4 zeros.
@
@

\<< -3 CF 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
\>>


[cdmackay]

(02-24-2019 07:54 AM)Joe Horn Wrote:  Unfortunately, the User RPL program above cannot be used on the HP 48, because it uses long integers,

thanks very much indeed, Joe.

Indeed, I'd read about the long integer issue, and then forgotten, sorry.

I didn't find those posts, I think because I was searching for "PDQ".

I'll add both HP's NEW2Q and your SysRPL D2F to my 48GX, thanks!