home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / forth / compiler / fpc / source / float4th.seq < prev    next >
Text File  |  1990-04-27  |  49KB  |  1,579 lines
  1. \\ High Level FLoating Point .  IEEE Format.     14:57 13Nov88RS)
  2.  
  3.           Copyright 1988 by Robert L. Smith
  4.       2300 St. Francis Drive, Palo Alto, CA 94303
  5.                  (415) 856-9321
  6.  
  7.  These routines may be freely used, provided only that the
  8.  copyright notice be displayed and preserved.
  9.  
  10. Please: If you wish to distribute any of your changes to this
  11.  package, please give your name, address and telephone number
  12.  so that any other users can contact you if they find problems.
  13.  
  14.       ( Approximate times on 4.77 MHz 8088 are: )
  15.       ( F+    10 milliseconds )
  16.       ( F/    20 milliseconds )
  17.       ( FLN   30 milliseconds )
  18.  
  19. \ Load Block                                    15:10 23Jul89RS)
  20.  
  21. : COPYRIGHT
  22.      CR ." Floating Forth  Version 1.2   Jan. 25, 1990 " CR
  23.      CR ." written by Robert L. Smith "
  24.      CR ." 2300 St. Francis Drive, Palo Alto, CA 94303 " CR ;
  25.  
  26. {
  27.  
  28. ( Tests for Intel or Motorola processor type.   07:24 31Oct88RS)
  29.  
  30. CREATE TRIBUF  0 , 0 , 0 ,
  31.  
  32. TRIBUF 1+  CONSTANT TRIBUF1
  33. TRIBUF 2+  CONSTANT TRIBUF2
  34. TRIBUF 3 + CONSTANT TRIBUF3
  35. TRIBUF 4 + CONSTANT TRIBUF4
  36.  
  37. }
  38. \ Byte-order dependency load screen.            20:41 23Oct88RS)
  39.  
  40. TRIBUF is an array used for shifting by 8 bits.  Its use
  41.    depends on the byte order of the computer and implementation
  42.    being used.  In order to determine how to perform the shift,
  43.    the byte order is determined by the function BIGENDIAN .
  44.    Depending on the results of that test, one screen or the
  45.    other is then loaded.
  46.  
  47. {
  48. ( Little Endian i.e., INTEL forms.              08:25 15Oct88RS)
  49.  
  50. ( For Little Endian byte order such as Intel uses: )
  51.  
  52. : DSHFT8        ( lo hi -- lowest middle high )
  53.                 TRIBUF3 !  TRIBUF1 !
  54.                 TRIBUF @ TRIBUF2 @ TRIBUF4 @ ;
  55.  
  56. }
  57. \ 8-bit shift functions
  58.  
  59. DSHFT8   Perform an 8 bit shift on a double precision number,
  60.    returning a triple word result (no loss of bits).  The
  61.    results can be considered as either right or left shifted.
  62.  
  63. {
  64. \  MA  MB  MC  XMD  Y*  FRACT*                 07:07 15Oct88RS)
  65.  
  66. 8 ARRAY MA
  67. MA 2 + CONSTANT MB
  68. MA 4 + CONSTANT MC
  69. MA 6 + CONSTANT XMD
  70.  
  71. : Y*   ( d1 d2 -- t )  ( High order 8 bits of d1 and d2 = 0 )
  72.         MA 2!  MC 2!  MB @ XMD @ UM* 0
  73.         MA @ XMD @ UM* D+  MB @ MC @ UM* D+
  74.         MA @ MC @ * + ;
  75.  
  76. : FRACT*   ( n1 n2 -- n3 )
  77.         UM* SWAP DROP ;
  78.  
  79. }
  80. \  Partial and fractional multiplication.       20:51 23Oct88RS)
  81.  
  82. MA  MB  MC  XMD  Variables used for multiplication.
  83.  
  84. Y*   Take the product of two 24-bit numbers, returning a triple
  85.    precision result.  The 24-bit multiplicands are double
  86.    numbers with the most significant 8 bits cleared to 0.
  87.  
  88. FRACT*  Take the product of two single numbers, treated as
  89.    fractions with the binary point at the left.
  90.  
  91. {
  92. \  F#BYTES  FPSTAT  FPSTAT  FPS0  FDEPTH        23:44 15Nov88RS)
  93.  
  94. 4 CONSTANT F#BYTES           ( Size of floating point number )
  95. CREATE FPSTAT 0 , 0 ,          ( Holds error information )
  96.  
  97.  
  98. 20 F#BYTES * CONSTANT FPSSIZE           \ Size of floating pt. stack
  99.  
  100. FPSSIZE 3 F#BYTES * + ARRAY FPSTACK     \ Floating Point Stack
  101.                                         \ with room for underflow
  102.  
  103. : FSP0          ( -- a1 )               \ Point to base of stack.
  104.                 FPSTACK FPSSIZE + ;
  105.  
  106. VARIABLE FSP                            \ Points to top of stack
  107.                                         \ initialized by FCLEAR
  108.  
  109. : FDEPTH   ( -- n )
  110.         FSP0 FSP @ - F#BYTES / ;
  111.  
  112. }
  113. \ Definition of the Floating Point stack.       21:36 23Oct88RS)
  114.  
  115. F#BYTES   Size of a floating point number in bytes.
  116. FPSTAT    A variable to hold the Floating Point Status.
  117.  
  118. FPSSIZE   The size of the floating point stack, in bytes.
  119.  
  120. FPSTACK   The floating point stack.
  121.  
  122. FSP0      The base of the floating point stack.
  123. FSP       Contains the pointer to the F.P. Stack top.
  124.  
  125. FDEPTH    The depth of the floating point stack, in FPSSIZE
  126.           units.
  127.  
  128. {
  129. \  XBIAS  FCLEAR  FDROP  FPUSH                  14:40 08Nov88RS)
  130.  
  131. $3F80 CONSTANT XBIAS             ( Exponent bias )
  132.  
  133. : FPUSH   ( F: -- r ; d -- )
  134.         -4 FSP +!  FSP @ 2! ;
  135.  
  136. : FCLEAR   ( -- )
  137.         FSP0 FSP ! ;
  138.  
  139. FCLEAR          \ initialize the floating point system
  140.  
  141. : FDROP   ( F: r -- )
  142.         F#BYTES FSP +! ;
  143.  
  144. : F2DROP  ( -- )
  145.         FDROP FDROP ;
  146.  
  147. }
  148. \  XBIAS  FCLEAR  FDROP                         14:17 08Nov88RS)
  149.  
  150. XBIAS   The exponent bias.
  151.  
  152. FCLEAR  Empties the floating point stack.
  153.  
  154. FDROP   Drop one floating point element from the F.P. stack.
  155.  
  156. FNSWAP  Exchange the n-th item on the F.P. stack with the
  157.         zeroth item.
  158.  
  159. {
  160. \ FPICK  FDUP  FOVER  FSWAP                     14:32 08Nov88RS)
  161.  
  162. : FPICK   ( F: rn rn-1 ... r0 -- fn rn-1 ... r0 rn ; n -- )
  163.         F#BYTES * FSP @ + 2@ FPUSH ;
  164.  
  165. : FDUP   ( F: r -- r r )
  166.         FSP @ 2@  F#BYTES NEGATE FSP +!
  167.         FSP @ 2! ;
  168.  
  169. : FOVER   ( F: r1 r2 -- r1 r2 r1 )
  170.         F#BYTES NEGATE FSP +!
  171.         FSP @ 8 + 2@ FSP @ 2! ;
  172.  
  173. : F2DUP   ( ? -- ? )
  174.         FOVER FOVER ;
  175.  
  176. : FSWAP   ( F: r1 r2 -- r2 r1 )
  177.         FSP @ 2@  FSP @ 4 + 2@
  178.         FSP @ 2!  FSP @ 4 + 2! ;
  179. }
  180. \ FDUP  FOVER  FSWAP                            14:35 08Nov88RS)
  181.  
  182. FPICK   Push a copy of the nth element of the F.P. stack
  183.         onto the F.P. stack.
  184.  
  185. FDUP    Duplicate the top element on the F.P. stack.
  186.  
  187. FOVER   Push a copy of the second element on the F.P. stack
  188.         onto the F.P. stack.
  189.  
  190. FSWAP   Interchange the top two elements on the F.P. stack.
  191.  
  192. {
  193. \  FNEGATE  FNIP  FROT  F-ROT                   21:40 23Oct88RS)
  194.  
  195. : FNEGATE   ( F: r1 -- r2 )
  196.         FSP @ DUP @ $8000 XOR SWAP ! ;
  197.  
  198. : FNIP   ( F: r1 r2 -- r2 )
  199.         FSP @ 2@  4 FSP +!  FSP @ 2! ;
  200.  
  201. : FROT   ( F: r1 r2 r3 -- r2 r3 r1 )
  202.         FSP @ 2@  FSP @ 4 + 2@  FSP @ 8 + 2@
  203.         FSP @ 2!  FSP @ 8 + 2!  FSP @ 4 + 2! ;
  204.  
  205. : F-ROT   ( F: r1 r2 r3 -- r3 r1 r2 )
  206.         FSP @ 2@  FSP @ 4 + 2@  FSP @ 8 + 2@
  207.         FSP @ 4 + 2!  FSP @ 2!  FSP @ 8 + 2! ;
  208.  
  209. }
  210. \  FNEGATE  FNIP  FROT  F-ROT                   21:44 23Oct88RS)
  211.  
  212. FNEGATE   Reverse the sign of the element at the top of the
  213.           F.P. stack.
  214.  
  215. FNIP      Drop the second item from the F.P. stack.
  216.  
  217. FROT      Rotate the top three items on the F.P. stack,
  218.           bring the third element to the top.
  219.  
  220. F-ROT     Rotate the top three items on the F.P. stack,
  221.           bringing the second item to the top, and moving the
  222.           former top item to the third position.
  223.  
  224. {
  225. \ FPOP  FPCOPY  F0=  FNSWAP              14:39 014:45 08Nov88RS)
  226.  
  227. : FPOP   ( F: r -- ; -- d )
  228.         FSP @ 2@ FDROP ;
  229.  
  230. : FPCOPY   ( F: r -- r ; -- d )
  231.         FSP @ 2@ ;
  232.  
  233. : F0=   ( F: r -- ; -- flag )
  234.         FPOP $7FFF AND OR 0= ;
  235.  
  236. : FPOP0=   ( F: r -- ; -- d flag )
  237.         FPOP 2DUP $7FFF AND OR 0= ;
  238.  
  239. }
  240. \ FPOP  FPUSH  FPCOPY  F0=                      23:02 23Oct88RS)
  241.  
  242. FPOP    Pop the top number from the F.P. stack, and push it
  243.         onto the parameter stack as a double number.
  244.  
  245. FPUSH   Pop the top double number and
  246.  
  247. FPCOPY  Get a copy of the top number on th F.P. stack and push
  248.         it on the parameter stack.
  249.  
  250. F0=     Pop the top member of the F.P. stack and test its value.
  251.         If the value is zero, push true, else push false onto
  252.         the parameter stack.
  253.  
  254. {
  255. \  F0<  F=                                      14:45 08Nov88RS)
  256.  
  257. : F0<   ( F: r -- ; -- flag )
  258.         FPOP DUP 0<
  259.         IF      $7FFF AND OR 0= 0=
  260.         ELSE    2DROP 0   THEN ;
  261.  
  262. : F=   ( F: r1 r2 -- ; -- flag )
  263.         FPOP0=
  264.         IF      2DROP FPOP0= SWAP DROP SWAP DROP
  265.         ELSE    FPOP D=   THEN ;
  266.  
  267. : FNSWAP   ( F: rn rn-1 ... r0 -- r0 rn-1 ... rn ; n -- )
  268.         F#BYTES * FSP @ + DUP >R 2@  FSP @ 2@
  269.         R> 2!  FSP @ 2! ;
  270.  
  271. }
  272. \  F0<  FPOP0=  F=                              23:13 23Oct88RS)
  273.  
  274. F0<     Pop and test the number at the top of the F.P. stack.
  275.         If the sign is negative and the value is non-zero, push
  276.         a true flag.  Otherwise push a false flag.
  277.  
  278. FPOP0=  Pop the top member of the F.P. stack, test, and push it
  279.         onto the parameter stack.  If the number has a value of
  280.         zero, also push a true flag; otherwise push a false
  281.         flag.
  282. F=      Pop the top two elements from the F.P. stack.  If the
  283.         two numbers are equal, push a true flag.  Otherwise,
  284.         push a false flag.
  285.  
  286. {
  287. \  F<                                           06:28 18Jul89RS)
  288.  
  289. : F<    ( F: r1 r2 -- ; -- flag )
  290.         FPOP FPOP DUP 0<
  291.         IF      DU<
  292.         ELSE    2SWAP D<
  293.         THEN ;
  294.  
  295. }
  296. \  F<                                           23:16 23Oct88RS)
  297.  
  298. F<      Pop the top two F.P. numbers from the F.P. stack and
  299.         compare them.  If the second number is arithmetically
  300.         less than the first number, push a true flag.
  301.         Otherwise push a false flag.
  302.  
  303. {
  304. \ F>                                            06:44 18Jul89RS)
  305.  
  306. : F>    ( F: r1 r2 -- ; -- flag )
  307.         FPOP FPOP DUP 0<
  308.         IF      2SWAP DU<
  309.         ELSE    D<
  310.         THEN ;
  311.  
  312. }
  313. \ F>                                            23:18 23Oct88RS)
  314.  
  315. F>      Pop the top two numbers from the F.P. stack and
  316.         compare them.  If the second number is arithmetically
  317.         greater than the top, push a true flag onto the
  318.         parameter stack.  Otherwise push a false flag.
  319.  
  320. {
  321. \  FABS  FMIN  FMAX                             06:31 18Jul89RS)
  322.  
  323.  
  324. : FABS  ( F: r1 -- r2 )
  325.         FSP @ DUP @ $7FFF AND SWAP ! ;
  326.  
  327. : FMIN   ( F: r1 r2 -- r3 )
  328.         FPOP FPOP 2OVER 2OVER DUP 0<
  329.         IF   DU<
  330.         ELSE 2SWAP D<
  331.         THEN
  332.         IF   2SWAP  THEN
  333.         2DROP FPUSH ;
  334.  
  335. : FMAX   ( F: r1 r2 -- r3 )
  336.         FPOP FPOP 2OVER 2OVER DUP 0<
  337.         IF      2SWAP DU<
  338.         ELSE    D<
  339.         THEN
  340.         IF  2SWAP  THEN
  341.         2DROP FPUSH ;
  342.  
  343. }
  344. \  FABS  FMIN  FMAX                             23:26 23Oct88RS)
  345.  
  346. FABS    Set the sign of the top of the F.P. stack to 0.
  347.  
  348. FMIN    Pop the top two members from the F.P. stack.  Push the
  349.         arithmetically smaller back onto the F.P. stack.
  350.  
  351. FMAX    Pop the top two members from the F.P. stack.  Push the
  352.         arithmetically larger back onto the F.P. stack.
  353.  
  354. {
  355. \ F@  F!  FCONSTANT  FVARIABLE                  08:03 20Oct88RS)
  356.  
  357. : F@           ( F: -- r ; addr -- )
  358.                2@ FPUSH ;
  359.  
  360. : F!           ( F: r -- ; addr -- )
  361.                >R FPOP R> 2! ;
  362.  
  363. : fconstant    ( F: r -- )  ( F: -- r )
  364.                CREATE fpop , ,
  365.                DOES>  F@ ;
  366.  
  367. : FVARIABLE    ( -- )  ( -- addr )
  368.                CREATE 4 ALLOT ;
  369.  
  370. }
  371. \ F@  F!  FCONSTANT  FVARIABLE                  00:03 23Oct88RS)
  372.  
  373. F@      Fetch the F.P. variable at the address specified by the
  374.         top of the parameter stack.  Push the contents of the
  375.         variable onto the F.P. stack.  Pop and discard the
  376.         address at the top of the parameter stack.
  377. F!      Store the number at the top of the F.P. stack into
  378.         memory at the address at the top of the parameter stack.
  379.         Pop the number from the F.P. stack, and pop the address
  380.         from the parameter stack.
  381. FCONSTANT  Create a F.P. constant with a value equal to the
  382.         number poped off the F.P. stack.
  383. FVARIABLE  Create a F.P. variable.
  384.  
  385. {
  386. \ Various FCONSTANTs: F1.0 PI F0.0 FLOG10E      07:20 15Oct88RS)
  387.  
  388. $0000 $3F80 fpush fconstant F1.0
  389. $0FDB $4049 fpush fconstant PI
  390. $0000 $0000 fpush fconstant F0.0
  391. $5BD9 $3EDE fpush fconstant FLOG10E
  392. $5D8E $4013 fpush fconstant FLN10.0
  393. $0000 $4120 fpush fconstant F10.0
  394. $0000 $3F00 fpush fconstant F0.5
  395.  
  396. }
  397. \ Various FCONSTANTs: F1.0 PI F0.0 FLOG10E      00:07 23Oct88RS)
  398.  
  399. F1.0      Floating point 1.
  400. PI        Floating point pi
  401. F0.0      Floating point 0.
  402. FLOG10E   Floating point log base 10 of e
  403. FLN10.0   Floating point natural log of 10
  404. F10.0     Floating point 10.
  405. F0.5      Floating point 0.5
  406.  
  407. {
  408. \ T2/  T2*  T>SHIFT                             07:21 15Oct88RS)
  409.  
  410. : T2/   ( t1 -- t2 )
  411.         >R D2/ R@ 1 AND
  412.         IF    $8000 OR   ELSE   $7FFF AND   THEN    R> 2/ ;
  413.  
  414. : T2*   ( t1 -- t2 )
  415.         >R DUP 0<
  416.         IF   D2* R> 2* 1 OR   ELSE   D2* R> 2*  THEN ;
  417.  
  418. : T>SHIFT   ( t1 n -- t2 )   ( n is multiple of 80 hex )
  419.         ?DUP
  420.         IF   0  DO  T2/  $80 +LOOP  THEN ;
  421.  
  422. }
  423. \ T2/  T2*  T>SHIFT                             00:10 23Oct88RS)
  424.  
  425. T2/     Shift right a triple precision number.
  426.  
  427. T2*     Shift left a triple precision number.
  428.  
  429. T>SHIFT  Shift right a triple precision number by a count equal
  430.          to the argument divided by 128 (80 hex).
  431.  
  432. {
  433. \ D>SHIFT  D<SHIFT  D-CY                        07:21 15Oct88RS)
  434.  
  435. : D>SHIFT   ( d1 n -- d2 )   ( n a multiple of 80 hex )
  436.         ?DUP  IF  0  DO D2/  $80 +LOOP  THEN ;
  437.  
  438. : D<SHIFT   ( d1 n -- d2 )  ( n is multiple of 80 hex )
  439.         ?DUP
  440.         IF      0 DO  D2*  $80 +LOOP
  441.         THEN ;
  442.  
  443. : D-CY   ( d1 d2 -- d3 n )   ( n is borrow or carry )
  444.         2OVER 2OVER DU<
  445.         IF      D- -1
  446.         ELSE    D- 0
  447.         THEN ;
  448.  
  449. }
  450. \ D>SHIFT  D<SHIFT  D-CY                        00:14 23Oct88RS)
  451.  
  452. D>SHIFT  Shift a double precision number right by a count equal
  453.          to the top parameter divided by 128 (80 hex).
  454.  
  455. D<SHIFT  Shift a double precision number left by a count equal
  456.          to the top parameter divided by 128 (80 hex).
  457.  
  458. D-CY     Subtract the top double number from the second.
  459.          Return the difference and a borrow flag.
  460.  
  461. {
  462. \  D+CY  ZSIGN  ZEXP  FIXGRS                    07:21 15Oct88RS)
  463.  
  464. : D+CY   ( d1 d2 -- d3 n )
  465.         DNEGATE D-CY ;
  466.  
  467. CREATE ZSIGN   0 ,
  468. CREATE ZEXP    0 ,
  469.  
  470. : FIXGRS   ( n1 -- n2 )
  471.         DUP $3FFF AND
  472.         IF      $C000 AND $0F OR
  473.         THEN ;
  474.  
  475. }
  476. \  D+CY  ZSIGN  ZEXP  FIXGRS                    00:18 23Oct88RS)
  477.  
  478. D+CY    Add the top two double numbers and return a carry flag.
  479.  
  480. ZSIGN   A variable to carry the sign of a result.
  481. ZEXP    A vaiable to carry the resultant exponent.
  482.  
  483. FIXGRS  Set "sticky" bit flags.
  484.  
  485. {
  486. \  UNNORMALIZE                                  07:21 15Oct88RS)
  487.  
  488. : UNNORMALIZE   ( d1 n -- d2 grs )  ( n is multiple of 80X )
  489.         >R R@ $480 <
  490.         IF      0 -ROT R> T>SHIFT ROT FIXGRS EXIT   THEN
  491.         R@ $880 <
  492.         IF      $800 R> - D<SHIFT SWAP 0 SWAP FIXGRS EXIT  THEN
  493.         R@ $D80 <
  494.         IF      SWAP DUP $3FF AND
  495.                 IF      $4000 OR   THEN  SWAP
  496.                 R> $800 - D>SHIFT 0 ROT FIXGRS EXIT
  497.         THEN
  498.         R> DROP OR
  499.         IF      0 0 $0F  ELSE   0 0 0  THEN ;
  500.  
  501. }
  502. \  UNNORMALIZE                                  00:21 23Oct88RS)
  503.  
  504. UNNORMALIZE   Unnormalize the double number by a count specified
  505.               at the top of the stack.  The number at the top of
  506.               the stack is a signed count multiplied by 128
  507.               (128 hex).
  508.  
  509. {
  510. \  >NORMALIZE                                   07:22 15Oct88RS)
  511.  
  512. : >NORMALIZE   ( d1 -- grs d2 n )
  513.         0 -ROT 0 $480 0
  514.         DO      DROP DUP $100 U<
  515.                 IF      I LEAVE
  516.                 ELSE    T2/ I
  517.                 THEN
  518.        $80 +LOOP ;
  519.  
  520. }
  521. \  >NORMALIZE                                   00:26 23Oct88RS)
  522.  
  523. >NORMALIZE   Normalize the double number at the top of the
  524.              parameter stack.  Return GRS (Guard, Round, and
  525.              Sticky) bits and a normalizing count multiplied
  526.              by 128.
  527.  
  528. {
  529. \ MROUND  EVROUND  DENORMALIZE1                 08:23 18Jul89RS)
  530.  
  531. : MROUND   ( d1 evenflg -- d2 )
  532.         IF      1 0 D+ SWAP $FFFE AND SWAP
  533.         ELSE    1 0 D+
  534.         THEN ;
  535.  
  536. : EVROUND   ( d1 grs -- d2 )
  537.         DUP 0<
  538.         IF      $8000 =   MROUND
  539.         ELSE    DROP
  540.         THEN ;
  541.  
  542. : DENORMALIZE1   ( d1 n -- d2 )
  543.                 NEGATE $80 + UNNORMALIZE EVROUND ;
  544.  
  545. }
  546. \ MROUND  EVROUND  DENORMALIZE1                 07:57 31Oct88RS)
  547.  
  548. MROUND  Round the double number on the stack.  If the flag at
  549.         the top of the stack is true, then round to even.
  550.  
  551. EVROUND  If the guard bit is set, round the double number on the
  552.          stack.
  553.  
  554. DENORMALIZE1  Denormalize the double number according to the
  555.         shifted count at the top of the stack.  Round toward
  556.         even.
  557. {
  558. \ DENORMALIZE2  1NORMALIZE                      07:56 31Oct88RS)
  559.  
  560. : DENORMALIZE2   ( grs d1 n -- d2 )
  561.         NEGATE $80 + UNNORMALIZE >R ROT R> OR EVROUND ;
  562.  
  563. : 1NORMALIZE    ( d1 -- d2 n )
  564.         $8000 $C00 0
  565.         DO      DROP DUP $7F >
  566.                 IF      $7F AND I NEGATE LEAVE
  567.                 ELSE    D2* $8000
  568.                 THEN
  569.        $80 +LOOP ;
  570.  
  571. }
  572.                                                 08:04 31Oct88RS)
  573.  
  574. DENORMALIZE2  Denormalize the combination of d1 and GRS
  575.               according to the shifted count n .
  576.  
  577. 1NORMALIZE    Shift the double number left until it is in a
  578.               normalized form.  A shifted form of the count is
  579.               left on the stack.  The shifted form is used to
  580.               speed up the conversion process.
  581.  
  582. {
  583. \ NORMALIZE  4NORMALIZE                         08:16 31Oct88RS)
  584.  
  585. : NORMALIZE   ( d1 -- d2 n )
  586.         2DUP D0=
  587.         IF      0 EXIT  THEN
  588.         1NORMALIZE ;
  589.  
  590. : 4NORMALIZE   ( d1 -- d2 )
  591.         2DUP OR 0=  IF   EXIT   THEN
  592.         $1000 0
  593.         DO      D2* DUP $7F >
  594.                 IF      $7F AND $3F00 I - OR LEAVE
  595.                 THEN
  596.        $80 +LOOP ;
  597.  
  598. }
  599.                                                 08:26 31Oct88RS)
  600.  
  601. NORMALIZE  Perform a normalization process on the double number.
  602.            Push the shifted count on the stack.
  603.  
  604. 4NORMALIZE  Normalize a double number by shifting left.  This is
  605.             used only by 2NORMALIZE (and ultimately by FLN).
  606.             The result is really a floating point number on the
  607.             parameter stack.
  608.  
  609. {
  610. \ 3NORMALIZE                                    08:18 31Oct88RS)
  611.  
  612. : 3NORMALIZE   ( d1 -- d2 )
  613.         $4480 $4000
  614.         DO      DUP $200 <
  615.                 IF      1 0 D+ D2/ DUP $0FF >
  616.                         IF      D2/
  617.                         ELSE    $7F AND
  618.                         THEN
  619.                         I + LEAVE
  620.                 ELSE    D2/
  621.                 THEN
  622.        $80 +LOOP ;
  623.  
  624. }
  625.                                                 08:25 31Oct88RS)
  626.  
  627. 3NORMALIZE  Normalize a double number by shifting right.  This
  628.             routine is used only by 2NORMALIZE (and ultimately
  629.             by FLN).  The result is really a floating point
  630.             number on the parameter stack.
  631.  
  632. {
  633. \ 2NORMALIZE                                    08:20 31Oct88RS)
  634.  
  635. : 2NORMALIZE   ( d1 -- d2 )
  636.         $7FFF AND DUP $0FF >
  637.         IF      3NORMALIZE
  638.         ELSE    DUP $80 <
  639.                 IF      4NORMALIZE
  640.                 ELSE
  641.                         $3F80 OR
  642.                 THEN
  643.         THEN ;
  644.  
  645. }
  646.                                                 08:27 31Oct88RS)
  647.  
  648. 2NORMALIZE  Normalize the double number, performing either a
  649.             left or right shift, as required.
  650.  
  651. {
  652. \ FLOAT                                         07:02 18Jul89RS)
  653.  
  654. : FLOAT         ( F: -- r ; d -- )
  655.                 2DUP OR 0= IF  FPUSH EXIT THEN
  656.                 DUP 0< $8000 AND >R DABS DUP $80 U<
  657.                 IF      NORMALIZE >R $7F AND $4B00 R> + OR
  658.                         R> OR FPUSH EXIT
  659.                 THEN    >NORMALIZE $4B00 + >R ROT DUP 0<
  660.                 IF      $8000 =  MROUND
  661.                         DUP $100 U< 0=
  662.                         IF      D2/ $7F AND R> $80 + OR
  663.                                 R> OR FPUSH EXIT
  664.                         THEN
  665.                 ELSE    DROP
  666.                 THEN    $7F AND R> OR
  667.                 R> OR FPUSH ;
  668.  
  669. }
  670.                                                 08:28 31Oct88RS)
  671.  
  672. FLOAT   Convert the double number on the parameter stack to a
  673.         floating point number on the floating point stack.
  674.  
  675. {
  676. \ DINTABS                                       05:49 01Nov88RS)
  677.  
  678. : DINTABS   ( F: r -- ; -- d flag )
  679.         FPOP DUP $7F80 AND DUP $3F80 <
  680.         IF      DROP 2DROP 0 0 0 EXIT  THEN
  681.         SWAP $7F AND $80 OR SWAP $4B00 - DUP 0<
  682.         IF      0 SWAP DO   D2/  $80 +LOOP  0
  683.         ELSE    $0400 MIN DUP
  684.                 IF      0 DO   D2*   $80 +LOOP
  685.                 ELSE    DROP
  686.                 THEN    DUP 0<
  687.         THEN  ;
  688.  
  689. }
  690.                                                 05:45 01Nov88RS)
  691.  
  692. DINTABS  Pop the top number from the floating point stack.
  693.          Take the absolute value and convert it to a double
  694.          number on the parameter stack.  If the resulting
  695.          number is positive, push a 0 onto the stack.
  696.          Otherwise, push a true flag (-1) onto the stack.
  697.  
  698. {
  699. \ INT   BMASK                                   07:23 15Oct88RS)
  700.  
  701. : INT   ( F: r -- ; -- d )
  702.         FDUP F0<
  703.         IF      DINTABS >R DNEGATE R>
  704.         ELSE    DINTABS
  705.         THEN    IF   ."  Out of range in INT " FCLEAR ABORT  THEN ;
  706.  
  707. CREATE BMASK    $0000 , $8000 , $C000 , $E000 , $F000 , $F800 ,
  708.                 $FC00 , $FE00 , $FF00 , $FF80 , $FFC0 , $FFE0 ,
  709.                 $FFF0 , $FFF8 , $FFFC , $FFFE , $FFFF ,
  710.  
  711.  
  712. }
  713.                                                 05:54 01Nov88RS)
  714.  
  715. INT     Pop a floating point number from the f.p. stack and
  716.         convert it to a double number.
  717.  
  718. BMASK   An array of bit masks used to obtain the integer part
  719.         of a floating point number.
  720.  
  721. {
  722. \ FINT                                          07:24 15Oct88RS)
  723.  
  724. : FINT   ( F: r1 -- r2 )
  725.         FPOP DUP $7F80 AND DUP $3F80 <
  726.         IF      DROP $8000 AND NIP 0 SWAP FPUSH EXIT  THEN
  727.         DUP $4B00 <
  728.         IF      $3B00 - 2* $100 / 2* DUP $20 <
  729.                 IF      BMASK + @ AND
  730.                 ELSE    $1F AND BMASK + @ ROT AND SWAP
  731.                 THEN
  732.         ELSE    DROP
  733.         THEN
  734.         FPUSH ;
  735.  
  736. }
  737.                                                 05:54 01Nov88RS)
  738.  
  739. FINT    Convert the number at the top of the f.p. stack to its
  740.         integer part represented as a floating point number.
  741.  
  742. {
  743. \ ROUND1                                        07:24 15Oct88RS)
  744.  
  745. : ROUND1   ( grs d1 n1 -- d2 n2 )
  746.         >R ROT DUP 0<
  747.         IF      $8000 = MROUND
  748.                 DUP $0FF >
  749.                 IF      D2/ R> $80 +
  750.                 ELSE    R>
  751.                 THEN
  752.         ELSE    DROP R>
  753.         THEN ;
  754.  
  755. }
  756.                                                 05:58 01Nov88RS)
  757.  
  758. ROUND1  The input parameters represent a floating point number
  759.         broken into an exponent at the top, a double number,
  760.         and the GRS (guard, round, sticky) bits.  Round the
  761.         number according to GRS.
  762.  
  763. {
  764. \ ROUND2                                        07:24 15Oct88RS)
  765.  
  766. : ROUND2   ( d1 grs -- d2 n )
  767.         $8000 =
  768.         IF      SWAP $FFFE AND SWAP
  769.         THEN
  770.         DUP $0FF >
  771.         IF      D2/ $80
  772.         ELSE    0
  773.         THEN ;
  774.  
  775. }
  776.                                                 05:58 01Nov88RS)
  777.  
  778. ROUND2  Round the double number according to the GRS bits at
  779.         the top of the stack.
  780.  
  781. {
  782. \ Aux for F+                                    10:11 18Jul89RLS
  783.  
  784. : (F-X1=X2)   ( d1 d2 -- d3 )  ( F- : Equal exponents )
  785.         D- DUP 0<
  786.         IF      ( mantissa1 < mantissa2 )
  787.                 DNEGATE ZSIGN @ $8000 XOR ZSIGN !
  788.         ELSE    2DUP D0=  IF   EXIT   THEN
  789.         THEN   ZEXP @ $80 >
  790.         IF      ( Equal exponents, normal r1' )
  791.                 NORMALIZE ZEXP @ + DUP $80 <
  792.                 IF      ( denormalize )
  793.                         DENORMALIZE1
  794.                 ELSE    SWAP $7F  AND OR
  795.                 THEN
  796.          THEN ( ZSIGN @ OR ) ;
  797.  
  798. }
  799.                                                 06:02 01Nov88RS)
  800.  
  801. (F-X1=X2)  Auxilliary floating point subtraction of magnitudes
  802.            function for the case of equal exponents for the two
  803.            operands.
  804.  
  805. {
  806. \ Auxilliary for F+                             10:10 18Jul89RLS
  807.  
  808. : (F+X1=X2)   ( d1 d2 -- d3 )
  809.         D+ DUP $0FF >
  810.         IF      ( normal case )  OVER 1 AND
  811.                 IF      ( Round to even case )
  812.                         1 0 D+ D2/ SWAP $FFFE AND SWAP
  813.                 ELSE    D2/
  814.                 THEN   ZEXP @ + DUP 0<
  815.                 IF   ."  Overflow in F+ " FCLEAR ABORT  THEN
  816.         ELSE    DUP $7F >
  817.                 IF      $7F AND ZEXP @ DUP 0=
  818.                         IF      DROP $80   THEN  OR
  819.                 THEN
  820.         THEN  ( ZSIGN @ OR ) ;
  821.  
  822. }
  823.                                                 06:01 01Nov88RS)
  824.  
  825. (F+X1=X2)  Auxillary floating point addition of magnitudes
  826.            function for the case of equal exponents of the
  827.            operands.
  828.  
  829. {
  830. \ Auxilliary for F+                             07:25 15Oct88RS)
  831.  
  832. : (1F-)   ( d1 d2 x1-x2 -- grs d3 n )
  833.         UNNORMALIZE NEGATE DUP >R
  834.         IF    1 0 D+  THEN   D- R> -ROT        ( grs d3 )
  835.         DUP $7F >  IF   0 EXIT  THEN
  836.         T2* DUP $7F >
  837.         IF      $FF80 EXIT  THEN
  838.         T2* 0 $C800 $0100
  839.         DO      DROP DUP $7F >
  840.                 IF      I NEGATE LEAVE THEN
  841.                 D2* I NEGATE $80
  842.         +LOOP ;
  843.  
  844. }
  845.                                                 06:03 01Nov88RS)
  846.  
  847. (1F-)   An auxillary function for floating point subtraction of
  848.         magnitudes.
  849.  
  850. {
  851. \ Aux for F+                                    07:25 15Oct88RS)
  852.  
  853. : (F+-AUX)   ( d2 sx2 -- d3 sx2 x1-x2 flag )
  854.         DUP $7F80 AND DUP 0=
  855.         IF      DROP >R $7F AND R> $80 OR $80   THEN
  856.         ZEXP @ SWAP - DUP 0< ;
  857.  
  858. }
  859.                                                 06:05 01Nov88RS)
  860.  
  861. (F+-AUX)  An auxillary function for floating point addition.
  862.  
  863. {
  864. \ Aux for F+                                    07:26 15Oct88RS)
  865.  
  866. : (F-)   ( d1 d2 sx2 -- d3 )        ( Subtract magnitudes )
  867.         (F+-AUX)
  868.         IF      ( x2 > x1 )
  869.                 NEGATE SWAP DUP $7F80 AND ZEXP !
  870.                 $8000 AND ZSIGN !  >R 2SWAP R>
  871.         ELSE    ( x2 <= x1 )  NIP DUP 0=
  872.                 IF      ( Exponents are equal )
  873.                         DROP (F-X1=X2) EXIT
  874.                 THEN
  875.         THEN    ( d1 d2 x1-x2 )
  876.         (1F-) ZEXP @ + DUP $80 <
  877.         IF   DENORMALIZE2  ELSE   ROUND1 SWAP $7F AND OR  THEN ;
  878.  
  879. }
  880.                                                 06:06 01Nov88RS)
  881.  
  882. (F-)    Auxilliary function for the subtraction of magnitudes
  883.         of floating point numbers.
  884.  
  885. {
  886. \ Auxilliary for F+                             08:46 20Oct88RS)
  887.  
  888. : (F+)   ( d1 d2 sx2 -- d3 )
  889.         (F+-AUX)
  890.         IF      ( x2 > x1 )  NEGATE SWAP DUP $7F80 AND ZEXP !
  891.                 $8000 AND ZSIGN !  >R 2SWAP R>
  892.         ELSE    ( x2 <= x1 )   NIP DUP 0=
  893.                 IF      ( x1 = x2 ) DROP (F+X1=X2) EXIT  THEN
  894.         THEN    ( d1 d2 x1-x2 ) UNNORMALIZE
  895.         >R D+ DUP $0FF >   IF  1 0 D+ R> 0=
  896.                 IF      SWAP $FFFC AND SWAP  THEN  D2/ $80
  897.         ELSE    R@ 0<
  898.                 IF   1 0 D+ R> ROUND2  ELSE   R> DROP 0  THEN
  899.         THEN  ZEXP @ + DUP 0<
  900.         IF      ."  Overflow in F+ " FCLEAR ABORT  THEN   DUP
  901.         IF      SWAP $7F AND SWAP  THEN  OR ;
  902. }
  903. (F+)                                            06:11 01Nov88RS)
  904.  
  905. (F+)    Auxillary function for addition of floating point
  906.         numbers having the same sign.
  907.  
  908. {
  909. \  F+                                           07:26 15Oct88RS)
  910.  
  911. : F+   ( F: r1 r2 -- r3 )
  912.         FSP @ 4 + 2@ DUP $8000 AND ZSIGN !
  913.         DUP $7F80 AND DUP ZEXP !
  914.         IF      $7F AND $80 OR
  915.         ELSE    $7F AND 2DUP OR 0=
  916.                 IF   2DROP FNIP EXIT   THEN
  917.         THEN
  918.         FPOP0=  IF    2DROP 2DROP EXIT   THEN
  919.         FDROP DUP $7F AND $80 OR SWAP $FF80 AND DUP
  920.         ZSIGN @ XOR 0<
  921.         IF      (F-)  ELSE    (F+)   THEN
  922.         ZSIGN @ XOR FPUSH ;
  923.  
  924. }
  925. F+                                              06:11 01Nov88RS)
  926.  
  927. F+      Floating point addition.
  928.  
  929. {
  930. \ F-  FIX  ZDEN  ZQUOT                          08:46 20Oct88RS)
  931.  
  932. : F-   ( F: r1 r2 -- r3 )
  933.         FNEGATE F+ ;
  934.  
  935. : FIX   ( F: r -- ; -- d )
  936.         FDUP F0<
  937.         IF      FABS F0.5 F+ DINTABS >R DNEGATE R>
  938.         ELSE    F0.5 F+ DINTABS
  939.         THEN
  940.         IF   ."  Out of range in FIX " FCLEAR ABORT  THEN  ;
  941.  
  942. 2VARIABLE ZDEN   2VARIABLE ZQUOT
  943.  
  944. }
  945. F-  FIX  ZDEN  ZQUOT                            06:10 01Nov88RS)
  946.  
  947. F-      Floating point subtraction.
  948.  
  949. FIX     Pop a number from the floating point stack, convert it
  950.         to a double number and push the result on the parameter
  951.         stack.  Issue an error message if the number cannot be
  952.         properly converted.
  953.  
  954. ZDEN    A variable for temporary results in f.p. division.
  955. ZQUOT   A variable for temporary results in f.p. division.
  956.  
  957. {
  958. \ Aux for F/                                    07:27 15Oct88RS)
  959.  
  960. : (1F/)   ( 0 d1 d2 -- t )
  961.         ( Set quotient exponent, possible num adjust )
  962.         FDROP DUP $7F AND $80 OR SWAP $7F80 AND DUP
  963.         IF      XBIAS - NEGATE
  964.         ELSE    DROP NORMALIZE
  965.         THEN
  966.         2/ ZEXP +!
  967.         2OVER 2OVER DU< 0=
  968.         IF      ZDEN 2!  T2/
  969.         ELSE    ZDEN 2!  $-40 ZEXP +!
  970.         THEN
  971.         ZDEN 2@ DSHFT8 DROP ZDEN 2! ;  ( Shift den left by 8 )
  972.  
  973. }
  974. \ Aux for F/                                    06:13 01Nov88RS)
  975.  
  976. (1F/)   An auxillary function for floating point division.
  977.         Set the quotient exponent, possibly adjust the
  978.         numerator.
  979.  
  980. {
  981. \ Aux for F/                                    06:15 01Nov88RS)
  982.  
  983. : (2F/)   ( t -- d )
  984.         ZDEN @ UM/MOD DUP ZQUOT !
  985.         ZDEN 2+ @ UM* D-CY
  986.         IF      ZDEN 2@ D+CY
  987.                 IF      ZDEN 2@ D+ -2
  988.                 ELSE    -1
  989.                 THEN
  990.                 ZQUOT +!
  991.         THEN ;
  992.  
  993. }
  994. \ Aux for F/                                    06:14 01Nov88RS)
  995.  
  996. (2F/)   An auxilliary function for floating point division.
  997.         In this routine we obtain the most significant part
  998.         of the quotient.
  999.  
  1000. {
  1001. \ Aux for F/                                    15:24 14Oct88RLS
  1002.  
  1003. : (3F/)   ( 0 d1 n -- d2 )  ( Usual generate low quotient )
  1004.         UM/MOD DUP ZQUOT 2+ !
  1005.         ZDEN 2+ @ UM* D-CY
  1006.         IF      ZDEN 2@ D+CY
  1007.                 IF      ZDEN 2@ D+ -2
  1008.                 ELSE    -1
  1009.                 THEN
  1010.                 ZQUOT 2+ +!
  1011.         THEN ;
  1012.  
  1013. }
  1014.                                                 06:16 01Nov88RS)
  1015.  
  1016. (3F/)   Auxillary function for floating point division.  This
  1017.         routine is normally called to obtain the low order part
  1018.         of the quotient.
  1019.  
  1020. {
  1021. \ Aux for F/                                    15:23 14Oct88RLS
  1022.  
  1023. : (4F/)   ( 0 d1 n -- d2 )
  1024.         ( Gen low quotient for hi num = hi den )
  1025.         DROP NIP ZDEN 2+ @ SWAP 0 ZDEN @ 0 D+
  1026.         ZDEN 2+ @ 0 D- NIP 0<           \ TJZ 01/25/90 17:36:44.78 ADDED NIP
  1027.         IF      ZDEN 2@ D+CY
  1028.                 IF      ZDEN 2@ D+ -3
  1029.                 ELSE    -2
  1030.                 THEN
  1031.         ELSE    -1
  1032.         THEN
  1033.         ZQUOT 2+ ! ;
  1034.  
  1035. }
  1036. \ Aux for F/                                    06:18 01Nov88RS)
  1037.  
  1038. (4F/)   Auxilliary function for floating point division.
  1039.         Generate the low part of the quotient for the case of
  1040.         the high part of the numerator equal to the high part
  1041.         of the denominator.
  1042.  
  1043. {
  1044. \ SROUND                                        07:27 15Oct88RS)
  1045.  
  1046. : SROUND        ( d1 -- d2 )
  1047.                 DUP 0<
  1048.                 IF      ( Round up )  2DROP ZQUOT 2@ 1 0 D+ EXIT
  1049.                 THEN    D2* 2DUP ZDEN 2@ DU< 0=
  1050.                 IF      ( Do we round to even? )
  1051.                         ZDEN 2@ D=
  1052.                         IF      ( Yes )
  1053.                                 ZQUOT 2@ 1 0 D+
  1054.                                 SWAP $FFFE AND SWAP
  1055.                         ELSE    ZQUOT 2@ 1 0 D+
  1056.                         THEN
  1057.                 ELSE    2DROP ZQUOT 2@
  1058.                 THEN    ;
  1059.  
  1060. }
  1061.                                                 06:20 01Nov88RS)
  1062.  
  1063. SROUND  A rounding function used in floating point division.
  1064.  
  1065. {
  1066. \ Auxilliary for F/                             07:27 15Oct88RS)
  1067.  
  1068. : (UF/)   ( d1 -- d2 )   ( Generate unnormalized quotient )
  1069.         ZEXP @ 2* NEGATE $80 + UNNORMALIZE DUP 0<
  1070.         IF      $8000 = MROUND
  1071.         ELSE    DROP
  1072.         THEN ;
  1073.  
  1074. VARIABLE ZTEMP
  1075.  
  1076. }
  1077.                                                 06:22 01Nov88RS)
  1078.  
  1079. (UF/)   An auxilliary function used in floating point division
  1080.         to generate an unnormalized quotient.
  1081.  
  1082. ZTEMP   Another variable for temporary storage.
  1083.  
  1084. {
  1085. \  F/                                           08:47 20Oct88RS)
  1086.  
  1087. : F/   ( F: r1 r2 -- r3 )
  1088.       0 0 FSP @ 4 + 2@ DUP FSP @ @ XOR $8000 AND ZSIGN !
  1089.       DUP $7F80 AND DUP 2/ ZEXP !
  1090.       IF      $7F AND $80 OR
  1091.       ELSE    $7F AND 2DUP OR 0=
  1092.               IF  2DROP FDROP FDROP 0 ZSIGN @ FPUSH EXIT
  1093.               THEN  NORMALIZE $80 + 2/ ZEXP !
  1094.       THEN  FPOP0=
  1095.       IF  ."  Floating Division by 0" FCLEAR ABORT  THEN  (1F/) (2F/)
  1096.       ZDEN @ 2DUP U<  IF   (3F/)  ELSE  (4F/)  THEN
  1097.       ZEXP @ $40 <    IF  (UF/)  ELSE  SROUND ZEXP @ $3FC0 >
  1098.               IF  ."  Overflow in F/ "  FCLEAR ABORT  THEN
  1099.               $7F AND ZEXP @ 2* OR
  1100.       THEN  ZSIGN @ OR FPUSH ;
  1101. }
  1102.                                                 06:31 01Nov88RS)
  1103.  
  1104. F/      The floating point division routine.  The rather sneaky
  1105.         technique used here is attributable to Roedy Green, the
  1106.         author of BBL/Abundance.  To do a division by a double
  1107.         number, shift the denominator until the m.s. bit is
  1108.         set.  Use the high order part to obtain the first
  1109.         approximation, along with its remainder.  If neccessary,
  1110.         make one or two stages of correction to obtain the
  1111.         exact high order part and temporary remainder.  Repeat
  1112.         the process to obtain the low order part.
  1113.  
  1114. {
  1115. \ (F*CLEANUP)    Aux for F*                     08:47 20Oct88RS)
  1116.  
  1117. : (F*CLEANUP)   ( t1 flag -- d2 )
  1118.                 IF    ZEXP @ 2* NEGATE $80 + UNNORMALIZE DUP 0<
  1119.                         IF    $8000 = ZTEMP @ 0= AND MROUND
  1120.                         ELSE  DROP
  1121.                         THEN  EXIT             \ leave here
  1122.                 THEN    ROT DUP 0<
  1123.                 IF      $8000 = ZTEMP @ 0= AND MROUND DUP $0FF >
  1124.                         IF      D2/ $40 ZEXP +!  ZEXP @ $3FC0 >
  1125.                                 IF  ."  Overflow in F* " FCLEAR ABORT THEN
  1126.                         THEN
  1127.                 ELSE    DROP
  1128.                 THEN  $7F AND ZEXP @ 2* OR ;
  1129.  
  1130. }
  1131. \ (F*CLEANUP)    Aux for F*                     06:32 01Nov88RS)
  1132.  
  1133. (F*CLEANUP)  Auxillary function for F*
  1134.  
  1135. {
  1136. \ 1F*  Aux for F*                               07:30 15Oct88RS)
  1137.  
  1138. : 1F*   ( n1 -- n1 )
  1139.      DUP $7F80 AND DUP 2/ $40 + ZEXP !
  1140.      IF      $7F AND $80 OR
  1141.      ELSE    $7F AND NORMALIZE $80 + 2/ ZEXP !
  1142.      THEN ;
  1143.  
  1144. }
  1145.                                                 06:32 01Nov88RS)
  1146.  
  1147. 1F*     Auxillary function for F*
  1148.  
  1149. {
  1150. \  F*  Floating point multiply                  08:47 20Oct88RS)
  1151.  
  1152. : F*   ( F: r1 r2 -- r3 )
  1153.      FPOP FSP @ @ OVER XOR $8000 AND ZSIGN !  $7FFF AND 2DUP D0=
  1154.      IF      2DROP FDROP 0 ZSIGN @ FPUSH EXIT  THEN
  1155.      1F* FPOP $7FFF AND 2DUP D0=
  1156.      IF      2DROP 2DROP 0 ZSIGN @ FPUSH EXIT  THEN
  1157.      DUP $7F80 AND DUP
  1158.      IF      XBIAS - 2/ ZEXP +! $7F AND $80 OR
  1159.      ELSE    $7F AND NORMALIZE $80 + XBIAS - 2/ ZEXP +!
  1160.      THEN   Y* DUP 0< 0=
  1161.      IF      T2* $-40 ZEXP +!  THEN
  1162.      ROT ZTEMP !  DSHFT8 ZEXP @ DUP $3FC0 >
  1163.      IF  ."  Overflow in F* " FCLEAR ABORT  THEN   $40 < (F*CLEANUP)
  1164.      ZSIGN @ OR FPUSH ;
  1165.  
  1166. }
  1167. F*                                              06:32 01Nov88RS)
  1168.  
  1169. F*      Floating point multiplication function.
  1170.  
  1171. {
  1172. \ F**+N   Raise fp to positive integer power.   07:30 15Oct88RS)
  1173.  
  1174. : F**+N     ( F: r1 -- r2 ; n -- )
  1175.         $7FFF AND >R F1.0
  1176.         BEGIN R@ 1 AND
  1177.               IF   FOVER F*  THEN
  1178.               R> 2/ DUP
  1179.         WHILE >R FSWAP FDUP F* FSWAP
  1180.         REPEAT
  1181.         DROP FNIP ;
  1182.  
  1183. }
  1184.                                                 06:34 01Nov88RS)
  1185.  
  1186. F**+N   Raise the floating point number at the top of the f.p.
  1187.         stack to the positive integer power at the top of the
  1188.         parameter stack.
  1189.  
  1190. {
  1191. \ F**N  F**N*  MF**2                            13:52 16Jul89RS)
  1192.  
  1193. : F**N     ( F: r1 -- r2 ; n -- )  ( r1^n )
  1194.         DUP 0<
  1195.         IF    ABS F**+N F1.0 FSWAP F/
  1196.         ELSE  F**+N
  1197.         THEN ;
  1198.  
  1199. : F**N*    ( F: r1 r2 -- r3 ; n -- ) ( r1 * [r2^n] )
  1200.         DUP 0< IF   ABS F**+N F/ ELSE    F**+N F*  THEN ;
  1201.  
  1202. : MF**2   ( xlo xhi -- x^2lo x^2hi )
  1203.         DUP >R UM* NIP 0 D2*
  1204.         R> DUP UM* D+ ;
  1205.  
  1206. }
  1207. \ F**N  F**N*  MF**2                            06:39 01Nov88RS)
  1208.  
  1209. F**N    Raise the number at the top of the f.p. stack to the
  1210.         power specified at the top of the parameter stack.
  1211.  
  1212. F**N    Raise the number at the top of the f.p. stack to the
  1213.         power specified at the top of the parameter stack,
  1214.         then multiply by the number second on the f.p. stack.
  1215. MF**2   Square the double number mantissa on the parameter
  1216.         stack.
  1217.  
  1218. {
  1219. \  D2**N                                        07:31 15Oct88RS)
  1220.  
  1221. CREATE 2**NTAB  $0001 , $0002 , $0004 , $0008 , $0010 , $0020 ,
  1222.                 $0040 , $0080 , $0100 , $0200 , $0400 , $0800 ,
  1223.                 $1000 , $2000 , $4000 , $8000 ,
  1224.  
  1225.  
  1226. : D2**N    ( n -- d )
  1227.         $1F AND DUP $10 <
  1228.         IF      2* 2**NTAB + @ 0
  1229.         ELSE    $0F AND 2* 2**NTAB + @ 0 SWAP
  1230.         THEN ;
  1231.  
  1232. 2VARIABLE DROOT
  1233.  
  1234. }
  1235.                                                 06:47 01Nov88RS)
  1236.  
  1237. 2**NTAB  A table of 2 raised to various powers.
  1238.  
  1239. D2**N    Return a double number representing 2 raised to the
  1240.          power specified.
  1241.  
  1242. DROOT    A double number variable for temporary use with
  1243.          square roots.
  1244.  
  1245. {
  1246. \ Auxilliaries for Square Root.                 14:18 14Oct88RLS
  1247.  
  1248. : SQRTSTEP   ( d1 n1 -- d2 n2 )
  1249.         2* >R D2* D2* DUP R@ >
  1250.         IF      R@ 1 OR -  R> 2 OR
  1251.         ELSE    R> THEN ;
  1252.  
  1253. : DSQRTSTEP   ( d1 d2 -- d3 d4 )
  1254.         D2* DROOT 2!  D2* D2* 2DUP DROOT 2@ 2SWAP D<
  1255.         IF      DROOT 2@ 1 0 D+ D-
  1256.                 DROOT 2@ 2 0 D+
  1257.         ELSE
  1258.                 DROOT 2@
  1259.         THEN ;
  1260.  
  1261. }
  1262.                                                 06:48 01Nov88RS)
  1263.  
  1264. SQRTSTEP  Auxillary function for FSQRT .
  1265.  
  1266. DSQRTSTEP  Another auxillary function for FSQRT .
  1267.  
  1268. {
  1269. \ Auxilliary Square Root function.              08:45 20Oct88RS)
  1270.  
  1271. : FSQRT1   ( F: r1 -- ; -- n1 n2 n3 )
  1272.         FPOP0=  IF  FPUSH EXIT  THEN  DUP 0<
  1273.         IF      ."  Negative argument for FSQRT " FCLEAR ABORT  THEN
  1274.         DUP $7F80 AND DUP ZEXP !
  1275.         IF      $7F AND $80 OR
  1276.         ELSE    $7F AND NORMALIZE $80 + ZEXP +!
  1277.         THEN   DSHFT8 DROP  ZEXP @ XBIAS - DUP $80 AND
  1278.         IF      $80 - 2/ XBIAS + ZEXP !  0 D2* D2*
  1279.         ELSE    2/ XBIAS + ZEXP !  0 D2*
  1280.         THEN
  1281.         1- 2  7 0
  1282.         DO      SQRTSTEP  LOOP ;
  1283.  
  1284. }
  1285.                                                 06:48 01Nov88RS)
  1286.  
  1287. FSQRT1  Yet another auxillary function for FSQRT .
  1288.  
  1289. {
  1290. \ FSQRT                                         07:32 15Oct88RS)
  1291.  
  1292. : FSQRT   ( F: r1 -- r2 )
  1293.         FSQRT1 ROT DROP  5 0  DO  SQRTSTEP  LOOP
  1294.         ROT DROP 0 SWAP 0  $0C 0
  1295.         DO      DSQRTSTEP  LOOP
  1296.         D2/ OVER 1 AND
  1297.         IF      D2/ 2SWAP D0=
  1298.                 IF      1 0 D+ SWAP $FFFE AND SWAP
  1299.                 ELSE    1 0 D+
  1300.                 THEN
  1301.         ELSE    D2/ 2SWAP 2DROP
  1302.         THEN
  1303.         DUP $0FF >
  1304.         IF      D2/ -1 ZEXP +!   THEN
  1305.         $7F AND ZEXP @ OR FPUSH ;
  1306. }
  1307.                                                 06:49 01Nov88RS)
  1308.  
  1309. FSQRT   Replace the number at the top of the f.p. stack with
  1310.         its square root.
  1311.  
  1312. {
  1313. \ Tables for logarithms.                        07:32 15Oct88RS)
  1314.  
  1315. CREATE LOGTAB1
  1316.         $0000 , $0000 ,   $00FC , $14D8 ,   $01F0 , $A30C ,
  1317.         $02DE , $1A51 ,   $03C4 , $E0EE ,   $04A5 , $54BE ,
  1318.         $057F , $CC1C ,   $0654 , $96A7 ,   $0723 , $FDF2 ,
  1319.         $07EE , $461B ,   $08B3 , $AE56 ,   $0974 , $715D ,
  1320.         $0A30 , $C5E1 ,   $0AE8 , $DEE0 ,   $0B9C , $EBFB ,
  1321.         $0C4D , $19C3 ,   $0CF9 , $91F6 ,   $0D22 , $7BBE ,
  1322.         $0E47 , $FBE4 ,
  1323.  
  1324.  
  1325. CREATE LOGTAB2
  1326.         $0000 , $0000 ,   $0208 , $2BB1 ,   $0421 , $662D ,
  1327.         $064C , $D797 ,   $088B , $C741 ,   $0ADF , $A036 ,
  1328.         $0D49 , $F69E ,   $0FCC , $8E36 ,
  1329.  
  1330.  
  1331. }
  1332.                                                 06:50 01Nov88RS)
  1333.  
  1334. LOGTAB1  Auxillary table used for logarithm functions.
  1335.  
  1336. LOGTAB2  Auxillary table used for logarithm functions.
  1337.  
  1338. {
  1339. \ Auxilliary functions                          08:09 31Oct88RS)
  1340.  
  1341. : YLN2*   ( n -- d )
  1342.         2* DUP >R $B172 UM* R> $17F8 UM*
  1343.         $8000 0 D+ NIP 0 D+ 2NORMALIZE
  1344.         DUP IF  $80 -  THEN ;
  1345.  
  1346. : XLN2*   ( n -- d )  ( Multiply a shifted exponent by ln 2 )
  1347.         DUP 0<
  1348.         IF      ABS YLN2* $8000 OR
  1349.         ELSE    YLN2*
  1350.         THEN ;
  1351.  
  1352. }
  1353.                                                 08:12 31Oct88RS)
  1354.  
  1355. YLN2*   Multiply the positive shifted exponent by the natural
  1356.         logarithm of 2.
  1357.  
  1358. XLN2*   Multiply the shifted exponent by the natural log of 2.
  1359.  
  1360. {
  1361. \ Auxilliary function for FLN                   07:32 15Oct88RS)
  1362.  
  1363. : FLN+   ( F: -- r ; d1 n -- )
  1364.         $3F80 - XLN2* FPUSH
  1365.         DUP $FC AND DUP $7C AND >R  $100 * >R      ( R: div 4J )
  1366.         3 AND >R 0 SWAP D2/ D2/ D2/ $1FFF AND
  1367.         0 R> D2/ D2/ D2/ DROP OR
  1368.         R@ UM/MOD SWAP 0 SWAP R> UM/MOD NIP SWAP
  1369.         DUP >R $20 R@ 0<  IF  1-  THEN
  1370.         R@ FRACT* $0555 SWAP - DUP 0=
  1371.         IF      DROP R@
  1372.         ELSE    NEGATE R@ FRACT*
  1373.         THEN    R> FRACT* 0 SWAP DSHFT8 ROT DROP D2* D2*
  1374.         D- DSHFT8 ROT DROP R> LOGTAB1 + 2@ D+
  1375.         2NORMALIZE DUP IF  $300 -  THEN
  1376.         FPUSH F+ ;
  1377. }
  1378.                                                 06:56 01Nov88RS)
  1379.  
  1380. FLN+    Auxillary function for the natural logarithm function
  1381.         used when the mantissa is less than or equal to 1.5625 .
  1382.  
  1383. {
  1384. \ Auxilliary function for Logarithm             13:37 20Jul89RLS
  1385.  
  1386. : FLN-   ( F: -- r ; d n -- )
  1387.         $3F00 - XLN2* FPUSH DUP $F8 AND 8 + $100 * >R
  1388.         0 $100 2SWAP D- DUP $F8 AND 2/ R> SWAP >R >R
  1389.         7 AND >R 0 SWAP D2/ D2/ D2/ D2/ $FFF AND
  1390.         0 R> D2/ D2/ D2/ D2/ DROP OR R@ DUP
  1391.         IF  UM/MOD SWAP 0 SWAP R> UM/MOD NIP SWAP
  1392.         ELSE    R> 2DROP
  1393.         THEN   DUP >R 2DUP MF**2 DUP R@ $033 FRACT* $400 +
  1394.         R@ FRACT* $5555 + R> FRACT* FRACT* 0 SWAP D2/ D2/ D2/
  1395.         D+ D2/ D2/ D2/ D2/ D2/ D+  D2/ $7FFF AND
  1396.         D2/ D2/ D2/ D2/ D2/ R> LOGTAB2 + 2@ D+ 2DUP D0= 0=
  1397.         IF      2NORMALIZE $380 - $8000 OR  THEN
  1398.         FPUSH F+ ;
  1399.  
  1400. }
  1401.                                                 06:56 01Nov88RS)
  1402.  
  1403. FLN-    Auxillary function for natural logarithm, used when
  1404.         the mantissa is greater than 1.5625 and less than 2.0 .
  1405.  
  1406. {
  1407. \ FLN   Natural Logarithm                       08:47 20Oct88RS)
  1408.  
  1409. : FLN   ( F: r1 -- r2 )
  1410.         FPOP0=
  1411.         IF      CR ."  Zero argument for FLN "
  1412.                 2DROP -1 -1 FPUSH EXIT
  1413.         THEN  DUP 0<
  1414.         IF    ."  Negative argument for FLN " FCLEAR ABORT  THEN
  1415.         DUP $7F80 AND DUP 0=
  1416.         IF      1NORMALIZE
  1417.         ELSE    SWAP $7F AND $80 OR SWAP
  1418.         THEN
  1419.         OVER $C8 >
  1420.         IF      FLN-   ELSE    FLN+   THEN ;
  1421.  
  1422. }
  1423.                                                 06:57 01Nov88RS)
  1424.  
  1425. FLN     Replace the number at the top of the f.p. stack with
  1426.         its natural logarithm.
  1427.  
  1428. {
  1429. \ FLOG ( Common Logarithm ) and FPARTS          07:41 18Jul89RS)
  1430.  
  1431. : FLOG   ( F: r1 -- r2 )
  1432.         FLN FLOG10E F* ;
  1433.  
  1434. : FPARTS   ( F: r1 -- ; n -- d exp sign )   ( Aux for E. )
  1435.         8 MIN 1 MAX FDUP F0=
  1436.         IF      FDROP DROP 0 0 0 0 EXIT  THEN
  1437.         FSP @ @ 0< >R
  1438.         FABS FDUP FLOG INT DROP
  1439.         2DUP - F10.0 F**N* F0.5 F+ FINT
  1440.         SWAP FDUP F10.0 F**N F< 0=
  1441.         IF      F10.0 F/ 1+  THEN
  1442.         >R INT R> R> ;
  1443.  
  1444. }
  1445.                                                 06:59 01Nov88RS)
  1446.  
  1447. FLOG    Replace the number at the top of the f.p. stack with
  1448.         its common (base 10) logarithm.
  1449.  
  1450. FPARTS  An auxillary function for  (E.) .
  1451.  
  1452. {
  1453. \ Numeric output E.                             15:08 23Jul89RS)
  1454.  
  1455. VARIABLE F#PLACES
  1456.  
  1457. : (E.)   ( F: r -- ; n -- addr cnt )
  1458.         F#PLACES @
  1459.         FPARTS BASE @ >R >R DECIMAL <# DUP ABS 0 # # 2DROP 0<
  1460.         IF     $2D   ELSE $2B  THEN     ( Send "-" or "+" )
  1461.         HOLD $45 HOLD                    ( Send the "E" )
  1462.         F#PLACES @ 0 DO  #  LOOP                 ( Send the fraction )
  1463.         $2E HOLD R> 0<                   ( Send "." Check sign )
  1464.         IF      $2D   ELSE  $20  THEN
  1465.         HOLD                            ( Send "-" or space )
  1466.         R> BASE ! #> ;                  ( Restore BASE )
  1467.  
  1468. : E.            ( F: r -- )
  1469.                 8 F#PLACES ! (E.) TYPE SPACE ;
  1470.  
  1471. : E.R           ( F: r -- ; places width -- )
  1472.                 >R 8 min 1 max F#PLACES ! (E.)   R> OVER - SPACES TYPE   ;
  1473.  
  1474. : .FS           ( -- )
  1475.                 FDEPTH  IF   CR FSP @ FDEPTH F#BYTES * 0 DO  DUP I + 2@
  1476.                 FPUSH E. F#BYTES  +LOOP  DROP  ELSE ."  Empty"  THEN ;
  1477.  
  1478. : (F.)          ( F: r -- ; n -- addr cnt )
  1479.                 F#PLACES @ FPARTS BASE @ >R >R DECIMAL <#
  1480.                 F#PLACES @ SWAP - 0max 0 ?DO # LOOP
  1481.                 $2E HOLD #S R> 0<               ( Send "." Check sign )
  1482.                 IF      $2D   ELSE  $20  THEN
  1483.                 HOLD                            ( Send "-" or space )
  1484.                 R> BASE ! #> ;                  ( Restore BASE )
  1485.  
  1486. : F.            ( F: r -- )
  1487.                 8 F#PLACES ! (F.) TYPE SPACE ;
  1488.  
  1489. : F.R           ( F: r -- ; places width -- )
  1490.                 >R 8 min 1 max F#PLACES ! (F.)   R> OVER - SPACES TYPE   ;
  1491.  
  1492. }
  1493.  
  1494. (E.)  E.                                        16:06 23Jul89RS)
  1495.  
  1496. (E.)    Auxillary function for E.
  1497.  
  1498. E.      The floating point output routine.
  1499.  
  1500. .FS     A utility for checking the contents of the F.P. stack.
  1501.  
  1502. {
  1503. \ Auxilliary finctions for numeric input.       14:05 14Oct88RLS
  1504.  
  1505. : Ee(?   ( n -- flag )
  1506.         DUP 40 =                     ( Check for left paren )
  1507.         IF    DROP -1
  1508.         ELSE  DUP 69 = SWAP 101 = OR  ( Check for "e" or "E" )
  1509.         THEN ;
  1510.  
  1511. : -?    ( addr1 -- addr2 flag )
  1512.         DUP 1+ C@ 45 =
  1513.         IF      1+ -1
  1514.         ELSE    0
  1515.         THEN ;
  1516.  
  1517. }
  1518.                                                 07:03 01Nov88RS)
  1519.  
  1520. E.(?    Auxillary function to test for exponential indicator
  1521.         in floating point input.  The indicator should be one
  1522.         of the following:   E  e  (
  1523.  
  1524. -?      Check the character at the specified address for
  1525.         ASCII - sign.  If found, increment the address pointer
  1526.         and return a true flag.  Otherwise, return a false flag.
  1527.  
  1528. {
  1529. \ Numeric Input Conversion                      11:27 07Nov88RS)
  1530.  
  1531. : (FNUMBER?)    ( a1 -- f1 ; F: -- r )  \ convert string a1 to floating point #
  1532.         0 0     ROT -? >R DUP 1+ C@ Ee(?
  1533.         IF      1+ ROT DROP 1 -ROT 0 >R   ( 1 to mantissa )
  1534.         ELSE    CONVERT DUP C@ 46 =          ( Check for "." )
  1535.                 IF      DUP >R CONVERT DUP R> - 1-
  1536.                 ELSE    0  THEN    >R
  1537.         THEN
  1538.         DUP C@ Ee(?
  1539.         IF      -? >R >R 0 0 R> CONVERT NIP
  1540.         ELSE    0 SWAP 0 >R THEN
  1541.         C@ DUP 0= SWAP DUP 32 = SWAP 41 = OR OR 0=
  1542.         IF      DROP 2DROP R>DROP R>DROP R>DROP F0.0 FALSE
  1543.         ELSE    R> IF   NEGATE  THEN  R> - BASE @ 0 FLOAT F**N
  1544.                 FLOAT F* R>  IF  FNEGATE  THEN
  1545.                 ( TRUE ) 1      \ return 1 to be compatible with SFLOAT
  1546.         THEN    ;
  1547.  
  1548. : $F#   ( a1 -- F: -- r )       \ convert string a1 to floating point #
  1549.         (FNUMBER?) 0= IF ."  Bad Floating Point Input" FCLEAR ABORT THEN ;
  1550.  
  1551. : F#    ( F: -- r )             \ convert string from input stream to
  1552.                                 \ a floating point nubmer
  1553.         $20 WORD $F# ;
  1554.  
  1555. : [F#]  ( F: -- r )             \ convert string in colon def to an
  1556.                                 \ inline literal floating point number
  1557.         F# FPOP SWAP
  1558.         [COMPILE] LITERAL [COMPILE] LITERAL
  1559.         COMPILE FPUSH ; IMMEDIATE
  1560.  
  1561. : F,    ( F: r -- )             \ compile a floating point number into
  1562.                                 \ the dictionary
  1563.         FPOP , , ;
  1564.  
  1565. }
  1566.                                                 07:07 01Nov88RS)
  1567.  
  1568. $F#     This function converts a counted string into a floating point
  1569.         number.
  1570.  
  1571. F#      This is the function used to get a floating point number
  1572.         from the input stream.  Usage examples follow:
  1573.  
  1574.                 F# -2.34
  1575.                 F# 34.5e6
  1576.                 F# -.1E-2
  1577.                 F# 2.34(5)
  1578.  
  1579.