\ Double Number Word Set \ by Luca Masini 2008) \ ---- logic and shift words for double cell operands ------------------------- : dand ( d d -- d ) rot and >r and r> ; : dinvert ( d -- d ) swap invert swap invert ; : dlshift ( d u -- d ) 0 ?do d2* loop ; : dor ( d d -- d ) rot or >r or r> ; : dnot ( d -- 0.|1. ) d0= abs s>d ; : drshift ( d u -- d ) 0 ?do d2/ loop dabs ; : dxor ( d d -- d ) rot xor >r xor r> ; \ ---- arithmetic words for unsigned double cell operands --------------------- : d* ( d d -- d ) 3 pick * >r tuck * >r um* r> + r> + ; : t* ( ud un -- ut ) dup rot um* 2>r um* 0 2r> d+ ; : t/ ( ut un -- ud ) >r r@ um/mod swap rot 0 r@ um/mod swap rot r> um/mod swap drop 0 2swap swap d+ ; : u*/ ( ud un un -- ud ) >r t* r> t/ ; \ initialize SUPERBASE (normally $100000000. on 32 bits machine, with fallback value of $10000.) s" max-u" environment? 0= [if] $10000. [then] 0 1. d+ 2constant SUPERBASE : d/ { D: u D: v -- ud_quot } v 0. d= if -10 throw then \ throw for division by zero v u du> if 0. exit then \ if v is bigger then 0. v u d= if 1. exit then \ if v is equal then 1. v 0= if >r u 1 r> u*/ exit then \ use mixed precision word drop v { v1 v0 } \ v = v0 * b + v1 v0 -1 = if 1 else SUPERBASE 1 v0 1+ u*/ drop then { d } \ d = b/(v0+1) v d 1 u*/ { w1 w0 } \ w = d * v = w0 * b + w1 u over 0 w1 w0 u*/ d- d w0 u*/ nip 0 ; : dmod { D: d1 D: d2 -- d } d1 2dup d2 d/ d2 d* d- ; : dumin ( d1 d2 -- d ) 2over 2over du> if 2swap then 2drop ; \ d is min(d1,d2) \ initialize UMAX (normally $ffffffffffffffff. on 32 bits machine, with fallback value of $ffffffff.) s" max-ud" environment? 0= [if] $ffffffff. [then] 2constant UMAX : d+mod { D: a D: b D: m -- d } \ addition modulo m a m dmod to a \ normalize a b m dmod to b \ normalize b a UMAX b d- d<= if \ no overflow.. a b d+ m dmod exit \ ..built-in computation then \ ---- go with the algorithm ;-) a m a d- dumin { D: aA } b m b d- dumin { D: bB } b bB d= ( -- f ) \ leave a flag on stack a aA d= if bB aA du> if if aA bB d+ m dmod else m bB aA d- d- then exit ( f -- ) \ ..consume the flag else >r aA bB r> if d+ else d- then m dmod exit ( f -- ) \ ..consume the flag then else if aA bB du> if m aA bB else m m bB aA d- then d- d- exit then ( f -- ) \ ..consume the flag then m aA bB d+ m dmod d- ; : d*mod { D: a D: b D: m -- d } \ multiplication modulo m a m dmod to a \ normalize a b m dmod to b \ normalize b a 1. d= if b exit then b 1. d= if a exit then a m a d- dumin { D: aA } b m b d- dumin { D: bB } aA d0= bB d0= or if 0. exit then aA 1. d= if m b d- exit then bB 1. d= if m a d- exit then aA a d= bB b d= and aA a d<> bB b d<> and or { pos } \ pos is True if positive, False otherwise aA UMAX bB d/ du<= if aA bB d* m dmod pos 0= if m 2swap d- m dmod then exit then aA d2/ { D: a0 } aA a0 d- { D: a1 } bB d2/ { D: b0 } bB b0 d- { D: b1 } a1 b1 m recurse { D: p4 } 0. 0. 0. { D: p1 D: p2 D: p3 } a0 a1 d= b0 b1 d= and if p4 to p1 p4 to p2 p4 to p3 else a0 a1 d= if p4 m a1 d- m dmod m d+mod 2dup to p3 to p1 p4 to p2 else p4 m b1 d- m dmod m d+mod to p2 b0 b1 d= if p2 to p1 p4 to p3 else p4 m a1 d- m dmod m b1 d- m dmod 1. m d+mod m d+mod m d+mod to p1 p4 m a1 d- m dmod m d+mod to p3 then then then p1 p2 p3 p4 m d+mod m d+mod m d+mod pos 0= if m 2swap d- m dmod then ; : d**mod { D: base D: power D: m -- d } \ exponentiation modulo m 1. { D: res } begin power 0. du> while 1. power dand drop if \ if power is odd res base m d*mod to res then base base m d*mod to base power 1 drshift to power repeat res ; .( double number wordset included) cr ( finis)