(     Title:  Convert among floating-point 32-, 64-, 80-, and
              128-bit raw, C99 mixed hex/decimal string, raw hex
              string, and Forth default formats
       File:  mixfloat.fs
  Test file:  mixfloat-test.fs
   Log file:  mixfloat.log
    Version:  0.9.0
    Revised:  December 23, 2010
     Author:  David N. Williams
    License:  LGPL

We believe any of this code derived from other authors to be
either in the public domain or otherwise compatible with the
LGPL.  For the sake of the LGPL, the rest is
)
\ Copyright (C) 2003, 2005, 2009, 2010  David N. Williams
(
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or at your option any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

If you take advantage of the option in the LGPL to put a
particular version of this library under the GPL, the author[s]
would regard it as polite if you would put any direct
modifications under the LGPL as well, and include this paragraph
in your license notice.  A direct modification is
one that enhances or extends the library in line with its
original concept, as opposed to developing a distinct
application or library which might use it.

The library converts among a number of floating-point formats
based on the binary part of the IEEE 754-2008 standard,
including mixed strings with hex mantissa digits and decimal
radix two exponent.

The central floating-point encoding formats are "raw" formats,
raw32, raw64, raw80, and raw128, which are defined on the data
stack.  The numbers in the names are the storage widths in bits.
Whether a cell is 32 or 64 bits, the raw formats occupy the
high bits of enough cells on the data stack to contain them,
with the most significant cell topmost, the next most
significant just below it, and so on, towards the bottom of the
stack.

The raw32, raw64, and raw128 formats correspond to IEEE 754-2008
interchange formats, binary32, binary64, and binary128, which
all have an implicit leading integer bit for normal and
subnormal numbers.  IEEE does not define a binary interchange
format for intel 80-bit floating-point numbers, and that seems
to be deliberate.  The raw80 format defined here has an explicit
leading integer bit, to simplify interaction with the intel
format.

Many of the conversions are independent of the DPANS94 default
float format.  Only near the end does the code depend on the
representation of floating point numbers in memory, and only
closer to the end does it depend on the floating-point stack.

The library works with both big- and little-endian systems, as
long as the memory representation is totally one or the other.

We have used this code as a tool for manipulating floating point
constants for high-precision algorithms, and as a tool for
checking floating point extensions to ANS Forth.

Upon inspection, neither here nor in mixfloat-test.fs do we see
words with mixed data/float inputs, outputs, or intermediate
calculations.  Everything should work whether or not the
floating-point stack is separate from the data stack, and that
has been checked with both options in pfe.  See the section
"FLOATING-POINT ENVIRONMENT" for how to select an integrated
stack in pfe.

There is an environmental dependence on lower case.

The library eventually refuses to load if cells are not 32 bits,
chars are not 8 bits and 1 address unit, or the arithmetic is
not two's complement.

The library currently does not allow nans in mixed input
strings.  It does produce them in mixed output strings,
including signs, but ignoring loads and quietness.  We expect
eventually to treat nan loads.

The test file mixfloat-test.fs contains many syntax examples,
and a dictionary of floating-point bit patterns.

Notation:
"s" or "something.s" in stack comments stands for an ANS
    Forth string [addr len] pair.
"mix" in a word name indicates an ANS Forth string
    representing a floating-point number with hex significand
    and radix two decimal power.
"hex" in a word name indicates an ANS Forth string
    representing a raw format as hex digits, including all
    leading zero digits.  A trailing "8", "16", "20", or "32"
    indicates the number of hex digits corresponding to raw32,
    raw60, raw80, and raw64, respectively.
"p", "emax", and "emin" are the characteristic parameters of an
    IEEE 754 binary format, where p is the number of bits in a
    normalized mantissa, including the leading integer bit,
    whether implicit or not, where emax is the maximum power of
    2 for a normalized number, and where emin = 1 - emax < 0 <
    emax.
)

\ *** USER WORDS

\ CONSTANTS
\   bits/au  bits/cell  bits/float

\ UTILITY
\   lshift>ud  rshift>ud
\   (qlshift)  (qrshift) qlshift  qrshift  qshift
\   qfaligned  qfalign
\   -spaces  hex>ud  hex>uq  sdec>n
 
\ MIX/RAW CONVERSIONS
\   mix>raw32  mix>raw64  mix>raw80  mix>raw128  mix>raw
\   raw32>mix  raw64>mix  raw80>mix  raw128>mix  raw>mix

\ HEX/RAW CONVERSIONS
\   hex8>raw32  hex16>raw64  hex20>raw80  hex32>raw128  hex>raw
\   raw32>hex8  raw64>hex16  raw80>hex20  raw128>hex32  raw>hex

\ HEX/MIX CONVERSIONS
\   hex8>mix  hex16>mix  hex20>mix  hex32>mix  hex>mix
\   mix>hex8  mix>hex16  mix>hex20  mix>hex32  mix>hex

\ RAW FETCH/STORE
\   raw32!  raw64!  raw80!  raw128!  raw!
\   raw32@  raw64@  raw80@  raw128@  raw@

\ RAW/FLOAT CONVERSIONS
\   raw32>f  raw64>f  raw80>f  raw128>f  raw>f
\   f>raw32  f>raw64  f>raw80  f>raw128  f>raw

\ HEX/FLOAT CONVERSIONS
\   hex8>f  hex16>f  hex20>f  hex32>f  hex>f
\   f>hex8  f>hex16  f>hex20  f>hex32  f>hex

\ MIX/FLOAT CONVERSIONS
\   mix>f  f>mix  fm.


\ *** GENERAL USE

decimal

[UNDEFINED] -rot [IF] : -rot   rot rot ; [THEN]
[UNDEFINED] 0>=  [IF] : 0>=    0< 0= ;   [THEN]
[UNDEFINED] 0<=  [IF] : 0<=    0> 0= ;   [THEN]

[UNDEFINED] unders+ [IF]
  : unders+  ( n1 n2 n3 -- n1+n2 n3 ) -rot + swap ; [THEN]
[UNDEFINED] rdrop [IF]
  : rdrop  ( r: a -- )  s" r> drop" evaluate ; immediate [THEN]
