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

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