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