(*
 * LANGUAGE    : ANS Forth 
 * PROJECT     : Number Theory
 * CATEGORY    : The infinite precision prime test for BIGNUMbers
 * AUTHOR      : Albert van der Horst, 911108
 * LAST CHANGE : ah 911108 cosmetic
 * LAST CHANGE : Marcel Hendrix, 19 March 2001 cosmetic
 * LAST CHANGE : Marcel Hendrix, 24 March 2001 bug fix in RABIN? ( CHOOSE )
 *)



	NEEDS -bignum

	REVISION -rabin "ÄÄÄ M.O Rabin's oracle  Version 2.01 ÄÄÄ"

	PRIVATES


MAX.DIGITS BIGNUM TheNumber

-- Just an alias "most significant bit is up"

: MSB 	EVAL" 0< " ;P IMMEDIATE


DOC RABIN
(*
	 RR{new} = RR{old}^ (2^32) * VV^U modulo MM
	 Perform BIT/CELL (typical 32) bits of an exponentiation modulo MM
	 of a Vnumber VV to the bits of a single cell power U
	 RR is the intermediate product.
	 Loop invariant escapes me but it IS correct.
	 (Something like RR*VV^(TOS>>I) modulo MM )
	 e.g. for V^S   starts with 1*VV^U, ends with RR*VV^0
	 Infinite loop whenever MM is not `REPACK'ed
	 Result in RR ( remainder )
	 IF THIS IS DIFFICULT TRY TO UNDERSTAND IN COMBINATION WITH V^S
*)
ENDDOC


-- Assumes RR is initialized!

: V^bits ( VV MM U -- )
	SWAP LOCAL MM
	SWAP LOCAL VV
	BITS/CELL  
	0 DO	RR RR VV*		\ Square: Result in PP
		PP MM VV/MOD		\ Result in RR
		DUP MSB IF		\ If bit up
		   	   RR VV VV*	\ Also multiply by VV : Result in PP
	 	   	   PP MM VV/MOD	\ Result in RR
	     	     ENDIF
	        1 LSHIFT
	LOOP DROP ;P


-- RR = VV^U modulo MM
-- Put a Vnumber VV to a single cell power U modulo MM
-- Infinite loop whenever MM is not `REPACK'ed
-- Result in RR ( remainder )

: V^S ( VV MM U -- )
	1 RR V! 
	V^bits ;


-- RR = VV^UU modulo MM
-- Put a Vnumber VV to the power UU modulo MM
-- Result in RR ( remainder )
-- Infinite loop whenever MM is not `REPACK'ed

: V^V ( VV UU MM -- )
	LOCAL MM
	LOCAL UU
	LOCAL VV
	1 RR V!
	UU SIZE@ 
	0 ?DO	VV  MM  I UU []HEAD @	\ Prepare parameters
		V^bits
	 LOOP ;


-- V{in} = V{out}*2^exponent
-- where V{out} modulo 2 is 0
-- Convert Vnum to yet another exponent, mantissa format
-- MH: a bit inefficient. Better: check for amount of zero bits (up to 32/loop)

: EVEN-NORMALIZE  		\ ( V -- exponent )
	0 LOCAL exp
	  LOCAL V
	V REPACK  V VS0= ABORT" even-normalize :: 0"
	BEGIN
	  V 'TAIL@ @ 		\ least significant cell might move ..
	  1 AND 0= 		\ .. with V>> (see REPACK) [bug 24 March 2001 -mhx]
	WHILE			\ l.s. bit 0
	  1 +TO exp
	  V 1 V>>		\ does a REPACK
	REPEAT 
        exp ;P


-- Rabin? tells us whether TheNumber **might be** prime
-- FALSE : certainly not
-- TRUE : probably prime

MAX.DIGITS BIGNUM PROBE  PRIVATE

: RABIN?  ( TheNumber -- flag )
	LOCAL TheNumber
	TheNumber TO SS
	SS 1 VS-		\ The random number put to this power
	SS EVEN-NORMALIZE 
	  LOCAL two's		\ Separate power of two
	TheNumber 
	SIZE@ 1 = IF TheNumber 'TAIL@ @ CHOOSE
	        ELSE RANDOM
	       ENDIF PROBE V!	\ Start PROBE up with a random number < TheNumber
	PROBE SS TheNumber V^V	\ Result in RR
	#DIGITS RR 1 =  
	0 TAIL  RR 1 = 
	AND IF			\ RR equal 1?
		TRUE EXIT	\ Prime!
	 ENDIF

	BEGIN
	  RR 1 VS+
	  RR TheNumber VV= 
	     IF
		TRUE EXIT	\ Result -1 mod N -> probably prime
    	  ENDIF
	  two's 
	WHILE
	  RR 1 VS-		\ Restore RR
	  RR RR VV*		\ Square : Result in PP
	  PP TheNumber VV/MOD	\ Result in RR
	  -1 +TO two's
	REPEAT
	FALSE ;			\ Not prime


:ABOUT	CR ." M.O. Rabin's prime probability tester."
	CR .~ Enter  "TheNumber RABIN? ."~
	CR ." The following numbers are prime and can be used to test the above:"
	CR ."     TheNumber VREAD 5704689200685129054721"
	CR ."     TheNumber VREAD 37866809061660057264219253397"
	CR ."     TheNumber VREAD 59649589127497217" ;

	DEPRIVE
	.ABOUT -rabin

				(* End of File *)