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))
  25 END-CODE
  26
  27 CODE f/
  28     b = reinterpret(Float64, popPS())
  29     a = reinterpret(Float64, popPS())
  30     pushPS(reinterpret(Int64, a/b))
  31 END-CODE
  32
  33 CODE f^
  34     b = reinterpret(Float64, popPS())
  35     a = reinterpret(Float64, popPS())
  36     pushPS(reinterpret(Int64, a^b))
  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))
  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 : f/mod
 149     2dup fmod -rot f/ ;
 150
 151 : 0.0
 152     [ 0 i->f ] literal ;
 153
 154 : 1.0
 155     [ 1 i->f ] literal ;
 156
 157 : -1.0
 158     [ -1 i->f ] literal ;
 159
 160 : 10.0
 161     [ 10 i->f ] literal ;
 162
 163 : flog10
 164     flog [ 10 i->f flog ] literal f/ ;
 165
 166 : floor
 167     dup 0.0 f>= if
 168         dup 1.0 fmod f-
 169     else
 170         dup 1.0 fmod dup 0.0 <> if
 171             f- 1.0 f-
 172         else
 173             drop
 174         then
 175     then
 176 ;
 177
 178 : fasin ( float -- res )
 179     dup
 180     dup f* 1.0 swap f- fsqrt
 181     f/
 182
 183     fatan
 184 ;
 185
 186 : facos ( float -- res )
 187     dup f* 1.0 swap f/ 1.0 f- fsqrt
 188     fatan
 189 ;
 190
 191 : fhead ( float -- )
 192     floor f->i
 193     0 .R  ;
 194
 195 : ftail ( float prec -- )
 196     dup 0<= if 2drop exit then
 197
 198     swap
 199
 200     1.0 fmod 10.0 f*
 201
 202     dup floor f->i 0 .R
 203
 204     1.0 fmod dup 0.0 f> if
 205         swap 1- recurse
 206     else
 207         2drop
 208     then
 209 ;
 210
 211 variable precision
 212 16 precision !
 213
 214 : f.plain ( float -- )
 215
 216     dup 0.0 = if
 217         ." 0.0"
 218         drop exit
 219     then
 220
 221     dup 0.0 f< if
 222         [char] - emit
 223         -1.0 f*
 224     then
 225
 226     dup fhead
 227
 228     [char] . emit
 229
 230     precision @ over flog10 floor f->i -
 231     ftail
 232 ;
 233
 234 : f.scientific ( float -- )
 235     dup 0.0 = if
 236         ." 0.0"
 237         drop exit
 238     then
 239
 240     dup 0.0 f< if
 241         [char] - emit
 242         -1.0 f*
 243     then
 244
 245     dup flog10 floor dup -rot
 246     10.0 swap f^ f/ f.plain
 247     ." e" f->i 0 .R
 248 ;
 249
 250 : f.nospace ( float -- )
 251     dup fabs dup 1000000 i->f f>= swap 1 i->f 10000 i->f f/ f< or if
 252         f.scientific
 253     else
 254         f.plain
 255     then
 256 ;
 257
 258 : f. ( float -- )
 259     f.nospace space ;