src/float.4th

   1 \ Floating point arithmetic
   2
   3 ( Cheating for now by using forth.jl CODE/END-CODE to
   4   access Julia's floating point support.  This isn't
   5   at all portable.  That said, the year is 2016 and
   6   I've only cheated for words that have corresponding
   7   x87 FPU instructions, so I don't feel too bad! )
   8
   9 CODE f+
  10     b = reinterpret(Float64, popPS())
  11     a = reinterpret(Float64, popPS())
  12     pushPS(reinterpret(Int64, a+b))
  13 END-CODE
  14
  15 CODE f-
  16     b = reinterpret(Float64, popPS())
  17     a = reinterpret(Float64, popPS())
  18     pushPS(reinterpret(Int64, a-b))
  19 END-CODE
  20
  21 CODE f*
  22     b = reinterpret(Float64, popPS())
  23     a = reinterpret(Float64, popPS())
  24     pushPS(reinterpret(Int64, a*b + 0.0))
  25 END-CODE
  26
  27 CODE f/
  28     b = reinterpret(Float64, popPS())
  29     a = reinterpret(Float64, popPS())
  30     pushPS(reinterpret(Int64, a/b + 0.0))
  31 END-CODE
  32
  33 CODE f^
  34     b = reinterpret(Float64, popPS())
  35     a = reinterpret(Float64, popPS())
  36     pushPS(reinterpret(Int64, a^b + 0.0))
  37 END-CODE
  38
  39 CODE fsqrt
  40     a = reinterpret(Float64, popPS())
  41     pushPS(reinterpret(Int64, sqrt(a)))
  42 END-CODE
  43
  44 CODE f>
  45     b = reinterpret(Float64, popPS())
  46     a = reinterpret(Float64, popPS())
  47     if a > b
  48         pushPS(-1)
  49     else
  50         pushPS(0)
  51     end
  52 END-CODE
  53
  54 CODE f<
  55     b = reinterpret(Float64, popPS())
  56     a = reinterpret(Float64, popPS())
  57     if a < b
  58         pushPS(-1)
  59     else
  60         pushPS(0)
  61     end
  62 END-CODE
  63
  64 CODE f=
  65     b = reinterpret(Float64, popPS())
  66     a = reinterpret(Float64, popPS())
  67     if a == b
  68         pushPS(-1)
  69     else
  70         pushPS(0)
  71     end
  72 END-CODE
  73
  74 : f<=
  75     f> invert ;
  76
  77 : f>=
  78     f< invert ;
  79
  80 CODE fmod
  81     b = reinterpret(Float64, popPS())
  82     a = reinterpret(Float64, popPS())
  83     pushPS(reinterpret(Int64, a%b + 0.0))
  84 END-CODE
  85
  86 CODE flog
  87     a = reinterpret(Float64, popPS())
  88     pushPS(reinterpret(Int64, log(a)))
  89 END-CODE
  90
  91 CODE fexp
  92     a = reinterpret(Float64, popPS())
  93     pushPS(reinterpret(Int64, exp(a)))
  94 END-CODE
  95
  96 CODE fsin
  97     a = reinterpret(Float64, popPS())
  98     pushPS(reinterpret(Int64, sin(a)))
  99 END-CODE
 100
 101 CODE fcos
 102     a = reinterpret(Float64, popPS())
 103     pushPS(reinterpret(Int64, cos(a)))
 104 END-CODE
 105
 106 CODE ftan
 107     a = reinterpret(Float64, popPS())
 108     pushPS(reinterpret(Int64, tan(a)))
 109 END-CODE
 110
 111 CODE fatan
 112     a = reinterpret(Float64, popPS())
 113     pushPS(reinterpret(Int64, atan(a)))
 114 END-CODE
 115
 116 CODE fnan?
 117     a = reinterpret(Float64, popPS())
 118     if isnan(a)
 119         pushPS(-1)
 120     else
 121         pushPS(0)
 122     end
 123 END-CODE
 124
 125 CODE finf?
 126     a = reinterpret(Float64, popPS())
 127     if isinf(a)
 128         pushPS(-1)
 129     else
 130         pushPS(0)
 131     end
 132 END-CODE
 133
 134 CODE i->f
 135     pushPS(reinterpret(Int64, Float64(popPS())))
 136 END-CODE
 137
 138 CODE f->i
 139     a = reinterpret(Float64, popPS())
 140     pushPS(Int64(round(a)))
 141 END-CODE
 142
 143 CODE fabs
 144     a = reinterpret(Float64, popPS())
 145     pushPS(reinterpret(Int64, abs(a)))
 146 END-CODE
 147
 148 CODE fround
 149     a = reinterpret(Float64, popPS())
 150     pushPS(reinterpret(Int64, round(a)))
 151 END-CODE
 152
 153 : f/mod
 154     2dup fmod -rot f/ ;
 155
 156 : 0.0
 157     [ 0 i->f ] literal ;
 158
 159 : 1.0
 160     [ 1 i->f ] literal ;
 161
 162 : -1.0
 163     [ -1 i->f ] literal ;
 164
 165 : 10.0
 166     [ 10 i->f ] literal ;
 167
 168 : flog10
 169     flog [ 10 i->f flog ] literal f/ ;
 170
 171 : truncate
 172     dup 1.0 fmod f- ;
 173
 174 : floor
 175     dup 0.0 f>= if
 176         truncate
 177     else
 178         dup 1.0 fmod dup 0.0 <> if
 179             f- 1.0 f-
 180         else
 181             drop
 182         then
 183     then
 184 ;
 185
 186 : ceiling
 187     dup 0.0 f>= if
 188         dup 1.0 fmod dup 0.0 <> if
 189             f- 1.0 f+
 190         else
 191             drop
 192         then
 193     else
 194         truncate
 195     then
 196 ;
 197
 198 : fasin ( float -- res )
 199     dup
 200     dup f* 1.0 swap f- fsqrt
 201     f/
 202
 203     fatan
 204 ;
 205
 206 : facos ( float -- res )
 207     dup f* 1.0 swap f/ 1.0 f- fsqrt
 208     fatan
 209 ;
 210
 211 : fhead ( float -- )
 212     floor f->i
 213     0 .R  ;
 214
 215 : ftail ( float prec -- )
 216     dup 0<= if 2drop exit then
 217
 218     swap
 219
 220     1.0 fmod 10.0 f*
 221
 222     dup floor f->i 0 .R
 223
 224     1.0 fmod dup 0.0 f> if
 225         swap 1- recurse
 226     else
 227         2drop
 228     then
 229 ;
 230
 231 variable precision
 232 16 precision !
 233
 234 : f.plain ( float -- )
 235
 236     dup 0.0 = if
 237         ." 0.0"
 238         drop exit
 239     then
 240
 241     dup 0.0 f< if
 242         [char] - emit
 243         -1.0 f*
 244     then
 245
 246     dup fhead
 247
 248     [char] . emit
 249
 250     precision @ over flog10 floor f->i -
 251     ftail
 252 ;
 253
 254 : f.scientific ( float -- )
 255     dup 0.0 = if
 256         ." 0.0"
 257         drop exit
 258     then
 259
 260     dup 0.0 f< if
 261         [char] - emit
 262         -1.0 f*
 263     then
 264
 265     dup flog10 floor dup -rot
 266     10.0 swap f^ f/ f.plain
 267     ." e" f->i 0 .R
 268 ;
 269
 270 : f.nospace ( float -- )
 271     dup fabs dup 1000000 i->f f>= swap 1 i->f 10000 i->f f/ f< or if
 272         f.scientific
 273     else
 274         f.plain
 275     then
 276 ;
 277
 278 : f. ( float -- )
 279     f.nospace space ;