git.alexw.nyc home about git garden
    1
    2
    3
    4
    5
    6
    7
    8
    9
   10
   11
   12
   13
   14
   15
   16
   17
   18
   19
   20
   21
   22
   23
   24
   25
   26
   27
   28
   29
   30
   31
   32
   33
   34
   35
   36
   37
   38
   39
   40
   41
   42
   43
   44
   45
   46
   47
   48
   49
   50
   51
   52
   53
   54
   55
   56
   57
   58
   59
   60
   61
\ double arithmetic
\ copied from https://gitlab.com/b2495/uf/-/blob/master/lib/dmath.f?ref_type=heads
\ somewhat sloppy
\ TODO deal with signed numbers the 'dusk os' way
\ TODO annotate properly

\ d == double ( 2 cells ) 

: abs dup 0< if -1 * then ;
: 2dup over over ;
: 2nip rot drop rot drop ;
: 2swap >r rot> r> rot> ;

: d=  ( n n n n -- f )
    >r rot xor swap r> xor or 0 = ;
: d+ ( augend . addend . -- sum . )
    rot + >r over +  dup rot < if  r> 1+ else r> then ;
: s>d dup 0<  if  -1  else  0  then ;
: dnegate
    ^ swap ^ swap 1 0 d+ ;
: dabs ( d -- ud ) dup 0< if dnegate then ;
: d- dnegate d+ ;

: d< ( d d -- f ) rot 2dup = if 2drop < else 2nip > then ;
: d0 ( d -- f ) = or 0 = ;
: d0< ( d -- f ) 0 0 d< ;
: d<< ( d -- d ) 2dup d+ ;
: d2/ ( d -- d ) dup 31 lshift >r >> swap >> r> or swap ;
: tf ( n -- f ) if  -1  else  0  then ;

0 value scratch
: um*  ( u1 u2 -- ud )
    to scratch 
    0 0 32 for 2dup d+ rot dup 0< if << rot> scratch 0 d+
    else << rot> then next rot drop ;
: abssgn    ( a b -- |a| |b| negf )
    2dup xor 0< tf >r abs swap abs swap r> ;
: m*  ( n1 n2 -- d ) abssgn >r um* r> if dnegate then ;
: divstep
    ( divisor dq hi )
    << over 0< if 1+ then
    swap << swap rot
    2dup >= if tuck - swap rot 1+            
    rot else rot> then ;

: um/mod ( ud u1 -- u2 u3 )
    rot> 32 for divstep next 
    rot drop swap ;
: sm/rem ( d n -- r q )  ( symmetric )
  over >r >r  dabs r@ abs um/mod
  r> r@ xor 0< if ^ then  r> 0< if >r ^ r> then ;
: */mod >r m* r> sm/rem ;
: */    */mod nip ;
: t<<  over >r >r d<< r> << r> 0< tf 1 and + ;

0 value divisor
: m*/mod ( d n n -- d ) 
    to divisor 
    tuck um* 2swap um* swap >r 0 d+ r>   rot>             
    64 for t<< dup divisor >= if divisor - rot 1+ rot> then next ;
: m*/  ( d1 n1 +n2 -- d2 ) m*/mod drop ;