# Forth-eV Wiki

### Webseiten-Werkzeuge

en:examples:project_euler_-_problem_p146

### Project Euler - Problem p146

Solution of the Problem and definition of forth words:

• double precision logic and shift words
• double precision math
• Miller-Rabin primality test
• random number generator 1)
```#! /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```
en/examples/project_euler_-_problem_p146.txt · Zuletzt geändert: 2013-06-06 21:26 (Externe Bearbeitung)