HP Forums
(50g) Fun with Farey sequences - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (50g) Fun with Farey sequences (/thread-4162.html)

(50g) Fun with Farey sequences - Tugdual - 06-15-2015 08:37 PM

A video to start with.
A wikipedia page for the algorithm.
And some code for the 50g
Code:
%%HP: T(3)A(R)F(,);
\<< 0 1 1 0 1 \-> N A B C D nb
  \<< N 'D' STO A B /
    WHILE C N <
    REPEAT 'nb' INCR DROP N B + D / IP R\->I DUP C * A - SWAP D * B - C 'A' STO D 'B' STO 'D' STO 'C' STO A B /
    END nb \->LIST
  \>>
\>>
Usage: enter order N and call the function; returns list of items for the Farey sequence order N.

RE: (50g) Fun with Farey sequences - Thomas Klemm - 06-15-2015 08:53 PM

More fun …

Cheers
Thomas

RE: (50g) Fun with Farey sequences - Tugdual - 06-15-2015 10:34 PM

(06-15-2015 08:53 PM)Thomas Klemm Wrote:  More fun …

Cheers
Thomas
Nice!
Note: Had no idea one could sum List with objects to push back that object in the list. The 48 code is very slim.

RE: (50g) Fun with Farey sequences - Didier Lachieze - 06-16-2015 07:09 PM

Thanks for the links. Great video !

I've adapted the Python program for the Prime:

Code:
#pragma mode( separator(.,;) integer(h32) )
EXPORT farey(n)
BEGIN
     LOCAL a:=0,b:=1,c:=1,d:=n,k,l;
     l:={ a/b };
     WHILE c < n DO
         k:=IP((n+b)/d);
         c:=−a+k*(a:=c); d:=-b+k*(b:=d);
         l:=CONCAT(l,a/b);
     END;  
     exact(l);
END;

The result is a list of reals, you need to press the "a b/c" key to display it as fractions. I don't know how to do it programmatically.
EDIT: edited to add exact() to the result as suggested by Joe below.

RE: (50g) Fun with Farey sequences - Joe Horn - 06-16-2015 10:50 PM

(06-16-2015 07:09 PM)Didier Lachieze Wrote:  The result is a list of reals, you need to press the "a b/c" key to display it as fractions. I don't know how to do it programmatically.

The "exact" function will serve your purposes. exact({1.2, 3.4}) returns {6/5, 17/5}, even in a Home program. Doesn't work on arrays in Home, just lists. Works on arrays and lists in CAS programs, though.

RE: (50g) Fun with Farey sequences - Didier Lachieze - 06-17-2015 05:20 AM

Thanks !

RE: (50g) Fun with Farey sequences - Didier Lachieze - 06-17-2015 08:04 AM

(06-15-2015 10:34 PM)Tugdual Wrote:  The 48 code is very slim.
yes and you can have the UserRPL program returning the same result as the SysRPL version (a list of algebraic fractions) by replacing: a b →V2 + by: a b / →Q +

RE: (50g) Fun with Farey sequences - Thomas Klemm - 06-18-2015 09:33 PM

(06-15-2015 10:34 PM)Tugdual Wrote:  The 48 code is very slim.

This program avoids the STO command:
Code:
\<< \-> n
    \<< { [ 0 1 ] } 0 1 DUP n
        WHILE OVER n <
        REPEAT \-> a b c d
            \<< c d \->V2 +
                c d
                n b + d / IP
                DUP c * a -
                SWAP d * b -
            \>>
        END
        DROP2 DROP2
    \>>
\>>

Cheers
Thomas

RE: (50g) Fun with Farey sequences - Gerald H - 06-23-2015 09:54 AM

I have this for the 49G.

Code:
::
  CK1&Dispatch
  BINT1
  ::
    COERCE
    BINT1
    FPTR2 ^2LAMBIND
    Z0_
    BINT0
    ONEONE
    2GETLAM
    BEGIN
    1GETLAM
    #1+
    1PUTLAM
    OVER
    FPTR2 ^#>Z
    OVER
    FPTR2 ^#>Z
    FPTR2 ^NDXFext
    5UNROLL
    4ROLL
    4ROLL
    2GETLAM
    OVER#+
    4PICK
    #/
    SWAPDROP
    5PICK
    OVER
    #*
    4ROLL
    #-
    3UNROLL
    4PICK
    #*
    SWAP#-
    2DUP#>
    UNTIL
    BINT4
    NDROP
    1GETABND
    {}N
  ;
;


RE: (50g) Fun with Farey sequences - Thomas Klemm - 06-24-2015 10:23 PM

(06-23-2015 09:54 AM)Gerald H Wrote:  
Code:
    5UNROLL         ( … a b c d )
    4ROLL           ( … b c d a )
    4ROLL           ( … c d a b )
    2GETLAM         ( … c d a b n )
    OVER#+          ( … c d a b n+b )
    4PICK           ( … c d a b n+b d )
    #/              ( … c d a b k )
    SWAPDROP        ( … c d a k )

Couldn't that be replaced by:
Code:
    5UNROLL         ( … a b c d )
    4ROLL           ( … b c d a )
    4ROLL           ( … c d a b )
    2GETLAM         ( … c d a b n )
    #+              ( … c d a n+b )
    3PICK           ( … c d a n+b d )
    #/              ( … c d a k )

Unfortunately I'm not familiar with the HP-49/50. In addition to that I wasn't able to compile your program successfully with Debug4x. It appears to be a problem with flash-pointers:

Entry FPTR2 does not exist.

Therefore I can't run the program and must admit that I don't fully understand what's going on in this section:
Code:
    OVER
    FPTR2 ^#>Z
    OVER
    FPTR2 ^#>Z
    FPTR2 ^NDXFext

Would you mind explaining it? Maybe even provide stack diagrams as comment?
Without these I would have had a hard time to understand how the HP-48 is calculating complex arccos and arcsin functions.

Kind regards
Thomas

RE: (50g) Fun with Farey sequences - Gerald H - 06-25-2015 05:46 AM

Dear Thomas,

The code replacement you suggest is not good.

#/ is bint division, stack diagram

Code:

bint a    bint remainder(a/b)
bint b -> bint quotient(a/b)

FPTR2 ^NDXFext does exist, stack diagram

Code:
zint a 
zint b -> 'zint a/zint b'

You may have a problem with Z0_ which is shorthand for ZINT 0

The section you don't fully understand does this

Code:

          bint a
bint a    bint b
bint b -> 'zint a/zint b'

I suggest you change Z0_ to ZINT 0 & try to compile again.

Thanks for your interest.

RE: (50g) Fun with Farey sequences - Thomas Klemm - 06-26-2015 11:19 PM

(06-25-2015 05:46 AM)Gerald H Wrote:  #/ is bint division, stack diagram

Code:

bint a    bint remainder(a/b)
bint b -> bint quotient(a/b)

From "Programming in System RPL" by Eduardo M Kalinowski and Carsten Dominik:
Code:
Addr.  Name   Description
303E9  %/     ( % %' → %/%' )
303D3  %%/    ( %% %%' → %%/%%' )
03EF7  #/     ( # #' → #r #q )

Thus #/ is what is called /MOD in Forth:
Code:
/MOD     ( n1 n2 -- rem quot )     Divides. Returns the remainder and quotient.

Who would think that this was a good idea to use #/ instead?

Thanks for the explanation of ^NDXFex. I assumed something like this but didn't know where to look it up.

Quote:Thanks for your interest.

I try my best to keep up with your production of programs. But I must admit that I'm not fluent in System RPL. Thus I'm a slow reader.

Cheers
Thomas