Forth-eV Wiki

\ 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)