[UNDEFINED] toupper [IF]
  : toupper ( char -- char' )
    dup [char] a [char] { within
    IF [ 1 5 lshift invert ] literal and THEN ; [THEN]

: hex-digit  ( char -- +n true | false )
  ( char) toupper
  dup [char] 0 [char] : within 
  IF   [char] 0 -
  ELSE dup [char] A [char] G within
       IF   [ char A 10 - ] literal -
       ELSE drop false EXIT THEN
  THEN true ;

[UNDEFINED] bits/cell [If]
s" ADDRESS-UNIT-BITS" environment? 0=
[IF] cr .( ***Can't determine ADDRESS-UNIT-BITS.) ABORT [THEN]
                  constant bits/au
bits/au 1 cells * constant bits/cell
[THEN]

[UNDEFINED] lshift>ud [IF]
: lshift>ud  ( u +n -- ud' )
(
Assume +n <= BITS/CELL.  Logical shift u to the left by +n bits
into an unsigned double number.
)
  LOCALS| n u |
  n 0=          IF u 0 EXIT THEN
  n bits/cell = IF 0 u EXIT THEN
  u n lshift u bits/cell n - rshift  ;
[THEN]

[UNDEFINED] rshift>ud [IF]
: rshift>ud  ( u +n -- ud )
(
Assume +n <= BITS/CELL.  Logical shift u to the right by +n bits
into an unsigned double number.
)
  LOCALS| n u |
  n 0=          IF 0 u EXIT THEN
  n bits/cell = IF u 0 EXIT THEN
  u bits/cell n - lshift u n rshift ;
[THEN]

[UNDEFINED] (qlshift) [IF]
: (qlshift)  ( q +n -- q' )
(
Assume +n <= BITS/CELL.  Logical shift the quadruple q to the
left by +n bits.
)
  LOCALS| n q0 q1 q2 q3 |
  n bits/cell = IF 0 ELSE q0 n lshift THEN ( q0'.hi)
  q1 n lshift>ud ( q0'.hi q1'.hi q0'.lo)
  rot + ( q0') to q0
  ( q1'.hi) q2 n lshift>ud ( q1'.hi q2'.hi q1'.lo)
  rot + ( q1') to q1
  ( q2'.hi) q3 n lshift>ud ( q2'.hi q3'.hi q2'.lo)
  rot + q1 q0 ;
[THEN]

[UNDEFINED] (qrshift) [IF]
: (qrshift)  ( q +n -- q' )
(
Assume +n <= BITS/CELL.  Logical shift the quadruple q to the
right by +n bits.
)
  LOCALS| n q0 q1 q2 q3 |
  n bits/cell = IF 0 ELSE q3 n rshift THEN ( q3'.lo)
  q2 n rshift>ud ( q3'.lo q3'.hi q2'.lo)
  unders+ ( q3' q2'.lo)
  ( q2'.lo) q1 n rshift>ud ( q2'.lo q2'.hi q1'.lo)
  unders+ ( q3' q2' q1'.lo)
  ( q1'.lo) q0 n rshift>ud ( q1'.lo q1'.hi q0'.lo)
  unders+ ;
[THEN]

: qlshift  ( g +n -- q' )
(
Logical shift the quadruple q to the left by +n bits.
)
  bits/cell /mod >r (qlshift) r>
  0 ?DO bits/cell (qlshift) LOOP ;

: qrshift  ( g +n -- q' )
(
Logical shift the quadruple q to the right by +n bits.
)
  bits/cell /mod >r (qrshift) r>
  0 ?DO bits/cell (qrshift) LOOP ;

: qshift  ( q n -- q' )
(
Logical shift the quadruple q to the left/right by |n| bits when
n is positive/negative.
)
  dup 0< IF negate qrshift ELSE qlshift THEN ;

[UNDEFINED] qfaligned [IF]
  : qfaligned  ( addr -- addr' )  127 + [ 127 invert ] literal and ;
[THEN]

[UNDEFINED] qfalign [IF]
  : qfalign  ( -- )  here qfaligned here - allot ;
[THEN]

[UNDEFINED] -spaces [IF]
: -spaces  ( s -- s' )
(
Remove all blanks from the ANS Forth string s and leave the
compacted string s'.  The new string body address is at PAD.
This works even if the input string s starts somewhere in PAD,
but then it may get partially written over.
)
  pad >r
  BEGIN
    ( len) dup 0>
  WHILE
    over c@ bl <>
    IF over c@ r@ c! r> 1+ >r THEN
    1 /string
  REPEAT 2drop
  pad r> pad - ;
[THEN]

[UNDEFINED] hex>ud [IF]
: hex>ud  ( s -- ud flag ) 
(
If the string contains only hexadecimal digits with no sign, or
is empty, leave true and the unsigned double number conversion
[0 if there are no digits].  Otherwise leave an undefined ud
result and false.
)
  base @ >r hex
  0 0 2swap >number ( ud rest.s)
  nip 0=   r> base ! ;

: ?hex>ud  ( s -- ud )  hex>ud 0= ABORT"  ***Not a hex string" ;

: hex>uq  ( s -- uq flag ) 
(
If the string contains only hexadecimal digits with no sign, or
is empty, leave true and the unsigned quadruple number
conversion [0 if there are no digits].  Otherwise leave an
undefined uq result and false.
)
  ( len) dup 16 <=
  IF over 0 2swap ELSE 16 - 2dup + 16 THEN ( hi.s lo.s)
  LOCALS| lo.len lo.addr hi.len hi.addr |
  base @ >r hex
  0 0 lo.addr lo.len >number
  IF
    ( lo.ud addr) 0 false
  ELSE  
    ( lo.ud addr) drop 0 0 hi.addr hi.len >number
    ( lo.ud hi.ud addr' len') nip 0=
  THEN  
  r> base ! ;

: ?hex>uq  ( s -- uq )  hex>uq 0= ABORT" ***Not a hex string" ;
[THEN]

: {}in  ( char &chars -- char|0 )
(
Leave char if it is in the null-terminated list of characters
starting at &chars.  Otherwise leave zero.
)
  ( &chars) >r
  BEGIN ( char)
    r@ c@ dup
  WHILE
    r> 1+ >r
    over =
  UNTIL
  ELSE
    2drop false
  THEN r> drop ;

: {}scan  ( s &chars -- s' )
(
Leave the substring of s starting at the first occurrence of one
of the characters in the null-terminated string starting at
&chars.  If no such character is found and s is not empty, s' is
the empty string, with address just after the end of s.  If s is
empty, s' is s.

This is a variation on SCAN from Wil Baden's Tool Belt.
)
  ( &chars) >r
  BEGIN
    ( len) dup
  WHILE
    over c@ r@ {}in 0=
  WHILE
    1 /string
  REPEAT
  THEN r> drop ;

: {}scan-past  ( s &chars -- s' char|0 )
(
Scan s for the first occurrence of any of the characters in the
null-terminated string starting at &chars.  If the character is
found, leave it with s' the part of s that follows it, or the
empty string with address just after the character if it is
found at the end of s.  If s is empty, leave s and zero.  If s
is nonempty and the character is not found, leave zero and s'
the empty string with address just after the end of s.
)
  {}scan dup
  IF over c@ >r 1 /string r>
  ELSE 0 THEN ;

: /split  ( a m a+i m-i -- a+i m-i a i )
(
Split the character string a m at the place given by a+i m-i. 
Wil Baden's "cut-split", from the Tool Belt.
)
  dup >r  2swap  r> - ;

: {}separate  ( s &chars -- after.s before.s char | s false )
(
Scan s for the first occurrence of any of the null-terminated
sequence of characters at &chars.  If a character is found,
leave it as char, and leave the strings after and before.  If
not found, leave zero instead and the original string.
)
  ( &chars) >r
  ( s) 2dup r>
  {}scan-past ( char|0) dup
  IF   ( char) >r ( s after.s) /split 1- r>
  ELSE -rot 2drop THEN ;

: skip  ( addr len char -- addr+i len-i )
(
Advance past leading characters.  From Wil Baden's Tool Belt.
)
  >r BEGIN dup WHILE over c@ r@ =
     WHILE 1 /string REPEAT THEN
  r> drop ;

: skip-back  ( addr len char -- addr len-i )
(
Trim trailing instances of char from the input string.  Based on
TRIM from Wil Baden's Tool Belt.
)
    ( char) >r
    BEGIN
      ( len) dup
    WHILE
      1-  2dup + c@ r@ <>
    UNTIL
      1+
    THEN r> drop ;

: scan  ( addr len char -- addr+i len-i )
(
Look for a particular character in the input string.  From Wil
Baden's Tool Belt.
)
    >r  BEGIN  dup WHILE  over c@ r@ -
        WHILE  1 /string  REPEAT THEN
    r> drop ;


\ *** MIXED TO FLOATING-POINT PARTS

: minus?  ( s -- s' flag )
(
If s is not empty and has a plus or minus sign as its first
character, s' is s skipped past the first character.  Otherwise
s' is s.  If the first character is a minus sign, flag is true.
Flag is false if the first character is a plus sign, not a sign,
or if s is empty.
)
  ( len) dup IF
    over c@ >r
    r@ [char] - =   r@ [char] + = or
    IF 1 /string THEN
    r> [char] - = 
  ELSE
    false
  THEN ;

create p|P  char p c,   char P c,   0 c,

: sdec>n  ( s -- n flag )
(
If the string consists of an optional leading plus or minus sign
followed only by optional decimal digits, leave true and the
signed single number conversion [0 if there are no digits].
Otherwise leave an undefined result and false.  There is no
check for overflow.
)
  base @ >r decimal
  minus? >r
  0 0 2swap >number nip nip             ( ud.lo len)
  0= swap                               ( valid? ud.lo)
  r> ( minus?) IF negate THEN
  swap                                  ( n valid?)
  r> base ! ;

: -first-dot ( s -- s' #before )
(
Delete the first "." [dot] from s.  If there is no dot, s' is
the same as s, and #before is a copy of its length.  If there is
a dot, any characters following it are moved up one position, so
s' has the same address as s and length reduced by one.  Then
#before is the number of characters before the dot in s.
)
  2dup [char] . scan
  ( len-i) dup 0= IF
    2drop dup
  ELSE  ( s dot.addr 1+#after)
    1- ( #after) >r
    ( dot.addr) dup 1+ swap ( after.addr dot.addr) r@ cmove
    ( s.len) 1- dup r> -
  THEN ;

: bin-defect  ( nonzero.nibble -- -#leading.bin.zeros )
  0 >r  BEGIN 2* dup 16 and 0=
        WHILE r> 1+ >r REPEAT drop r> negate ;

1 bits/cell 1- lshift constant fsign-mask

: mix>fparts  ( s -- uqmant #mantbits binexp signbit )
(
The input string is a mixed hexadecimal number with radix two
exponent written in decimal, in the style of the C99 standard
but without a leading "0x".

Blanks are allowed anywhere but are ignored.

The hex part on the left and the power part on the right must be
separated by "p" or "P".

A plus or minus sign may precede either part, and the hex part
can include at most one hex point.

The absence of any digits in either part makes it zero.

The quadruple number uqmant is the integer that results from
the hex part by dropping the hex point, if present, and then
dropping trailing hex zeros.  It is considered to have an
implicit binary point just to the right of its leading nonzero
bit, called the "normalized position", unless it is zero.  There
is an abort if the hex mantissa is too large to fit into four
cells.  If the hex part is zero, including having no digits,
uqmant is zero.

The number of binary digits in uqmant, starting with the leading
nonzero bit, is #mantbits.  Although a nonzero uqmant has no
trailing zero hex digits, it may have up to three trailing zero
bits; uqmant is zero iff #mantbits is zero.

When uqmant is nonzero, binexp is the signed power of two
required to move the implicit binary point to the correct
position in nonexponential notation, taking into account both
the power of two part and any hex point in the hex part.  It is
said to be normalized and is assumed to fit into a single cell.

When uqmant is zero, binexp is zero.

The leading-bit flag signbit determines the sign associated with
the mantissa. That can be minus even when uqmant is zero.
)
  -spaces
  minus? IF fsign-mask ELSE 0 THEN 0 0
  LOCALS| binexp #mantbits signbit |

  p|P {}separate
  0= ABORT"  ***Missing 'p' or 'P'."    ( exp.s hex.s)

  \ Remove leading 0's from the hex part.
  [char] 0 skip

  \ Remove any hex point and get #intdigs.  Note that
  \ if the hex point is the first char, #intdigs is 0. 
  -first-dot ( #intdigs) >r

  \ Remove trailing 0's from the hex part.
  [char] 0 skip-back

  \ Remove and negatively count any leading frac 0's
  \ revealed if the hex point was the first char.
  ( len) dup >r   [char] 0 skip         ( exp.s mant.s)
      ( len') dup r> ( len) -
  
  r> ( -#fracdigs #intdigs)  \ one of these *has* to be zero
      + ( hexexp) 4 * to binexp  \ maybe zero
  \ At this point the implicit binary point is just left of
  \ the leading hex digit, if there is one.

  ( #hexdigs=len) dup 8 cells > ABORT" ***Hex mantissa too large."
  dup 4 * to #mantbits  \ includes leading hex digit

  \ Adjust for leading 0 bits in the first hex digit, and move
  \ the implicit binary point just right of the leading 1 bit.
  ( #hexdigs) dup
  IF
    over c@ hex-digit 0= ABORT" ***Number not hexadecimal."                    
    ( nibble) bin-defect
        dup #mantbits + to #mantbits
        binexp + 1- to binexp  \ binary point moved 1 bit to right
  THEN
   
  ( exp.s mant.s) 2swap
  ( exp.s) sdec>n   0= ABORT" ***Improper decimal power of 2."
  #mantbits IF ( power) binexp + ELSE ( power) drop 0 THEN
  to binexp

  ( mant.s) hex>uq  0= ABORT" ***Number not hexadecimal."
  #mantbits binexp signbit ;

: >p-bits  ( mant.uq #bits #p -- mant.p.uq )
(
Assume that the number of precison bits #p is no larger than the
number of bits in a quadruple number.  If #bits is zero, assume
that the quadruple mant.uq is zero, and leave mant.p.uq as
mant.uq.  Else assume that mant.uq is nonzero and that the
leading nonzero bit and any trailing bits occupy its least
significant #bits.  Do a logical shift of mant.uq so that its
leading nonzero bit and as many trailing zero bits as necessary
occupy the rightmost p bits of the quadruple, and leave that as
mant.p.uq.  If #bits is more than p, mant.p.uq is the truncated
result.
)
  ( p) >r ( #bits) ?dup
  IF ( #bits) r@ - negate qshift THEN 
  rdrop ;


\ *** MIX/RAW CONVERSIONS
(
From here on, we assume that a cell is 32 bits, and two's
complement arithmetic.

We do not assume anything about endianness, nor do we use the
floating-point stack, until further notice.
)

bits/cell 32 <> [IF] cr
  .( ***bits/cell is ) bits/cell . .( bits, not 32) ABORT [THEN]

hex FFFFFFFF decimal -1 <> [IF] cr
  .( ***arithmetic is not two's complement) ABORT [THEN]

\ These masks are only needed for the high 32 bits of the raw
\ format.
decimal     \ assume 2's complement
fsign-mask invert constant fbody-mask
1 23 lshift 1-    constant smantissa-mask
1 20 lshift 1-    constant dmantissa-mask
1 16 lshift 1-    constant qmantissa-mask

\ The low 16 bits of raw80 are in the high 16 bits of the third
\ cell:
qmantissa-mask invert constant xmantissa-lo-mask

fbody-mask smantissa-mask xor constant sexp-mask
fbody-mask dmantissa-mask xor constant dexp-mask
fbody-mask qmantissa-mask xor constant qexp-mask

hex
7F800000          constant  Inf-raw32       \  8-bit exp
00000000 7FF00000 2constant Inf-raw64       \ 11-bit exp
00000000 7FFF8000 2constant Inf-raw80.2hi   \ 15-bit exp
00000000 7FFF0000 2constant Inf-raw128.2hi  \ 15-bit exp
decimal

\ *** MIXED TO RAW

\ Note that the least significant 16 bits of mant.64 are
\ unspecified, and the leading cell is zero.
: >24-bits  ( mant.uq #bits -- mant.24 )   24 >p-bits 2drop drop ;
: >53-bits  ( mant.uq #bits -- mant.53 )   53 >p-bits 2drop ;
: >113-bits ( mant.uq #bits -- mant.113 ) 113 >p-bits ;
: >64-bits  ( mant.uq #bits -- mant.64 )   80 >p-bits drop ;

\ MIX>RAW32, MIX>RAW64, MIX>RAW80, MIX>RAW128
(
Convert a string in mixed hex with radix two decimal exponent
floating point notation to the indicated raw format.  Truncate
any mantissa bits that do not fit.

Abort if the string is not legally formatted.

If the number is nonzero and the normalized binary exponent is
larger than emax, leave the raw Inf with the appropriate sign.

If the number is nonzero and the normalized binary exponent is
less than emin, leave the corresponding raw subnormal number. 
Truncate any mantissa bits that do no fit, and leave zero with
the appropriate sign when none fit.
)

: mix>raw32  ( s -- raw32 )
  mix>fparts
  LOCALS| signbit exp #bdigs mant0 mant1 mant2 mant3 |
 
  #bdigs IF
    mant3 mant2 mant1 mant0 #bdigs >24-bits to mant0

    exp 127 > IF Inf-raw32 signbit or EXIT THEN

    exp 127 + to exp
    exp 0<=    \ subnormal?
    IF exp -22 <
      IF 0 ( hi)
      ELSE mant0 1 exp - rshift
      THEN to mant0 0 to exp
    THEN

    mant0 smantissa-mask and exp 23 lshift or

  ELSE 0 THEN signbit or ;

: mix>raw64  ( s -- raw64 )
  mix>fparts
  LOCALS| signbit exp #bdigs mant0 mant1 mant2 mant3 |

  #bdigs IF
    exp 1023 > IF Inf-raw64
    ELSE
      mant3 mant2 mant1 mant0 #bdigs >53-bits to mant0 to mant1
      exp 1023 + to exp
      exp 0<=  \ subnormal?
      IF exp -51 <
        IF   0 to mant1 0 to mant0  \ zero
        ELSE 0 0 ( mant3 mant2) mant1 mant0 exp 1- qshift
             to mant0 to mant1 2drop
        THEN 0 to exp
      THEN
      mant1 mant0 dmantissa-mask and exp 20 lshift or
    THEN
  ELSE 0 0 THEN signbit or ;

: mix>raw128  ( s -- raw128 )
  mix>fparts
  LOCALS| signbit exp #bdigs mant0 mant1 mant2 mant3 |

  #bdigs IF
    exp 16383 > IF [ 0 0 ] 2literal Inf-raw128.2hi
    ELSE
      mant3 mant2 mant1 mant0 #bdigs >113-bits
          to mant0 to mant1 to mant2 to mant3
      exp 16383 + to exp
      exp 0<=  \ subnormal?
      IF exp -111 <
        IF   0 to mant1 0 to mant0
             0 to mant2 0 to mant3 \ zero
        ELSE mant3 mant2 mant1 mant0 exp 1- qshift
             to mant0 to mant1 to mant2 to mant3
        THEN 0 to exp
      THEN
      mant3 mant2 mant1 mant0 qmantissa-mask and
      exp 16 lshift or
    THEN
 ELSE [ 0 0 ] 2literal [ 0 0 ] 2literal THEN signbit or ;

: mix>raw80  ( s -- raw80 )
  mix>fparts
  LOCALS| signbit exp #bdigs mant0 mant1 mant2 mant3 |

  #bdigs IF
    exp 16383 > IF 0 Inf-raw80.2hi
    ELSE
      mant3 mant2 mant1 mant0 #bdigs >64-bits
          to mant0 to mant1 to mant2
      exp 16383 + to exp
      exp 0<=  \ subnormal?
      IF exp -62 <
        IF   0 to mant1 0 to mant0 0 to mant2  \ zero
        ELSE 0 ( mant3) mant2 mant1 mant0 exp 1- qshift
             to mant0 to mant1 to mant2 drop
        THEN 0 to exp
      THEN
      mant2 xmantissa-lo-mask and
          mant1 mant0 qmantissa-mask and
      exp 16 lshift or
    THEN
 ELSE 0 [ 0 0 ] 2literal THEN signbit or ;

\ *** RAW TO MIXED

: num>s  ( n -- n.s )
(
Convert the number according to the current BASE to a Forth
string in transient memory.  Traditional.
)
  <# dup >r abs s>d #s r> sign #> ;

(
Because char-counted strings are deprecated, we use "measured
strings" [mstrings], which are strings with an implementation-
defined count-field size.  Our library defaults to cell-size.

In stack comments, a $ suffix indicates the address of an
mstring.
)

: m$>s  ( m$ -- s )  \ COUNT for cell count-field
  >r r@ cell+ ( addr) r> @  ( len) ;

: s-m$+  ( s m$ -- )
(
Append the Forth string to the mstring.  Transliteration of Wil
Baden's APPEND.
)
  2dup 2>r   m$>s +   swap move 2r> +! ;

hex
0000FFFF constant lo16-mask
FFFF0000 constant hi16-mask
decimal

create mix$ 128 chars allot

: 0mix$  ( -- )    0 mix$ ! ;
: mix$>s ( -- )    mix$ m$>s ;
: mix+   ( s -- )  mix$ s-m$+ ;
: bl+    ( -- )    s"  " mix+ ;
: num+   ( n -- )  num>s mix+ ;

: fsign+  ( fraw.hi32 -- )
  fsign-mask and IF s" -" ELSE s" +" THEN mix+ ;

: mix+>s  ( s -- mix.s )  mix+ mix$>s ;

: "Inf"+>s  ( -- )  s" Inf" mix+>s ;
: "NaN"+>s  ( -- )  s" NaN" mix+>s ;

: digs+  ( +num #digs -- )
  ( #digs) >r num>s
  r> over ( len) - dup 0< ABORT" ***Too many digits."
  ( #digs-len) 0 ?DO s" 0" mix+ LOOP
  ( num.s) mix+ ;

: bl-4digs+  ( +4.digit.num -- )  bl+ 4 digs+ ;

\ MIXED OUTPUT SPECIFICATION
(
There is a different mixed output format for each of the raw
formats, the point being to display the mantissa with hex digits
that correspond directly to the raw bits.

RAW32>MIX
normal:     <"+"|"-">.hh hhhh p<"+"|"-">d..d
subnormal:  <"+"|"-">.hh hhhh p-125
leading bit of leading "h" is 1/0 for normal/subnormal numbers

RAW64>MIX
normal:     <"+"|"-">1.h hhhh hhhh hhhh p<"+"|"-">d..d
subnormal:  <"+"|"-">0.h hhhh hhhh hhhh p-1022

RAW80>MIX
normal:     <"+"|"-">.hhhh hhhh hhhh hhhh p<"+"|"-">d..d
subnormal:  <"+"|"-">.hhhh hhhh hhhh hhhh p-16381
leading bit of leading "h" is 1/0 for normal/subnormal numbers

RAW128>MIX
normal:    <"+"|"-">1.hhhh hhhh hhhh hhhh hhhh hhhh hhhh p<"+"|"-">d..d
subnormal: <"+"|"-">0.hhhh hhhh hhhh hhhh hhhh hhhh hhhh p-16382

Otherwise the output is "+Inf", "-Inf", "+NaN", or "-NaN"
ignoring the arbitrary bits of a raw NaN.
)

\ RAW32-FNUM+, RAW6-FNUM+, RAW128-FNUM+, RAW80-FNUM+  
(
These words are the main factors in the corresponding raw format
to mixed output words.  They assume that HEX$ is already the
sign of the mantissa, which they ignore, and that the raw format
corresponds to a normal or subnormal number.  They convert and
append the mantissa and exponent to HEX$ as a mixed hex/decimal
string, using the mixed output specification corresponding to
the raw format.

The implementations convert NaN and infinity to something
unspecified.  They convert zero correctly, but that is neither
used nor guaranteed.
)

: raw32-fnum+  ( raw32 -- )
  base @ >r

  s" ." mix+
  dup sexp-mask and >r

  hex
  ( raw) dup smantissa-mask and lo16-mask xor 16 rshift
  r@ ( sexp) IF 128 ( 0x80) or THEN
  2 digs+
  ( raw) lo16-mask and bl-4digs+

  decimal
  s"  p" mix+
  r> ( sexp)
  dup 0= IF ( 0) drop 1 ELSE 23 rshift THEN
  126 - dup 0>= IF s" +" mix+ THEN num+

  r> base ! ;

: raw64-fnum+  ( raw64 -- )
  base @ >r

  dup dexp-mask and IF s" 1." ELSE s" 0." THEN mix+

  hex  \ mantissa after "1." or "0."

  \ single digit
  ( raw.hi32) dup dmantissa-mask and
      dup lo16-mask xor 16 rshift num+

  \ remaining 4 digits of high mantissa
  ( raw.hi32) lo16-mask and bl-4digs+

  \ high 4 digits of low mantissa
  swap ( mantissa.low32) dup hi16-mask and
      16 rshift bl-4digs+

  \ low 4 digits of low mantissa
  ( mantissa.low32) lo16-mask and bl-4digs+

  decimal  \ binary exponent

  s"  p" mix+
  ( raw.hi32) dexp-mask and 20 rshift
  dup 0= IF 1+ THEN  \ adjust when subnormal
  1023 - dup 0>= IF s" +" mix+ THEN num+

  r> base ! ;

: raw128-fnum+  ( raw128 -- )
  base @ >r

  ( raw.hi32) >r
  r@ qexp-mask and IF s" 1." ELSE s" 0." THEN mix+

  hex  \ mantissa frac

  \ frac.hi digits
  r@ ( raw.hi32)
  qmantissa-mask and 4 digs+

  3 0 DO
    dup hi16-mask and 16 rshift bl-4digs+
    lo16-mask and bl-4digs+
  LOOP

  decimal  \ binary exponent

  s"  p" mix+
  r> ( raw.hi32) qexp-mask and 16 rshift
  dup 0= IF 1+ THEN  \ adjust when subnormal
  16383 - dup 0>= IF s" +" mix+ THEN num+

  r> base ! ;

: raw80-fnum+  ( raw80 -- )
  base @ >r

  s" ." mix+

  hex  \ mantissa

  ( raw.hi32) >r
  r@ ( raw.hi32)
  qmantissa-mask and 4 digs+
  dup hi16-mask and 16 rshift bl-4digs+
  lo16-mask and bl-4digs+
  hi16-mask and 16 rshift bl-4digs+

  decimal  \ binary exponent

  s"  p" mix+
  r> ( raw.hi32) qexp-mask and 16 rshift
  dup 0= IF 1+ THEN  \ adjust for subnormal
  1+                 \ adjust for leading decimal point
  16383 - dup 0>= IF s" +" mix+ THEN num+

  r> base ! ;

\ RAW32>MIX, RAW64>MIX, RAW128>MIX, RAW80>MIX  
(
Convert the indicated raw formats to a mixed hex/decimal string
according to the mixed output specification.
)

: raw32>mix  ( raw32 -- mix.s )
  0mix$
  dup fsign+

  fbody-mask and Inf-raw32 over <
  IF drop "NaN"+>s EXIT THEN

  dup Inf-raw32 =
  IF drop "Inf"+>s EXIT THEN

  dup 0=
  IF drop s" .00 0000 p+0" mix+>s EXIT THEN

  raw32-fnum+ mix$>s ;

: raw64>mix  ( raw64 -- mix.s )
  0mix$
  dup fsign+

  fbody-mask and Inf-raw64 2over d<
  IF 2drop "NaN"+>s EXIT THEN

  2dup Inf-raw64 d=
  IF 2drop "Inf"+>s EXIT THEN

  2dup d0=
  IF 2drop s" 0.0 0000 0000 0000 p+0" mix+>s EXIT THEN

  raw64-fnum+ mix$>s ;

: raw128>mix  ( raw128 -- mix.s )
  0mix$
  dup fsign+ fbody-mask and

  2over d0= >r
  2dup Inf-raw128.2hi d= r> and
  IF 2drop 2drop "Inf"+>s EXIT THEN

  dup qexp-mask and qexp-mask =
  IF 2drop 2drop "NaN"+>s EXIT THEN 

  2dup d0= >r
  2over d0= r> and
  IF 2drop 2drop s" 0.0000 0000 0000 0000 0000 0000 0000 p+0"
     mix+>s EXIT THEN

  raw128-fnum+ mix$>s ;

: raw80>mix  ( raw80 -- mix.s )
  0mix$
  dup fsign+ fbody-mask and

  2 pick 0= >r
  2dup Inf-raw80.2hi d= r> and
  IF 2drop drop "Inf"+>s EXIT THEN

  dup qexp-mask and qexp-mask =
  IF 2drop drop "NaN"+>s EXIT THEN 
  
  2 pick 0= >r
  2dup d0= r> and
  IF 2drop drop s" .0000 0000 0000 0000 p+0" mix+>s EXIT THEN

  raw80-fnum+ mix$>s ;


\ *** HEX/RAW CONVERSIONS

\ HEX8>RAW32, HEX16>RAW64, HEX20>RAW80, HEX32>RAW128
(
Translate a Forth string representation of a raw floating-point
datum to the corresponding raw format on the data stack.  These
words are for convenience in expressing raw data without concern
for right to left stack notation, and with embedded blanks for
readability.  They all ABORT if the string does not contain
exactly the required number of hex digits.
)

: hex8>raw32  ( hex8.s -- raw32 )
  -spaces dup 8 <> ABORT" ***Not 8 hex digits."
  ?hex>ud ( hi) drop ;

: hex16>raw64  ( hex16.s -- raw64 )
  -spaces dup 16 <> ABORT" ***Not 16 hex digits."
  ?hex>ud ;

: hex20>raw80  ( hex20.s -- raw80 )
  -spaces dup 20 <> ABORT" ***Not 20 hex digits."
  ?hex>uq 16 (qlshift) drop ;

: hex32>raw128  ( hex32.s -- raw128 )
  -spaces dup 32 <> ABORT" ***Not 32 hex digits."
  ?hex>uq ;

\ RAW32>HEX8, RAW64>HEX16, RAW80>HEX20, RAW128>HEX32
(
Translate the raw formats to hex strings of hex digits formatted
in 4-digit, blank-separated groups.
)

: u>hex+  ( u -- )
(
Translate u to a string of 8 hex digits, in two 4-digit,
blank-separated groups, and append it to hex$.
)
  base @ >r hex dup hi16-mask and 16 rshift 4 digs+
  lo16-mask and bl-4digs+ r> base ! ;

: raw32>hex8   ( raw32 -- hex8.s )   0mix$ u>hex+ mix$>s ;
: raw64>hex16  ( raw64 -- hex16.s )  0mix$ u>hex+ bl+ u>hex+ mix$>s ;

: raw80>hex20  ( raw80 -- hex20.s )
  0mix$ u>hex+ bl+ u>hex+ bl+ u>hex+ mix$>s 5 - ;

: raw128>hex32  ( raw128 -- hex32.s )
  0mix$ u>hex+ 3 0 DO bl+ u>hex+ LOOP mix$>s ;


\ *** HEX/MIX CONVERSIONS

\ HEX8>MIX, HEX16>MIX, HEX20>MIX, HEX32>MIX,
\ MIX>HEX8, MIX>HEX16, MIX>HEX20, MIX>HEX32, 
(
Translate among the hex raw string and mixed string formats.
)

: hex8>mix   ( hex8.s -- mix.s )   hex8>raw32   raw32>mix ;
: hex16>mix  ( hex16.s -- mix.s )  hex16>raw64  raw64>mix ;
: hex20>mix  ( hex20.s -- mix.s )  hex20>raw80  raw80>mix ;
: hex32>mix  ( hex32.s -- mix.s )  hex32>raw128 raw128>mix ;

: mix>hex8   ( mix.s -- hex8.s )   mix>raw32  raw32>hex8 ;
: mix>hex16  ( mix.s -- hex16.s )  mix>raw64  raw64>hex16 ;
: mix>hex20  ( mix.s -- hex20.s )  mix>raw80  raw80>hex20 ;
: mix>hex32  ( mix.s -- hex32.s )  mix>raw128 raw128>hex32 ;


\ *** FLOATING-POINT ENVIRONMENT
(
This code will abort if a Forth char is not 8 bits, and if the
default float format is not 32, 64, 80, or 128 bits.  It further
assumes the default float format to be either totally big- or
little-endian in memory, and that when floats are little-endian,
so are cells.  The tests in mixfloat-test.fs should catch the
failure of the latter assumption.
)
decimal

bits/au 1 chars * 8 <> [IF] cr
  .( ***1 CHARS is not 8 bits.) ABORT [THEN]

\ Use this instead for pfe's integrated fp stack option:
\ s" FPNOSTACK-EXT"
s" FLOATING-EXT"
environment? [IF] ( flag) drop
[ELSE] cr .( ***Floating point not available.) ABORT [THEN]

1 floats 16 >
[IF] cr .( ***1 FLOATS is too big.) ABORT [THEN]

[UNDEFINED] qfpad [IF]
qfalign here 5 cells allot constant qfpad  \ 1 extra cell
[THEN]

MARKER -FIND-FLOAT-FORMAT
qfpad 4 cells 1+ char U fill

: #fbytes  ( -- #bytes )
  0 4 cells 1+ 0
  DO drop i qfpad i + c@ [char] U = IF LEAVE THEN LOOP ;

1E qfpad f! #fbytes ( #bytes) bits/au * ( bits/float)
[UNDEFINED] bits/float

qfpad c@ 0= ( lendian?)
[UNDEFINED] [LITTLE-ENDIAN]

-FIND-FLOAT-FORMAT ( bits/float undef? lendian? undef?)
[IF] constant [LITTLE-ENDIAN] immediate [ELSE] drop [THEN]
[IF] constant bits/float                [ELSE] drop [THEN]

bits/float 32 =     bits/float 64 = or
bits/float 80 = or  bits/float 128 = or
0= [IF] cr
  .( ***Unrecognized float format, bits/float = )
  bits/float . ABORT [THEN]


\ *** MORE MIX/RAW CONVERSIONS

: mix>raw  ( s -- raw.default )
  [ bits/float 32 = ]  [IF] mix>raw32 [ELSE]
  [ bits/float 64 = ]  [IF] mix>raw64 [ELSE]
  [ bits/float 80 = ]  [IF] mix>raw80 [ELSE]
  [ bits/float 128 = ] [IF] mix>raw128 [THEN] [THEN] [THEN] [THEN] ;

: raw>mix  ( raw.default -- mix.s )
  [ bits/float 32 = ]  [IF] raw32>mix [ELSE]
  [ bits/float 64 = ]  [IF] raw64>mix [ELSE]
  [ bits/float 80 = ]  [IF] raw80>mix [ELSE]
  [ bits/float 128 = ] [IF] raw128>mix [THEN] [THEN] [THEN] [THEN] ;


\ *** MORE HEX/RAW CONVERSIONS

: hex>raw  ( hex.s -- raw.default )
  [ bits/float 32 = ]  [IF] hex8>raw32 [ELSE]
  [ bits/float 64 = ]  [IF] hex16>raw64 [ELSE]
  [ bits/float 80 = ]  [IF] hex20>raw80 [ELSE]
  [ bits/float 128 = ] [IF] hex64>raw128 [THEN] [THEN] [THEN] [THEN] ;

: raw>hex  ( raw64.default -- hex.s )
  [ bits/float 32 = ]  [IF] raw32>hex8 [ELSE]
  [ bits/float 64 = ]  [IF] raw64>hex16 [ELSE]
  [ bits/float 80 = ]  [IF] raw80>hex20 [ELSE]
  [ bits/float 128 = ] [IF] raw128>hex64 [THEN] [THEN] [THEN] [THEN] ;


\ *** MORE HEX/MIX CONVERSIONS

: hex>mix  ( hex.s -- mix.s )  hex>raw raw>mix ;
: mix>hex  ( mix.s -- hex.s )  mix>raw raw>hex ;


[UNDEFINED] raw! [IF] \ BEGIN rawfloat.fs not loaded
\ *** RAW FETCH/STORE

: raw32!  ( raw32 sf-addr -- )  ! ;
: raw64!  ( raw64 df-addr -- )  [LITTLE-ENDIAN] [IF] >r swap r> [THEN] 2! ;

: raw80!  ( raw80 xf-addr -- )
  dup 10 + [LITTLE-ENDIAN]
  [IF]   1-   DO 0 8 (qlshift) i c! -1 +LOOP
  [ELSE] swap DO 0 8 (qlshift) i c!     LOOP [THEN] drop 2drop ;

: raw128!  ( raw128 qf-addr -- )
  >r [LITTLE-ENDIAN]
  [IF]   2swap swap r@ 2! swap r> 8 + 2!
  [ELSE] r@ 2! r> 8 + 2! [THEN] ;

: raw32@  ( sf-addr -- raw32 )  @ ;
: raw64@  ( df-addr -- raw64 )  2@ [LITTLE-ENDIAN] [IF] swap [THEN] ;

: raw80@  ( xf-addr -- raw80 )
  >r 0 0 0 r>
  dup 10 + [LITTLE-ENDIAN] 
  [IF]   swap DO i c@ 8 (qrshift) drop     LOOP
  [ELSE] 1-   DO i c@ 8 (qrshift) drop -1 +LOOP [THEN] ;

: raw128@  ( qf-addr -- raw128 )
  >r r@ 2@ [LITTLE-ENDIAN]
  [IF]   swap r> 8 + 2@ swap
  [ELSE] r> 8 + 2@ 2swap [THEN] ;

: raw!  ( raw f-addr -- )
  [ bits/float 32 = ]  [IF] raw32! [ELSE]
  [ bits/float 64 = ]  [IF] raw64! [ELSE]
  [ bits/float 80 = ]  [IF] raw80! [ELSE]
  [ bits/float 128 = ] [IF] raw128! [THEN] [THEN] [THEN] [THEN] ;

: raw@  ( f-addr -- raw )
  [ bits/float 32 = ]  [IF] raw32@ [ELSE]
  [ bits/float 64 = ]  [IF] raw64@ [ELSE]
  [ bits/float 80 = ]  [IF] raw80@ [ELSE]
  [ bits/float 128 = ] [IF] raw128@ [THEN] [THEN] [THEN] [THEN] ;


\ *** RAW/FLOAT CONVERSIONS

: raw32>f  ( raw32 -- f: r )      qfpad raw32! qfpad sf@ ;
: raw64>f  ( raw64 -- f: r )      qfpad raw64! qfpad df@ ;
: f>raw32  ( f: r -- s: raw32 )   qfpad sf! qfpad raw32@ ;
: f>raw64  ( f: r -- s: raw64 )   qfpad df! qfpad raw64@ ;

[DEFINED] XF@ [IF]  \ assume XF! defined, too
: raw80>f  ( raw80 -- f: r )      qfpad raw80! qfpad xf@ ;
: f>raw80  ( f: r -- s: raw80 )   qfpad xf! qfpad raw80@ ;
[THEN]

[DEFINED] QF@ [IF]  \ assume QF! defined, too
: raw128>f  ( raw128 -- f: r )    qfpad raw128! qfpad qf@ ;
: f>raw128  ( f: r -- s: raw128 ) qfpad qf! qfpad raw128@ ;
[THEN]

: raw>f  ( raw.default -- f: r )
  [ bits/float  32 = ]  [IF] raw32>f [ELSE]
  [ bits/float  64 = ]  [IF] raw64>f [ELSE]
  [ bits/float  80 = ]  [IF] raw80>f [ELSE]
  [ bits/float  128 = ] [IF] raw128>f [THEN] [THEN] [THEN] [THEN] ;

: f>raw  ( f: r -- s: raw.default )
  [ bits/float  32 = ]  [IF] f>raw32 [ELSE]
  [ bits/float  64 = ]  [IF] f>raw64 [ELSE]
  [ bits/float  80 = ]  [IF] f>raw80 [ELSE]
  [ bits/float  128 = ] [IF] f>raw128 [THEN] [THEN] [THEN] [THEN] ;

[THEN]  \ END rawfloat.fs not loaded


\ *** HEX/FLOAT CONVERSIONS

: hex8>f   ( hex8.s -- f: r )     hex8>raw32  qfpad raw32! qfpad sf@ ;
: hex16>f  ( hex16.s -- f: r )    hex16>raw64 qfpad raw64! qfpad df@ ;
: f>hex8   ( f: r -- s: hex8.s )  qfpad sf! qfpad raw32@ raw32>hex8 ;
: f>hex16  ( f: r -- s: hex16.s ) qfpad df! qfpad raw64@ raw64>hex16 ;

[DEFINED] XF@ [IF]  \ assume XF! defined, too
: hex20>f  ( hex20.s -- f: r )    hex20>raw80 qfpad raw80! qfpad xf@ ;
: f>hex20  ( f: r -- s: hex20.s ) qfpad xf! qfpad raw80@ raw80>hex20 ;
[THEN]

[DEFINED] QF@ [IF]  \ assume QF! defined, too
: hex32>f  ( hex32.s -- f: r )    hex32>raw128 qfpad raw128! qfpad qf@ ;
: f>hex32  ( f: r -- s: hex32.s ) qfpad qf! qfpad raw128@ raw128>hex32 ;
[THEN]

: hex>f  ( hex.default.s -- f: r )
  [ bits/float  32 = ]  [IF] hex8>f [ELSE]
  [ bits/float  64 = ]  [IF] hex16>f [ELSE]
  [ bits/float  80 = ]  [IF] hex20>f [ELSE]
  [ bits/float  128 = ] [IF] hex32>f [THEN] [THEN] [THEN] [THEN] ;


: f>hex  ( f: r -- s: hex.default.s )
  [ bits/float  32 = ]  [IF] f>hex8 [ELSE]
  [ bits/float  64 = ]  [IF] f>hex16 [ELSE]
  [ bits/float  80 = ]  [IF] f>hex20 [ELSE]
  [ bits/float  128 = ] [IF] f>hex32 [THEN] [THEN] [THEN] [THEN] ;


\ *** MIX/FLOAT CONVERSIONS

: mix>f  ( mix.s -- f: r )     mix>raw raw>f ;
: f>mix  ( f: r -- s: mix.s )  f>raw raw>mix ;

: fm.  ( f: r -- )  f>mix type ;