==== Project Euler - Problem 146 ==== Lösung des Problems und Definitionen von Forth Worten zu den Stichworten: * Doppelte Genauigkeit (double precision) * logic and shift words for double precision operands * double precision math * Miller Rabin Primzahltest * random number generator (([[http://www.taygeta.com/fsl/sciforth.html|Scientific Forth Library]])) Die Beschreibung findest du im Forth Magazin "Vierte Dimension" Heft 1/2008 im download Bereich der www.forth-ev.de #! /usr/bin/gforth \ ---- 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> ; \ ---- ran4 : a random number generator --------------------------------------- : (FuncG) ( d dc1 dc2 -- d ) 2>r dxor 2dup um* 2swap dup um* dinvert rot dup um* d+ swap 2r> dxor d+ ; : (PseudoDes) ( d d -- d d ) 2swap 2over $BAA96887E34C383B. $4B0F3B583D02B5F8. (FuncG) dxor 2swap 2over $1E17D32C39F74033. $E874F0C39226BF1A. (FuncG) dxor 2swap 2over $03BCDC3C60B43DA7. $6955C5A61D38CD47. (FuncG) dxor 2swap 2over $0F33D1B265E9215B. $55A7CA46F358B432. (FuncG) dxor ; 2variable Counter 2variable Sequence# : start-sequence ( dcounter dseq# -- ) Sequence# 2! Counter 2! ; : ran4 ( -- d ) Sequence# 2@ Counter 2@ (PseudoDes) 2swap 2drop Counter 2@ 1. d+ Counter 2! ; \ ---- 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 ; : miller_rabin { D: n W: rounds -- f } n 1. du<= if false exit then n 19. du<= if n drop \ to use single cell in the following tests dup 2 = if drop true exit then dup 3 = if drop true exit then dup 5 = if drop true exit then dup 7 = if drop true exit then dup 11 = if drop true exit then dup 13 = if drop true exit then dup 17 = if drop true exit then 19 = if true exit then \ do not dup in the last test false exit else n 2. dmod 0. d= if false exit then n 3. dmod 0. d= if false exit then n 5. dmod 0. d= if false exit then n 7. dmod 0. d= if false exit then n 11. dmod 0. d= if false exit then n 13. dmod 0. d= if false exit then n 17. dmod 0. d= if false exit then n 19. dmod 0. d= if false exit then then n 1. d- { D: d } 0 { W: s } begin d 1. dand 0. d= while d 1 drshift to d s 1+ to s repeat 0. 0. 0. 0 { D: a D: tst D: Rbig W: Rsmall } begin rounds 1- dup to rounds 0> while ran4 n 1. d- dmod 1. d+ to a a d n d**mod 1. du<> if 1. to Rbig 0 to Rsmall begin Rsmall s u< while a Rbig d d* n d**mod to tst \ tst = powmod( a, R * d, n ) n 1. d- tst du= if s to Rsmall \ cause a break from the loop else Rbig Rbig d+ to Rbig \ double Rbig Rsmall 1+ to Rsmall then repeat n 1. d- tst du<> if false exit \ return false (composite number) then then repeat true \ return true (maybe prime) ; \ ---- main and its internal words -------------------------------------------- : (main-inc-p) ( u -- u ) dup 2 + 3 mod 0= if 4 else 2 then + ; : (main-ba) { W: n -- n' } 11 { W: p } 0 { W: nd } n s>d 2dup d* { D: n2 } begin n2 p s>d dmod d>s to nd nd 1+ p mod 0= if 0 exit then nd 3 + p mod 0= if 0 exit then nd 7 + p mod 0= if 0 exit then nd 9 + p mod 0= if 0 exit then nd 13 + p mod 0= if 0 exit then nd 27 + p mod 0= if 0 exit then p (main-inc-p) to p \ inc. p by 2 or by 4 p n 1+ > if n2 15. d+ 10 miller_rabin if 0 exit then n2 19. d+ 10 miller_rabin if 0 exit then n2 21. d+ 10 miller_rabin if 0 exit then n2 25. d+ 10 miller_rabin if 0 exit then ." Found n=" n . cr n exit then again ; : main { W: limit -- d } 0. { D: sum } limit 1+ 10 do \ i must be divisible by 10, loop only through multiples of 10 i 3 mod 0<> if \ i can't be divisible by 3 i 4 + 7 mod 1 <= if \ i % 7 must be either 3 or 4 i (main-ba) s>d sum d+ to sum then then 10 +loop ." sum=" sum ud. cr ; 150000000 main bye