\ Solver for Ramanujan taxi problem.  Method: first generate sums
\ pairs of numbers, store sum of cubes together with one of nambers
\ in table.  Then sort to find duplicate sum.  Use sum and one
\ of numbers to reconstruct the other.

\ quicksort from http://www.complang.tuwien.ac.at/forth/programs/sort.fs
\ replaced original heapsort code

\ : ]L ] postpone literal ;
\ : cell- [ -1 cells ]L + ;

: exchange ( addr1 addr2 -- )
    over @ over @ >r swap ! r> swap ! ;

: partition ( al ah -- al1 ah1 al2 ah2 )
    \ assert( dup aligned dup = )
    over 2/ over 2/ + -1 cells and @ >r
    2dup cell- begin ( al ah l r )
	2dup u<= while
	    swap begin
		dup @ r@ < while
		    cell+ repeat
	    swap begin
		dup @ r@ > while
		    cell- repeat
	    2dup u<= if
		2dup exchange swap cell+ swap cell- then
    repeat
    r> drop
    cell+ swap rot ;

: sort2 ( al ah -- )
    begin
	over cell+ over u< while
	    partition recurse
    repeat
    2drop ;

: sortm ( addr u -- )
    cells over + sort2 ;

: heapsort ( addr u -- addr )
    \ compatibility word for the rest of the taxi_hs.fs code
    >r dup r> sortm ;

\ end of quicksort code

\ Numbers are limited to 16 bit and sum of qubes is limited to 48 bit.
\ So we can pack one number and sum in 64 bits.
: pack_pair ( n n n -- n )
    DROP SWAP 
    16 LSHIFT OR
;

: 3pow ( n -- n ) DUP DUP * * ;

: pack_taxi ( n n -- n)
    DUP >R 3pow SWAP
    DUP >R 3pow +
    R> R> pack_pair
;

VARIABLE t_cnt
0 t_cnt !
VARIABLE t_a_addr

: store_taxi ( n -- )
    t_a_addr @ t_cnt DUP >R @ DUP >R CELLS + !
    R> 1+ R> !
;
  
: store_taxis ( n -- )
    2 DO
        I 1 DO I J pack_taxi store_taxi LOOP
    LOOP
;

: sqm ( n -- n ) DUP * ;

: sqm1 ( n r -- n r1 ) sqm OVER * ;

: pow10923 ( n -- n )
    65535 AND DUP
    sqm sqm1 sqm sqm1 sqm sqm1 sqm sqm1 sqm sqm1 sqm sqm1 sqm1
    65535 AND SWAP DROP
;

: 3root ( n -- n )
    0
    OVER 16777215 AND 0= IF
        8 +
        SWAP 24 RSHIFT SWAP
    THEN
    OVER 4095 AND 0= IF
        4 +
        SWAP 12 RSHIFT SWAP
    THEN
    OVER 63 AND 0= IF
        2 +
        SWAP 6 RSHIFT SWAP
    THEN
    OVER 7 AND 0= IF
        1+
        SWAP 3 RSHIFT SWAP
    THEN
    SWAP
    pow10923
    SWAP LSHIFT
;

: print_taxi ( n -- )
    DUP 65535 AND >R R@ .
    16 RSHIFT DUP R> 3pow - 3root .
    . CR
;

: print_taxis ( ra n -- )
    0 DO DUP @ print_taxi CELL+ LOOP DROP
;

: get_taxi ( n -- n ) CELLS t_a_addr @ + @ ;

: remove_single ( n -- )
    0 t_cnt !
    0 get_taxi SWAP
    1 DO
        I get_taxi OVER 16 RSHIFT OVER 16 RSHIFT
        = IF
            DUP ROT
            store_taxi store_taxi
        ELSE
            SWAP DROP
        THEN
    LOOP
;

: rama_taxis ( n -- )
    DUP 1- DUP 1- * 2/
    HERE t_a_addr ! CELLS ALLOT 
    store_taxis
    t_a_addr @ t_cnt @ heapsort
    t_a_addr @ t_cnt @ remove_single
    t_a_addr @ t_cnt @ print_taxis
    CR
;

9000 rama_taxis