git.gag.com Git - debian/pforth/blob - fth/floats.fth

   1 \ @(#) floats.fth 98/02/26 1.4 17:51:40
   2 \ High Level Forth support for Floating Point
   3 \
   4 \ Author: Phil Burk and Darren Gibbs
   5 \ Copyright 1994 3DO, Phil Burk, Larry Polansky, David Rosenboom
   6 \
   7 \ Permission to use, copy, modify, and/or distribute this
   8 \ software for any purpose with or without fee is hereby granted.
   9 \
  10 \ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
  11 \ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
  12 \ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL
  13 \ THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
  14 \ CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
  15 \ FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF
  16 \ CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  17 \ OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  18 \
  19 \ 19970702 PLB Drop 0.0 in REPRESENT to fix  0.0 F.
  20 \ 19980220 PLB Added FG. , fixed up large and small formatting
  21 \ 19980812 PLB Now don't drop 0.0 in REPRESENT to fix  0.0 F.  (!!!)
  22 \              Fixed F~ by using (F.EXACTLY)
  23
  24 ANEW TASK-FLOATS.FTH
  25
  26 : FALIGNED  ( addr -- a-addr )
  27     1 floats 1- +
  28     1 floats /
  29     1 floats *
  30 ;
  31
  32 : FALIGN    ( -- , align DP )
  33     dp @ faligned dp !
  34 ;
  35
  36 \ account for size of create when aligning floats
  37 here
  38 create fp-create-size
  39 fp-create-size swap - constant CREATE_SIZE
  40
  41 : FALIGN.CREATE  ( -- , align DP for float after CREATE )
  42     dp @
  43     CREATE_SIZE +
  44     faligned
  45     CREATE_SIZE -
  46     dp !
  47 ;
  48
  49 : FCREATE  ( <name> -- , create with float aligned data )
  50     falign.create
  51     CREATE
  52 ;
  53
  54 : FVARIABLE ( <name> -- ) ( F: -- )
  55     FCREATE 1 floats allot
  56 ;
  57
  58 : FCONSTANT
  59     FCREATE here   1 floats allot   f!
  60     DOES> f@
  61 ;
  62
  63 : F0SP ( -- ) ( F: ? -- )
  64     fdepth 0 max  0 ?DO fdrop LOOP
  65 ;
  66
  67 \ Floating point structure member.
  68 : FFIELD:   ( n1 "name" -- n2 ; addr1 -- addr2 )
  69     FALIGNED 1 FLOATS +FIELD
  70 ;
  71
  72 \ Convert between single precision integer and floating point
  73 : S>F ( s -- ) ( F: -- r )
  74     s>d d>f
  75 ;
  76 : F>S ( -- s ) ( F: r -- )
  77     f>d d>s
  78 ;
  79
  80 : (F.EXACTLY) ( r1 r2 -f- flag , return true if encoded equally ) { | caddr1 caddr2 fsize fcells }
  81     1 floats -> fsize
  82     fsize cell 1- + cell 1- invert and  \ round up to nearest multiple of stack size
  83     cell / -> fcells ( number of cells per float )
  84 \ make room on data stack for floats data
  85     fcells 0 ?DO 0 LOOP
  86     sp@ -> caddr1
  87     fcells 0 ?DO 0 LOOP
  88     sp@ -> caddr2
  89 \ compare bit representation
  90     caddr1 f!
  91     caddr2 f!
  92     caddr1 fsize caddr2 fsize compare 0=
  93     >r fcells 2* 0 ?DO drop LOOP r>  \ drop float bits
  94 ;
  95
  96 : F~ ( -0- flag ) ( r1 r2 r3 -f- )
  97     fdup F0<
  98     IF
  99         frot frot  ( -- r3 r1 r2 )
 100         fover fover ( -- r3 r1 r2 r1 r2 )
 101         f- fabs    ( -- r3 r1 r2 |r1-r2| )
 102         frot frot  ( -- r3  |r1-r2| r1 r2 )
 103         fabs fswap fabs f+ ( -- r3 |r1-r2|  |r1|+|r2| )
 104         frot fabs f* ( -- |r1-r2|  |r1|+|r2|*|r3| )
 105         f<
 106     ELSE
 107         fdup f0=
 108         IF
 109             fdrop
 110             (f.exactly)  \ f- f0=  \ 19980812 Used to cheat. Now actually compares bit patterns.
 111         ELSE
 112             frot frot  ( -- r3 r1 r2 )
 113             f- fabs    ( -- r3 |r1-r2| )
 114             fswap f<
 115         THEN
 116     THEN
 117 ;
 118
 119 \ FP Output --------------------------------------------------------
 120 fvariable FVAR-REP  \ scratch var for represent
 121 : REPRESENT { c-addr u | n flag1 flag2 --  n flag1 flag2 , FLOATING } ( F: r -- )
 122     TRUE -> flag2   \ FIXME - need to check range
 123     fvar-rep f!
 124 \
 125     fvar-rep f@ f0<
 126     IF
 127         -1 -> flag1
 128         fvar-rep f@ fabs fvar-rep f!   \ absolute value
 129     ELSE
 130         0 -> flag1
 131     THEN
 132 \
 133     fvar-rep f@ f0=
 134     IF
 135 \       fdrop \ 19970702 \ 19980812 Remove FDROP to fix "0.0 F."
 136         c-addr u [char] 0 fill
 137         0 -> n
 138     ELSE
 139         fvar-rep f@
 140         flog
 141         fdup f0< not
 142         IF
 143             1 s>f f+ \ round up exponent
 144         THEN
 145         f>s -> n
 146 \ ." REP - n = " n . cr
 147 \ normalize r to u digits
 148         fvar-rep f@
 149         10 s>f u n - s>f f** f*
 150         1 s>f 2 s>f f/ f+   \ round result
 151 \
 152 \ convert float to double_int then convert to text
 153         f>d
 154 \ ." REP - d = " over . dup . cr
 155         <#  u 1- 0 ?DO # loop #s #>  \ ( -- addr cnt )
 156 \ Adjust exponent if rounding caused number of digits to increase.
 157 \ For example from 9999 to 10000.
 158         u - +-> n
 159         c-addr u move
 160     THEN
 161 \
 162     n flag1 flag2
 163 ;
 164
 165 variable FP-PRECISION
 166
 167 \ Set maximum digits that are meaningful for the precision that we use.
 168 1 FLOATS 4 / 7 * constant FP_PRECISION_MAX
 169
 170 : PRECISION ( -- u )
 171     fp-precision @
 172 ;
 173 : SET-PRECISION ( u -- )
 174     fp_precision_max min
 175     fp-precision !
 176 ;
 177 7 set-precision
 178
 179 32 constant FP_REPRESENT_SIZE
 180 64 constant FP_OUTPUT_SIZE
 181
 182 create FP-REPRESENT-PAD FP_REPRESENT_SIZE allot  \ used with REPRESENT
 183 create FP-OUTPUT-PAD FP_OUTPUT_SIZE allot     \ used to assemble final output
 184 variable FP-OUTPUT-PTR            \ points into FP-OUTPUT-PAD
 185
 186 : FP.HOLD ( char -- , add char to output )
 187     fp-output-ptr @ fp-output-pad 64 + <
 188     IF
 189         fp-output-ptr @ tuck c!
 190         1+ fp-output-ptr !
 191     ELSE
 192         drop
 193     THEN
 194 ;
 195 : FP.APPEND { addr cnt -- , add string to output }
 196     cnt 0 max   0
 197     ?DO
 198         addr i + c@ fp.hold
 199     LOOP
 200 ;
 201
 202 : FP.STRIP.TRAILING.ZEROS ( -- , remove trailing zeros from fp output )
 203     BEGIN
 204         fp-output-ptr @ fp-output-pad u>
 205         fp-output-ptr @ 1- c@ [char] 0 =
 206         and
 207     WHILE
 208         -1 fp-output-ptr +!
 209     REPEAT
 210 ;
 211
 212 : FP.APPEND.ZEROS ( numZeros -- )
 213     0 max   0
 214     ?DO [char] 0 fp.hold
 215     LOOP
 216 ;
 217
 218 : FP.MOVE.DECIMAL   { n prec -- , append with decimal point shifted }
 219     fp-represent-pad n prec min fp.append
 220     n prec - fp.append.zeros
 221     [char] . fp.hold
 222     fp-represent-pad n +
 223     prec n - 0 max fp.append
 224 ;
 225
 226 : (EXP.) ( n -- addr cnt , convert exponent to two digit value )
 227     dup abs 0
 228     <# # #s
 229     rot 0<
 230     IF [char] - HOLD
 231     ELSE [char] + hold
 232     THEN
 233     #>
 234 ;
 235
 236 : FP.REPRESENT ( -- n flag1 flag2 ) ( r -f- )
 237 ;
 238
 239 : (FS.)  ( -- addr cnt ) ( F: r -- , scientific notation )
 240     fp-output-pad fp-output-ptr !  \ setup pointer
 241     fp-represent-pad   precision  represent
 242 \ ." (FS.) - represent " fp-represent-pad precision type cr
 243     ( -- n flag1 flag2 )
 244     IF
 245         IF [char] - fp.hold
 246         THEN
 247         1 precision fp.move.decimal
 248         [char] e fp.hold
 249         1- (exp.) fp.append \ n
 250     ELSE
 251         2drop
 252         s" <FP-OUT-OF-RANGE>" fp.append
 253     THEN
 254     fp-output-pad fp-output-ptr @ over -
 255 ;
 256
 257 : FS.  ( F: r -- , scientific notation )
 258     (fs.) type space
 259 ;
 260
 261 : (FE.)  ( -- addr cnt ) ( F: r -- , engineering notation ) { | n n3 -- }
 262     fp-output-pad fp-output-ptr !  \ setup pointer
 263     fp-represent-pad precision represent
 264     ( -- n flag1 flag2 )
 265     IF
 266         IF [char] - fp.hold
 267         THEN
 268 \ convert exponent to multiple of three
 269         -> n
 270         n 1- s>d 3 fm/mod \ use floored divide
 271         3 * -> n3
 272         1+ precision fp.move.decimal \ amount to move decimal point
 273         [char] e fp.hold
 274         n3 (exp.) fp.append \ n
 275     ELSE
 276         2drop
 277         s" <FP-OUT-OF-RANGE>" fp.append
 278     THEN
 279     fp-output-pad fp-output-ptr @ over -
 280 ;
 281
 282 : FE.  ( F: r -- , engineering notation )
 283     (FE.) type space
 284 ;
 285
 286 : (FG.)  ( F: r -- , normal or scientific ) { | n n3 ndiff -- }
 287     fp-output-pad fp-output-ptr !  \ setup pointer
 288     fp-represent-pad precision represent
 289     ( -- n flag1 flag2 )
 290     IF
 291         IF [char] - fp.hold
 292         THEN
 293 \ compare n with precision to see whether we do scientific display
 294         dup precision >
 295         over -3 < OR
 296         IF  \ use exponential notation
 297             1 precision fp.move.decimal
 298             fp.strip.trailing.zeros
 299             [char] e fp.hold
 300             1- (exp.) fp.append \ n
 301         ELSE
 302             dup 0>
 303             IF
 304 \ POSITIVE EXPONENT - place decimal point in middle
 305                 precision fp.move.decimal
 306             ELSE
 307 \ NEGATIVE EXPONENT - use 0.000????
 308                 s" 0." fp.append
 309 \ output leading zeros
 310                 negate fp.append.zeros
 311                 fp-represent-pad precision fp.append
 312             THEN
 313             fp.strip.trailing.zeros
 314         THEN
 315     ELSE
 316         2drop
 317         s" <FP-OUT-OF-RANGE>" fp.append
 318     THEN
 319     fp-output-pad fp-output-ptr @ over -
 320 ;
 321
 322 : FG.  ( F: r -- )
 323     (fg.) type space
 324 ;
 325
 326 : (F.)  ( F: r -- , normal or scientific ) { | n n3 ndiff prec' -- }
 327     fp-output-pad fp-output-ptr !  \ setup pointer
 328     fp-represent-pad  \ place to put number
 329     fdup flog 1 s>f f+ f>s precision max
 330     fp_precision_max min dup -> prec'
 331     represent
 332     ( -- n flag1 flag2 )
 333     IF
 334 \ add '-' sign if negative
 335         IF [char] - fp.hold
 336         THEN
 337 \ compare n with precision to see whether we must do scientific display
 338         dup fp_precision_max >
 339         IF  \ use exponential notation
 340             1 precision fp.move.decimal
 341             fp.strip.trailing.zeros
 342             [char] e fp.hold
 343             1- (exp.) fp.append \ n
 344         ELSE
 345             dup 0>
 346             IF
 347     \ POSITIVE EXPONENT - place decimal point in middle
 348                 prec' fp.move.decimal
 349             ELSE
 350     \ NEGATIVE EXPONENT - use 0.000????
 351                 s" 0." fp.append
 352     \ output leading zeros
 353                 dup negate precision min
 354                 fp.append.zeros
 355                 fp-represent-pad precision rot + fp.append
 356             THEN
 357         THEN
 358     ELSE
 359         2drop
 360         s" <FP-OUT-OF-RANGE>" fp.append
 361     THEN
 362     fp-output-pad fp-output-ptr @ over -
 363 ;
 364
 365 : F.  ( F: r -- )
 366     (f.) type space
 367 ;
 368
 369 : F.S  ( -- , print FP stack )
 370     ." FP> "
 371     fdepth 0>
 372     IF
 373         fdepth 0
 374         DO
 375             cr?
 376             fdepth i - 1-  \ index of next float
 377             fpick f. cr?
 378         LOOP
 379     ELSE
 380         ." empty"
 381     THEN
 382     cr
 383 ;
 384
 385 \ FP Input ----------------------------------------------------------
 386 variable FP-REQUIRE-E   \ must we put an E in FP numbers?
 387 false fp-require-e !   \ violate ANSI !!
 388
 389 : >FLOAT { c-addr u | dlo dhi u' fsign flag nshift -- flag }
 390     u 0= IF false exit THEN
 391     false -> flag
 392     0 -> nshift
 393 \
 394 \ check for minus sign
 395     c-addr c@ [char] - =     dup -> fsign
 396     c-addr c@ [char] + = OR
 397     IF   1 +-> c-addr   -1 +-> u   \ skip char
 398     THEN
 399 \
 400 \ convert first set of digits
 401     0 0 c-addr u >number -> u' -> c-addr -> dhi -> dlo
 402     u' 0>
 403     IF
 404 \ convert optional second set of digits
 405         c-addr c@ [char] . =
 406         IF
 407             dlo dhi c-addr 1+ u' 1- dup -> nshift >number
 408             dup nshift - -> nshift
 409             -> u' -> c-addr -> dhi -> dlo
 410         THEN
 411 \ convert exponent
 412         u' 0>
 413         IF
 414             c-addr c@ [char] E =
 415             c-addr c@ [char] e =  OR
 416             IF
 417                 1 +-> c-addr   -1 +-> u'   \ skip E char
 418                 u' 0>
 419                 IF
 420                     c-addr c@ [char] + = \ ignore + on exponent
 421                     IF
 422                         1 +-> c-addr   -1 +-> u'   \ skip char
 423                     THEN
 424                     c-addr u' ((number?))
 425                     num_type_single =
 426                     IF
 427                        nshift + -> nshift
 428                        true -> flag
 429                     THEN
 430                 ELSE
 431                     true -> flag   \ allow "1E"
 432                 THEN
 433             THEN
 434         ELSE
 435 \ only require E field if this variable is true
 436             fp-require-e @ not -> flag
 437         THEN
 438     THEN
 439 \ convert double precision int to float
 440     flag
 441     IF
 442         dlo dhi d>f
 443         10 s>f nshift s>f f** f*   \ apply exponent
 444         fsign
 445         IF
 446             fnegate
 447         THEN
 448     THEN
 449     flag
 450 ;
 451
 452 3 constant NUM_TYPE_FLOAT   \ possible return type for NUMBER?
 453
 454 : (FP.NUMBER?)   ( $addr -- 0 | n 1 | d 2 | r 3 , convert string to number )
 455 \ check to see if it is a valid float, if not use old (NUMBER?)
 456     dup count >float
 457     IF
 458         drop NUM_TYPE_FLOAT
 459     ELSE
 460         (number?)
 461     THEN
 462 ;
 463
 464 defer fp.old.number?
 465 variable FP-IF-INIT
 466
 467 : FP.TERM    ( -- , deinstall fp conversion )
 468     fp-if-init @
 469     IF
 470         what's  fp.old.number? is number?
 471         fp-if-init off
 472     THEN
 473 ;
 474
 475 : FP.INIT  ( -- , install FP converion )
 476     fp.term
 477     what's number? is fp.old.number?
 478     ['] (fp.number?) is number?
 479     fp-if-init on
 480     ." Floating point numeric conversion installed." cr
 481 ;
 482
 483 FP.INIT
 484 if.forgotten fp.term
 485
 486
 487 0 [IF]
 488
 489 23.8e-9 fconstant fsmall
 490 1.0 fsmall f- fconstant falmost1
 491 ." Should be 1.0 = " falmost1 f. cr
 492
 493 : TSEGF  ( r -f- , print in all formats )
 494 ." --------------------------------" cr
 495     34 0
 496     DO
 497         fdup fs. 4 spaces  fdup fe. 4 spaces
 498         fdup fg. 4 spaces  fdup f.  cr
 499         10.0 f/
 500     LOOP
 501     fdrop
 502 ;
 503
 504 : TFP
 505     1.234e+22 tsegf
 506     1.23456789e+22 tsegf
 507     0.927 fsin 1.234e+22 f* tsegf
 508 ;
 509
 510 [THEN]