\ 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.

\ Classic heapsort
: sift
    ( rra ra ir ii )
    DUP 2 * 1+  ( rra ra ir ii jj )
    BEGIN 
        >R OVER R@ SWAP R> OVER OVER ( rra ra ir ii jj ir jj ir jj )
    > WHILE ( rra ra ir ii jj ir jj )
        1+ > IF ( rra ra ir ii jj )
            >R >R OVER R> SWAP R> SWAP OVER ( rra ra ir ii jj ra jj )
            CELLS + DUP @ SWAP CELL+ @ < IF ( rra ra ir ii jj )
                1 +
            THEN
        THEN ( rra ra ir ii jj )
        >R >R >R OVER OVER R> ROT ROT R> ROT ROT R@
            ROT ROT R> ( rra ra ir ii jj rra ra jj )
        CELLS OVER + @ ROT ( rra ra ir ii jj ra rr_jj rra )
        OVER < IF ( rra ra ir ii jj ra rr_jj )
            >R >R OVER CELLS R> + R> SWAP !
            SWAP DROP DUP DUP + 1+
        ELSE
            DROP DROP DROP ( rra ra ir ii )
            OVER 1 +
        THEN ( rra ra ir ii jj )
    REPEAT
    DROP DROP DROP SWAP DROP ( rra ra ii )
    CELLS + !
;

: heapsort ( ra n )
    0 OVER 1 - 2 / DO
        OVER OVER OVER ( ra n ra n ra )
        I CELLS + @  ROT ROT I ( ra n rra ra n I )
        sift
    -1 +LOOP
    1 - 1 SWAP DO ( ra )
        DUP DUP I CELLS + ( ra ra ra+I )
        DUP @ ( ra ra ra+I rra )
        SWAP ROT ( ra rra ra+I ra )
        @ SWAP ! ( ra rra )
        I 1 = IF
            OVER !
        ELSE
            OVER I 0 ( rra ra I 0 )
            sift
        THEN
    -1 +LOOP
;


\ 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