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