Forth-eV Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
==== 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 (([[http://www.taygeta.com/fsl/sciforth.html|Scientific Forth Library]])) 

<code forth>
#! /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
</code>