( FACTORF - prime factors with HP71 FORTH)
( Tapani Tarvainen 85/12/08)
( ?FORGET TIME .TIMED  from FLIB)

?FORGET (UD.)  HERE  DECIMAL
: (UD.) ( ud)   <# #S #> TYPE ;
: (U.) ( u)   0 (UD.) ;
: UD. ( ud)   (UD.) SPACE ;
: D0= ( d -- ?)   OR 0= ;
: D= ( d1 d2 -- ?)   D- D0= ;
: D<> ( d1 d2 -- ?)   D= NOT ;
: INPUTD ( -- d)   QUERY BL WORD NUMBER ; ( Beware floating-point)
: UD->U ( ud -- u)   IF DROP -1 THEN ; ( overflow to 2^20-1)
: UM/ ( ud1 u -- ud2)   M/MOD ROT DROP ;
: UM2/ ( ud -- u)   2 UM/ UD->U ;
: UMIN ( u1 u2 -- min{u1,u2})   2DUP SWAP U<  IF SWAP THEN  DROP ;
: SQITERATE ( ud u1 -- u2)   DUP >R UM/  R@ 0 D+ UM2/  R> UMIN ;
: SQRTLOOP ( ud guess -- u=sqrt{ud})
   >R BEGIN 2DUP R@ SQITERATE  R> SWAP  DUP >R = UNTIL 2DROP R> ;
: UDSQRT ( ud -- u=sqrt{ud})   2DUP UD->U  SQRTLOOP ;

: FOUND ( n -- u ;  return n+I as result of SEARCH')
   R> DROP R> +  R> DROP ;
: SEARCH' ( ud u1 u2 -- ud u3)
   SWAP 1+ SWAP DO
   2DUP I       UM/MOD DROP 0= IF  0 FOUND THEN 
   2DUP I 2+ 2+ UM/MOD DROP 0= IF  4 FOUND THEN 
   2DUP I 5+ 1+ UM/MOD DROP 0= IF  6 FOUND THEN 
   2DUP I 5+ 5+ UM/MOD DROP 0= IF 10 FOUND THEN 
   2DUP I 12 +  UM/MOD DROP 0= IF 12 FOUND THEN 
   2DUP I 16 +  UM/MOD DROP 0= IF 16 FOUND THEN 
   2DUP I 22 +  UM/MOD DROP 0= IF 22 FOUND THEN 
   2DUP I 24 +  UM/MOD DROP 0= IF 24 FOUND THEN 
   30 +LOOP 1 ;
: SEARCH ( ud limit start -- ud u ;  u=factor or 1)
   OVER 0< IF   OVER >R  SEARCH'  DUP 1 <>  IF R> DROP EXIT THEN
   DROP R> 524317 THEN   SEARCH' ;
: NEXTSMALL ( u1 -- u2)   DUP 2 U< SWAP 3 U< 2* + 5+ ;
5 CONSTANT LASTSMALL
: INIT-SEARCH ( u1 -- u2)    DUP  0  30 UM/MOD DROP  -  7 + ;
: FINDBIG ( ud u1 -- ud u2)   >R  2DUP UDSQRT  R>
   INIT-SEARCH  2DUP U<  IF 2DROP 1  ELSE SEARCH THEN ;
: TRY ( ud u -- ud u ;  returns past one level if u is factor)
   DUP 2OVER ROT  UM/MOD DROP  0= IF  R> DROP THEN ;
: NEXTFACTOR ( ud u1 -- ud u2 ;  u1=previous factor, u2=next)
   BEGIN  DUP LASTSMALL U< WHILE  NEXTSMALL TRY  REPEAT  FINDBIG ;

: EXPONENT? ( ud1 u -- ud2 u n ;  ud2=ud1/u^n)
   0 >R BEGIN   >R  2DUP R@ M/MOD ROT
   0= WHILE  2SWAP 2DROP R>  R> 1+ >R  REPEAT   2DROP R> R> ;
: FACTOR ( ud1 u1 -- ud2 n ud3 ;  ud1 = ud2*ud3^n)
   NEXTFACTOR  DUP 1 = IF  0 2SWAP  1 ROT ROT
   ELSE  EXPONENT? SWAP 0  THEN ;
: .EXPONENT ( u)   DUP 1 <> IF ." ^" (U.) ELSE DROP THEN ;
: .FACTORS ( ud)   0  BEGIN FACTOR  2DUP (UD.)  DROP >R  .EXPONENT
   2DUP 1, D<> WHILE ." *" R> REPEAT  2DROP R> DROP ;
: T.F ( ud)   TIME  2DUP (UD.) ." ="  .FACTORS  SPACE .TIMED ;
: BEEP    " BEEP" BASICX ;
: FACTORS   ." number=" INPUTD  2DUP  D0= ABORT" 0 ?!?"
   T.F  BEEP QUIT ;
: PRIME? ( ud -- ?)   0 NEXTFACTOR  DUP 1 = >R  0 D=  R> XOR ;
: DOWNTOPRIME ( ud1 -- ud2)   BEGIN  1, D-  2DUP PRIME? UNTIL ;
: FINDPRIME   ." down from:" INPUTD  DOWNTOPRIME  BEEP UD. ;
HERE SWAP - . .( nibs )