RfD: IEEE-FP 0.5.2
David N. Williams
old one has become unwieldy. That doesn't mean that ongoing
conversations in the old thread shouldn't continue there.
SUMMARY OF CHANGES
* "number" is no longer used for infinities and nans -- "IEEE
datum" is used instead.
* no nan load i/o
* two environmental queries instead of one
* binary format parameter definition moved to a new section,
details moved to an informative appendix
* the default format (fp stack format) can be any IEEE binary
format, and nothing in the document depends on encoding specifics
* nans and infinities removed from decimal text
interpretation, and added as Forth constants, with the bnf
restored to *nearly* the same as DPANS94
* hexadecimal input removed
* >FLOAT alternatives suggested
* F>D alternatives suggested
* IEEE comparisons reworked
* F~ reworked
* exceptions and rounding put on hold pending discussion
As for exceptions and rounding, it doesn't seem to me
controversial that the basic IEEE functionality should be
implemented, but ideas for how haven't gelled yet.
-- David
-----------------------------------------
PROPOSAL FOR AN OPTIONAL IEEE 754, BINARY
FLOATING-POINT WORD SET
version 0.5.2
dnw 07-Jul-09
TABLE of CONTENTS
1 INTRODUCTION
2 TERMINOLOGY AND NOTATION
2a IEEE BINARY FLOATING-POINT FORMATS
3 IMPLEMENTATION
4 DATA TYPES
5 ENVIRONMENTAL QUERIES
6 TEXT INPUT
6.1 Constants
6.2 Decimal input
6.3 Hexadecimal input
7 GLOSSARY
7.1 Conversion
7.2 Output
7.3 Comparison
7.4 Classification
7.5 Arithmetic
7.6 Math functions
7.7 Sign bit operations
7.8 Nearest integer functions
7.9 Number manipulation
7.10 Exceptions
7.11 Rounding modes
8 REFERENCES
A.2a IEEE FLOATING-POINT FORMATS
A.6.1 NaN signs and loads
A.7.3 Comparison (informative rationale)
1 INTRODUCTION
This is a proposal for an optional Forth 200x word set, called
the "IEEE floating-point word set", that supports the binary
part of the IEEE 754-2008 standard for floating-point arithmetic
[1]. The most recent, freely available, but less comprehensive
version is IEEE 754 draft 1.2.9, January 27, 2007 [2]. There is
also a Wikipedia summary [3].
The standard [1] is hereafter referred to as "IEEE 754-2008",
with section numbers indicated by IEEE <section number>.
This specification requires that ISO Forth [4,5] floating-point
and floating-point extension words in the optional
floating-point word set, when present *with the IEEE
floating-point word set*, satisfy additional IEEE 754-2008
requirements. Words in that word set and this that correspond
to mathematical, including logical, operations or functions in
IEEE 754-2008 adopt the behavior required or recommended there
by reference, as far as that is possible and makes sense, unless
otherwise stated.
The specification is compatible with, rather than conformant to,
IEEE 754-2008, because it includes only a subset of the IEEE
requirements.
Reference [4], the final draft of "ANSI X3.215-1994, American
National Standard for Information Systems--Programming
Languages--Forth", is hereafter referred to as "DPANS94". It is
believed to be the same as the published version, ISO/IEC
15145:1997 [5]. This document adopts the official terminology
of DPANS94 unless otherwise stated. Section numbers in that
document are indicated by DPANS94 <section number>.
When it refers to the IEEE floating-point word set, the term
"required" in this document is to be understood in the context
of the following two paragraphs from DPANS94 A.1.3.1, which
discuss the meaning of "optional word set":
The basic requirement is that if the implementor claims to
have a particular optional word set the entire required
portion of that word set must be available. If the
implementor wishes to offer only part of an optional word set,
it is acceptable to say, for example, "This system offers
portions of the [named] word set", particularly if the
selected or excluded words are itemized clearly.
and
Optional word sets may be offered in source form or otherwise
factored so that the user may selectively load them.
The current C99 standard [6-8], ISO/IEC 9899:1999, has a
comprehensive treatment of IEEE 754-1985, which offers a route
to implementation, for those Forth systems that can call C
libraries. Section numbers from reference [7] are indicated by
C99 <section number>.
[**Bracketed statements like this are for editorial questions
and comments, eventually to be removed.]
2 TERMINOLOGY AND NOTATION
"fp" or "bfp": Short for "binary floating point". The Forth
floating-point stack is called the "fp stack".
"IEEE special datum", or an "IEEE special": Signed zero, a
quiet or signaling signed nan, or signed infinity.
"full IEEE set": For an IEEE binary format, the set of normal
and subnormal numbers plus special data that it represents.
"IEEE datum": Any member of a full IEEE set.
"IEEE arithmetic": Arithmetic defined by IEEE 754-2008 for IEEE
data.
"affinely extended reals": Finite real numbers and +/-infinity,
with -infinity < {every finite number} < +infinity.
"nan load" or "nan payload": The value of the fractional bits
in the binary format of a nan, excluding the quiet bit,
considered as a positive integer. The smallest signaling load
is unity, and the smallest quiet load is zero.
"qnan", resp., "snan": A quiet or signaling nan, respectively,
of any sign or load.
"single": In the context of Forth fp, an IEEE 754-2008 32-bit
interchange format.
"double": In the context of Forth fp, an IEEE 754-2008 64-bit
interchange format.
"default": In the context of Forth fp, the float format for
data that can appear on the fp stack.
"exception": Used in the sense of IEEE 2.1.18.
[**QUOTE
An event that occurs when an operation on some particular
operands has no outcome suitable for every reasonable
application. That operation might signal one or more
exceptions by invoking the default or, if explicitly
requested, a language-defined alternate handling. Note that
"event", "exception", and "signal" are defined in diverse
ways in different programming environments.
]
2a IEEE BINARY FLOATING-POINT FORMATS
Each IEEE bfp format has two fixed parameters, p > 0 (precision)
and emax > 0 (maximum exponent), and defines emin = 1 - emax
(minimum exponent). Each such format represents all real
numbers of the form
r = (-1)^s * 2^e * b_0.b_1 ... b_{p-1}
where
s = 0 or 1, emin <= e <= emax,
b_i = 0 or 1, p = #significand bits.
See Sec. A.2a for more information about IEEE bfp formats.
3 IMPLEMENTATION
According to DPANS94, Section 3, "Usage requirements":
A system shall provide all of the words defined in 6.1 Core
Words. It may also provide any words defined in the optional
word sets and extension word sets.
The DPANS94 Floating-Point word set is an optional word set,
and so is the word set described by this document.
The word "shall" in the remainder of this document states a
requirement when the environmental query for IEEE-FP returns
true. "Should" means "strongly recommended".
The internal fp representation of default fp data, i.e., data
that can appear on the fp stack, shall correspond to one of the
IEEE basic, or extended, full binary formats.
4 DATA TYPES
For the purpose of this document, the DPANS94 r type is extended
to include all IEEEE data for the default fp format.
5 ENVIRONMENTAL QUERIES
[**CHANGES
There are two queries instead of one. It seemed essential that
there be a query returning the fp format parameters, even when
not all words are present.
]
Value
String Data Type Constant? Meaning
------------------------------------------------------------------
IEEE-FP flag no IEEE and DPANS94 floating-point word
sets present
IEEE-FP-FORMAT d no in usual stack notation, the default
format has IEEE parameters ( emax p )
------------------------------------------------------------------
A true result for the IEEE-FP environmental query (not the data
value) shall mean that any words that are present from the
DPANS94 floating-point word set, the DPANS94 floating-point
extensions word set, or the IEEE floating-point word set shall
obey the specifications of this document.
A false data value means that only a subset of the IEEE and
DPANS94 word sets is present.
A true result for the IEEE-FP-FORMAT query shall mean that the
DPANS94 MAX-FLOAT query shall return true and the largest,
finite IEEE number in the default format.
Nothing in this document depends on the encoding of the format
corresponding to emax and p.
6 TEXT INPUT
6.1 Constants
-------------
+INF ( f: -- +Inf )
-INF ( f: -- -Inf )
+NAN ( f: -- +NaN )
-NAN ( f: -- -NaN )
These words return, respectively, IEEE signed infinity and the quiet
signed nan with zero load, in the default format.
See Sec. A.6.1 for more information about the encoding of NaN.
6.2 Decimal input
-----------------
IEEE requires that conversion between text and binary fp formats
shall include signed zero, signed infinity, and signed nans,
with and without loads. See IEEE 5.4.2, "Conversion operations
for floating-point formats and decimal character sequences", and
IEEE 5.12, "Details of conversion between floating-point data
and external character sequences".
Conversion of nan loads is not included in this specification.
Signed infinity and signed, unloaded nans are covered by the
constants defined in Sec. 6.1. Signed zero is already included
in the syntax specification in DPANS94 12.3.7, "Text input
number conversion". When IEEE-FP is present, that specification
shall be replaced by
Convertible string := <significand><exponent>
<significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
[**Note the extra .<digits> option.]
<exponent> := <e-char>[<sign>]<digits0>
[**Note that a sign with no digits is still
recognized, as in DPANS94.]
<digits> := <digit><digits0>
<digits0> := <digit>*
<sign> := { + | - }
<e-char> := { E | e }
[**Note the extra "e" option.]
<digit> := { 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 }
Interpretation shall convert to the default IEEE format, with
roundTiesToEven for real numbers, and subject to overflow and
underflow exceptions. See IEEE 7.4 for overflow and IEEE 7.5
for underflow.
6.3 Hexadecimal input
---------------------
IEEE requires a text format for numbers with a hexadecimal
significand, and decimal radix two exponent, with exact
conversion to and from binary fp formats where possible. See
IEEE 5.12.3, "External hexadecimal-significand character
sequences representing finite numbers".
Conversion of that format is not included in this specification.
7 GLOSSARY
Unless otherwise stated, all fp words that do computations or
comparisons shall obey the requirements and recommendations of
IEEE 5 and IEEE 6, for binary formats.
7.1 Conversion
--------------
D>F ( d -- f: r )
The DPANS94 specification is amended to require that when d
cannot be represented precisely in the default fp format, r
shall be the roundTiesToEven value.
[**IEEE and C99 background:
IEEE does not seem to specify an equivalent. C99-n1256 6.3.1.4,
paragraph 2 says:
When a value of integer type is converted to a real floating
type, if the value being converted can be represented exactly
in the new type, it is unchanged. If the value being
converted is in the range of values that can be represented
but cannot be represented exactly, the result is either the
nearest higher or nearest lower representable value, chosen in
an implementation-defined manner. If the value being
converted is outside the range of values that can be
represented, the behavior is undefined.
Note that d is always in range for any of the formats
binary32, binary64, binary80, and binary128.
]
F>D ( f: r -- s: d )
[**Suggested specificaton #1:
The DPANS94 specification is amended to state that an
ambiguous condition exists, not only when the integer part of
r is not representable by a signed, double number, but also
when r is a nan or infinity.
]
[**Suggested specificaton #2:
The DPANS94 specification is amended to require that a Forth
exception be thrown when r is a nan or infinity, or has
integer part out of range for a signed, double number.
]
[**Suggested specificaton #3:
The DPANS94 specification is amended to require that the
invalid operation exception be quietly signaled when r is a
nan or infinity, or has integer part out of range for a
signed, double number. In that case the value of d is
undefined.
]
[**IEEE background:
The DPANS94 version corresponds to
convertToIntegerTowardZero(); see IEEE 5.8, "Details of
conversions from floating-point to integer formats". The IEEE
version requires that the invalid operation exception be
signaled when r is a nan or infinity or out of range of the
destination format.
IEEE also requires
convertToIntegerTiesToEven()
convertToIntegerTowardPositive()
convertToIntegerTowardNegative()
convertToIntegerTiesToAway()
plus versions of the five conversions that signal inexact when
appropriate.
]
[**Candidate for rationale section:
Note that some of the other conversions in IEEE 5.8, "Details
of conversions from floating-point to integer formats", are
equivalent to Forth phrases such as "FROUND F>D".
]
[** ALTERNATIVE 1:
>FLOAT ( c-addr u -- [r: r s: true]|[false] )
The DPANS94 specification is extended to include IEEE
specials, with modified syntax. Conversion shall be governed
by IEEE 5.12, "Details of conversion between floating-point
data and external character sequences", except for paragraph 4
on hexadeximal conversions. If the string represents a valid
IEEE datum in the syntax below the datum r and true are
returned. Otherwise only false is returned.
A string of blanks [**or the empty string?] shall be treated
as +0.
Syntax of a convertible string := { <significand>[exponent]
| <special> }
<significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
<exponent> := <marker><digits0>
<marker> := { <e-form> | <sign> }
<e-form> := <e-char>[<sign>]
<sign> := { + | - }
<e-char> := { D | d | E | e }
<special> := [<sign>]{ <inf> | <nan> }
<inf> := { Inf | inf | INF | infinity | Infinity }
<nan> := { NaN | nan | NAN }
[**mh prefers a mandatory <sign> in <special>. dnw favors
keeping it optional, in line with the omnivorous philosophy of
DPANS94 A.12.6.1.0558.]
]
[**ALTERNATIVE 2:
>FLOAT ( c-addr u -- [r: r s: true]|[false] )
The expression "the string represents a valid floating-point
number" in the DPANS94 specification shall be interpreted to
mean that it represents a finite number in the range of the
default format.
>IEEE-FLOAT ( c-addr u -- [r: r s: true]|[false] )
The specification is exactly that of >FLOAT in alternative 1.
]
[**I prefer alternative 1. I mean, really, a "float" in the
IEEE context is any IEEE default datum. But I would accept a
variant of 2 if 1 should prove to be a show stopper.]
REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
The DPANS94 specification is extended to include IEEE
specials. The valid-result flag, flag2, is false if r is
infinity or a nan.
[**IMHO this is still unsettled, still needing discussion.
1. Should the current rounding mode be used, instead of
roundTiesToEven? I'll have to check to what extent that
has already been discussed.
2. Marcel says "It should be possible to format/output +Inf or
+NaN, also with the existing REPRESENT." The DPANS94 spec
for false flag2, that n and flag1 are implementation
defined, does seem (perhaps intentionally?) well-suited to
the spec that flag2 continues to be the negative sign flag,
while n distinguishes among nans and infinity. For
example, n = 0 means infinity, and n <> 0 means nan with
implementation defined values related to quietness and
load.
]
SF! ( f: r s: sf-addr -- )
SF@ ( sf-addr -- f: f )
DF! ( f: r s: df-addr -- )
DF@ ( sf-addr -- f: f )
The specification for these DPANS94 words is amended to
explicitly require conversion to or from the respective IEEE
binary32 or binary64 interchange formats, with exact
conversion in either direction for signed zero, signed
infinity, and real numbers to a wider format, and with
the current rounding mode for conversion of real numbers to a
narrower format (see IEEE Sec. 5.4.2, formatOf-convertFormat).
The conversion of nans shall preserve the sign and signaling
bit, shall not signal an exception, and should treat payloads
according to IEEE 6.2.3.
The addresses sf-addr and df-addr are aligned according to
SFALIGNED and DFALIGNED. The memory representation of the
corresponding binary interchange formats is implementation
defined [**AFAIK that's the common interpretation of IEEE].
[**CHANGE: current rounding mode is used instead of the
former roundTiesToEven.]
[**I think an RfD for XF@ and XF!, for binary80, is in order.
And perhaps names like QF@ and QF!, for binary128, should be
reserved.]
7.2 Output
----------
F. ( f: r -- )
FE. ( f: r -- )
FS. ( f: r -- )
The DPANS94 specification is extended to include IEEE
specials, with output text of the appropriate form below, with
implementation-dependent case sensitivity:
[**Default rounding? Or use the current rounding mode?]
[<sign>]0{ E | e }<space>
[<sign>]{ Inf | INF | inf }<space>
[<sign>]{ NaN | NAN | nan }<space>
[**CHANGE: Removed nan load.]
7.3 Comparison
--------------
IEEE has twenty-two required comparisons which apply to the full
set of IEEE data. Twelve of these are quiet, and ten are
signaling. See IEEE 5.6.1, "Comparisons", and IEEE 5.11,
"Details of comparison predicates".
This proposal requires only quiet comparisons, which do not
raise exceptions, and of those, only a subset of five, which is
sufficient for expressing all twelve.
See Sec. A.7 for rationale and more information about the
remaining comparisons, with high-level implementation examples.
IEEE identifies four fundamental, mutually exclusive
comparisons: less than ("<"), equal ("="), greater than (">"),
and unordered (see rule 3 below). Each of these is true iff
each of the others is false.
The basic rules are the following:
1. The sign of zero is ignored.
2. The sign of infinity is not ignored, and is treated in the
natural way for the "ordinary" comparisons with real numbers
or infinity, namely <, =, and >. In particular, either
signed infinity is equal to itself.
3. The unordered comparison is true iff at least one of its two
arguments is a nan. That implies that any of the other
three, "ordinary" comparisons involving a nan is false.
The five required comparisons are "<", ">", "=", "<=", and ">=",
where "<=" and ">=" stand for the usual phrases "less than or
equal" and "greater than or equal". Note that familiar
identities for real numbers are generally not satisfied by IEEE
comparisons. For example, the negation of "<" is not the same
as ">=". See Sec. A.7.
F< ( f: r1 r2 -- s: [r1<r2]? )
F= ( f: r1 r2 -- s: [r1=r2]? )
F> ( f: r1 r2 -- s: [r1>r2]? )
F<= ( f: r1 r2 -- s: [r1<=r2]? )
F>= ( f: r1 r2 -- s: [r1>=r2]? )
F0< ( f: r1 r2 -- s: [r<0]? )
F0= ( f: r1 r2 -- s: [r=0]? )
F0> ( f: r1 r2 -- s: [r>0]? )
F0<= ( f: r1 r2 -- s: [r<=0]? )
F0>= ( f: r1 r2 -- s: [r>=0]? )
The data stack outputs are DPANS94 flags corresponding to the
indicated IEEE predicates. In particular, the specifications
for the existing DPANS94 words F<, F0<, and F0= are extended
to include IEEE specials.
[**
All of the F0< family of words could be replaced by simple
phrases like "0E F<", so maybe we shouldn't have them at all.
An argument for keeping them is that the principle of least
surprise demands the extension of the already exsisting DPANS94
F0= to IEEE, just like F<. Given that, the same principle
demands that the names and meanings of the other five be
reserved.
An alternative might be to declare ambiguous conditions for F0=
and F0< evaluated at IEEE specials, and remove all five from the
IEEE-FP word set in favor of the corresponding phrases. That
violates the principle of least surprise, because F< has no
ambiguity.
]
F~ ( f: r1 r2 r3 -- s: flag )
If the r3 is positive and not a nan or zero, flag is true iff
the absolute value of r1 minus r2 is less than r3, taking into
account IEEE arithmetic and comparison rules. [**+Inf is
regarded as positive.]
If r3 is signed zero, flag is true iff r1 and r2 have
identical formats.
If r3 is negative and not a nan or zero, flag is true iff the
absolute value of r1 minus r2 is less than the absolute value
of r3 times the sum of the absolute values of r1 and r2,
taking into account IEEE arithmetic and comparison rules.
[**-Inf is regarded as negative.]
If r3 is a nan, flag is false.
[**
There is a lot of feeling that F~ is a nasty word. It has
three distinct functions:
1. Do a comparison of the form:
|r1 - r2| < |r3|
2. Do a comparison of the form:
|r1 - r2| < |r3| * (|r1| + |r2|)
3. Test whether r1 and r2 have identical fp encodings.
In the IEEE context, the specification above that flag be
false when r3 is a nan is natural because of comparisons 1 and 2.
In the IEEE context, something equivalent to comparison 3 is a
convenience for comparing to signed zero, since F= ignores the
sign of zero. It can be used to distinguish nan types and
loads, although better tools might be designed for that.
Basically, the above specification for r3 equal to signed zero
seems the natural extension of the DPANS94 spec.
Any alternative would have to specify some other behavior when
any of r1, r2, or r3 is and IEEE special.
Either way, IMHO, a separate RfD should be submitted to
declare F~ obsolescent and implement its functionality with
two or three new words.
]
7.4 Classification
------------------
IEEE 5.7.2, "General operations", requires a large number of
classification operations. This documents defines only those
corresponding to:
isSignMinus
isNormal
isFinite
isZero
isSubnormal
isInfinite
isNaN
isSignaling
Actually isSignMinus corresponds to FSIGNBIT, and isZero
corresponds to F0=, which leaves the following:
FINITE? ( r: r -- s: [normal|subnormal]? )
FNORMAL? ( r: r -- s: normal? )
FSUBNORMAL? ( r: r -- s: subnormal? )
FINFINITE? ( r: r -- s: [+|-]Inf? )
FNAN? ( r: r -- s: nan? )
FSIGNALING? ( r: r -- s: snan? )
7.5 Arithmetic
--------------
See IEEE 5.4.1, "Arithmetic operations".
F* ( f: r1 r2 -- r1*r2 )
F*+ ( f: r1 r2 r3 -- [r2*r3]+r1 ) [**new]
F+ ( f: r1 r2 -- r1+r2 )
F- ( f: r1 r2 -- r1-r2 )
F/ ( f: r1 r2 -- r1/r2 )
FSQRT ( f: r -- sqrt[r] )
The DPANS94 specification is extended to IEEE arithmetic. See
IEEE 5.1, "Overview" for precision, rounding, special data
treatment, and exceptions. See IEEE 5.4.1, "Arithmetic
operations", for the arithmetic words.
[** F*+ is the Forth name for the IEEE required
fusedMultiplyAdd operation.]
[**Note that IEEE requires these to be "correctly rounded".]
7.6 Math functions
-------------------
The Forth words FABS, FMAX, FMIN, and FSQRT are covered
elsewhere.
The DPANS94 specification for the following words is extended to
adopt the corresponding IEEE behavior. See IEEE 9.2, "Recommended
correctly rounded functions", and 9.2.1, "Special values".
F** FACOS FACOSH FALOG FASIN FASINH FATAN FATAN2
FATANH FCOS FCOSH FEXP FEXPM1 FLN FLNP1 FLOG FSIN
FSINCOS FSINH FSQRT FTAN FTANH
IEEE recommends additional functions, whose recommended Forth
names would be:
FEXP2 FEXP2M1 FEXP10 FEXP10M1
FLOG2 FLOGP1 FLOG2P1 FHYPOT 1/FSQRT
FCOMPOUND FROOTN F**N |F|**
FSINPI FCOSPI FATANPI FATAN2PI
[**Separate RfD for the new names?]
7.7 Sign bit operations
-----------------------
FSIGNBIT ( f: r -- s: minus? )
This word corresponds to isSignMinus in IEEE 5.7.2, "General
operations". The name is based on C99.
The following are all required by IEEE. See IEEE 5.5.1, "Sign
bit operations". The IEEE copy() function is superfluous in
Forth [**IIUC].
FNEGATE ( f: r -- -r )
FABS ( f: r -- |r| )
The DPANS94 specification is extended to IEEE specials.
FCOPYSIGN ( f: r1 r2 -- r3 )
The output r3 is r1 with its sign bit replaced by that of r2.
7.8 Nearest integer functions
-----------------------------
FCEIL ( f: r1 -- r2 ) [**new]
FLOOR ( f: r1 -- r2 )
FROUND ( f: r1 -- r2 )
FTRUNC ( f: r1 -- r2 )
These words correspond to the respective IEEE required
operations:
roundToIntegralTowardPositive
roundToIntegralTowardNegative
roundToIntegralTiesToEven
roundToIntegralTowardZero
See IEEE 5.3.1, "General operations" and 5.9, "Details of
operations to round a floating-point datum to integral value".
The names are based on C99. No word is defined for IEEE
roundToIntegralTiesToAway.
[**Separate RfD for FCEIL ?]
FNEARBYINT ( f: r1 -- r2 )
This word corresponds to the IEEE required operation:
roundToIntegralExact
It performs the function of whichever of the other four
corresponds to the current rounding mode, and shall be
provided only if the rounding mode words are provided.
7.9 Number manipulation
------------------------
FMAX ( f: r1 r2 -- r3 )
FMIN ( f: r1 r2 -- r3 )
The DPANS94 specification is extended to IEEE specials. See
minNum and maxNum in IEEE 5.3.1, "General operations" and 6.2,
"Operations with NaNs".
FNEXTUP ( f: r1 -- r2 )
When r1 is a nonzero real number, FNEXTUP returns the next
affinely extended real in the default format that compares
larger than r1. See IEEE 5.3.1, "General operations" for the
behavior when r1 is an IEEE special. FNEXTDOWN is not
defined. According to IEEE, nextDown(x) is -nextUp(-x).
FSCALBN ( f: r s: n -- f: r*2^n )
The output is efficiently scaled by 2^n. See IEEE 5.3.3,
"logBFormat operations".
FLOGB ( f: r -- e )
Leave the radix-two exponent e of the fp representation as an
fp integer. If r is subnormal, the exponent is computed as if
r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
operations" for treatment of IEEE specials.
7.10 Exceptions
---------------
[**UNDER CONSTRUCTION
This is an important section, because the IEEE treatment of
exceptions is central to the philosophy of its treatment of
nans, infinities, and subnormal numbers. AFAIU, that was
strongly influenced by William Kahan.
However, it doesn't follow that it needs to be in this document.
It certainly needs discussion.
]
7.11 Rounding modes
-------------------
[**UNDER CONSTRUCTION
This section is currently under discussion in comp.lang.forth.
See IEEE 9.3, "Operations on dynamic modes for attributes". Only
words corresponding to 9.3.1, "Operations on individual dynamic
modes", are expected to be implemented, and among those,
roundTiesToAway is not expected to be implemented.
]
8 REFERENCES
[1] "IEEE Standard for Floating-Point Arithmetic", approved June
12, 2008 as IEEE Std 754-2008:
http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4610935
[2] "DRAFT Standard for Floating-Point Arithmetic P754", IEEE
754 draft 1.2.9, January 27, 2007:
http://www.validlab.com/754R/nonabelian.com/754/comments/Q754.129.pdf
[3] Wikipedia, "IEEE 754-2008":
http://en.wikipedia.org/wiki/IEEE_754
[4] ANSI X3.215-1994 final draft:
http://www.taygeta.com/forth/dpans.html
[5] ISO/IEC 15145:1997:
http://webstore.ansi.org/RecordDetail.aspx?sku=ISO%2fIEC+15145%3a1997
http://www.iso.org/iso/catalogue_detail.htm?csnumber=26479
[6] ISO/IEC 9899:1999 (December 1, 1999),
ISO/IEC 9899:1999 Cor. 1:2001(E),
ISO/IEC 9899:1999 Cor. 2:2004(E),
ISO/IEC 9899:1999 Cor. 3:2007(E):
http://www.open-std.org/jtc1/sc22/wg14/www/standards.html#9899
[7] C99 + TC1 + TC2 + TC3 is included in the freely available
WG14/N1256, September 7, 2007 [**thanks to David Thompson]:
http://www.open-std.org/jtc1/sc22/WG14/www/docs/n1256.pdf
[8] Single UNIX 3, AFAICS duplicates the C99 library spec, with
some things pinned down more tightly [**thanks to David
Thompson]:
http://www.unix.org/single_unix_specification/
A.2a IEEE BINARY FLOATING-POINT FORMATS
IEEE 754-2008 defines three basic binary fp formats, binary32,
binary64, and binary128, plus three corresponding extended
binary formats, whose parameters are shown in Tables 1 and 2
below. It also defines the four binary interchange formats
shown in Table 3, plus those with storage widths of more than
128 bits that are a multiple of 32 bits.
Table 1: Parameters for IEEE 754-2008
basic binary formats.
binary32 binary64 binary128
---------------------------------------
p = 24 53 113
emax = 127 1023 16383
Table 2: Parameters for IEEE 754-2008
extended binary formats.
binary32 binary64 binary128
---------------------------------------
p >= 32 64 128
emax >= 1023 16383 65535
Table 3: Parameters for IEEE 754-2008
binary interchange formats
(k is the storage width in
bits).
binary 16 binary32 binary64 binary128
---------------------------------------------------
k = 16 32 64 128
p = 11 24 53 113
emax = 15 127 1023 16383
Note that the intel 80-bit format corresponds to one of the
extended binary64 formats, with p = 64 and emax = 16383. Its
precision is greater than that of basic binary64 and less than
that of basic binary128, with exponent range the same as basic
binary128. Although it is not defined as a basic IEEE binary
format, it may be called the "binary80" basic format. Its
implementation normally differs from that of the other basic
formats by having an explicit leading bit for normal and
subnormal numbers.
Binary interchange formats are logical formats, with unspecified
memory layout. They all have an implicit leading bit for normal
and subnormal numbers.
Note that the binary128 interchange format is the only one in
Table 3 that can contain the binary80 basic format. IEEE does
not define a binary80 interchange format.
A.6.1 NAN SIGNS AND LOADS
IEEE allows the load for nan results like 0E 0E F/ to be
anything, so the following does not necessarily give a zero
load:
0E 0E F/ FABS FCONSTANT +NAN
Aside from the ambiguous load, the FABS (extended to nan) is
necessary here, because not only does IEEE not specify it, but
both pfe and gforth actually give opposite signs for 0E 0E F/
under ppc (+) vs. intel (-) Mac OS X. They do both give quiet
nans with zero load.
As a matter of fact, the intel QNaN "floating-point indefinite"
is the qnan with zero load and negative sign, according to
"Intel(R) 64 and IA-32 Architectures Software Developer's
Manual, Volume 1: Basic Architecture", Table 4-3:
http://www.intel.com/Assets/PDF/manual/253665.pdf
A.7.3 COMPARISON (INFORMATIVE RATIONALE)
[**Any high-level definitions in this section assume a separate
floating-point stack.]
The twelve IEEE required comparisons are the following, where
"N" means logical negation, "?" stands for "unordered", and "<?"
and ">?" stand for "less than or unordered" and "greater than or
unordered":
< = > ? N< N= N> N? <= >= <? >?
Unfortunately, the common notation "?" for the unordered
predicate clashes with Forth practice, where "?" usually
[**always?] indicates a flag or test. The ? notation is used
here as a convenience for IEEE predicates, and does not appear
in any corresponding Forth names.
The <= and >= comparisons are no longer simple negations of >
and <, but are rather the AND's of those negations with N?. It
can be shown that <? is N(>=), and >? is N(<=). See IEEE Table
5.3, "Required unordered-quiet predicates and negations".
It is possible to implement all of the IEEE comparisons via
high-level definitions in terms of a few low-level words, even
fewer than the five required for this word set. The following
remarks are offered as a guide to possibilities for the choice
of low-level words.
DPANS94 has only F< and F~. Since IEEE "=" is semantically
different from 0E F~, low-level implementation of F= seems
inevitable. The two words F< and F= are probably a minimum set
for high-level implementation of the rest.
For example, IEEE ">" is not expressible in terms of "<" and "="
plus logical operations, but F> can be defined as:
: F> ( f: r1 r2 -- s: [r1>r2]? ) FSWAP F< ;
FUNORDERED is not required by this document; in particular, the
name is not reserved; but it can be defined as
: FUNORDERED ( f: r1 r2 -- s: [r1?r2]? )
FDUP F= FDUP F= AND 0= ;
The logical negations N<, N=, N>, and N? can be expressed with
the Forth phrases "F< 0=", etc.
Forth words for the <= and >= predicates can be defined as
: F2DUP ( f: r1 r2 -- r1 r2 r1 r2 ) FOVER FOVER ;
: F<= ( f: r1 r2 -- s: [r1<=r2]? ) F2DUP F< F= OR ;
: F>= ( f: r1 r2 -- s: [r1>=r2]? ) F2DUP F> F= OR ;
The <? and >? predicates can be expressed by the Forth
phrases "F>= 0=" and "F<= 0=".
Thus all twelve IEEE predicates can be expressed with the five
required words, and as few as two.
[**Treat the following as a footnote.]
It can be shown that the closure of the four fundamental IEEE
comparison predicates under AND, OR, and negation consists of
sixteen independent relations, including the twelve that IEEE
requires. Two additional elements are the trivial,
identically true and false relations, and the other two are
"less than or greater than" and its negation, "unordered or
equal". The five words of this word set, plus their five
negations, implement the only nontrivial, transitive relations
among the sixteen.
On the other hand, several current cpu's have efficient
operations for all four of the fundamental IEEE comparisons, <,
>, =, and ?. Low-level implementations of at least F<, F>, and
F=, would be natural for such systems.
All of the words in the F0< family can be defined by analogy to
: F0< ( f: r -- s: [r<0]? ) 0E F< ;
Krishna Myneni
> Here's the current draft. I started a new thread because the
> old one has become unwieldy. That doesn't mean that ongoing
> conversations in the old thread shouldn't continue there.
It will take some time to digest all of this. Thanks for the effort.
I have not been following the earlier thread too closely, but I recall
mention of some test code for compliance with this optional
floating-point word set.
Krishna
David N. Williams
> [...]
>
> I have not been following the earlier thread too closely, but I recall
> mention of some test code for compliance with this optional
> floating-point word set.
There are three, so far -- actually I forgot to post the third
until now. Ah! Now I see why. Here are the first two:
http://www-personal.umich.edu/~williams/archive/forth/utilities/fatan2-test.fs
http://www-personal.umich.edu/~williams/archive/forth/utilities/ieee-fprox-test.fs
I believe that fatan2-test.fs needs to be revised a bit in view
of what I've learned about nans.
The fprox tests assume the proposed extension to IEEE-FP.
The third is now posted:
http://www-personal.umich.edu/~williams/archive/forth/utilities/ieee-arith-test.fs
Aside from not testing the nonexistent F*+, it was incomplete
for FSQRT, namely, there was a test
t{ -1e fsqrt fabs -> +nan }t
which failed with Mac OS X ppc. The reason is that the nan
produced is loaded, whereas +NAN isn't. Demonstrations for my
ppc and intel systems are copied at the end of this message.
The correct way to do the test is:
t{ -1e fsqrt fnan? -> true }t
FNAN? doesn't exist yet, so I put together a huge kludge to
define and partially test it. The problem was (and is) how to
define it with DPANS94 words most likely to be IEEE compliant.
If you're interested in making some IEEE special tests for the
other fp functions, you could have a look at IEEE 754-2008,
Sec. 9.2.1 (reference [1] in the draft). If you don't have a
copy, let me know and we'll try to figure something out.
-- David
-------------------------------------------------------------
Gforth gives the same result as pfe for both of the following
systems.
---------------
[jost:~/utilities] dnwills% pfe hexfloat.fs
\ Portable Forth Environment 0.33.61 (Aug 21 2008 11:00:57)
Copyright (C) Dirk Uwe Zoller 1993 - 1995.
Copyright (C) Tektronix, Inc. 1998 - 2003.
ANS/ffa ITC Forth - Please enter LICENSE and WARRANTY.
"[VOID]" is redefined
Running on powerpc darwin9.4.0 - to quit say BYE. ok
-1e fsqrt f>raw hex ud. 7FF8000000000001 ok
----------------
That's the quiet, positive nan with load 1.
------------------
schwartz:utilities dnwills$ pfe hexfloat.fs
\ Portable Forth Environment 0.33.70 (Mar 15 2009 12:12:20)
Copyright (C) Dirk Uwe Zoller 1993 - 1995.
Copyright (C) Tektronix, Inc. 1998 - 2003.
Copyright (C) Guido U. Draheim 2005 - 2008.
ANS/ffa ITC Forth - Please enter LICENSE and WARRANTY.
"[VOID]" is redefined
Running on i386 darwin9.6.2 - to quit say BYE. ok
-1e fsqrt f>raw hex ud. FFF8000000000000 ok
------------------
That's the quiet, negative nan with load 0.
So Mac OS X ppc produces a loaded, positive nan for sqrt(-1),
while Mac OS X intel produces an unloaded, negative nan.
Marcel Hendrix
> Anton Ertl wrote:
>> "David N. Williams" <will...@umich.edu> writes:
[..]
>> It seems that these are supposed to be words (and not just values
>> recognized by the FP input conversion), right. If so, define them as
>> words. I don't think that two names for each infinity is going to
>> fly, though.
> Right, they're just Forth words. Okay, how about this:
> ----------------------------
> +INF ( f: -- +Inf )
> -INF ( f: -- -Inf )
> +NAN ( f: -- +NaN )
> -NAN ( f: -- -NaN )
> These words return, respectively, IEEE +/-infinity and the quiet
> +/-nan with zero load in the default format.
> -----------------------------
> IIRC Marcel wanted to have the names +INFINITY and -INFINITY for
> the infinities. The three-letter version does seem more
> IEEE-like to me, but I don't feel strongly about it.
I wanted +INFINITY and -INFINITY if they were text recognized
by >FLOAT (i.e. user input). +Inf and -Inf are OK for programmers.
>> However, Gforth does not recognize, e.g., "1e-", and I have never
>> heard a complaint about that, so maybe there are no such standard
>> programs.
> Okay, I'll leave this out and mark it for a future RfD.
Just for reference: iForth recognizes "1e-". Why did you want to
remove it?
Is there a full current version of the proposal somewhere?
-marcel
Marcel Hendrix
[..]
> That's the quiet, negative nan with load 0.
> So Mac OS X ppc produces a loaded, positive nan for sqrt(-1),
> while Mac OS X intel produces an unloaded, negative nan.
iForth version 3.0.3, generated 10:38:09, July 6, 2009.
x86_32 binary, native floating-point, double precision.
Copyright 1996 - 2009 Marcel Hendrix.
FORTH> -1e fsqrt pad df! pad D@ HEX UD. DECIMAL FFF8000000000000 ok
FORTH> -1e fsqrt pad xf! pad 10 dump
$005B5588 00 00 00 00 00 00 00 C0--FF FF
iForth version 3.0.3, generated 14:01:56, June 28, 2009.
x86_64 binary, native floating-point, double precision.
Copyright 1996 - 2009 Marcel Hendrix.
FORTH> -1e fsqrt pad df! pad @ H. $FFF8000000000000 ok
FORTH> -1e fsqrt pad xf! pad 10 dump
$011BB020 00 00 00 00 00 00 00 C0--FF FF
ok
FORTH> -1e fsqrt pad xf! pad d@ hex ud. FFFFC000000000000000 ok
FORTH>
-marcel
Ed
> Here's the current draft. I started a new thread because the
> old one has become unwieldy. That doesn't mean that ongoing
> conversations in the old thread shouldn't continue there.
> ... When IEEE-FP is present, that specification
> shall be replaced by
>
> Convertible string := <significand><exponent>
>
> <significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
> [**Note the extra .<digits> option.]
This requires users to remember that when they're writing non-IEEE
f/p code that they must insert a digit before the decimal point.
Either this change is made for all f/p users, or one sticks to the
'94 convention. I use internal syntax checking so I'm never caught
out with leading decimal points and it enforces consistent forth
source style. F. requires a leading digit. >FLOAT permits
leading decimal point, IMO that's sufficient.
David N. Williams
Thanks! I especially appreciate seeing explicit examples of how
xf! works.
The conclusion for intel seems pretty clear -- the hardware
FSQRT is being used by gcc, as I assume it is by iForth.
According to intel documentation, for nonzero negative arguments
it returns the standard intel "floating-point indefinite" QNaN,
which is the qnan with zero load and negative sign.
For ppc, there is no hardware FSQRT (at least not on earlier
cpu's), but presumably there is a trap to system firmware. I
just haven't bothered to check the assembly language output for
gcc's sqrt() ...
-- David
David N. Williams
> David N. Williams wrote:
>> Here's the current draft. I started a new thread because the
>> old one has become unwieldy. That doesn't mean that ongoing
>> conversations in the old thread shouldn't continue there.
>
>> ... When IEEE-FP is present, that specification
>> shall be replaced by
>>
>> Convertible string := <significand><exponent>
>>
>> <significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
>> [**Note the extra .<digits> option.]
>
> This requires users to remember that when they're writing non-IEEE
> f/p code that they must insert a digit before the decimal point.
They're required to remember it already, which I think is not a
good thing.
> Either this change is made for all f/p users, or one sticks to the
> '94 convention.
I do agree with this point. So consider the extra .<digits>
option removed. I hope somebody will submit an RfD to put it
back in the non-IEEE part of Forth 200x.
That leaves only the extra "e" option for <e-char>.
-- David
Ed
> ...
> [**ALTERNATIVE 2:
>
> >FLOAT ( c-addr u -- [r: r s: true]|[false] )
>
> The expression "the string represents a valid floating-point
> number" in the DPANS94 specification shall be interpreted to
> mean that it represents a finite number in the range of the
> default format.
>
> >IEEE-FLOAT ( c-addr u -- [r: r s: true]|[false] )
>
> The specification is exactly that of >FLOAT in alternative 1.
> ]
>
> [**I prefer alternative 1. I mean, really, a "float" in the
> IEEE context is any IEEE default datum. But I would accept a
> variant of 2 if 1 should prove to be a show stopper.]
In the event the inconceivable happens, here's some PD code to
get you going. Works with any f/p stack model. Tested briefly.
--
Thinking about this system in more general terms ...
A generic name for this function should suffice e.g. >XFLOAT .
It represents an 'extended' float input - not just for IEEE, but for
any f/p system which has "non-number" data that must be input.
Unlike >FLOAT which was intended to be portable, >XFLOAT
is specific to the f/p system which defines it. Thus an IEEE
>XFLOAT will support Nans/Infs, but a >XFLOAT in another
f/p system would support non-numbers peculiar to that system.
There would be no name/function conflict as these words would
exist mutually exclusive to each other.
This isn't a proposal - just seeing where it leads.
--
\ NAN/INF constants for text interpreter
0 1 D>F 0 0 D>F F/ ( +INF )
FDUP FCONSTANT +INF
FNEGATE FCONSTANT -INF
0 0 D>F 0 0 D>F F/ ( -NAN )
FDUP FCONSTANT -NAN
FNEGATE FCONSTANT +NAN
\ >XFLOAT implementation for IEEE
\ CAPS-COMPARE is a case-insensitive COMPARE.
\ getc gets are support functions from my PD >FLOAT
\ package:
: getc ( a u -- a' u' c )
1 /STRING OVER 1- C@ ;
\ get sign
: gets ( a u -- a' u' n|0 )
DUP IF
getc DUP [CHAR] - = IF EXIT THEN
[CHAR] + <> /STRING
THEN 0 ;
: >XFLOAT ( c-addr u -- r true | false )
2DUP 2>R gets >R
2DUP S" NAN" CAPS-COMPARE IF
2DUP S" INF" CAPS-COMPARE IF
2DUP S" INFINITY" CAPS-COMPARE IF
2DROP R> DROP 2R> >FLOAT EXIT
THEN
THEN 2DROP +INF
ELSE 2DROP +NAN THEN
R> IF FNEGATE THEN TRUE 2R> 2DROP ;
\ Test compliance
: CHK ( addr len flag )
>R CR [CHAR] " EMIT 2DUP TYPE [CHAR] " EMIT
10 OVER - SPACES >XFLOAT DUP >R IF FDROP THEN R>
." --> " DUP IF ." TRUE " ELSE ." FALSE" THEN
R> - IF ." *fail* " ELSE ." pass " THEN ;
: TEST ( -- )
CR ." Checking >XFLOAT IEEE compliance ..." CR
\ S" ." FALSE CHK
\ S" E" FALSE CHK
\ S" .E" FALSE CHK
\ S" .E-" FALSE CHK
\ S" +" FALSE CHK
\ S" -" FALSE CHK
\ S" 9" FALSE CHK
\ S" 9 " FALSE CHK
\ S" " TRUE CHK
\ S" " TRUE CHK
\ S" 1+1" TRUE CHK
\ S" 1-1" TRUE CHK
S" 9" TRUE CHK
S" 9." TRUE CHK
S" .9" TRUE CHK
S" 9E" TRUE CHK
\ S" 9e+" TRUE CHK
\ S" 9d-" TRUE CHK
S" +NaN" TRUE CHK
S" -NAN" TRUE CHK
S" nan" TRUE CHK
S" NAN " FALSE CHK
S" +Inf" TRUE CHK
S" -INF" TRUE CHK
S" infinity" TRUE CHK
S" -Infinity" TRUE CHK
;
TEST
Marcel Hendrix
> Here's the current draft. I started a new thread because the
> old one has become unwieldy. That doesn't mean that ongoing
> conversations in the old thread shouldn't continue there.
Some preliminary remarks, found when trying to implement the
new stuff. I think (I have got them all by now.)
> FSIGNALING? ( r: r -- s: snan? )
How to test this? SNaNs are not produced by (Intel) hardware.
When loading an SNaN bitpattern, the hardware seems to convert
it to a QNaN immediately.
> FNEARBYINT ( f: r1 -- r2 )
Isn't this FROUND ?
> FNEXTUP ( f: r1 -- r2 )
Isn't this (the new) FCEIL ?
> FLOGB ( f: r -- e )
> Leave the radix-two exponent e of the fp representation as an
> fp integer. If r is subnormal, the exponent is computed as if
> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
> operations" for treatment of IEEE specials.
An fp "integer?" Why not ( F: r -- ) ( -- e ) then?
Is "8.2222e FLOGB F." resulting in 3.039524e or in 3e ?
-marcel
Ed
> ...
> > FNEARBYINT ( f: r1 -- r2 )
>
> Isn't this FROUND ?
If only it were.
IIUC, FNEARBYINT performs according to the currently
set rounding mode.
As perhaps you've found, the name can be misleading. E.g.
2.9e would produce 2.0e when truncate was in force - not
exactly what most would associate with "nearby". However
Forth has no need to adopt ambiguous names used by other
languages!
In Swiftforth, it is FROUND that behaves according to the
current rounding mode. The name is well-suited to the task
since it is unspecific.
ANS never said FROUND should be fixed at RNE for all time
- only that generic f/p programs can expect it and that systems
should be able to supply it. IEEE now gives users more options.
Marcel Hendrix
> The third is now posted:
> http://www-personal.umich.edu/~williams/archive/forth/utilities/ieee-arith-test.fs
[..]
This is the result of running ieee-arith~test.fs.
FORTH> in ieee-arith-test.frt
Creating John Hopkins TESTER Version 2.00
TESTING F+
INCORRECT FP RESULT: t{ -nan +nan f+ -> -nan }t
TESTING F-
INCORRECT FP RESULT: t{ -nan +nan f- -> -nan }t
TESTING F*
INCORRECT FP RESULT: t{ -nan +nan f* -> -nan }t
TESTING F/
INCORRECT FP RESULT: t{ -nan +nan f/ -> -nan }t
TESTING FSQRT
NOT TESTING F*+
#ERRORS: 4
ok
The actual result of these expressions is:
FORTH> -nan +nan f+ f. +NAN ok
FORTH> -nan +nan f- f. +NAN ok
FORTH> -nan +nan f* f. +NAN ok
FORTH> -nan +nan f/ f. +NAN ok
This is the result of the Intel hardware. Maybe these four
tests ( with -nan +nan ) should not be counted as errors,
or the test could t{ -nan +nan <op> fabs -> +nan }t
(see also the remarks in ieee-arith~test.fs).
-marcel
David N. Williams
Good to know! The remarks were like this:
\ These are for information only. IEEE only requires the output
\ to be one of the two input nans.
t{ +nan +nan f+ -> +nan }t
t{ -nan +nan f+ -> -nan }t
t{ +nan -nan f+ -> +nan }t
t{ -nan -nan f+ -> -nan }t
The tests are fixed as you suggested. See version 0.5.1:
http://www-personal.umich.edu/~williams/archive/forth/utilities/ieee-arith-test.fs
I was in a hurry when I uploaded the first version, and somehow
missed the obvious test. At least the old one did elicit information!
-- David
David N. Williams
> Marcel Hendrix wrote:
>> ...
>>> FNEARBYINT ( f: r1 -- r2 )
>> Isn't this FROUND ?
>
> If only it were.
>
> IIUC, FNEARBYINT performs according to the currently
> set rounding mode.
Right. From the draft:
-------------------------------
FNEARBYINT ( f: r1 -- r2 )
This word corresponds to the IEEE required operation:
roundToIntegralExact
It performs the function of whichever of the other four
corresponds to the current rounding mode, and shall be
provided only if the rounding mode words are provided.
-------------------------------
The draft should have mentioned that the name comes from C99's
nearbyint(), in the same style as the four specific rounding
words.
> As perhaps you've found, the name can be misleading. E.g.
> 2.9e would produce 2.0e when truncate was in force - not
> exactly what most would associate with "nearby". However
> Forth has no need to adopt ambiguous names used by other
> languages!
>
> In Swiftforth, it is FROUND that behaves according to the
> current rounding mode. The name is well-suited to the task
> since it is unspecific.
The name FNEARBYINT is also unspecific, although I agree FROUND
would have been a better name. But DPANS94 is specific about
requiring it to use "round to nearest".
> ANS never said FROUND should be fixed at RNE for all time
> - only that generic f/p programs can expect it and that systems
> should be able to supply it. IEEE now gives users more options.
I'm really puzzled why you think this is okay for FROUND and not
for >FLOAT?
-- David
David N. Williams
> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
>
>> Here's the current draft. I started a new thread because the
>> old one has become unwieldy. That doesn't mean that ongoing
>> conversations in the old thread shouldn't continue there.
>
> Some preliminary remarks, found when trying to implement the
> new stuff. I think (I have got them all by now.)
>
>> FSIGNALING? ( r: r -- s: snan? )
>
> How to test this? SNaNs are not produced by (Intel) hardware.
> When loading an SNaN bitpattern, the hardware seems to convert
> it to a QNaN immediately.
Uh-oh!
Let's replace it by
FQUIET? ( r: r -- s: snan? )
Then FSIGNALING? is equivalent to 0= FQUIET?, and doesn't need
to be mentioned. Presumably there are some intel exceptions
that generate a signalling nan, so the word isn't superfluous,
even for intel.
>> FNEARBYINT ( f: r1 -- r2 )
>
> Isn't this FROUND ?
No, see Ed's reply, and my reply to that.
>> FNEXTUP ( f: r1 -- r2 )
>
> Isn't this (the new) FCEIL ?
No. FCEIL produces an integer, whereas FNEXTUP produces the
next representable affine real. It's the same as C99's
nextafter(x,y), but with the argument y wired to be in the
positive direction from x.
>> FLOGB ( f: r -- e )
>
>> Leave the radix-two exponent e of the fp representation as an
>> fp integer. If r is subnormal, the exponent is computed as if
>> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
>> operations" for treatment of IEEE specials.
>
> An fp "integer?" Why not ( F: r -- ) ( -- e ) then?
> Is "8.2222e FLOGB F." resulting in 3.039524e or in 3e ?
I'll reply to this separately, because there's a complication.
-- David
David N. Williams
>> FLOGB ( f: r -- e )
>
>> Leave the radix-two exponent e of the fp representation as an
>> fp integer. If r is subnormal, the exponent is computed as if
>> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
>> operations" for treatment of IEEE specials.
>
> An fp "integer?" Why not ( F: r -- ) ( -- e ) then?
> Is "8.2222e FLOGB F." resulting in 3.039524e or in 3e ?
It's just the exponent. In mixed hex/decimal notation,
8.2222e taken as a double is 1.0 71C4 32CA 57A8 p+3, so the
answer is 3e . :-)
Thanks for picking up on this -- it's indeed a problem.
The earliest draft had both FLOGB and FILOGB, but it was
suggested that IEEE requires logB but not ilogB. IEEE 754-2008
is maybe more abstract. It requires that logB and scaleB use
the same logBFormat, which is allowed to be either an ordinary
integer or a floating-point integer.
That's already inconsistent in the draft because FSCALBN uses an
integer logBFormat while FLOGB uses a floating point logBFormat.
I remember now why I did it that way. The integer format for
FSCALBN seems really all around better, but IEEE 5.3.3,
"logBFormat operations", requires the following for logB:
When logBFormat is an integer format, then logB(NaN), logB(Inf),
and logB(0) return language-defined values outside the range
+/-2 x (emax + p - 1) and signal the invalid operation
exception.
That's a bit complicated, but one probably *would* usually want
a non-fp integer when there's no exception.
How about the following replacement for FLOGB?
---------------------------------
FILOGB ( f: r -- s: n )
Leave the radix-two exponent n of the fp representation as an
integer. If r is subnormal, the exponent is computed as if r
were normalized, with e < emin. When r is signed zero,
infinity, or nan, the invalid operation exception is signaled.
In those cases, n is the DPANS94 MAX-N for +infinity or a
positive nan, and minus MAX-N for -infinity, a negative nan,
or zero with either sign.
----------------------------------
What "signaled" is going to mean remains to be seen.
The names FILOGB and FSCALBN are based on C99.
What do people think?
-- David
David N. Williams
> Marcel Hendrix wrote:
> > "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
> >
> >> Here's the current draft. I started a new thread because the
> >> old one has become unwieldy. That doesn't mean that ongoing
> >> conversations in the old thread shouldn't continue there.
> >
> > Some preliminary remarks, found when trying to implement the
> > new stuff. I think (I have got them all by now.)
> >
> >> FSIGNALING? ( r: r -- s: snan? )
> >
> > How to test this? SNaNs are not produced by (Intel) hardware.
> > When loading an SNaN bitpattern, the hardware seems to convert
> > it to a QNaN immediately.
>
> Uh-oh!
>
> Let's replace it by
>
> FQUIET? ( r: r -- s: snan? )
>
> Then FSIGNALING? is equivalent to 0= FQUIET?, and doesn't need
> to be mentioned. Presumably there are some intel exceptions
> that generate a signalling nan, so the word isn't superfluous,
> even for intel.
Damn! FSIGNALING? is *not* equivalent to 0= FQUIET?, because
the latter is true for ordinary numbers, etc., while FSIGNALING?
is true only for nans that are signaling.
So both are needed.
Then, as you say, how to test FSIGNALING? Let me rummage around
in IEEE to see if there's something that really is required to
produce a signaling nan. Or maybe somebody already knows such a
case for intel?
-- David
Marcel Hendrix
> Marcel Hendrix wrote:
>> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
[..]
>>> FSIGNALING? ( r: r -- s: snan? )
>>
>> How to test this? SNaNs are not produced by (Intel) hardware.
>> When loading an SNaN bitpattern, the hardware seems to convert
>> it to a QNaN immediately.
> Uh-oh!
> Let's replace it by
> FQUIET? ( r: r -- s: snan? )
> Then FSIGNALING? is equivalent to 0= FQUIET?, and doesn't need
> to be mentioned. Presumably there are some intel exceptions
> that generate a signalling nan, so the word isn't superfluous,
> even for intel.
I don't understand your solution ... ?
The problem I have is in *testing*, I don't see another way than really
generating a SNaN, and don't know how to do it. Loading a pattern with
DF@ is not working: on load the SNaN is converted to a QNaN.
I guess I wouldn't miss FSIGNALING? much.
>>> FNEARBYINT ( f: r1 -- r2 )
>> Isn't this FROUND ?
> No, see Ed's reply, and my reply to that.
This is clear now.
>>> FNEXTUP ( f: r1 -- r2 )
>> Isn't this (the new) FCEIL ?
> No. FCEIL produces an integer, whereas FNEXTUP produces the
> next representable affine real. It's the same as C99's
> nextafter(x,y), but with the argument y wired to be in the
> positive direction from x.
Could you comment on the (ugly) implementation below?
>>> FLOGB ( f: r -- e )
>>> Leave the radix-two exponent e of the fp representation as an
>>> fp integer. If r is subnormal, the exponent is computed as if
>>> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
>>> operations" for treatment of IEEE specials.
>> An fp "integer?" Why not ( F: r -- ) ( -- e ) then?
>> Is "8.2222e FLOGB F." resulting in 3.039524e or in 3e ?
> I'll reply to this separately, because there's a complication.
[..]
> When logBFormat is an integer format, then logB(NaN), logB(Inf),
> and logB(0) return language-defined values outside the range
> +/-2 x (emax + p - 1) and signal the invalid operation
> exception.
> That's a bit complicated, but one probably *would* usually want
> a non-fp integer when there's no exception.
> How about the following replacement for FLOGB?
No thanks, you explained quite well that the logic of FLOGB is much
simpler wrt NaN and such, so I have no issue anymore. I also note
that FLOGB generates better code than FILOGB :-)
-marcel
-- ------------
ANEW -nextafter
DOC
(*
FORTH> 1.23e FDUP 1e FNEXTAFTER F< . 0 ok
FORTH> 1.23e FDUP 1e FNEXTAFTER F> . -1 ok
FORTH> 1.23e FDUP 2e FNEXTAFTER F> . 0 ok
FORTH> 1.23e FDUP 2e FNEXTAFTER F< . -1 ok
FORTH> +nan +nan FNEXTAFTER f. +NAN ok
FORTH> -nan +nan FNEXTAFTER f. +NAN ok
FORTH> +nan -nan FNEXTAFTER f. +NAN ok
FORTH> -nan -nan FNEXTAFTER f. -NAN ok
FORTH> 1e FDUP 1e FNEXTAFTER F= . -1 ok
FORTH> 0e +1e FNEXTAFTER f. 0.000000 ok
FORTH> 0e -1e FNEXTAFTER f. 0.000000 ok
FORTH> +Inf +Inf FNEXTAFTER f. +INF ok
FORTH> +Inf 0e FNEXTAFTER f. 1.797693e308 ok
FORTH> -Inf 0e FNEXTAFTER f. -1.797693e308 ok
FORTH> -Inf -Inf FNEXTAFTER f. -INF ok
*)
ENDDOC
-- Only for doubles. Even uglier than F~
-- Needs FLOCALS that you can take the address of.
: FNEXTAFTER ( F: x y -- z )
FLOCALS| y x |
'OF x 4 + @ LOCAL hx \ HI(x)
'OF x @ LOCAL lx \ LO(x)
'OF y 4 + @ LOCAL hy \ HI(y)
'OF y @ LOCAL ly \ LO(y)
hx $7FFFFFFF AND LOCAL ix \ |x|
hy $7FFFFFFF AND LOCAL iy \ |y|
ix $7FF00000 >= ix $7FF00000 - lx OR 0<> AND \ x is nan
iy $7FF00000 >= iy $7FF00000 - ly OR 0<> AND \ y is nan
OR IF x y F+ EXIT ENDIF ( both NaN )
x y F= IF x EXIT ENDIF
ix lx OR 0= IF hy $80000000 AND 'OF x 4 + ! \ HI(x)
1 'OF x ! \ LO(x)
x FSQR TO y \ return +-minsubnormal
y x F= IF y ELSE x ENDIF EXIT
ENDIF
hx 0>=
IF hx hy = lx ly > AND hx hy >
OR IF lx 0= IF -1 +TO hx ENDIF -1 +TO lx \ x > y, x -= ulp
ELSE 1 +TO lx lx 0= IF 1 +TO hx ENDIF \ x < y, x += ulp
ENDIF
ELSE hx hy = lx ly > AND hx hy > OR hy 0>=
OR IF lx 0= IF -1 +TO hx ENDIF -1 +TO lx \ x > y, x -= ulp
ELSE 1 +TO lx lx 0= IF 1 +TO hx ENDIF \ x < y, x += ulp
ENDIF
ENDIF
hx $7FF00000 AND TO hy
hy $7FF00000 >= IF x FSQR EXIT ENDIF \ overflow
hy $00100000 < IF x FSQR TO y \ underflow
x y F<> IF hx 'OF y 4 + !
lx 'OF y ! y EXIT
ENDIF
ENDIF
hx 'OF x 4 + !
lx 'OF x ! x ;
: FNEXTUP ( F: x -- y ) FDUP F1+ FNEXTAFTER ;
David N. Williams
> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
>
>> Marcel Hendrix wrote:
>>> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
> [..]
> >>> FSIGNALING? ( r: r -- s: snan? )
> >>
> >> How to test this? SNaNs are not produced by (Intel) hardware.
> >> When loading an SNaN bitpattern, the hardware seems to convert
> >> it to a QNaN immediately.
>
>> Uh-oh!
>
>> Let's replace it by
>
>> FQUIET? ( r: r -- s: snan? )
>
>> Then FSIGNALING? is equivalent to 0= FQUIET?, and doesn't need
>> to be mentioned. Presumably there are some intel exceptions
>> that generate a signalling nan, so the word isn't superfluous,
>> even for intel.
>
> I don't understand your solution ... ?
>
> The problem I have is in *testing*, I don't see another way than really
> generating a SNaN, and don't know how to do it. Loading a pattern with
> DF@ is not working: on load the SNaN is converted to a QNaN.
>
> I guess I wouldn't miss FSIGNALING? much.
My response to myself about this apparently crossed this in the
ether.
>>>> FNEARBYINT ( f: r1 -- r2 )
>
>>>[...]
>
> This is clear now.
Good!
>>>> FNEXTUP ( f: r1 -- r2 )
>
>>> Isn't this (the new) FCEIL ?
>
>> No. FCEIL produces an integer, whereas FNEXTUP produces the
>> next representable affine real. It's the same as C99's
>> nextafter(x,y), but with the argument y wired to be in the
>> positive direction from x.
>
> Could you comment on the (ugly) implementation below?
I plan to do so, but need a little time to digest it.
>>>> FLOGB ( f: r -- e )
>
>>>> Leave the radix-two exponent e of the fp representation as an
>>>> fp integer. If r is subnormal, the exponent is computed as if
>>>> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
>>>> operations" for treatment of IEEE specials.
>
>>> An fp "integer?" Why not ( F: r -- ) ( -- e ) then?
>>> Is "8.2222e FLOGB F." resulting in 3.039524e or in 3e ?
>
>> I'll reply to this separately, because there's a complication.
>
> [..]
>
>> When logBFormat is an integer format, then logB(NaN), logB(Inf),
>> and logB(0) return language-defined values outside the range
>> +/-2 x (emax + p - 1) and signal the invalid operation
>> exception.
>
>> That's a bit complicated, but one probably *would* usually want
>> a non-fp integer when there's no exception.
>
>> How about the following replacement for FLOGB?
>
> No thanks, you explained quite well that the logic of FLOGB is much
> simpler wrt NaN and such, so I have no issue anymore. I also note
> that FLOGB generates better code than FILOGB :-)
I'm not so sure it was that clear, but you discerned my exact point
about simplicity wrt NaN, etc. :-)
So unless there's an objection, I propose to forget about
FILOGB, keep both FSCALBN and FLOGB, and eschew IEEE logBFormat
consistency.
-- David
David N. Williams
> David N. Williams wrote:
> > Marcel Hendrix wrote:
>
> Damn! FSIGNALING? is *not* equivalent to 0= FQUIET?, because
> the latter is true for ordinary numbers, etc., while FSIGNALING?
> is true only for nans that are signaling.
>
> So both are needed.
>
> Then, as you say, how to test FSIGNALING? Let me rummage around
> in IEEE to see if there's something that really is required to
> produce a signaling nan. Or maybe somebody already knows such a
> case for intel?
I didn't see anything illuminating in IEEE. I found the following in
"Intel(R) 64 and IA-32 Architectures Software Developer's
Manual, Volume 1: Basic Architecture":
http://www.intel.com/Assets/PDF/manual/253665.pdf
SNaNs are typically used to trap or invoke an exception
handler. They must be inserted by software; that is, the
processor never generates an SNaN as a result of a
floating-point operation.
By unmasking the invalid operation exception, the programmer
can use signaling NaNs to trap to the exception handler. The
generality of this approach and the large number of NaN values
that are available provide the sophisticated programmer with a
tool that can be applied to a variety of special situations.
I'm in way over my head here, but I'm wondering what effect
unmasking would have on DF@.
IIUC, you found something like the following. Assuming an
aligned PAD initially containing a double sNaN,
PAD 2@ PAD DF@ PAD DF! PAD 2@ D=
returns false, because PAD now contains qNaN. Would that still be
the case after unmasking? After unmasking, is DF@ going to trap?
As I say, it's over my head.
-- David
Andreas
>
> I didn't see anything illuminating in IEEE. I found the following in
> "Intel(R) 64 and IA-32 Architectures Software Developer's
> Manual, Volume 1: Basic Architecture":
>
> http://www.intel.com/Assets/PDF/manual/253665.pdf
>
> 4.8.3.4:
>
> SNaNs are typically used to trap or invoke an exception
> handler. They must be inserted by software; that is, the
> processor never generates an SNaN as a result of a
> floating-point operation.
>
>
> 4.8.3.6:
>
> By unmasking the invalid operation exception, the programmer
> can use signaling NaNs to trap to the exception handler. The
> generality of this approach and the large number of NaN values
> that are available provide the sophisticated programmer with a
> tool that can be applied to a variety of special situations.
>
> I'm in way over my head here, but I'm wondering what effect
> unmasking would have on DF@.
>
> IIUC, you found something like the following. Assuming an
> aligned PAD initially containing a double sNaN,
>
> PAD 2@ PAD DF@ PAD DF! PAD 2@ D=
>
> returns false, because PAD now contains qNaN. Would that still be
> the case after unmasking? After unmasking, is DF@ going to trap?
>
> As I say, it's over my head.
>
> -- David
There is one thing that I can see as useable for applications: QNANs
propagate through instructions. So you can use the available QNAN range
for encoding additional information.
For example whether a fp-number shall be treated in reality as a 32 bit
address pointer, or an integer - similar to unions in C. IIRC
Strongforth uses a separate type info stack; perhaps one could get rid
of it by exploiting QNANs.
Just a green idea...
Andreas
Marcel Hendrix
> David N. Williams wrote:
>> David N. Williams wrote:
>> > Marcel Hendrix wrote:
[...]
> IIUC, you found something like the following. Assuming an
> aligned PAD initially containing a double sNaN,
> PAD 2@ PAD DF@ PAD DF! PAD 2@ D=
> returns false, because PAD now contains qNaN. Would that still be
> the case after unmasking? After unmasking, is DF@ going to trap?
I was indeed trying that
Unfortunately, SNaN detection is muxed with FPU stack under- and
overflow. When I unmask the invalid operation exception (IM) in the
controlword, I get an infinite loop of FLOAT STACK exceptions. Either
there is a pending float stack exception, or iForth does not handle the
FPU stack well. The infinite loop is because the iForth stack handler
does not reset the proper bits (this handler normally never runs and
apparently it is wrong). For the moment, I can't provide more
information.
However, it is clear that if exceptions/traps are to be part of the
IEEE FLOAT extension, the host Forth must allow to write and install
proper trap handlers. That might be a bit overdone, considering the
size of the perceived problem.
-marcel
Andrew Haley
> ---------------------------------
> FILOGB ( f: r -- s: n )
> Leave the radix-two exponent n of the fp representation as an
> integer. If r is subnormal, the exponent is computed as if r
> were normalized, with e < emin. When r is signed zero,
> infinity, or nan, the invalid operation exception is signaled.
> In those cases, n is the DPANS94 MAX-N for +infinity or a
> positive nan, and minus MAX-N for -infinity, a negative nan,
> or zero with either sign.
I don't think there's a great deal to be gained by treasting +/- NaN
specially. Returning minus MAX-N for -infinity is definitely wrong.
IMO, this function should ignore the sign of the argument in all
cases.
All these gotchas suggest to me that it might be much easier to leave
the result on the fp stack!
Andrew.
Ed
> ...
>
> The DPANS94 specification is extended to include IEEE
> specials. The valid-result flag, flag2, is false if r is
> infinity or a nan.
>
> [**IMHO this is still unsettled, still needing discussion.
>
> 1. Should the current rounding mode be used, instead of
> roundTiesToEven? I'll have to check to what extent that
> has already been discussed.
The point of having a "change current rounding mode"
function is so that functions which "round" will change.
ANS defines which functions "round". In '94, they are:
DF! DF@ SF! SF@ FROUND REPRESENT and by
implication F. FS. FE.
To exclude certain functions from this list would be
entirely arbitrary. I don't know what criteria one could
possibly use.
If one wants to temporarily change the behaviour of one
of these functions then: set the system to the rounding mode
wanted, execute the function, and restore it afterwards.
If one accepts that a mode change should encompass all
functions which "round", it then leaves one without a
permanent RNE function since FROUND will follow the
current mode. An explicit RNE function would need to
be provided.
Ed
> Ed wrote:
> > ...
> > In Swiftforth, it is FROUND that behaves according to the
> > current rounding mode. The name is well-suited to the task
> > since it is unspecific.
>
> The name FNEARBYINT is also unspecific, although I agree FROUND
> would have been a better name. But DPANS94 is specific about
> requiring it to use "round to nearest".
>
> > ANS never said FROUND should be fixed at RNE for all time
> > - only that generic f/p programs can expect it and that systems
> > should be able to supply it. IEEE now gives users more options.
>
> I'm really puzzled why you think this is okay for FROUND and not
> for >FLOAT?
Look at the consequences.
Nothing changes for FROUND. When one sets the system rounding
mode to RNE, all one's portable ANS programs run as before.
However when one expands >FLOAT to include nan/infs, it opens
a back-door whereby a portable program may fail e.g. the infinite
loop example.
Whereas a programmer can guarantee any portable program will
run with FROUND (by ensuring the system is set to the proper
rounding mode beforehand) - he can't give the same guarantee with
an expanded >FLOAT because he may have no control over the
data input to that program.
I really don't mind what happens. Should I ever need it, I can
always define:
: >FLOAT >FLOATX ;
See my other post for more FROUND related stuff.
Marcel Hendrix
[..]
> ANS defines which functions "round". In '94, they are:
> DF! DF@ SF! SF@ FROUND REPRESENT and by
> implication F. FS. FE.
FROUND rounds explicitly to nearest.
> To exclude certain functions from this list would be
> entirely arbitrary. I don't know what criteria one could
> possibly use.
> If one wants to temporarily change the behaviour of one
> of these functions then: set the system to the rounding mode
> wanted, execute the function, and restore it afterwards.
> If one accepts that a mode change should encompass all
> functions which "round", it then leaves one without a
> permanent RNE function since FROUND will follow the
> current mode. An explicit RNE function would need to
> be provided.
FROUND rounds to nearest, FNEARBYINT uses the rounding mode.
No problem.
-marcel
Andrew Haley
> David N. Williams wrote:
> > ...
> > REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
> >
> > The DPANS94 specification is extended to include IEEE
> > specials. The valid-result flag, flag2, is false if r is
> > infinity or a nan.
> >
> > [**IMHO this is still unsettled, still needing discussion.
> >
> > 1. Should the current rounding mode be used, instead of
> > roundTiesToEven? I'll have to check to what extent that
> > has already been discussed.
> The point of having a "change current rounding mode" function is so
> that functions which "round" will change.
> ANS defines which functions "round". In '94, they are: DF! DF@ SF!
> SF@ FROUND REPRESENT and by implication F. FS. FE.
And F+, F-, etc, etc.
> To exclude certain functions from this list would be entirely
> arbitrary. I don't know what criteria one could possibly use.
Well, it's not really arbitrary: it all depends on how hard it would
be to implement and how useful the result would be.
Insisting that all functions follow the current rounding mode makes a
bunch of optimizations impossible or very difficult. For example, an
optimizer can't precompute 0.5E FSIN if it doesn't know the mode in
which it will be executed. We should be *extremely* careful about
this, since it potentially has a great many knock-on effects and
imposes great burdens on implementers for doubtful benefits.
I'd be in favour of making REPRESENT in any mode other than
roundTiesToEven undefined.
> If one wants to temporarily change the behaviour of one of these
> functions then: set the system to the rounding mode wanted, execute
> the function, and restore it afterwards.
Right.
> If one accepts that a mode change should encompass all functions
> which "round",
I certainly don't: see above.
> it then leaves one without a permanent RNE function since FROUND
> will follow the current mode. An explicit RNE function would need
> to be provided.
I can't see the point of a permanent RNE function. You can always set
the mode to RNE and call FROUND.
Andrew.
Andrew Haley
> I can't see the point of a permanent RNE function. You can always
> set the mode to RNE and call FROUND.
Actually, I withdraw that: a permanent RNE function might be a
convenience.
Andrew.
Marcel Hendrix
> What's the timeframe for shipments to the US?
Your CD went in the mail yesterday (Thursday). Delivery is
unpredictable, shortest time I've seen is a few days.
Thanks again for your support!
-marcel
Marcel Hendrix
> Andrew Haley <andr...@littlepinkcloud.invalid> wrote:
Given exception problems and unreadability resulting from a global
setting (e.g. calling OS or library functions), wouldn't it be better
to not allow setting the rounding mode at all? Just explicitly use
FROUND, FTRUNC, FLOOR and FCEIL where needed, assuming the system is
always in RNE (or default set by OS) mode. FNEARBYINT would then just
be a (much) faster default rounder (probably FROUND).
-marcel
Ed
> "Ed" <nos...@invalid.com> writes Re: IEEE-FP 0.5.2 - Rounding
> [..]
> > ANS defines which functions "round". In '94, they are:
> > DF! DF@ SF! SF@ FROUND REPRESENT and by
> > implication F. FS. FE.
>
> FROUND rounds explicitly to nearest.
As do all those functions. Read their definitions.
You have no problem making the rounding of these other
functions user-defined. But you want to leave FROUND
out. It's arbitrary and I see no rational reason or need
for it.
> FROUND rounds to nearest, FNEARBYINT uses the rounding mode.
> No problem.
You have it backwards. Consistency demands FROUND
follows the rounding mode of the system as will the other
functions. FNEAR can be the permanent RNE function.
Marcel Hendrix
> Marcel Hendrix wrote:
>> "Ed" <nos...@invalid.com> writes Re: IEEE-FP 0.5.2 - Rounding
[..]
>> > ANS defines which functions "round". In '94, they are:
>>> DF! DF@ SF! SF@ FROUND REPRESENT and by
>>> implication F. FS. FE.
>> FROUND rounds explicitly to nearest.
> As do all those functions. Read their definitions.
e.g.:
12.6.2.1203 DF! "d-f-store" FLOATING EXT
( df-addr -- ) ( F: r -- ) or ( r df-addr -- )
Store the floating-point number r as a 64-bit IEEE double-precision
number at df-addr. If the significand of the internal representation
of r has more precision than the IEEE double-precision format, it will
be rounded using the "round to nearest" rule. An ambiguous condition
exists if the exponent of r is too large to be accommodated in IEEE
double-precision format.
It is unclear what they mean by this, but it certainly doesn't
involve rounding in the sense of 'executing FROUND/FTRUNC/FCEIL/FLOOR'.
(I guess this is where one would need FNEXTAFTER/FNEXTUP or worse.)
> You have no problem making the rounding of these other
> functions user-defined. But you want to leave FROUND
> out. It's arbitrary and I see no rational reason or need
> for it.
I hope I cleared it up.
>> FROUND rounds to nearest, FNEARBYINT uses the rounding mode.
>> No problem.
This statement holds. It has nothing to do with the details of
DF! DF! SF! SF@. How REPRESENT, F. FS., and FE. would be affected
I don't know. I expect that for numbers representable in
IEEE double-precision there is no change. What a system whose
'significand of the internal representation of r has more
precision than the IEEE double-precision format' does with
numbers more precise than double-precision (or what it seem to
report in MAX-FLOAT) I don't know. But I'm pretty sure it doesn't
involve FROUND/FTRUNC/FCEIL/ FLOORing them.
-marcel
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes Re: IEEE-FP 0.5.2 - Rounding
> > Andrew Haley <andr...@littlepinkcloud.invalid> wrote:
> >> I can't see the point of a permanent RNE function. You can always
> >> set the mode to RNE and call FROUND.
> > Actually, I withdraw that: a permanent RNE function might be a
> > convenience.
> Given exception problems and unreadability resulting from a global
> setting (e.g. calling OS or library functions), wouldn't it be
> better to not allow setting the rounding mode at all?
I think not, because:
* Some systems already have the ability to set the rounding mode, and
it's useful to standardize the syntax to make programs portable
between such systems.
* I don't believe that there are any huge problems with exceptions,
and I'm not aware of any readbility problems at all.
* Mode control is useful for a bunch of things, particularly interval
arithmetic. Floored floating-point division is useful for
floating-point output, for example.
> Just explicitly use FROUND, FTRUNC, FLOOR and FCEIL where needed,
> assuming the system is always in RNE (or default set by OS)
> mode.
That doesn't get you floor and ceiling arithmetic.
> FNEARBYINT would then just be a (much) faster default rounder
> (probably FROUND).
Why would it be faster? On x86 FNEARBYINT is just a single FRNDINT
instruction, surely.
Andrew.
Marcel Hendrix
> Marcel Hendrix <m...@iae.nl> wrote:
>> Andrew Haley <andr...@littlepinkcloud.invalid> writes Re: IEEE-FP 0.5.2 - Rounding
[..]
>> Given exception problems and unreadability resulting from a global
>> setting (e.g. calling OS or library functions), wouldn't it be
>> better to not allow setting the rounding mode at all?
> I think not, because:
> * Some systems already have the ability to set the rounding mode, and
> it's useful to standardize the syntax to make programs portable
> between such systems.
Ok.
> * I don't believe that there are any huge problems with exceptions,
> and I'm not aware of any readbility problems at all.
Won't it need extra throw codes and words to catch such exceptions (e.g.
when sse2 is being used).
Readability: With a new rounding mode in effect, old and tested words may
work different from before (like using BASE for some programs). Multitasking
systems will need to set/restore rounding mode too, as external library calls.
> * Mode control is useful for a bunch of things, particularly interval
> arithmetic. Floored floating-point division is useful for
> floating-point output, for example.
I hope indeed somebody will use this. I am not aware having seen any
aiming-to-be-portable Forth code that needs it.
>> Just explicitly use FROUND, FTRUNC, FLOOR and FCEIL where needed,
>> assuming the system is always in RNE (or default set by OS)
>> mode.
> That doesn't get you floor and ceiling arithmetic.
At least it will need to be formulated what this means exactly.
>> FNEARBYINT would then just be a (much) faster default rounder
>> (probably FROUND).
> Why would it be faster? On x86 FNEARBYINT is just a single FRNDINT
> instruction, surely.
Indeed, FNEARBYINT will be much faster than the others (who need
to read and write the FPU control word, a nasty operation).
-marcel
-- ------------
VARIABLE y
y VALUE x
-- "Ref" approximates the integer work of the other rounders.
: RTEST CR ." RTEST -- ref " TIMER-RESET 1e #100000000 0 DO x @ DROP x @ DROP LOOP FDROP .ELAPSED
CR ." RTEST -- fnearbyint " TIMER-RESET 1e #100000000 0 DO FNEARBYINT LOOP FDROP .ELAPSED
CR ." RTEST -- fround " TIMER-RESET 1e #100000000 0 DO FROUND LOOP FDROP .ELAPSED
CR ." RTEST -- floor " TIMER-RESET 1e #100000000 0 DO FLOOR LOOP FDROP .ELAPSED
CR ." RTEST -- fceil " TIMER-RESET 1e #100000000 0 DO FCEIL LOOP FDROP .ELAPSED ;
FORTH> rtest
RTEST -- ref 0.207 seconds elapsed.
RTEST -- fnearbyint 0.787 seconds elapsed.
RTEST -- fround 1.168 seconds elapsed.
RTEST -- floor 1.167 seconds elapsed.
RTEST -- fceil 1.167 seconds elapsed. ok
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes Re: IEEE-FP 0.5.2 - Rounding
> > Marcel Hendrix <m...@iae.nl> wrote:
> >> Andrew Haley <andr...@littlepinkcloud.invalid> writes Re: IEEE-FP 0.5.2 - Rounding
> [..]
> >> Given exception problems and unreadability resulting from a global
> >> setting (e.g. calling OS or library functions), wouldn't it be
> >> better to not allow setting the rounding mode at all?
> > I think not, because:
> > * Some systems already have the ability to set the rounding mode, and
> > it's useful to standardize the syntax to make programs portable
> > between such systems.
> Ok.
> > * I don't believe that there are any huge problems with exceptions,
> > and I'm not aware of any readbility problems at all.
> Won't it need extra throw codes and words to catch such exceptions
> (e.g. when sse2 is being used).
No. I don't think that IEEE floating-point should have any
dependencies on the exception handling wordset. IEEE floating-point
in C doesn't -- C doesn't have exceptions at all -- and I don't see
why we need them either.
As we discussed before, a bracketed wordset like
trunc MODE{ stuff }MODE
is a solution for handling rounding modes. I'm still looking at a
nice way to do this, BTW.
An alternative is is ' stuff trunc MODE-EXECUTE which is probably
easier and a better candidate for standardization.
> Readability: With a new rounding mode in effect, old and tested
> words may work different from before (like using BASE for some
> programs).
Yes. On the other hand, old and tested words now have potentially
useful new behaviours. It's a metter of inspection and testing.
> Multitasking systems will need to set/restore rounding mode too, as
> external library calls.
Yes. But multitasking systems have to save and restore the hardware
floating-point status word anyway.
Both of these seem like very minor issues.
> > * Mode control is useful for a bunch of things, particularly interval
> > arithmetic. Floored floating-point division is useful for
> > floating-point output, for example.
> I hope indeed somebody will use this. I am not aware having seen any
> aiming-to-be-portable Forth code that needs it.
This seems like a circular argument. Clearly we can't do this
portably now, so any aiming-to-be-portable Forth code doesn't use it.
Being able to use the rounding modes portably means that some code in
the floating-point wordset support package that currently uses
system-dependent words will be able to use standard ones.
> >> Just explicitly use FROUND, FTRUNC, FLOOR and FCEIL where needed,
> >> assuming the system is always in RNE (or default set by OS)
> >> mode.
> > That doesn't get you floor and ceiling arithmetic.
> At least it will need to be formulated what this means exactly.
It's formulated by IEEE-754.
> >> FNEARBYINT would then just be a (much) faster default rounder
> >> (probably FROUND).
> > Why would it be faster? On x86 FNEARBYINT is just a single
> > FRNDINT instruction, surely.
> Indeed, FNEARBYINT will be much faster than the others (who need
> to read and write the FPU control word, a nasty operation).
Right. As I said before, I prefer simply to redefine FROUND, but I
understand the argument against that.
Andrew.
Ed
> "Ed" <nos...@invalid.com> writes Re: RfD: IEEE-FP 0.5.2
>
> > Marcel Hendrix wrote:
> >> "Ed" <nos...@invalid.com> writes Re: IEEE-FP 0.5.2 - Rounding
> [..]
> >> > ANS defines which functions "round". In '94, they are:
> >>> DF! DF@ SF! SF@ FROUND REPRESENT and by
> >>> implication F. FS. FE.
>
> >> FROUND rounds explicitly to nearest.
>
> > As do all those functions. Read their definitions.
>
> e.g.:
>
> 12.6.2.1203 DF! "d-f-store" FLOATING EXT
> ( df-addr -- ) ( F: r -- ) or ( r df-addr -- )
> Store the floating-point number r as a 64-bit IEEE double-precision
> number at df-addr. If the significand of the internal representation
> of r has more precision than the IEEE double-precision format, it will
> be rounded using the "round to nearest" rule. An ambiguous condition
> exists if the exponent of r is too large to be accommodated in IEEE
> double-precision format.
>
> It is unclear what they mean by this, but it certainly doesn't
> involve rounding in the sense of 'executing FROUND/FTRUNC/FCEIL/FLOOR'.
> (I guess this is where one would need FNEXTAFTER/FNEXTUP or worse.)
"round to nearest" doesn't necessarily mean rounding to an
integer. It's rounding in the same sense that, say, F+ must
round-off the extra bits used during the calculation.
In the case of DF! SF! the FPU memory store function
automatically rounds as needed to fit the target, according
to the current rounding mode.
> How REPRESENT, F. FS., and FE. would be affected
> I don't know.
> ...
F. FS. FE. use REPRESENT
REPRESENT uses FROUND
FROUND is one of the '94 functions specified to "round".
FROUND is typically defined as:
code FROUND ( r1 -- r2 )
sp di mov qword 0 [di] fld frndint
qword 0 [di] fstp wait next
end-code
FRNDINT rounds according to the current rounding mode.
When you want "94 rounding", set the FPU to RNE.
IEEE lets you choose other rounding methods (trunc,
ceil, floor etc).
Where's the problem?
Anton Ertl
>Insisting that all functions follow the current rounding mode makes a
>bunch of optimizations impossible or very difficult. For example, an
>optimizer can't precompute 0.5E FSIN if it doesn't know the mode in
>which it will be executed.
So what? This is Forth. If I want to optimize that, I can write
[ 0.5e fsin ] fliteral
or alternatively write sin(1/2), with that word defined as follows:
0.5e fsin fconstant sin(1/2)
This gives the programmer full control of which rounding I want and
how much optimization I get. We don't need the "make the language
unpredictable to allow optimization" meme in Forth.
>We should be *extremely* careful about
>this, since it potentially has a great many knock-on effects and
>imposes great burdens on implementers for doubtful benefits.
We should be extremely careful about this, since it has a great many
knock-on effects and imposes great burdens on *programmers* for
doubtful benefits.
- anton
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: http://www.forth200x.org/forth200x.html
EuroForth 2009: http://www.euroforth.org/ef09/
Anton Ertl
>Andrew Haley <andr...@littlepinkcloud.invalid> writes Re: IEEE-FP 0.5.2 - Rounding
>
>> Andrew Haley <andr...@littlepinkcloud.invalid> wrote:
>
>>> I can't see the point of a permanent RNE function. You can always
>>> set the mode to RNE and call FROUND.
>
>> Actually, I withdraw that: a permanent RNE function might be a
>> convenience.
>
>Given exception problems and unreadability resulting from a global
>setting (e.g. calling OS or library functions), wouldn't it be better
>to not allow setting the rounding mode at all? Just explicitly use
>FROUND, FTRUNC, FLOOR and FCEIL where needed, assuming the system is
>always in RNE (or default set by OS) mode.
The main point of the rounding mode is not to affect FROUND, on the
contrary. The main point of the rounding mode is to determine what
happens for an ordinary operation like F* when the result cannot be
represented exactly in an FP number (i.e., if there is an inexact
IEEE-exception). Which way should the result be rounded?
The main use of other rounding modes than RNE seems to be for getting
an idea how strongly the result is affected by precision.
When a programmer writes FROUND, he explicitly requests rounding to an
integer. That seems to me quite different from the inevitable and
implicit roundings on simple FP operations, and there the programmer
can and should explicitly state whether he wants RNE or rounding
according to the current rounding mode. And my guess is that in most
uses of FROUND RNE is more appropriate, so it should not be affected
by the current rounding mode.
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >Insisting that all functions follow the current rounding mode makes
> >a bunch of optimizations impossible or very difficult. For
> >example, an optimizer can't precompute 0.5E FSIN if it doesn't know
> >the mode in which it will be executed.
> So what? This is Forth. If I want to optimize that, I can write
> [ 0.5e fsin ] fliteral
> or alternatively write sin(1/2), with that word defined as follows:
> 0.5e fsin fconstant sin(1/2)
OK. It's an example. The point is that we have to be aware of what
we give up when making these decisions.
> This gives the programmer full control of which rounding I want and
> how much optimization I get.
I'm all in favour of that. What I'm trying not to do is to tie the
implementor's hands behind thir back.
> We don't need the "make the language unpredictable to allow
> optimization" meme in Forth.
Huh? I don't understand this. No-one was suggesting unpredictability
as a goal. Quite the reverse, in fact.
> >We should be *extremely* careful about this, since it potentially
> >has a great many knock-on effects and imposes great burdens on
> >implementers for doubtful benefits.
> We should be extremely careful about this, since it has a great many
> knock-on effects and imposes great burdens on *programmers* for
> doubtful benefits.
I don't understand what point you're trying to make here. I was
saying that making every word that takes a floating-point argument
honour the current rounding mode is a bad idea. It potentially has a
great bloating effect on the system, and it makes some optimizations
hard. There is IMO very small benefit for potentially great cost.
The non-default rounding modes are useful in some cases for a few
primitive arithmetic words, but extending that to all words that take
floating-point args would be a massive undertaking that would all but
guarntee that no-one implements this proposal.
Andrew.
Bernd Paysan
> I don't understand what point you're trying to make here. I was
> saying that making every word that takes a floating-point argument
> honour the current rounding mode is a bad idea. It potentially has a
> great bloating effect on the system, and it makes some optimizations
> hard. There is IMO very small benefit for potentially great cost.
>
> The non-default rounding modes are useful in some cases for a few
> primitive arithmetic words, but extending that to all words that take
> floating-point args would be a massive undertaking that would all but
> guarntee that no-one implements this proposal.
On x86 hardware, it is fairly trivial to make all operations honor the
current rounding mode, because the hardware has a current rounding mode. On
architectures where the rounding mode is hard-coded into the instruction, it
is a lot more difficult.
--
Bernd Paysan
"If you want it done right, you have to do it yourself"
http://www.jwdt.com/~paysan/
Andrew Haley
> Andrew Haley wrote:
> > I don't understand what point you're trying to make here. I was
> > saying that making every word that takes a floating-point argument
> > honour the current rounding mode is a bad idea. It potentially has a
> > great bloating effect on the system, and it makes some optimizations
> > hard. There is IMO very small benefit for potentially great cost.
> > The non-default rounding modes are useful in some cases for a few
> > primitive arithmetic words, but extending that to all words that take
> > floating-point args would be a massive undertaking that would all but
> > guarntee that no-one implements this proposal.
> On x86 hardware, it is fairly trivial to make all operations honor
> the current rounding mode, because the hardware has a current
> rounding mode.
It's easy to make the core instructions that the hardware performs
honour the current rounding mode, but everything else has to be
specially written to do it.
Andrew.
Bernd Paysan
>> On x86 hardware, it is fairly trivial to make all operations honor
>> the current rounding mode, because the hardware has a current
>> rounding mode.
>
> It's easy to make the core instructions that the hardware performs
> honour the current rounding mode, but everything else has to be
> specially written to do it.
Well, on x86, "everything else" is also implemented in hardware, or at least
mostly so ;-). I suppose we should go the way the current C standard goes:
Make that a documentation issue, i.e. whether the non-primitive functions
(not F+, F-, F*, F/) honor the current rounding mode or not needs
documentation, but is not mandated by the standard. Draft here:
http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf
The support of IEEE-FP is already optional, this would add another option.
Options help the implementer, in terms of "quality of implementation" this
is certainly not so good for the user. My impression is that if there is a
current rounding mode, honoring it is a higher quality of implementation
than not honoring it; but if it is too difficult to honor it, round to
nearest (even) is ok.
Anton Ertl
>Anton Ertl <an...@mips.complang.tuwien.ac.at> wrote:
>> We don't need the "make the language unpredictable to allow
>> optimization" meme in Forth.
>
>Huh? I don't understand this. No-one was suggesting unpredictability
>as a goal. Quite the reverse, in fact.
>I don't understand what point you're trying to make here. I was
>saying that making every word that takes a floating-point argument
>honour the current rounding mode is a bad idea.
The question is: What should the words do instead? Do you require RNE
for these words? Or do you make the rounding (or maybe even more)
undefined for these words? The latter is what I meant with an
unpredictable language. You did propose an undefined and thus
unpredictable result for REPRESENT.
BTW, if you require RNE for words like FSIN, this will cost extra on
systems that use the 387 instructions to implement them.
>The non-default rounding modes are useful in some cases for a few
>primitive arithmetic words, but extending that to all words that take
>floating-point args would be a massive undertaking that would all but
>guarntee that no-one implements this proposal.
We should design it such that it is useful. My impression until now
is that rounding-mode control is used by a few people in a few
circumstances (and performance is not the be-all-end-all in these
circumstances), but then it should work predictably for all FP
operations involved (because it is used on FP routines that can use
all FP operations).
If that is the case, then the following might be a good design:
0) Without rounding mode control the inexact IEEE-exception is handled
by RNE.
1) Rounding-mode control is a separate extension from the rest of
IEEE-FP so that some systems may implement IEEE-FP but not
rounding-mode control.
2) Rounding-mode control applies to most FP words (probably not
FROUND, though).
3) Many systems probably want to implement rounding-mode control by
loading an extension that redefines some words in a slower way that
deals correctly with non-RNE rounding (and normally has RNE-only
implementations of these words). Therefore programs must ask for
rounding-mode control through an extension query.
Anton Ertl
>Well, on x86, "everything else" is also implemented in hardware, or at least
>mostly so ;-). I suppose we should go the way the current C standard goes:
>Make that a documentation issue, i.e. whether the non-primitive functions
>(not F+, F-, F*, F/) honor the current rounding mode or not needs
>documentation, but is not mandated by the standard.
Unpredictability is the C way of pseudo-standardization. If we
don't really want to standardize non-RNE rounding modes, we should be
honest about it and just don't do it.
>http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf
Is this a draft for the next C standard? Certainly the 2007 date
makes it look like it's not a draft for C99.
>The support of IEEE-FP is already optional, this would add another option.
>Options help the implementer, in terms of "quality of implementation" this
>is certainly not so good for the user. My impression is that if there is a
>current rounding mode, honoring it is a higher quality of implementation
>than not honoring it; but if it is too difficult to honor it, round to
>nearest (even) is ok.
No, we should define a predictable result. That's something I expect
the user of a non-default rounding mode to want (otherwise why would
they use a non-default rounding mode?). If it's too difficult to
implement that, the implementation should not pretend to support the
extension.
Anton Ertl
>On
>architectures where the rounding mode is hard-coded into the instruction, it
>is a lot more difficult.
Which hardware would that be? Alpha can code the rounding mode into
the instruction, but one of the rounding modes it can encode is to use
the one from the mode flags.
I guess one would use the static rounding modes for implementing
library functions such as FSIN, where most operations should not be
affected by the current rounding mode.
- anton
Andrew Haley
> Andrew Haley wrote:
> >> On x86 hardware, it is fairly trivial to make all operations honor
> >> the current rounding mode, because the hardware has a current
> >> rounding mode.
> >
> > It's easy to make the core instructions that the hardware performs
> > honour the current rounding mode, but everything else has to be
> > specially written to do it.
> Well, on x86, "everything else" is also implemented in hardware, or
> at least mostly so ;-). I suppose we should go the way the current C
> standard goes: Make that a documentation issue, i.e. whether the
> non-primitive functions (not F+, F-, F*, F/) honor the current
> rounding mode or not needs documentation, but is not mandated by the
> standard.
That's a good model, and is what is currently proposed for Forth.
> The support of IEEE-FP is already optional, this would add another
> option. Options help the implementer, in terms of "quality of
> implementation" this is certainly not so good for the user. My
> impression is that if there is a current rounding mode, honoring it
> is a higher quality of implementation than not honoring it; but if
> it is too difficult to honor it, round to nearest (even) is ok.
That is the line that IEEE 754 itself takes: correctly rounded
transcendental functions is a "should", not a "shall". The x86
doesn't do correctly rounded transcendentals.
Andrew.
Andrew Haley
> >The support of IEEE-FP is already optional, this would add another
> >option. Options help the implementer, in terms of "quality of
> >implementation" this is certainly not so good for the user. My
> >impression is that if there is a current rounding mode, honoring it
> >is a higher quality of implementation than not honoring it; but if
> >it is too difficult to honor it, round to nearest (even) is ok.
> No, we should define a predictable result. That's something I
> expect the user of a non-default rounding mode to want (otherwise
> why would they use a non-default rounding mode?).
The rounding modes for the basic arithmetic functions make it possible
to do interval arithmetic. They are also useful for decimal - float
conversions, allowing words that currently non-portable to be written
in standard Forth. They allow Forth programmers to control the
hardware in a reasonably portable way.
While it might be nice to have, say, a correctly rounded cosine, it
would be a lot of work and we wouldn't be able to use the versions
provided by the C library. Also, the functions would have difficult
to predict execution time.
We would be doing a huge amound of work for the sake of pedantry.
This is not, IMO, the Forth way.
> If it's too difficult to implement that, the implementation should
> not pretend to support the extension.
Even IEEE 754 is not so unrealistic as to require correct rounding for
all functions. IMO it's crazy to suggest that we should require
*more* than even the latest IEEE 754 does.
Andrew.
Bernd Paysan
> Bernd Paysan <bernd....@gmx.de> writes:
>>On
>>architectures where the rounding mode is hard-coded into the instruction,
>>it is a lot more difficult.
>
> Which hardware would that be? Alpha can code the rounding mode into
> the instruction, but one of the rounding modes it can encode is to use
> the one from the mode flags.
Since IEEE FP (854) defines a rounding mode, I suppose architectures without
such a current rounding modes are those which predate IEEE FP, and therefore
not relevant within this discussion:
http://754r.ucbtest.org/standards/854.html#rounding
Note that the rounding mode only affects "all arithmetic operations (except
remainder)". It does not talk about complex functions, which are not part of
IEEE 854, anyway.
> I guess one would use the static rounding modes for implementing
> library functions such as FSIN, where most operations should not be
> affected by the current rounding mode.
Yes. The main point here is that FSIN is likely to be implemented only using
those, and correct rounding to the current rounding mode is not easy to do.
Even generating a to-the-bit accurate result is not done very often, as it
requires more accurate intermediate results.
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >Anton Ertl <an...@mips.complang.tuwien.ac.at> wrote:
> >> We don't need the "make the language unpredictable to allow
> >> optimization" meme in Forth.
> >
> >Huh? I don't understand this. No-one was suggesting unpredictability
> >as a goal. Quite the reverse, in fact.
> ...
> >I don't understand what point you're trying to make here. I was
> >saying that making every word that takes a floating-point argument
> >honour the current rounding mode is a bad idea.
> The question is: What should the words do instead? Do you require
> RNE for these words? Or do you make the rounding (or maybe even
> more) undefined for these words? The latter is what I meant with an
> unpredictable language. You did propose an undefined and thus
> unpredictable result for REPRESENT.
I don't think it matters, since the use of REPRESENT in any mode other
than RNE is pretty unlikely. Forcing RNE is fine by me.
> BTW, if you require RNE for words like FSIN, this will cost extra on
> systems that use the 387 instructions to implement them.
It would, yes. Happily, IEEE 754 doesn't require any partiular
behaviour for FSIN.
> >The non-default rounding modes are useful in some cases for a few
> >primitive arithmetic words, but extending that to all words that
> >take floating-point args would be a massive undertaking that would
> >all but guarntee that no-one implements this proposal.
> We should design it such that it is useful. My impression until now
> is that rounding-mode control is used by a few people in a few
> circumstances (and performance is not the be-all-end-all in these
> circumstances), but then it should work predictably for all FP
> operations involved (because it is used on FP routines that can use
> all FP operations).
The problem is that you are stuck with some unpredictability: if you
insist on correct rounding for all FP operations, then occasionally
some opearations might take 100 times as long to execute. There is
more than one axis to predictability.
The consequence of your suggestion would be that anyone who used
rounding-mode control would not also have the benefit of fast
transcendental functions. Since the hardware is perfectly capable of
delivering both these things, its is absurd to suggest that Forth
should refuse it for the sake of sheer pedantry.
> If that is the case, then the following might be a good design:
> 0) Without rounding mode control the inexact IEEE-exception is handled
> by RNE.
> 1) Rounding-mode control is a separate extension from the rest of
> IEEE-FP so that some systems may implement IEEE-FP but not
> rounding-mode control.
Yes. That's what C does.
> 2) Rounding-mode control applies to most FP words
No; see above.
> (probably not FROUND, though).
> 3) Many systems probably want to implement rounding-mode control by
> loading an extension that redefines some words in a slower way that
> deals correctly with non-RNE rounding (and normally has RNE-only
> implementations of these words). Therefore programs must ask for
> rounding-mode control through an extension query.
As I said above, you're suggesting that anyone who wants to use
rounding-mode control has to put up with slow transcendentals, and
that is absurd. Another consequence that anyone wanting to allow
rounding-mode control would also have to use a multiple-precision
floating-point library. Since IEEE 754 doesn't require this, I don't
understand why you want it so much.
Andrew.
Ed
> Ed <nos...@invalid.com> wrote:
> ...
> > ANS defines which functions "round". In '94, they are: DF! DF@ SF!
> > SF@ FROUND REPRESENT and by implication F. FS. FE.
>
> And F+, F-, etc, etc.
Only the functions listed are specified as "round to nearest".
FLOOR is required to floor. Everything else is implementation
defined. (12.3.1.2 , 12.4.1.1)
> > To exclude certain functions from this list would be entirely
> > arbitrary. I don't know what criteria one could possibly use.
>
> Well, it's not really arbitrary: it all depends on how hard it would
> be to implement and how useful the result would be.
>
> Insisting that all functions follow the current rounding mode makes a
> bunch of optimizations impossible or very difficult.
> ...
The SwiftForth default is to truncate. Consequently all the
forth "round-to-nearest" functions truncate. You could ask
them what "great burdens" and "doubtful benefits" they've
encountered as a consequence.
You claim to support IEEE. IEEE allows changing the
default rounding mode. SwiftForth did that.
> > If one wants to temporarily change the behaviour of one of these
> > functions then: set the system to the rounding mode wanted, execute
> > the function, and restore it afterwards.
>
> Right.
Nevertheless you want Forth to nobble the feature by preventing
it from working on certain functions, or the routines which use
them. It's "knock-on" restrictions like that which don't impress me.
I don't need Forth dictating to me what rounding mode a function
should have - I tell it what I want. When I choose to change the
system default rounding, I expect all "round to nearest" functions
to follow suit. I can change it back whenever I want. That's freedom.
Andrew Haley
> Andrew Haley wrote:
> > Ed <nos...@invalid.com> wrote:
> > ...
> > > ANS defines which functions "round". In '94, they are: DF! DF@ SF!
> > > SF@ FROUND REPRESENT and by implication F. FS. FE.
> >
> > And F+, F-, etc, etc.
> Only the functions listed are specified as "round to nearest".
> FLOOR is required to floor. Everything else is implementation
> defined. (12.3.1.2 , 12.4.1.1)
OK.
> > > To exclude certain functions from this list would be entirely
> > > arbitrary. I don't know what criteria one could possibly use.
> >
> > Well, it's not really arbitrary: it all depends on how hard it would
> > be to implement and how useful the result would be.
> >
> > Insisting that all functions follow the current rounding mode makes a
> > bunch of optimizations impossible or very difficult.
> > ...
> The SwiftForth default is to truncate. Consequently all the forth
> "round-to-nearest" functions truncate.
This is, IMO, a bug: it decreases floating-point accuracy in almost
all cases.
> > > If one wants to temporarily change the behaviour of one of these
> > > functions then: set the system to the rounding mode wanted, execute
> > > the function, and restore it afterwards.
> >
> > Right.
> Nevertheless you want Forth to nobble the feature by preventing it
> from working on certain functions, or the routines which use them.
> It's "knock-on" restrictions like that which don't impress me.
No, I don't want Forth to nobble the feature at all.
> I don't need Forth dictating to me what rounding mode a function
> should have - I tell it what I want. When I choose to change the
> system default rounding, I expect all "round to nearest" functions
> to follow suit. I can change it back whenever I want. That's
> freedom.
I think you've misunderstood what I'm proposing. IEEE-754 does not
require that all functions be correctly rounded, and as far as I'm
aware no programming language that supports IEEE-754 does either. I'm
proposing that Forth does not do so.
Andrew.
Anton Ertl
>While it might be nice to have, say, a correctly rounded cosine, it
>would be a lot of work and we wouldn't be able to use the versions
>provided by the C library. Also, the functions would have difficult
>to predict execution time.
Oh, I see what you are talking about. No, I was not thinking about
bit-accurate results (which are called "correctly rounded" by IEEE 754
for some reason).
What I was thinking about was that the rounding in the final step(s)
of words like FCOS should be done according to the current roundiung
mode. That should be helpful in achieving the kind of error estimate
I expect that non-RNE rounding would be used for.
It would probably be even better if there was an implementation of
these words that guaranteed that the result would be in the right
direction from the exact result for the floored, ceiling, and
towards-zero rounding modes.
Looking at Figure G-5 to G-8 of my Pentium Processor's User's manual,
the Pentium (and AFAIK the later Intel processors use the same
algorithms and produce the same results) produces results within
+/-0.5ulp of the exact result for FCOS for most inputs in the range
they looked at. With that kind of accuracy, even the weaker thing I
had in mind should be useful.
OTOH, I wonder if we should standardize non-RNE rounding modes at the
moment at all. How many potential users are there? They should tell
us (possibly through the system implementors) what they need. If
nobody cares, there's no need to work on that.
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >While it might be nice to have, say, a correctly rounded cosine, it
> >would be a lot of work and we wouldn't be able to use the versions
> >provided by the C library. Also, the functions would have difficult
> >to predict execution time.
> Oh, I see what you are talking about. No, I was not thinking about
> bit-accurate results (which are called "correctly rounded" by IEEE
> 754 for some reason).
> What I was thinking about was that the rounding in the final step(s)
> of words like FCOS should be done according to the current roundiung
> mode. That should be helpful in achieving the kind of error
> estimate I expect that non-RNE rounding would be used for.
I have no problem with that, but since FCOS isn't going to be
correctly rounded I can see no way to specify such behaviour.
> It would probably be even better if there was an implementation of
> these words that guaranteed that the result would be in the right
> direction from the exact result for the floored, ceiling, and
> towards-zero rounding modes.
Maybe, but that requires 0.5ulp precision, which in turn requires
multiple-precision arithmetic. I don't think we should go there: as I
said, IEEE 754 doesn't make this requirement of implementations.
> Looking at Figure G-5 to G-8 of my Pentium Processor's User's
> manual, the Pentium (and AFAIK the later Intel processors use the
> same algorithms and produce the same results) produces results
> within +/-0.5ulp of the exact result for FCOS for most inputs in the
> range they looked at. With that kind of accuracy, even the weaker
> thing I had in mind should be useful.
Yes, but not all implementations will be as good as that. And how do
you specify "most inputs" ? :-)
> OTOH, I wonder if we should standardize non-RNE rounding modes at
> the moment at all. How many potential users are there? They should
> tell us (possibly through the system implementors) what they need.
> If nobody cares, there's no need to work on that.
It's not just exotic cases. For example, if you want to do exact
integer arithmetic without transferring operands into integer
registers, then you need to set roundToZero mode. This is useful, for
example, for decimal output.
Andrew.
David N. Williams
> [...]
>
> Could you comment on the (ugly) implementation below?
> [...]
Okay, a number of things intervened, but here are my comments.
I ANSI-fied it, and added ttester tests with the help of the
mixed hex mantissa, decimal radix two conversion from my just
updated hexfloat.fs. I'll describe that in a separate
announcement.
The tests work on Mac OS X ppc and intel, and on linux.
I must say, I don't find it ugly. At least the algorithm is
nice, albeit specialized to doubles.
A good way to do this would seem to be with the FSCALBN and
FLOGB. This example raises a question in my mind about what
FLOGB (or FILOGB) should really return, but that's another
discussion.
Your FNEXTAFTER performs exactly as I expected, and I *think*
that means I translated it correctly.
All three systems report the same glitches for FNEXTUP:
----------
INCORRECT FP RESULT: t{ -Inf FNEXTUP -> -max-+n }t
INCORRECT FP RESULT: t{ max-+n-down FNEXTUP -> max-+n }t
INCORRECT FP RESULT: t{ -max-+n FNEXTUP -> -max-+n-up }t
#ERRORS: 3
-Inf FNEXTUP f>mix type -Infok
max-+n-down FNEXTUP f>mix type +1.F FFFF FFFF FFFE p+1023ok
-max-+n FNEXTUP f>mix type -1.F FFFF FFFF FFFF p+1023ok
----------
-- David
Watch out for bad line wraps!
------------------------------------------------
\ fnextup-mh.fs
\ Marcel Hendrix, comp.lang.forth, July 11, 2009
\ ANSI-fied: dnw, July 16, 2009
\ Needs ttester.fs
(
http://www-personal.umich.edu/~williams/archive/forth/utilities/hexfloat.fs
)
s" hexfloat.fs" included \ fails if bits/cell <> 32
decimal
bits/float 64 <>
[IF] cr .( ***This code requires 64 bits/float.) ABORT [THEN]
[UNDEFINED] \\ [IF] \ for degugging
: \\ ( -- ) -1 parse 2drop BEGIN refill 0= UNTIL ; [THEN]
\ *** BADEN-STYLE FLOAT LOCALS
(
Wil Baden's cheap, nestable locals written for floats.
)
: framef ( -- addr ) here falign 2 floats allot ;
: unframe ( addr -- ) here - allot ;
: fpframe>r ( -- ) postpone framef postpone >r ; immediate
: r>frame ( -- ) postpone r> postpone unframe ; immediate
: EXIT-frame ( -- ) postpone r> postpone unframe postpone EXIT ; immediate
: 'x.hi here [ 2 floats [LITTLE-ENDIAN] [IF] 4 - [THEN] ] literal - ;
: 'x.lo here [ 2 floats [LITTLE-ENDIAN] 0= [IF] 4 - [THEN] ] literal - ;
: 'y.hi here [ 1 floats [LITTLE-ENDIAN] [IF] 4 - [THEN] ] literal - ;
: 'y.lo here [ 1 floats [LITTLE-ENDIAN] 0= [IF] 4 - [THEN] ] literal - ;
: x here [ 2 floats ] literal - f@ ;
: y here [ 1 floats ] literal - f@ ;
: TO-x here [ 2 floats ] literal - f! ;
: TO-y here [ 1 floats ] literal - f! ;
\ Locals don't work well with the >R in EXIT-FRAME.
0 value hx 0 value lx 0 value hy 0 value ly
0 value ix 0 value iy
: FSQR ( f: r -- r*r ) fdup f* ;
hex
: FNEXTAFTER ( F: x y -- z )
fpframe>r TO-y TO-x
'x.hi @ TO hx \ HI(x)
'x.lo @ TO lx \ LO(x)
'y.hi @ TO hy \ HI(y)
'y.lo @ TO ly \ LO(y)
hx 7FFFFFFF AND TO ix \ |x|
hy 7FFFFFFF AND TO iy \ |y|
ix 7FF00000 >= ix 7FF00000 - lx OR 0<> AND \ x is nan
iy 7FF00000 >= iy 7FF00000 - ly OR 0<> AND \ y is nan
OR IF x y F+ EXIT-frame ENDIF ( both NaN )
x y F= IF x EXIT-frame ENDIF
ix lx OR 0= IF hy 80000000 AND 'x.hi ! \ HI(x)
1 'x.lo ! \ LO(x)
x FSQR TO-y \ return +-minsubnormal
y x F= IF y ELSE x ENDIF EXIT-frame
ENDIF
hx 0>=
IF hx hy = lx ly > AND hx hy >
OR IF lx 0= IF hx 1- TO hx ENDIF lx 1- TO lx \ x > y,
x -= ulp
ELSE lx 1+ TO lx lx 0= IF hx 1+ TO hx ENDIF \ x < y,
x += ulp
ENDIF
ELSE hx hy = lx ly > AND hx hy > OR hy 0>=
OR IF lx 0= IF hx 1- TO hx ENDIF lx 1- TO lx \ x > y,
x -= ulp
ELSE lx 1+ TO lx lx 0= IF hx 1+ TO hx ENDIF \ x < y,
x += ulp
ENDIF
ENDIF
hx 7FF00000 AND TO hy
hy 7FF00000 >= IF x FSQR EXIT-frame ENDIF \ overflow
hy 00100000 < IF x FSQR TO-y \ underflow
x y F<> IF hx 'y.hi !
lx 'y.lo ! y EXIT-frame
ENDIF
ENDIF
hx 'x.hi !
lx 'x.lo ! x
r>frame ;
decimal
: FNEXTUP ( F: x -- y ) FDUP 1e F+ FNEXTAFTER ;
s" ttester.fs" included
true verbose !
set-exact
decimal
variable #errors 0 #errors !
:noname ( c-addr u -- )
(
Display an error message followed by the line that had the
error.
)
1 #errors +! error1 ; error-xt !
0e 0e f/ fabs fconstant +NAN
+NAN fnegate fconstant -NAN
1e 0e f/ fconstant +INF
+INF fnegate fconstant -INF
s" 1.F FFFF FFFF FFFF p+1023" mix>f fconstant max-+n
s" 1.0 0000 0000 0000 p-1022" mix>f fconstant min-+n
s" 0.F FFFF FFFF FFFF p-1022" mix>f fconstant max-+subn
s" 0.0 0000 0000 0001 p-1022" mix>f fconstant min-+subn
s" -0.0 0000 0000 0001 p-1022" mix>f fconstant -min-+subn
s" -0.F FFFF FFFF FFFF p-1022" mix>f fconstant -max-+subn
s" -1.0 0000 0000 0000 p-1022" mix>f fconstant -min-+n
s" -1.F FFFF FFFF FFFF p+1023" mix>f fconstant -max-+n
s" 1.0 0000 0000 0001 p0" mix>f fconstant 1-up
s" 1.F FFFF FFFF FFFE p+1023" mix>f fconstant max-+n-down
s" -1.F FFFF FFFF FFFF p-1" mix>f fconstant -1-up
s" -1.F FFFF FFFF FFFE p+1023" mix>f fconstant -max-+n-up
t{ 1.23e FDUP 1e FNEXTAFTER F< -> 0 }t
t{ 1.23e FDUP 1e FNEXTAFTER F> -> -1 }t
t{ 1.23e FDUP 2e FNEXTAFTER F> -> 0 }t
t{ 1.23e FDUP 2e FNEXTAFTER F< -> -1 }t
t{ +nan +nan FNEXTAFTER -> +NAN }t
t{ -nan +nan FNEXTAFTER fabs -> +NAN }t
t{ +nan -nan FNEXTAFTER fabs -> +NAN }t
t{ -nan -nan FNEXTAFTER -> -NAN }t
t{ 1e FDUP 1e FNEXTAFTER F= -> -1 }t
t{ 0e +1e FNEXTAFTER -> min-+subn }t
t{ 0e -1e FNEXTAFTER -> -min-+subn }t
t{ +Inf +Inf FNEXTAFTER -> +INF }t
t{ +Inf 0e FNEXTAFTER -> max-+n }t
t{ -Inf 0e FNEXTAFTER -> -max-+n }t
t{ -Inf -Inf FNEXTAFTER -> -INF }t
t{ 1e 2e FNEXTAFTER -> 1-up }t
t{ 0e 1e FNEXTAFTER -> min-+subn }t
t{ min-+subn 1e FNEXTAFTER -> min-+subn fdup f+ }t
t{ max-+subn 1e FNEXTAFTER -> min-+n }t
t{ max-+n-down max-+n FNEXTAFTER -> max-+n }t
t{ +nan FNEXTUP -> +NAN }t
t{ -nan FNEXTUP -> -NAN }t
t{ +nan FNEXTUP -> +NAN }t
t{ -nan FNEXTUP -> -NAN }t
t{ 0e FNEXTUP -> min-+subn }t
t{ -0e FNEXTUP -> min-+subn }t
t{ -min-+subn FNEXTUP -> -0e }t
t{ +Inf FNEXTUP -> +INF }t
t{ -Inf FNEXTUP -> -max-+n }t
t{ 1e FNEXTUP -> 1-up }t
t{ -1e FNEXTUP -> -1-up }t
t{ min-+subn FNEXTUP -> min-+subn fdup f+ }t
t{ max-+subn FNEXTUP -> min-+n }t
t{ max-+n-down FNEXTUP -> max-+n }t
t{ -max-+n FNEXTUP -> -max-+n-up }t
verbose @ [IF]
cr .( #ERRORS: ) #errors @ . cr
[THEN]
Andrew Haley
> Here's the current draft.
Sorry for the delay in replying.
> I started a new thread because the old one has become unwieldy.
> That doesn't mean that ongoing conversations in the old thread
> shouldn't continue there.
> SUMMARY OF CHANGES
> * "number" is no longer used for infinities and nans -- "IEEE
> datum" is used instead.
That's good: it is what 754 does.
> "exception": Used in the sense of IEEE 2.1.18.
> [**QUOTE
> An event that occurs when an operation on some particular
> operands has no outcome suitable for every reasonable
> application. That operation might signal one or more
> exceptions by invoking the default or, if explicitly
> requested, a language-defined alternate handling. Note that
> "event", "exception", and "signal" are defined in diverse
> ways in different programming environments.
> ]
I think that you should not use the word "exception" in a Forth
standard to refer to anything other than a Forth exception. Maybe
call them "FP exceptions"?
> The internal fp representation of default fp data, i.e., data that
> can appear on the fp stack, shall correspond to one of the IEEE
> basic, or extended, full binary formats.
I'm not sure that this standard has any business talking about the
"internal representation" of anything. It's not possible to address
into data on the fp stack, so it can't possibly affect any ANS
program. If you just mean that the system should behave as though the
representation were one of the 754 basic, or extended, full binary
formats, then you can lose the word "internal".
I think there's a need for something like
F! ( f-addr -- ) ( F: r -- )
or ( r f-addr -- )
Store r at f-addr. The data at f-addr should be in the same format as
the fp stack.
> 6 TEXT INPUT
> 6.1 Constants
> -------------
> +INF ( f: -- +Inf )
> -INF ( f: -- -Inf )
> +NAN ( f: -- +NaN )
> -NAN ( f: -- -NaN )
> These words return, respectively, IEEE signed infinity and the quiet
> signed nan with zero load, in the default format.
> See Sec. A.6.1 for more information about the encoding of NaN.
> 6.2 Decimal input
> -----------------
> IEEE requires that conversion between text and binary fp formats
> shall include signed zero, signed infinity, and signed nans, with
> and without loads. See IEEE 5.4.2, "Conversion operations for
> floating-point formats and decimal character sequences", and IEEE
> 5.12, "Details of conversion between floating-point data and
> external character sequences".
> Conversion of nan loads is not included in this specification.
> Signed infinity and signed, unloaded nans are covered by the
> constants defined in Sec. 6.1. Signed zero is already included
> in the syntax specification in DPANS94 12.3.7, "Text input
> number conversion". When IEEE-FP is present, that specification
> shall be replaced by
> Convertible string := <significand><exponent>
> <significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
> [**Note the extra .<digits> option.]
> <exponent> := <e-char>[<sign>]<digits0>
> [**Note that a sign with no digits is still
> recognized, as in DPANS94.]
> <digits> := <digit><digits0>
> <digits0> := <digit>*
> <sign> := { + | - }
> <e-char> := { E | e }
> [**Note the extra "e" option.]
> <digit> := { 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 }
> Interpretation shall convert to the default IEEE format, with
> roundTiesToEven for real numbers, and subject to overflow and
> underflow exceptions. See IEEE 7.4 for overflow and IEEE 7.5
> for underflow.
I think this, if strictly interpreted, might be too hard for some
implementations. For example, correct roundTiesToEven rounding of
"4503599627370496.50000000000000000000000000001e0" requires
multiple-precision arithmetic. The question is whether we wish to
require every Forth that supports this extension to carry aroud a
multiple-precision package for such a corner case.
> 7 GLOSSARY
> Unless otherwise stated, all fp words that do computations or
> comparisons shall obey the requirements and recommendations of
> IEEE 5 and IEEE 6, for binary formats.
> 7.1 Conversion
> --------------
> D>F ( d -- f: r )
> The DPANS94 specification is amended to require that when d
> cannot be represented precisely in the default fp format, r
> shall be the roundTiesToEven value.
This is a curious decision to make: AFAIK all fp hardware uses the
current rounding mode.
> [**IEEE and C99 background:
> IEEE does not seem to specify an equivalent. C99-n1256 6.3.1.4,
> paragraph 2 says:
> When a value of integer type is converted to a real floating
> type, if the value being converted can be represented exactly
> in the new type, it is unchanged. If the value being
> converted is in the range of values that can be represented
> but cannot be represented exactly, the result is either the
> nearest higher or nearest lower representable value, chosen in
> an implementation-defined manner. If the value being
> converted is outside the range of values that can be
> represented, the behavior is undefined.
> Note that d is always in range for any of the formats
> binary32, binary64, binary80, and binary128.
> ]
> F>D ( f: r -- s: d )
> [**Suggested specificaton #1:
> The DPANS94 specification is amended to state that an
> ambiguous condition exists, not only when the integer part of
> r is not representable by a signed, double number, but also
> when r is a nan or infinity.
> ]
> [**Suggested specificaton #2:
> The DPANS94 specification is amended to require that a Forth
> exception be thrown when r is a nan or infinity, or has
> integer part out of range for a signed, double number.
> ]
IMO definitely not: it would be Very Bad Indeed to depend on the
exception handling wordset.
> [**Suggested specificaton #3:
> The DPANS94 specification is amended to require that the
> invalid operation exception be quietly signaled when r is a
> nan or infinity, or has integer part out of range for a
> signed, double number. In that case the value of d is
> undefined.
> ]
How would anyone know that an invalid operation exception has been
quietly signaled?
> [**I prefer alternative 1. I mean, really, a "float" in the IEEE
> context is any IEEE default datum. But I would accept a variant of
> 2 if 1 should prove to be a show stopper.]
Me too. As far as I can see this is a political issue, not a
technical one.
> REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
> The DPANS94 specification is extended to include IEEE
> specials. The valid-result flag, flag2, is false if r is
> infinity or a nan.
> [**IMHO this is still unsettled, still needing discussion.
> 1. Should the current rounding mode be used, instead of
> roundTiesToEven? I'll have to check to what extent that
> has already been discussed.
I don't think there's any real need for this, but I think it would be
no great burden. I'm unsure whether this might also turn out to
require multiple-precision arithmetic if we're not careful. (Watch
out, here be dragons.)
> 2. Marcel says "It should be possible to format/output +Inf or
> +NaN, also with the existing REPRESENT." The DPANS94 spec
> for false flag2, that n and flag1 are implementation
> defined, does seem (perhaps intentionally?) well-suited to
> the spec that flag2 continues to be the negative sign flag,
> while n distinguishes among nans and infinity. For
> example, n = 0 means infinity, and n <> 0 means nan with
> implementation defined values related to quietness and
> load.
> ]
That's fine.
> [**I think an RfD for XF@ and XF!, for binary80, is in order.
Really? This is a non-standard format.
> And perhaps names like QF@ and QF!, for binary128, should be
> reserved.]
Seems sensible.
> The five required comparisons are "<", ">", "=", "<=", and ">=",
> where "<=" and ">=" stand for the usual phrases "less than or equal"
> and "greater than or equal". Note that familiar identities for real
> numbers are generally not satisfied by IEEE comparisons. For
> example, the negation of "<" is not the same as ">=". See Sec. A.7.
> F< ( f: r1 r2 -- s: [r1<r2]? )
> F= ( f: r1 r2 -- s: [r1=r2]? )
> F> ( f: r1 r2 -- s: [r1>r2]? )
> F<= ( f: r1 r2 -- s: [r1<=r2]? )
> F>= ( f: r1 r2 -- s: [r1>=r2]? )
> F0< ( f: r1 r2 -- s: [r<0]? )
> F0= ( f: r1 r2 -- s: [r=0]? )
> F0> ( f: r1 r2 -- s: [r>0]? )
> F0<= ( f: r1 r2 -- s: [r<=0]? )
> F0>= ( f: r1 r2 -- s: [r>=0]? )
OK. I think I've lost the will to fight for a smaller set of
comparisons. :-)
> 7.4 Classification
> ------------------
> IEEE 5.7.2, "General operations", requires a large number of
> classification operations. This documents defines only those
> corresponding to:
> isSignMinus
> isNormal
> isFinite
> isZero
> isSubnormal
> isInfinite
> isNaN
> isSignaling
> Actually isSignMinus corresponds to FSIGNBIT, and isZero
> corresponds to F0=, which leaves the following:
> FINITE? ( r: r -- s: [normal|subnormal]? )
> FNORMAL? ( r: r -- s: normal? )
> FSUBNORMAL? ( r: r -- s: subnormal? )
> FINFINITE? ( r: r -- s: [+|-]Inf? )
> FNAN? ( r: r -- s: nan? )
> FSIGNALING? ( r: r -- s: snan? )
I think FSIGNALING? has to go, given that the default x87 fixup for
loading a signaling NaN from memory is silently to convert it to a
quiet NaN. This makes all distinctions between quiet and signaling
NaNs useless.
> The Forth words FABS, FMAX, FMIN, and FSQRT are covered
> elsewhere.
> The DPANS94 specification for the following words is extended to
> adopt the corresponding IEEE behavior. See IEEE 9.2, "Recommended
> correctly rounded functions", and 9.2.1, "Special values".
> F** FACOS FACOSH FALOG FASIN FASINH FATAN FATAN2
> FATANH FCOS FCOSH FEXP FEXPM1 FLN FLNP1 FLOG FSIN
> FSINCOS FSINH FSQRT FTAN FTANH
No, definitely not, for reasons I explained elsewhere. Correct
rounding for these functions is fabulously expensive, and 754 doesn't
require it.
Andrew.
David N. Williams
>> SUMMARY OF CHANGES
>> * "number" is no longer used for infinities and nans -- "IEEE
>> datum" is used instead.
>
> That's good: it is what 754 does.
That's why I used it. :-)
>> "exception": Used in the sense of IEEE 2.1.18.
>> [**QUOTE
>> An event that occurs when an operation on some particular
>> operands has no outcome suitable for every reasonable
>> application. That operation might signal one or more
>> exceptions by invoking the default or, if explicitly
>> requested, a language-defined alternate handling. Note that
>> "event", "exception", and "signal" are defined in diverse
>> ways in different programming environments.
>> ]
>
> I think that you should not use the word "exception" in a Forth
> standard to refer to anything other than a Forth exception. Maybe
> call them "FP exceptions"?
I'm certainly in favor if disambiguating language. I'll try!
>> The internal fp representation of default fp data, i.e., data that
>> can appear on the fp stack, shall correspond to one of the IEEE
>> basic, or extended, full binary formats.
>
> I'm not sure that this standard has any business talking about the
> "internal representation" of anything. It's not possible to address
> into data on the fp stack, so it can't possibly affect any ANS
> program. If you just mean that the system should behave as though the
> representation were one of the 754 basic, or extended, full binary
> formats, then you can lose the word "internal".
Not such good language on my part. That's what I intended by
"correspond to". Okay, that's replaced by
Data that can appear on the fp stack shall correspond to one
of the IEEE basic, or extended, full binary formats.
> I think there's a need for something like
>
> F! ( f-addr -- ) ( F: r -- )
> or ( r f-addr -- )
>
> Store r at f-addr. The data at f-addr should be in the same format as
> the fp stack.
Maybe with an explicit disclaimer that the memory layout is not
specified. As far as I understand, all IEEE formats are logical
formats, with unspecified memory layout.
So should we make it a "should"?
>> 7 GLOSSARY
>
>> Unless otherwise stated, all fp words that do computations or
>> comparisons shall obey the requirements and recommendations of
>> IEEE 5 and IEEE 6, for binary formats.
>
>> 7.1 Conversion
>> --------------
>
>> D>F ( d -- f: r )
>
>> The DPANS94 specification is amended to require that when d
>> cannot be represented precisely in the default fp format, r
>> shall be the roundTiesToEven value.
>
> This is a curious decision to make: AFAIK all fp hardware uses the
> current rounding mode.
Okay with me to make it the current rounding mode.
They'd have to check whether an fp exception flag was raised. But
let's hope we can do alternative 1.
>> [**I prefer alternative 1. I mean, really, a "float" in the IEEE
>> context is any IEEE default datum. But I would accept a variant of
>> 2 if 1 should prove to be a show stopper.]
>
> Me too. As far as I can see this is a political issue, not a
> technical one.
I'll have to check. I don't recall that this was controversial...
>> REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
>
>> The DPANS94 specification is extended to include IEEE
>> specials. The valid-result flag, flag2, is false if r is
>> infinity or a nan.
>
>> [**IMHO this is still unsettled, still needing discussion.
>
>> 1. Should the current rounding mode be used, instead of
>> roundTiesToEven? I'll have to check to what extent that
>> has already been discussed.
>
> I don't think there's any real need for this, but I think it would be
> no great burden. I'm unsure whether this might also turn out to
> require multiple-precision arithmetic if we're not careful. (Watch
> out, here be dragons.)
I have no deep reason to prefer one rounding spec over the other.
It does seem easier to remember if the current rounding mode is
used except in special situations.
>> 2. Marcel says "It should be possible to format/output +Inf or
>> +NaN, also with the existing REPRESENT." The DPANS94 spec
>> for false flag2, that n and flag1 are implementation
>> defined, does seem (perhaps intentionally?) well-suited to
>> the spec that flag2 continues to be the negative sign flag,
>> while n distinguishes among nans and infinity. For
>> example, n = 0 means infinity, and n <> 0 means nan with
>> implementation defined values related to quietness and
>> load.
>> ]
>
> That's fine.
Maybe that makes three of us...
>> [**I think an RfD for XF@ and XF!, for binary80, is in order.
>
> Really? This is a non-standard format.
The binary format with p = 64 and emax = 13683 is an IEEE
extended format, so I wouldn't call it non-standard. It's
certainly common. The names are natural, and used by iForth. I
get the impression that it's the elephant in the room, as far as
the IEEE standard is concerned.
>> And perhaps names like QF@ and QF!, for binary128, should be
>> reserved.]
>
> Seems sensible.
Good!
>> The five required comparisons are "<", ">", "=", "<=", and ">=",
>> where "<=" and ">=" stand for the usual phrases "less than or equal"
>> and "greater than or equal". Note that familiar identities for real
>> numbers are generally not satisfied by IEEE comparisons. For
>> example, the negation of "<" is not the same as ">=". See Sec. A.7.
>
>> F< ( f: r1 r2 -- s: [r1<r2]? )
>> F= ( f: r1 r2 -- s: [r1=r2]? )
>> F> ( f: r1 r2 -- s: [r1>r2]? )
>> F<= ( f: r1 r2 -- s: [r1<=r2]? )
>> F>= ( f: r1 r2 -- s: [r1>=r2]? )
>> F0< ( f: r1 r2 -- s: [r<0]? )
>> F0= ( f: r1 r2 -- s: [r=0]? )
>> F0> ( f: r1 r2 -- s: [r>0]? )
>> F0<= ( f: r1 r2 -- s: [r<=0]? )
>> F0>= ( f: r1 r2 -- s: [r>=0]? )
>
> OK. I think I've lost the will to fight for a smaller set of
> comparisons. :-)
It comes back to, how to reserve obvious names without coercing
implementations?
There has to be *some* way to realize signaling nans in
software, as intel says. But I've suggested elsewhere that
there should be an FQUIET? . I don't object to leaving both
out.
>> The Forth words FABS, FMAX, FMIN, and FSQRT are covered
>> elsewhere.
>
>> The DPANS94 specification for the following words is extended to
>> adopt the corresponding IEEE behavior. See IEEE 9.2, "Recommended
>> correctly rounded functions", and 9.2.1, "Special values".
>
>> F** FACOS FACOSH FALOG FASIN FASINH FATAN FATAN2
>> FATANH FCOS FCOSH FEXP FEXPM1 FLN FLNP1 FLOG FSIN
>> FSINCOS FSINH FSQRT FTAN FTANH
>
> No, definitely not, for reasons I explained elsewhere. Correct
> rounding for these functions is fabulously expensive, and 754 doesn't
> require it.
For IEEE, it's only recommended. But I don't mind leaving out the
"correctly rounded" clause.
-- David
Marcel Hendrix
> Marcel Hendrix wrote:
>> [...]
>> Could you comment on the (ugly) implementation below?
[...]
> All three systems report the same glitches for FNEXTUP:
Which shows that EVERYTHING needs to be tested, always!
> ----------
> INCORRECT FP RESULT: t{ -Inf FNEXTUP -> -max-+n }t
> INCORRECT FP RESULT: t{ max-+n-down FNEXTUP -> max-+n }t
> INCORRECT FP RESULT: t{ -max-+n FNEXTUP -> -max-+n-up }t
> #ERRORS: 3
> -Inf FNEXTUP f>mix type -Infok
> max-+n-down FNEXTUP f>mix type +1.F FFFF FFFF FFFE p+1023ok
> -max-+n FNEXTUP f>mix type -1.F FFFF FFFF FFFF p+1023ok
> ----------
Here's a fix for FNEXTUP (again, only for double precision):
CREATE -max-+n PRIVATE $FF C, $FF C, $FF C, $FF C, $FF C, $FF C, $EF C, $FF C,
: FNEXTUP ( F: x -- y )
FDUP FINFINITE? IF F0< IF -max-+n DF@ ELSE +INF ENDIF EXIT ENDIF
FDUP 1e292 F+ FNEXTAFTER ;
-marcel
---
iForth version 3.0.3, generated 19:24:55, July 16, 2009.
x86_32 binary, native floating-point, double precision.
Copyright 1996 - 2009 Marcel Hendrix.
FORTH> in /idfwforth/examples/internet/idata/ieee-arith-test.frt
Redefining skip
Redefining scan
Creating --- John Hopkins TESTER Version 2.00 ---
TESTING equality of floating-point encoding
TESTING absolute tolerance
TESTING relative tolerance
TESTING F+
TESTING F-
TESTING F*
TESTING F/
TESTING FSQRT
NOT TESTING F*+
TESTING FNEXTAFTER FNEXTUP
#ERRORS: 0
ok
FORTH>
Ed
> ...
> >> 2. Marcel says "It should be possible to format/output +Inf or
> >> +NaN, also with the existing REPRESENT." The DPANS94 spec
> >> for false flag2, that n and flag1 are implementation
> >> defined, does seem (perhaps intentionally?) well-suited to
> >> the spec that flag2 continues to be the negative sign flag,
> >> while n distinguishes among nans and infinity. For
> >> example, n = 0 means infinity, and n <> 0 means nan with
> >> implementation defined values related to quietness and
> >> load.
> >> ]
> >
> > That's fine.
>
> Maybe that makes three of us...
Just to keep it balanced, here's a dissenting view :)
If you have defined the string in the buffer to be signed e.g.
"+INF" then why waste flag1 if you're never going to use it.
Similarly if you're going to use REPRESENT to display
signalling NaN's etc then the appropriate message is
already be in the buffer, and the code in 'n' is superfluous.
IMO until one has demonstrated a real use for 'n' and 'flag1'
(I don't think it has) it would be better to define them as
"reserved" for future use.
David N. Williams
> [...]
>
> Here's a fix for FNEXTUP (again, only for double precision):
>
> CREATE -max-+n PRIVATE $FF C, $FF C, $FF C, $FF C, $FF C, $FF C, $EF
C, $FF C,
>
> : FNEXTUP ( F: x -- y )
> FDUP FINFINITE? IF F0< IF -max-+n DF@ ELSE +INF ENDIF EXIT
ENDIF
> FDUP 1e292 F+ FNEXTAFTER ;
Yep, that works for me with both big- and little-endian ppc and
intel, after the following portability mods, which assume
hexfloat.fs to be loaded:
[UNDEFINED] +INF [IF]
1e 0e f/ fconstant +INF [THEN]
[UNDEFINED] FINFINITE? [IF]
: FINFINITE? ( f: r -- s: [+|-]Inf? ) FABS +INF F= ; [THEN]
\ CREATE -max-+n PRIVATE $FF C, $FF C, $FF C, $FF C,
\ $FF C, $FF C, $EF C, $FF C,
2variable -max-+n-dr
s" FF EF FF FF FF FF FF FF" hex>raw -max-+n-dr raw!
: FNEXTUP ( F: x -- y )
FDUP FINFINITE? IF F0< IF -max-+n-dr DF@ ELSE +INF ENDIF EXIT ENDIF
FDUP 1e292 F+ FNEXTAFTER ;
The HEX>RAW ... RAW! business is overkill, but it allows me not
to have to think about endianness.
I've added one FNEXTUP test, which I think makes those tests
complete:
t{ max-+n FNEXTUP -> +INF }t
Here's the full list, with the same fp constants as before:
t{ +nan FNEXTUP -> +NAN }t
t{ -nan FNEXTUP -> -NAN }t
t{ 0e FNEXTUP -> min-+subn }t
t{ -0e FNEXTUP -> min-+subn }t
t{ -min-+subn FNEXTUP -> -0e }t
t{ +Inf FNEXTUP -> +INF }t
t{ -Inf FNEXTUP -> -max-+n }t
t{ max-+n FNEXTUP -> +INF }t
t{ 1e FNEXTUP -> 1-up }t
t{ -1e FNEXTUP -> -1-up }t
t{ min-+subn FNEXTUP -> min-+subn fdup f+ }t
t{ max-+subn FNEXTUP -> min-+n }t
t{ max-+n-down FNEXTUP -> max-+n }t
t{ -max-+n FNEXTUP -> -max-+n-up }t
Humn! I'm thinking 0e and -0e in the above should be replaced
by explicit constants, just to remove any possible ambiguity
about the interpretation of signed zero.
-- David
Ed
> > Convertible string := <significand><exponent>
>
> > <significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
> > [**Note the extra .<digits> option.]
> > <exponent> := <e-char>[<sign>]<digits0>
> > [**Note that a sign with no digits is still
> > recognized, as in DPANS94.]
> > <digits> := <digit><digits0>
> > <digits0> := <digit>*
> > <sign> := { + | - }
> > <e-char> := { E | e }
> > [**Note the extra "e" option.]
> > <digit> := { 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 }
>
> > Interpretation shall convert to the default IEEE format, with
> > roundTiesToEven for real numbers, and subject to overflow and
> > underflow exceptions. See IEEE 7.4 for overflow and IEEE 7.5
> > for underflow.
>
> I think this, if strictly interpreted, might be too hard for some
> implementations. For example, correct roundTiesToEven rounding of
> "4503599627370496.50000000000000000000000000001e0" requires
> multiple-precision arithmetic. The question is whether we wish to
> require every Forth that supports this extension to carry aroud a
> multiple-precision package for such a corner case.
For some reason I had never considered such an interpretation.
Rounding, if any, ought to be according to the current system
mode. It's easiest to implement and promotes consistency.
> ...
> > [**I prefer alternative 1. I mean, really, a "float" in the IEEE
> > context is any IEEE default datum. But I would accept a variant of
> > 2 if 1 should prove to be a show stopper.]
>
> Me too. As far as I can see this is a political issue, not a
> technical one.
Sadly, yes.
> ...
> > REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
>
> > The DPANS94 specification is extended to include IEEE
> > specials. The valid-result flag, flag2, is false if r is
> > infinity or a nan.
>
> > [**IMHO this is still unsettled, still needing discussion.
>
> > 1. Should the current rounding mode be used, instead of
> > roundTiesToEven? I'll have to check to what extent that
> > has already been discussed.
>
> I don't think there's any real need for this, but I think it would be
> no great burden. I'm unsure whether this might also turn out to
> require multiple-precision arithmetic if we're not careful. (Watch
> out, here be dragons.)
IMO there is a need, and no, it's not difficult. In fact it
takes more effort to keep functions RNE and the results
would be bizarre.
Here's the output from my x87 based forth:
1 set-precision ok
set-near ok
1.5e f. 2. ok
-1.5e f. -2. ok
set-trunc ok
1.5e f. 1. ok
-1.5e f. -1. ok
set-floor ok
1.5e f. 1. ok
-1.5e f. -2. ok
set-ceil ok
1.5e f. 2. ok
-1.5e f. -1. ok
I can't imagine I'd want a forth where the results were
any different. It would be utterly confusing.
Ed
> ...
> Signed infinity and signed, unloaded nans are covered by the
> constants defined in Sec. 6.1. Signed zero is already included
> in the syntax specification in DPANS94 12.3.7, "Text input
> number conversion". When IEEE-FP is present, that specification
> shall be replaced by
>
>
> <significand> := [<sign>]{ <digits>[.<digits0>] | .<digits> }
> [**Note the extra .<digits> option.]
> <exponent> := <e-char>[<sign>]<digits0>
> [**Note that a sign with no digits is still
> recognized, as in DPANS94.]
> <digits> := <digit><digits0>
> <digits0> := <digit>*
> <sign> := { + | - }
> <e-char> := { E | e }
> [**Note the extra "e" option.]
> <digit> := { 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 }
>
> Interpretation shall convert to the default IEEE format, with
> roundTiesToEven for real numbers, and subject to overflow and
> underflow exceptions. See IEEE 7.4 for overflow and IEEE 7.5
> for underflow.
Not sure which IEEE standard document you are using as your
reference but I assume they all say much the same.
The IEEE strategy of allowing one to accept "overflow" input
and then substitute it with various "good results" is disconcerting
as it weakens the ability of programs to discriminate. I assume
the same 'substitution strategy' applies to your IEEE >FLOAT
string input function (?)
The more I look at IEEE, the less I'm liking it ...
Marcel Hendrix
> Andrew Haley wrote:
>> David N. Williams <will...@umich.edu> wrote:
[...]
>> > REPRESENT ( f: r s: c-addr u -- n flag1 flag2 )
>> > The DPANS94 specification is extended to include IEEE
>> > specials. The valid-result flag, flag2, is false if r is
>> > infinity or a nan.
>> > [**IMHO this is still unsettled, still needing discussion.
>> > 1. Should the current rounding mode be used, instead of
>> > roundTiesToEven? I'll have to check to what extent that
>> > has already been discussed.
>> I don't think there's any real need for this, but I think it would be
>> no great burden. I'm unsure whether this might also turn out to
>> require multiple-precision arithmetic if we're not careful. (Watch
>> out, here be dragons.)
> IMO there is a need, and no, it's not difficult. In fact it
> takes more effort to keep functions RNE and the results
> would be bizarre.
Well, the current standard already explicitly specifies
round-to-nearest for REPRESENT.
> Here's the output from my x87 based forth:
> 1 set-precision ok
> set-near ok
> 1.5e f. 2. ok
[..]
Very good. But what happens if I boot your Forth and do
1 set-precision ok
1.5e f. ( 2. ? )
I mean, does it set the rounding mode to RNE on boot? (I guess it
has to, in order to provide consistent results and be ANS standard.)
Under what conditions does the system reset the rounding mode? ( ABORT,
QUIT, THROW, CATCH, -1 @, system call, callback, task start/stop/pause
etc..) I'd expect the system to document what it does under these
conditions, especially for the "etc." cases.
A more philosophical / tangential question is how .S should show
the FP stack contents (should it respect PRECISION and rounding
mode)? set-trunc PI F. does not mean the value of PI is 3.
-marcel
Ed
> > Here's the output from my x87 based forth:
>
> > 1 set-precision ok
> > set-near ok
> > 1.5e f. 2. ok
> [..]
>
> Very good. But what happens if I boot your Forth and do
>
> 1 set-precision ok
> 1.5e f. ( 2. ? )
>
> I mean, does it set the rounding mode to RNE on boot? (I guess it
> has to, in order to provide consistent results and be ANS standard.)
Yes. IIRC, the x87 FINIT instruction sets up RNE mode and I think
it's an IEEE boot-up requirement.
SwiftForth lets you boot up in any mode as it saves the rounding
mode in the registry.
> Under what conditions does the system reset the rounding mode? ( ABORT,
> QUIT, THROW, CATCH, -1 @, system call, callback, task start/stop/pause
> etc..) I'd expect the system to document what it does under these
> conditions, especially for the "etc." cases.
SwiftForth has a deferred command /NDP which resets rounding
to the boot-up default. Looking at the source code, it appears to
be executed under various abort conditions.
I wouldn't know whether this is optimum as I don't have any
experience in the area. SwiftForth is Windows, whereas I use
MS-DOS. The needs may be different.
> A more philosophical / tangential question is how .S should show
> the FP stack contents (should it respect PRECISION and rounding
> mode)? set-trunc PI F. does not mean the value of PI is 3.
Like the F. example, I would feel it confusing not to use the
same rounding mode. However it depends on your setup and
personal preference.
On my system the rounding isn't reset by ABORT, so it won't
change "behind my back". My .S uses G. (from FPOUT). It
displays up to PRECISION digits using the current rounding.
David N. Williams
FLOGB and FSCALBN. I've also included Marcel's version, which
uses FNEXTAFTER.
All three have been tested with an experimental pfe
implementation of part of the proposed extension, called
ieeefp-ext.c. I'll post it on the web later, after some
cleanup.
Default floats are doubles. The fp stack is separate.
The first new implementation uses IEEE data-type testing words
from the current draft.
The second shows that FCLASSIFY can be helpful, IMHO. FCLASSIFY
is not in the draft, but I'm wondering if it shouldn't be put
back in. It seems intended for case structures, which I
understand are heavily optimized in C (also in Forth?), whereas
the individual tests would seem better for just a few cases.
Should we reconsider having both?
In the pfe module, the FCLASSIFY enumerations are built-in, but
they can be defined with FCLASSIFY by phrases like:
+nan fclassify constant FP-NAN
FNEXTAFTER is included in the pfe module from its former
incarnation in fpieee-ext.c
-- David
P.S. I see there's a two's-complement dependency in the use of
FLOGB...
--- fnextup.fs -------------------------------------------------
( Title: FNEXTUP tryouts with pfe ieeefp-ext.c
File: fnextup.fs
Version: 0.1.0
Started: July 18, 2009
Revised: July 19, 2009
Authors: Marcel Hendrix, David N. Williams
License: Public Domain?
)
loadm ieeefp \ pfe experimental module
s" hexfloat.fs" included \ fails if bits/cell <> 32
decimal
bits/float 64 <>
[IF] cr .( ***This code requires 64 bits/float.) ABORT [THEN]
[UNDEFINED] \\ [IF] \ for degugging
: \\ ( -- ) -1 parse 2drop BEGIN refill 0= UNTIL ; [THEN]
0e fabs fconstant +0
+0 fnegate fconstant -0
0e 0e f/ fabs fconstant +nan
+nan fnegate fconstant -nan
1e +0 f/ fconstant +inf
+inf fnegate fconstant -inf
s" 1.F FFFF FFFF FFFF p+1023" mix>f fconstant max-+n
s" 1.0 0000 0000 0000 p-1022" mix>f fconstant min-+n
s" 0.F FFFF FFFF FFFF p-1022" mix>f fconstant max-+subn
s" 0.0 0000 0000 0001 p-1022" mix>f fconstant min-+subn
s" -0.0 0000 0000 0001 p-1022" mix>f fconstant -min-+subn
s" -0.F FFFF FFFF FFFF p-1022" mix>f fconstant -max-+subn
s" -1.0 0000 0000 0000 p-1022" mix>f fconstant -min-+n
s" -1.F FFFF FFFF FFFF p+1023" mix>f fconstant -max-+n
s" 1.0 0000 0000 0001 p0" mix>f fconstant 1-up
s" -1.F FFFF FFFF FFFF p-1" mix>f fconstant -1-up
s" 1.F FFFF FFFF FFFE p+1023" mix>f fconstant max-+n-down
s" -1.F FFFF FFFF FFFE p+1023" mix>f fconstant -max-+n-up
s" 1p-52" mix>f fconstant feps
s" 1p-53" mix>f fconstant fepsneg
true [IF]
2variable -max-+n-dr
s" FF EF FF FF FF FF FF FF" hex>raw -max-+n-dr raw!
: FNEXTUP ( F: x -- y ) \ mh 18Jul09, slightly revised dnw 18Jul09
FDUP FINFINITE?
IF F0<
IF -max-+n-dr DF@ ELSE +INF ENDIF EXIT
ENDIF
\ FDUP 1e292 F+ FNEXTAFTER ;
+INF FNEXTAFTER ;
[THEN]
false [IF]
: FNEXTUP ( F: x -- y ) \ dnw 18Jul09
fdup f0= IF min-+subn f+ EXIT THEN
fdup finfinite? IF f0< IF -max-+n ELSE +inf THEN EXIT THEN
fdup fnan? IF EXIT THEN
fdup -min-+subn f= IF fdrop -0 EXIT THEN
fdup fsubnormal? IF min-+subn f+ EXIT THEN
fdup flogb f>d drop
fdup f0< IF fepsneg ELSE feps THEN fscalbn f+ ;
[THEN]
false [IF]
: FNEXTUP ( F: x -- y ) \ dnw 18Jul09
fdup fclassify
CASE
FP-ZERO OF fdrop min-+subn ENDOF
FP-NORMAL OF fdup flogb f>d drop
fdup fsignbit
IF fepsneg ELSE feps THEN fscalbn f+ ENDOF
FP-INFINITE OF fsignbit IF -max-+n ELSE +inf THEN ENDOF
FP-SUBNORMAL OF fdup -min-+subn f=
IF fdrop -0 ELSE min-+subn f+ THEN ENDOF
( FP-NAN)
ENDCASE ;
[THEN]
s" ttester.fs" included
true verbose !
set-exact
decimal
variable #errors 0 #errors !
:noname ( c-addr u -- )
(
Display an error message followed by the line that had the
error.
)
1 #errors +! error1 ; error-xt !
t{ 1.23e FDUP 1e FNEXTAFTER F< -> 0 }t
t{ 1.23e FDUP 1e FNEXTAFTER F> -> -1 }t
t{ 1.23e FDUP 2e FNEXTAFTER F> -> 0 }t
t{ 1.23e FDUP 2e FNEXTAFTER F< -> -1 }t
t{ +nan +nan FNEXTAFTER -> +nan }t
t{ -nan +nan FNEXTAFTER fabs -> +nan }t
t{ +nan -nan FNEXTAFTER fabs -> +nan }t
t{ -nan -nan FNEXTAFTER -> -nan }t
t{ 1e FDUP 1e FNEXTAFTER F= -> -1 }t
t{ +0 +1e FNEXTAFTER -> min-+subn }t
t{ +0 -1e FNEXTAFTER -> -min-+subn }t
t{ +Inf +Inf FNEXTAFTER -> +INF }t
t{ +Inf +0 FNEXTAFTER -> max-+n }t
t{ -Inf +0 FNEXTAFTER -> -max-+n }t
t{ -Inf -Inf FNEXTAFTER -> -INF }t
t{ 1e 2e FNEXTAFTER -> 1-up }t
t{ +0 1e FNEXTAFTER -> min-+subn }t
t{ min-+subn 1e FNEXTAFTER -> min-+subn fdup f+ }t
t{ max-+subn 1e FNEXTAFTER -> min-+n }t
t{ max-+n-down max-+n FNEXTAFTER -> max-+n }t
t{ +nan FNEXTUP -> +nan }t
t{ -nan FNEXTUP -> -nan }t
t{ +0 FNEXTUP -> min-+subn }t
t{ -0 FNEXTUP -> min-+subn }t
t{ -min-+subn FNEXTUP -> -0 }t
t{ +Inf FNEXTUP -> +INF }t
t{ -Inf FNEXTUP -> -max-+n }t
t{ max-+n FNEXTUP -> +INF }t
t{ 1e FNEXTUP -> 1-up }t
t{ -1e FNEXTUP -> -1-up }t
t{ min-+subn FNEXTUP -> min-+subn fdup f+ }t
t{ max-+subn FNEXTUP -> min-+n }t
t{ max-+n-down FNEXTUP -> max-+n }t
t{ -max-+n FNEXTUP -> -max-+n-up }t
verbose @ [IF]
Andrew Haley
> Andrew Haley wrote:
>> David N. Williams <will...@umich.edu> wrote:
>> ...
>> I think this, if strictly interpreted, might be too hard for some
>> implementations. For example, correct roundTiesToEven rounding of
>> "4503599627370496.50000000000000000000000000001e0" requires
>> multiple-precision arithmetic. The question is whether we wish to
>> require every Forth that supports this extension to carry aroud a
>> multiple-precision package for such a corner case.
> For some reason I had never considered such an interpretation.
>
> Rounding, if any, ought to be according to the current system
> mode.
That doesn't really answer the question, which is whether decimal
floating-point input conversion must require prefectly correct
rounding in all cases. If it's not required, then the direction
doesn't much matter, since the difference is at most 1 ulp. If
correct rounding is required, then rounding correctly at the end is a
relatively small problem.
I'm concerned that this specification not end up requiring
multiple-precision arithmetic. I think that would be an unfortunate
(and, I suspect, unintended) side-effect. Maybe the reference
implementation will show that it's not a huge burden in this case,
though.
> It's easiest to implement and promotes consistency.
I don't think it's easier to implement: we can't use the hardware's
rounding mode when shortening decimal input to the precision of the
FPU.
Andrew.
Andrew Haley
> ------------------------
> FMAX ( f: r1 r2 -- r3 )
> FMIN ( f: r1 r2 -- r3 )
> The DPANS94 specification is extended to IEEE specials. See
> minNum and maxNum in IEEE 5.3.1, "General operations" and 6.2,
> "Operations with NaNs".
> FNEXTUP ( f: r1 -- r2 )
> When r1 is a nonzero real number, FNEXTUP returns the next
> affinely extended real in the default format that compares
> larger than r1. See IEEE 5.3.1, "General operations" for the
> behavior when r1 is an IEEE special. FNEXTDOWN is not
> defined. According to IEEE, nextDown(x) is -nextUp(-x).
> FSCALBN ( f: r s: n -- f: r*2^n )
> The output is efficiently scaled by 2^n. See IEEE 5.3.3,
> "logBFormat operations".
> FLOGB ( f: r -- e )
> Leave the radix-two exponent e of the fp representation as an
> fp integer. If r is subnormal, the exponent is computed as if
> r were normalized, with e < emin. See IEEE 5.3.3, "logBFormat
> operations" for treatment of IEEE specials.
We still seem to be missing FREMAINDER:
FREMAINDER ( f: x y -- r q )
When y is not 0, the remainder r = fremainder(x, y) is defined for
finite x and y regardless of the rounding-direction attribute by the
mathematical relation r = x - y * q , where q is the integer nearest
the exact number x/y ; whenever | n - x/y | = 1/2 , then n is even.
Thus, the remainder is always exact. If r = 0, its sign shall be that
of x. fremainder(x, inf) is x for finite x.
Although many floating-point applications don't use fremainder, it's a
useful factor of several standard floating-point functions.
Andrew.
Andrew Haley
The IEEE approach is to give the programmer control over the way that
exceptional conditions are handled but provide default fixups in all
cases. It always makes a best effort to do something sensible. This
is a major philosophical difference from the "if anything untoward
happens, bail out and return an error" approach.
> I assume the same 'substitution strategy' applies to your IEEE
> >FLOAT string input function (?)
> The more I look at IEEE, the less I'm liking it ...
Maybe it's not right for you.
Andrew.
Andrew Haley
Sorry, typo'd:
When y is not 0, the remainder r = fremainder(x, y) is defined for
finite x and y regardless of the rounding-direction attribute by the
mathematical relation r = x - y * q , where q is the integer nearest
the exact number x/y ; whenever | q - x/y | = 1/2 , then q is even.
Thus, the remainder is always exact. If r = 0, its sign shall be that
of x. fremainder(x, inf) is x for finite x.
Andrew.
Anton Ertl
>I think that you should not use the word "exception" in a Forth
>standard to refer to anything other than a Forth exception. Maybe
>call them "FP exceptions"?
I think that's still too easily confused. The first thing I think of
when I read this is what happens when something throws one of the
varlues -55, -54, -46..-41.
Either make an acronym of it (e.g., IFENFE for "IEEE FP exception, not
Forth exception"), or use a different word, e.g, condition.
>For example, correct roundTiesToEven rounding of
>"4503599627370496.50000000000000000000000000001e0" requires
>multiple-precision arithmetic.
I don't think so. You just need to remember whether you are sitting
on the fence. If there's anythink more, that pushes you on one side
of the fence, otherwise you just RNE the result. The same happens
with e.g., multiplication.
Anton Ertl
>Anton Ertl <an...@mips.complang.tuwien.ac.at> wrote:
>> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
>> >While it might be nice to have, say, a correctly rounded cosine, it
>> >would be a lot of work and we wouldn't be able to use the versions
>> >provided by the C library. Also, the functions would have difficult
>> >to predict execution time.
>
>> Oh, I see what you are talking about. No, I was not thinking about
>> bit-accurate results (which are called "correctly rounded" by IEEE
>> 754 for some reason).
>
>> What I was thinking about was that the rounding in the final step(s)
>> of words like FCOS should be done according to the current roundiung
>> mode. That should be helpful in achieving the kind of error
>> estimate I expect that non-RNE rounding would be used for.
>
>I have no problem with that, but since FCOS isn't going to be
>correctly rounded I can see no way to specify such behaviour.
Good point.
>> It would probably be even better if there was an implementation of
>> these words that guaranteed that the result would be in the right
>> direction from the exact result for the floored, ceiling, and
>> towards-zero rounding modes.
>
>Maybe, but that requires 0.5ulp precision
No. Why do you think so?
The result can be off by more then 1 ulp (say, 5 ulp), it just has to
be off in the right direction. You don't need 0.5ulp precision for
that, just an algorithm that's deviates in the right direction, or
that knows whether it is above or below the exact result.
>> OTOH, I wonder if we should standardize non-RNE rounding modes at
>> the moment at all. How many potential users are there? They should
>> tell us (possibly through the system implementors) what they need.
>> If nobody cares, there's no need to work on that.
>
>It's not just exotic cases. For example, if you want to do exact
>integer arithmetic without transferring operands into integer
>registers, then you need to set roundToZero mode. This is useful, for
>example, for decimal output.
Ok, that's a different use, and for that you don't need the
transcendental functions. So I guess the rounding modes are useful
even if we don't get the transcendental functions right. Getting the
transcendentals right is quite a burden on the implementation, so with
that use in mind, we should not require it.
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >I think that you should not use the word "exception" in a Forth
> >standard to refer to anything other than a Forth exception. Maybe
> >call them "FP exceptions"?
> I think that's still too easily confused. The first thing I think of
> when I read this is what happens when something throws one of the
> varlues -55, -54, -46..-41.
> Either make an acronym of it (e.g., IFENFE for "IEEE FP exception,
> not Forth exception"), or use a different word, e.g, condition.
Whatever works, but we also need to make things reasonably sane for
people more familiar with IEEE arithmetic than with Forth.
> >For example, correct roundTiesToEven rounding of
> >"4503599627370496.50000000000000000000000000001e0" requires
> >multiple-precision arithmetic.
> I don't think so. You just need to remember whether you are sitting
> on the fence. If there's anythink more, that pushes you on one side
> of the fence, otherwise you just RNE the result.
Maybe, but I'm not convinced it's even the easiest way. I don't think
you can do the RNE until you have multiplied the decimal fraction by
an appropriate power of ten, and I think that you have to maintain
full precision until you have done that. I'm not going to insist that
what you suggest is actually impossible.
> The same happens with e.g., multiplication.
It's not quite the same, since multiplication requires at most 2N bits
of an intermediate result. But yes, quite a few operations do require
extra intermediate precision.
Andrew.
Anton Ertl
>Anton Ertl <an...@mips.complang.tuwien.ac.at> wrote:
>> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
>> >For example, correct roundTiesToEven rounding of
>> >"4503599627370496.50000000000000000000000000001e0" requires
>> >multiple-precision arithmetic.
>
>> I don't think so. You just need to remember whether you are sitting
>> on the fence. If there's anythink more, that pushes you on one side
>> of the fence, otherwise you just RNE the result.
>
>Maybe, but I'm not convinced it's even the easiest way. I don't think
>you can do the RNE until you have multiplied the decimal fraction by
>an appropriate power of ten, and I think that you have to maintain
>full precision until you have done that. I'm not going to insist that
>what you suggest is actually impossible.
Hmm, I think your example was bad. For your example it's easy, the
hardest case is when you have a 53+-digit number (for 53-bit
mantissae), and falling on either side of the fence is possible until
the 53rd digit. If you are still on the fence at that point, the case
degenerates to the one above.
Doesn't "How to read floating point numbers accurately" solve this
problem? (I had a look at it some time ago, and it looked quite
messy).
>> The same happens with e.g., multiplication.
>
>It's not quite the same, since multiplication requires at most 2N bits
>of an intermediate result.
AFAIK three extra bits are sufficient for RNE, and that's what the
8087 has (i.e., a 67-bit multiplier for 64-bit mantissae, with special
handling of the last three bits).
Andrew Haley
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >Anton Ertl <an...@mips.complang.tuwien.ac.at> wrote:
> >> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
> >> >For example, correct roundTiesToEven rounding of
> >> >"4503599627370496.50000000000000000000000000001e0" requires
> >> >multiple-precision arithmetic.
> >
> >> I don't think so. You just need to remember whether you are sitting
> >> on the fence. If there's anythink more, that pushes you on one side
> >> of the fence, otherwise you just RNE the result.
> >
> >Maybe, but I'm not convinced it's even the easiest way. I don't
> >think you can do the RNE until you have multiplied the decimal
> >fraction by an appropriate power of ten, and I think that you have
> >to maintain full precision until you have done that. I'm not going
> >to insist that what you suggest is actually impossible.
> Hmm, I think your example was bad.
It's an example that illustrates a case in which simply accumulating
the digits won't do; there are many other such numbers.
> For your example it's easy, the hardest case is when you have a
> 53+-digit number (for 53-bit mantissae), and falling on either side
> of the fence is possible until the 53rd digit. If you are still on
> the fence at that point, the case degenerates to the one above.
> Doesn't "How to read floating point numbers accurately" solve this
> problem? (I had a look at it some time ago, and it looked quite
> messy).
That paper shows that it's necessary to use a little extra precision
in most cases and a lot more extra precision in a few cases. "When
using 64-bit ... to compute IEEE double results ... the algorithm
avoids higher-precision arithmetic 99% of the time."
> >> The same happens with e.g., multiplication.
> >
> >It's not quite the same, since multiplication requires at most 2N
> >bits of an intermediate result.
> AFAIK three extra bits are sufficient for RNE, and that's what the
> 8087 has (i.e., a 67-bit multiplier for 64-bit mantissae, with special
> handling of the last three bits).
Hence "at most".
But how do you get those three bits? It requires a bunch of special
logic. For RNE, the sticky bit is the OR of all the additional bits
in the result. One technique for calculating it observes that the
number of trailing 0s in the result is the sum of the number of
trailing 0s in the operands. So, the sticky bit is set to 1 iff that
sum is less than the number of bits that are to be discarded.
Counting strings of bits is fairly hard in software, so you might well
choose to generate 2N bits of the result.
Andrew.
Ed
Rounding happens automatically by virtue of the mathematic
functions F* F+ etc. used in the numeric conversion process.
It's hard to imagine IEEE would require anything that demanded
multi-precision arithmetic as it's too resource hungry.
Ed
I believe the programmer ought to stay in control.
> > I assume the same 'substitution strategy' applies to your IEEE
> > >FLOAT string input function (?)
>
> > The more I look at IEEE, the less I'm liking it ...
>
> Maybe it's not right for you.
Maybe.
I may not be personally enamoured by some of IEEE's more
strange dictums, but that doesn't mean an IEEE Forth needs
to be broken. i.e. it should still be possible to run '94 Standard
(non-IEEE) programs and applications without compromise.
AFAIK, IEEE defines certain behaviours but doesn't dictate
how Forth shall implement them. Those mistakes are reserved
for those defining this IEEE Forth extension :)
Ed
> ...
> A more philosophical / tangential question is how .S should show
> the FP stack contents (should it respect PRECISION and rounding
> mode)? set-trunc PI F. does not mean the value of PI is 3.
I've since checked SwiftForth. It's .S displays f/p numbers with a
precision of 8, regardless of PRECISION.
Rounding is currently broken on SwiftForth output display functions.
But were it not, .S would use the current rounding mode.
I personally don't like the idea that rounding can be reset by ABORT.
I recall an LMI forth that reset to DECIMAL on error. It drove me
nuts because I have the habit of using junk input to clear the stack :(
David N. Williams
> Here are two more implementations of FNEXTUP, both of which use
> FLOGB and FSCALBN. I've also included Marcel's version, which
> uses FNEXTAFTER.
>
> All three have been tested with an experimental pfe
> implementation of part of the proposed extension, called
> ieeefp-ext.c. I'll post it on the web later, after some
> cleanup.
[...]
Unfortunately, there's a serious bug in my versions of FNEXTUP,
and I'm sorry to say I found it by accident, not by systematic
testing. :-(
Now that I realize what's going on, it rings a bell -- I think
I've actually heard of the problem before.
Here's a new version that passes tests covering more corners.
: FNEXTUP ( f: x -- y )
(
When x is an affine extended real, y is the next affine extended
real that compares larger than x.
For normals, a loss of precision at negative numbers that are
powers of two requires special treatment. For -MIN-+N, rounding
by FSCALBN also has to be avoided.
)
fdup fclassify
CASE
FP-ZERO OF fdrop min-+subn ENDOF
FP-NORMAL OF fdup flogb f>d drop \ 2's complement dependency
( n) dup negate fscalbn
fdup -1e f=
IF fdrop ( n) dup -1022 = ( [x = -MIN-+N]?)
IF ( n) drop -max-+subn ELSE [-1]-up fscalbn THEN
ELSE
feps f+ fscalbn
THEN ENDOF
FP-INFINITE OF fsignbit IF -max-+n ELSE +inf THEN ENDOF
FP-SUBNORMAL OF fdup -min-+subn f=
IF fdrop -0 ELSE min-+subn f+ THEN ENDOF
( FP-NAN)
ENDCASE ;
This has been tested with pfe.
The more complete set of tests is included in the following
experimental file, with portable versions of Marcel's FNEXTAFTER
and FNEXTUP:
http://www-personal.umich.edu/~williams/archive/temp/fnextup-mh.fs
I don't promise the tests are definitive, because I need to
recheck my pencil and paper analysis.
There's a trivial fix in FNEXTAFTER, to make it pass these:
\ x y y
t{ +0 -0 FNEXTAFTER -> -0 }t
t{ -0 +0 FNEXTAFTER -> +0 }t
The point is that C99 requires y to be returned when x = y. The
nextafter function isn't in IEEE 754-2008 -- it's replaced by
nextup and nextdown, an improvement IMHO.
One of the other new tests for FNEXTAFTER fails:
t{ -max-+n -inf FNEXTAFTER -> -inf }t
I couldn't quickly see how to fix it. And it doesn't matter for
FNEXTUP.
I promise I'll get back to editing the draft... :-)
-- David
Ed
It may not be as bad as that.
Re-reading the section on IEEE overflows etc again, I note that
as well as providing default results, one is supposed to "signal"
the exception. How will an IEEE Forth signal exceptions?
In '94, I/O exceptions are signalled with an "ior" result from the
relevant function. The text interpreter doesn't have that - it
simply works on a pass/fail result. So if I enter the following
in immediate mode:
1e9000
should it signal the overflow exception - or quietly return an INF
(or whatever) ? Similarly should calculations entered in immediate
mode which result in divide-by-zero signal the exception?
IMO before one can define IEEE I/O functions, one must first
know how IEEE exceptions are to be handled. If one is going to
use ior values, then that will inform the stack effect of IEEE I/O
functions. If one plans to use e.g. >FLOAT for IEEE I/O then a
different method of handling exceptions will be needed.
Marcel Hendrix
[..]
> The more complete set of tests is included in the following
> experimental file, with portable versions of Marcel's FNEXTAFTER
> and FNEXTUP:
> http://www-personal.umich.edu/~williams/archive/temp/fnextup-mh.fs
> I don't promise the tests are definitive, because I need to
> recheck my pencil and paper analysis.
> There's a trivial fix in FNEXTAFTER, to make it pass these:
> \ x y y
> t{ +0 -0 FNEXTAFTER -> -0 }t
> t{ -0 +0 FNEXTAFTER -> +0 }t
> The point is that C99 requires y to be returned when x = y. The
> nextafter function isn't in IEEE 754-2008 -- it's replaced by
> nextup and nextdown, an improvement IMHO.
It says to "return x" in http://www.netlib.org/fdlibm/s_nextafter.c ...
> One of the other new tests for FNEXTAFTER fails:
> t{ -max-+n -inf FNEXTAFTER -> -inf }t
> I couldn't quickly see how to fix it. And it doesn't matter for
> FNEXTUP.
Fixed by special casing.
-marcel
-- ----------------------------------------------
FORTH> needs -ieee-arith-test
Creating --- Hexfloat Version 0.15 ---
Creating --- John Hopkins TESTER Version 2.00 ---
TESTING equality of floating-point encoding
TESTING absolute tolerance
TESTING relative tolerance
TESTING F+
TESTING F-
TESTING F*
TESTING F/
TESTING FSQRT
NOT TESTING F*+
TESTING FNEXTAFTER (legacy tests)
testing FNEXTAFTER nan propagation
testing FNEXTAFTER equality
testing FNEXTAFTER downward
testing FNEXTAFTER upward
testing FNEXTUP
#ERRORS: 0
ok
Andrew Haley
By setting the floating-point status, I think.
> In '94, I/O exceptions are signalled with an "ior" result from the
> relevant function. The text interpreter doesn't have that - it
> simply works on a pass/fail result. So if I enter the following
> in immediate mode:
>
> 1e9000
>
> should it signal the overflow exception
Interesting question. I think so, yes. But that requires the text
interpreter to leave the floating-point status untouched untili a user
reads it.
> - or quietly return an INF (or whatever) ? Similarly should
> calculations entered in immediate mode which result in
> divide-by-zero signal the exception?
Again, I think so.
> IMO before one can define IEEE I/O functions, one must first know
> how IEEE exceptions are to be handled. If one is going to use ior
> values, then that will inform the stack effect of IEEE I/O
> functions. If one plans to use e.g. >FLOAT for IEEE I/O then a
> different method of handling exceptions will be needed.
Different from what, exactly?
Andrew.
Ed
Since >FLOAT only returns a flag - not an ior - it would need to
update a status variable instead. Not that I'm recommending this.
David N. Williams
> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP
0.5.2 -- FNEXTUP
> [..]
>
>> There's a trivial fix in FNEXTAFTER, to make it pass these:
>
>> \ x y y
>> t{ +0 -0 FNEXTAFTER -> -0 }t
>> t{ -0 +0 FNEXTAFTER -> +0 }t
>
>> The point is that C99 requires y to be returned when x = y. The
>> nextafter function isn't in IEEE 754-2008 -- it's replaced by
>> nextup and nextdown, an improvement IMHO.
>
> It says to "return x" in http://www.netlib.org/fdlibm/s_nextafter.c ...
Indeed it does. It's good to know that resource! So Sun has
(or had) a discrepancy from what I understand to be official
C99, quoted at the end. Actually gcc on Mac OS agrees with C99,
at least for ppc.
>> One of the other new tests for FNEXTAFTER fails:
>
>> t{ -max-+n -inf FNEXTAFTER -> -inf }t
>
>> I couldn't quickly see how to fix it. And it doesn't matter for
>> FNEXTUP.
>
> Fixed by special casing.
Is your fix posted somewhere at your website?
-- David
REFERENCES
[6] ISO/IEC 9899:1999 (December 1, 1999),
ISO/IEC 9899:1999 Cor. 1:2001(E),
ISO/IEC 9899:1999 Cor. 2:2004(E),
ISO/IEC 9899:1999 Cor. 3:2007(E):
http://www.open-std.org/jtc1/sc22/wg14/www/standards.html#9899
[7] C99 + TC1 + TC2 + TC3 is included in the freely available
WG14/N1256, September 7, 2007 [**thanks to David Thompson]:
http://www.open-std.org/jtc1/sc22/WG14/www/docs/n1256.pdf
-----------
WG14/N1256 Committee Draft � Septermber 7, 2007 ISO/IEC 9899:TC3
7.12.11.3 The nextafter functions
Synopsis
#include <math.h>
double nextafter(double x, double y);
float nextafterf(float x, float y);
long double nextafterl(long double x, long double y);
Description
The nextafter functions determine the next representable value,
in the type of the function, after x in the direction of y,
where x and y are first converted to the type of the
function.211) The nextafter functions return y if x equals y. A
range error may occur if the magnitude of x is the largest
finite value representable in the type and the result is
infinite or not representable in the type.
----------
Andrew Haley
> FLOGB and FSCALBN. I've also included Marcel's version, which
> uses FNEXTAFTER.
>
> All three have been tested with an experimental pfe
> implementation of part of the proposed extension, called
> ieeefp-ext.c. I'll post it on the web later, after some
> cleanup.
>
> Default floats are doubles. The fp stack is separate.
>
> The first new implementation uses IEEE data-type testing words
> from the current draft.
>
> The second shows that FCLASSIFY can be helpful, IMHO. FCLASSIFY
> is not in the draft, but I'm wondering if it shouldn't be put
> back in. It seems intended for case structures, which I
> understand are heavily optimized in C (also in Forth?), whereas
> the individual tests would seem better for just a few cases.
> Should we reconsider having both?
This all looks like fantastically hard work! If we can get at the
bits of a floating-point number it all becomes a lot easier. This is
an example on a little-endian 32-bit Forth:
: nan? ( - t) ( f) fdup f= 0= ;
2variable foo
\ little-endian 2@ and 2!
: d! >r swap r> 2! ;
: d@ 2@ swap ;
: fnextup ( f: x - y)
fdup nan? if exit then
0.0e f+ \ Convert -0 to +0
foo df!
foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
foo df@ ;
I may have missed a corner case, but you get the idea.
Andrew.
David N. Williams
>> FLOGB and FSCALBN. I've also included Marcel's version, which
>> uses FNEXTAFTER.
>
> This all looks like fantastically hard work! If we can get at the
> bits of a floating-point number it all becomes a lot easier. This is
> an example on a little-endian 32-bit Forth:
>
> : nan? ( - t) ( f) fdup f= 0= ;
> 2variable foo
>
> \ little-endian 2@ and 2!
> : d! >r swap r> 2! ;
> : d@ 2@ swap ;
>
> : fnextup ( f: x - y)
> fdup nan? if exit then
> 0.0e f+ \ Convert -0 to +0
> foo df!
> foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
> foo df@ ;
>
> I may have missed a corner case, but you get the idea.
Well, I was skeptical, so I put it to the test:
--------------
[jost:~/forth/IEEE-FP] dnwills% pfe fnextup-ah.fs
\ Portable Forth Environment 0.33.61 (Aug 21 2008 11:00:57)
Copyright (C) Dirk Uwe Zoller 1993 - 1995.
Copyright (C) Tektronix, Inc. 1998 - 2003.
ANS/ffa ITC Forth - Please enter LICENSE and WARRANTY.
testing FNEXTUP
INCORRECT FP RESULT: t{ +Inf FNEXTUP -> +INF }t
#ERRORS: 1
--------------
Same result with gforth, and with big- or little-endian
Mac OS X.
Wonderful! Now I have to teach myself why any carries and
borrows into the exponent field are exactly what's needed.
And that it'll still work for binary80.
Unless my tests are incomplete. Either way I have something to
understand.
I modified your code to work with either endianness, with the
help of hexfloat.fs, which does give the kind of access to fp
bits that you used. The code and the tests are below.
-- David
---------------
\ Andrew Haley nntp://comp.lang.forth 23Jul09
\ Modified for portability by dnw, 23Jul09
s" hexfloat.fs" included
bits/float 64 <>
[IF] cr .( ***This code requires 64 bits/float.) ABORT [THEN]
: nan? ( - t) ( f) fdup f= 0= ;
0 [IF]
2variable foo
\ little-endian 2@ and 2!
: d! >r swap r> 2! ;
: d@ 2@ swap ;
: fnextup ( f: x - y)
fdup nan? if exit then
0.0e f+ \ Convert -0 to +0
foo df!
foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
foo df@ ;
[THEN]
: FNEXTUP ( f: x - y)
fdup nan? if exit then
0.0e f+ \ Convert -0 to +0
f>raw dup 0< if 1 0 d- else 1 0 d+ then raw>f ;
s" ttester.fs" included
true verbose !
set-exact
decimal
variable #errors 0 #errors !
:noname ( c-addr u -- )
(
Display an error message followed by the line that had the
error.
)
1 #errors +! error1 ; error-xt !
0e fabs fconstant +0
+0 fnegate fconstant -0
0e 0e f/ fabs fconstant +nan
+nan fnegate fconstant -nan
1e +0 f/ fconstant +inf
+inf fnegate fconstant -inf
s" 1.F FFFF FFFF FFFF p+1023" mix>f fconstant max-+n
s" 1.0 0000 0000 0000 p-1022" mix>f fconstant min-+n
s" 0.F FFFF FFFF FFFF p-1022" mix>f fconstant max-+subn
s" 0.0 0000 0000 0001 p-1022" mix>f fconstant min-+subn
s" -0.0 0000 0000 0001 p-1022" mix>f fconstant -min-+subn
s" -0.F FFFF FFFF FFFF p-1022" mix>f fconstant -max-+subn
s" -1.0 0000 0000 0000 p-1022" mix>f fconstant -min-+n
s" -1.F FFFF FFFF FFFF p+1023" mix>f fconstant -max-+n
s" 1.0 0000 0000 0001 p-1022" mix>f fconstant min-+n-up
s" -1.0 0000 0000 0001 p-1022" mix>f fconstant [-min-+n]-down
s" 1.0 0000 0000 0001 p0" mix>f fconstant 1-up
s" 1.F FFFF FFFF FFFF p-1" mix>f fconstant 1-down
s" -1.F FFFF FFFF FFFF p-1" mix>f fconstant [-1]-up
s" -1.0 0000 0000 0001 p0" mix>f fconstant [-1]-down
s" 1.0 0000 0000 0001 p1" mix>f fconstant 2-up
s" 1.F FFFF FFFF FFFF p0" mix>f fconstant 2-down
s" -1.F FFFF FFFF FFFF p0" mix>f fconstant [-2]-up
s" -1.0 0000 0000 0001 p1" mix>f fconstant [-2]-down
s" 1.F FFFF FFFF FFFE p+1023" mix>f fconstant max-+n-down
s" -1.F FFFF FFFF FFFE p+1023" mix>f fconstant [-max-+n]-up
cr
testing FNEXTUP
\ nan propagation
t{ +nan FNEXTUP -> +nan }t
t{ -nan FNEXTUP -> -nan }t
\ max-+n boundary
t{ +Inf FNEXTUP -> +INF }t
t{ max-+n FNEXTUP -> +INF }t
t{ max-+n-down FNEXTUP -> max-+n }t
\ +normal integer boundaries
t{ 2e FNEXTUP -> 2-up }t
t{ 2-down FNEXTUP -> 2e }t
t{ 1e FNEXTUP -> 1-up }t
t{ 1-down FNEXTUP -> 1e }t
\ +subnormal boundaries
t{ min-+subn FNEXTUP -> min-+subn f2* }t
t{ max-+subn FNEXTUP -> min-+n }t
t{ min-+n FNEXTUP -> min-+n-up }t
\ +/-0 boundaries
t{ +0 FNEXTUP -> min-+subn }t
t{ -0 FNEXTUP -> min-+subn }t
t{ -min-+subn FNEXTUP -> -0 }t
\ -subnormal boundaries
t{ -min-+subn f2* FNEXTUP -> -min-+subn }t
t{ [-min-+n]-down FNEXTUP -> -min-+n }t
t{ -min-+n FNEXTUP -> -max-+subn }t
\ -normal integer boundaries
t{ -1e FNEXTUP -> [-1]-up }t
t{ [-1]-down FNEXTUP -> -1e }t
t{ -2e FNEXTUP -> [-2]-up }t
t{ [-2]-down FNEXTUP -> -2e }t
\ -max-+n boundary
t{ -max-+n FNEXTUP -> [-max-+n]-up }t
t{ -Inf FNEXTUP -> -max-+n }t
verbose @ [IF]
Marcel Hendrix
> Andrew Haley wrote:
[..]
>> foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
>> foo df@ ;
[..]
> Wonderful! Now I have to teach myself why any carries and
> borrows into the exponent field are exactly what's needed.
> And that it'll still work for binary80.
> Unless my tests are incomplete. Either way I have something to
> understand.
I go for the incomplete test hypothesis, that "d+" doesn't look
good! I assume that the code will fail when testing some part of
the upper range.
BTW, the reason I didn't post the correction to FNEXTAFTER is because
also there I assume it tests OK but is still wrong.
-marcel
-- ----------------
CREATE -max-+n PRIVATE $FF C, $FF C, $FF C, $FF C, $FF C, $FF C, $EF C, $FF C,
-- Only for doubles
: FNEXTAFTER ( F: x y -- z )
DFLOCALS| y x |
'OF x 4 + 32B@ LOCAL hx \ HI(x)
'OF x 32B@ LOCAL lx \ LO(x)
'OF y 4 + 32B@ LOCAL hy \ HI(y)
'OF y 32B@ LOCAL ly \ LO(y)
hx $7FFFFFFF AND LOCAL ix \ |x|
hy $7FFFFFFF AND LOCAL iy \ |y|
ix $7FF00000 >= ix $7FF00000 - lx OR 0<> AND \ x is nan
iy $7FF00000 >= iy $7FF00000 - ly OR 0<> AND \ y is nan
OR IF x y F+ EXIT ENDIF ( both NaN )
x y F= IF y EXIT ENDIF
ix lx OR 0= IF hy $80000000 AND 'OF x 4 + 32B! \ HI(x)
1 'OF x 32B! \ LO(x)
x FSQR TO y \ return +-minsubnormal
y x F= IF y ELSE x ENDIF EXIT
ENDIF
hx 0>=
IF hx hy = lx ly > AND hx hy >
OR IF lx 0= IF -1 +TO hx ENDIF -1 +TO lx \ x > y, x -= ulp
ELSE 1 +TO lx lx 0= IF 1 +TO hx ENDIF \ x < y, x += ulp
ENDIF
ELSE hx hy = lx ly > AND hx hy > OR hy 0>=
OR IF lx 0= IF -1 +TO hx ENDIF -1 +TO lx \ x > y, x -= ulp
ELSE 1 +TO lx lx 0= IF 1 +TO hx ENDIF \ x < y, x += ulp
ENDIF
ENDIF
hx $7FF00000 AND TO hy
( >> ) hy $7FF00000 >= IF x -max-+n DF@ F= y -INF F= AND IF -INF EXIT ENDIF
x FSQR EXIT
ENDIF \ overflow
hy $00100000 < IF x FSQR TO y \ underflow
x y F<> IF hx 'OF y 4 + 32B!
lx 'OF y 32B! y EXIT
ENDIF
ENDIF
hx 'OF x 4 + 32B!
lx 'OF x 32B! x ;
: FNEXTUP ( F: x -- y )
FDUP FINFINITE? IF F0< IF -max-+n DF@ ELSE +INF ENDIF EXIT ENDIF
+INF FNEXTAFTER ;
David N. Williams
Marcel Hendrix wrote:
> "David N. Williams" <will...@umich.edu> writes Re: RfD: IEEE-FP 0.5.2
>
>> Andrew Haley wrote:
> [..]
>>> foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
>>> foo df@ ;
> [..]
>> Wonderful! Now I have to teach myself why any carries and
>> borrows into the exponent field are exactly what's needed.
>> And that it'll still work for binary80.
>
>> Unless my tests are incomplete. Either way I have something to
>> understand.
>
> I go for the incomplete test hypothesis, that "d+" doesn't look
> good! I assume that the code will fail when testing some part of
> the upper range.
Okay, I *think* I can prove it works when the leading integer
bit of the mantissa is implicit (unlike binary80), and integer
arithmetic is two's complement. It takes a bit of writing, so
I'll try that later.
Basically, you consider four cases where the number is nonzero,
namely, positive and negative numbers that are or are not powers
of two (zero or nonzero fractional part of mantissa). Then add
all binary one bits in the negative number case and just 1 in
the positive number case. It works out that the sign bit is
unchanged, and the exponent is modified by +/-1 exactly when it
should be.
This is always tricky stuff, so I do need to check it, but maybe
somebody knows a theorem about it?
As far as I can see, it's going to fail for positive numbers
when the leading integer bit is explicit, because the entire
mantissa gets turned into zero in the cases where there's a
carry out of the fraction, and the exponent gets bumped up by
one, giving a noncanonical representation of zero.
Or am I wrong?
> BTW, the reason I didn't post the correction to FNEXTAFTER is because
> also there I assume it tests OK but is still wrong.
Thanks for the copy. I'll have a look later.
-- David
David N. Williams
> [...]
> Okay, I *think* I can prove it works when the leading integer
> bit of the mantissa is implicit (unlike binary80), and integer
> arithmetic is two's complement. It takes a bit of writing, so
> I'll try that later.
I get it! Floating point data between +0 and +infinity
inclusive, considered as unsigned integers, is ordered in the
same way as the numbers they represent, with no gaps when the
integer bit is implicit. So adding 1 moves up one
representation step. The data with 1 in the sign bit,
considered as unsigned integers, has the same order as the
corresponding 0 sign-bit data, with no gaps, but the numbers
they represent have the opposite order, so *subtracting* one
moves up one representation step.
When the integer bit is explicit, there are gaps in the data,
requiring special treatment.
-- David
Andrew Haley
> Andrew Haley wrote:
> > David N. Williams <will...@umich.edu> wrote:
> >> Here are two more implementations of FNEXTUP, both of which use
> >> FLOGB and FSCALBN. I've also included Marcel's version, which
> >> uses FNEXTAFTER.
> >> [...]
> >
> > This all looks like fantastically hard work! If we can get at the
> > bits of a floating-point number it all becomes a lot easier. This is
> > an example on a little-endian 32-bit Forth:
> >
> > : nan? ( - t) ( f) fdup f= 0= ;
> > 2variable foo
> >
> > \ little-endian 2@ and 2!
> > : d! >r swap r> 2! ;
> > : d@ 2@ swap ;
> >
> > : fnextup ( f: x - y)
> > fdup nan? if exit then
> > 0.0e f+ \ Convert -0 to +0
> > foo df!
> > foo d@ dup 0< if 1 0 d- else 1 0 d+ then foo d!
> > foo df@ ;
> >
> > I may have missed a corner case, but you get the idea.
>
> Well, I was skeptical, so I put it to the test:
>
> --------------
Are you sure this is an error? I was assuming that the next float
after inf was NaN.
> Same result with gforth, and with big- or little-endian Mac OS X.
>
> Wonderful! Now I have to teach myself why any carries and borrows
> into the exponent field are exactly what's needed. And that it'll
> still work for binary80.
It probably won't work for binary80, because AFIAK binary80 does not
suppress the leading 1 bit in the significand. So, there has to be
special case code to handle that.
Andrew.
Andrew Haley
> > [...]
> > Okay, I *think* I can prove it works when the leading integer bit
> > of the mantissa is implicit (unlike binary80), and integer
> > arithmetic is two's complement. It takes a bit of writing, so
> > I'll try that later.
>
> I get it! Floating point data between +0 and +infinity inclusive,
> considered as unsigned integers, is ordered in the same way as the
> numbers they represent, with no gaps when the integer bit is
> implicit. So adding 1 moves up one representation step. The data
> with 1 in the sign bit, considered as unsigned integers, has the
> same order as the corresponding 0 sign-bit data, with no gaps, but
> the numbers they represent have the opposite order, so *subtracting*
> one moves up one representation step.
Exactly.
> When the integer bit is explicit, there are gaps in the data,
> requiring special treatment.
Yes; in essence, the leading 1 bit in binary80 is redundant, so
carries must propagate past it.
Andrew.
David N. Williams
>> Well, I was skeptical, so I put it to the test:
>>
>> #ERRORS: 1
>> --------------
>
> Are you sure this is an error? I was assuming that the next float
> after inf was NaN.
Indeed, the first datum after +inf is the first signalling nan.
But I just double checked. IEEE requires:
"...nextUp(+inf) is +inf,..."
The tests do cover every corner case that I could imagine would
be relevant. Of course one can never be sure what unnecessary
cases an implementation might include, like mine. An inherent
limitation of this sort of testing.
>> [...]
>
> It probably won't work for binary80, because AFIAK binary80 does not
> suppress the leading 1 bit in the significand. So, there has to be
> special case code to handle that.
Yep!
-- David
Andrew Haley wrote:
> David N. Williams <will...@umich.edu> wrote:
>> [...]
>> Well, I was skeptical, so I put it to the test:
>>
>> #ERRORS: 1
>> --------------
>
> Are you sure this is an error? I was assuming that the next float
> after inf was NaN.
Indeed, the first datum after +inf is the first signalling nan.
But I just double checked. IEEE requires:
"...nextUp(+inf) is +inf,..."
The tests do cover every corner case that I could imagine would
be relevant. Of course one can never be sure what unnecessary
cases an implementation might include, like mine. An inherent
limitation of this sort of testing.
Maybe one should encourage implementations that will be properly
tested by the tests. :-)
>> [...]
>
> It probably won't work for binary80, because AFIAK binary80 does not
> suppress the leading 1 bit in the significand. So, there has to be
> special case code to handle that.
Yep! I'll have to rethink the significance of this for the
tests. Some of the boundaries included seem relevant.
-- David
David N. Williams
Andrew Haley wrote:
> David N. Williams <will...@umich.edu> wrote:
>> [...]
>> Well, I was skeptical, so I put it to the test:
>>
>> #ERRORS: 1
>> --------------
>
> Are you sure this is an error? I was assuming that the next float
> after inf was NaN.
Indeed, the first datum after +inf is the first signalling nan.
But I just double checked. IEEE requires:
"...nextUp(+inf) is +inf,..."
The tests do cover every corner case that I could imagine would
be relevant. Of course one can never be sure what unnecessary
cases an implementation might include, like mine. An inherent
limitation of this sort of testing.
Maybe one should encourage implementations that will be properly
tested by the tests. :-)
>> [...]
>
> It probably won't work for binary80, because AFIAK binary80 does not
> suppress the leading 1 bit in the significand. So, there has to be
> special case code to handle that.
Yep! I'll have to rethink the significance of this for the
tests. Some of the boundaries included do seem relevant.
-- David
Andrew Haley
> Andrew Haley wrote:
> > David N. Williams <will...@umich.edu> wrote:
> >> [...]
> >> Well, I was skeptical, so I put it to the test:
> >>
> >> #ERRORS: 1
> >> --------------
> >
> > Are you sure this is an error? I was assuming that the next float
> > after inf was NaN.
>
> Indeed, the first datum after +inf is the first signalling nan.
> But I just double checked. IEEE requires:
>
> "...nextUp(+inf) is +inf,..."
Hmmm, I'm a bit surprised by that, but OK. I would have thought that
NaN would make more sense, given that the next real after infinity
*is* not a Number. But what the hell... :-)
> > It probably won't work for binary80, because AFIAK binary80 does not
> > suppress the leading 1 bit in the significand. So, there has to be
> > special case code to handle that.
>
> Yep!
Sure. The actual implementation should be trivial: wait for bit 63 to
overflow from 1 to 0 (and therefore propagate a carry to the exponent)
then set bit 63 to 1. Thereafter, bit 63 will always be 1. Err, I
think. :-)
> I'll have to rethink the significance of this for the tests. Some
> of the boundaries included do seem relevant.
Did you see http://www.math.utah.edu/~beebe/software/ieee/ ?
Andrew.
Thomas Pornin
> I get it! Floating point data between +0 and +infinity
> inclusive, considered as unsigned integers, is ordered in the
> same way as the numbers they represent, with no gaps when the
> integer bit is implicit.
And more generally, floating point data between -inf and +inf have the
same ordering than when their bit patterns are interpreted as integers,
with a sign bit + mantissa representation of integers. This is useful
for nextUp / nextDown, but also for comparisons. You still have to
add a specific test for NaN.
If you write an emulation of FP hardware with integer arithmetics, on a
CPU which uses two's complement for negative values, then you must also
take care of the sign and negative zeroes. There is such an emulator
written by D. Knuth himself, as part of his MMIX architecture:
http://www-cs-faculty.stanford.edu/~knuth/mmixware.html
It is a CWEB thing, i.e. both a book and a program. Released in public
domain, and quite instructive. Knuth uses a "classification function"
which returns 0, 1, 2 or 3, depending on whether the operand is a zero,
an infinite, a NaN, or something else. If the classification function is
called fcat(), then, for two operands a and b, the code usually begins
with the computation of fcat(a)+4*fcat(b), resulting in a value from 0
to 15. A simple switch statement (in a world of C) then handles the 15
special situations (0 to 14); afterwards, you know both operands are
normal, finite, non-zero numbers, and things are simpler.
--Thomas Pornin
David N. Williams
>> But I just double checked. IEEE requires:
>>
>> "...nextUp(+inf) is +inf,..."
>
> Hmmm, I'm a bit surprised by that, but OK. I would have thought that
> NaN would make more sense, given that the next real after infinity
> *is* not a Number. But what the hell... :-)
Depends whether you interpret the essence of nextup to be the
next representable affine real, or the next fp datum. :-)
>> I'll have to rethink the significance of this for the tests. Some
>> of the boundaries included do seem relevant.
>
> Did you see http://www.math.utah.edu/~beebe/software/ieee/ ?
Thanks. I didn't remember it, but did have it bookmarked. :-(
More useful stuff that'll eventually have to be waded through...
-- David
David Thompson
(Anton Ertl) wrote:
> Bernd Paysan <bernd....@gmx.de> writes:
> >http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf
>
> Is this a draft for the next C standard? Certainly the 2007 date
> makes it look like it's not a draft for C99.
>
It's the 'starting' draft for C1X; it consists of C99 + TC[123] and
doesn't have any of the proposed new stuff. (Later drafts do.)
This makes it a convenient (free, all in one file) reference, although
formally it's 'only' a draft and not 'the standard'.
n869 was the last draft of (what became) C99, if you want that.
And it's the last version produced in .txt (and .ps) as well as .pdf,
and there are some anti-PDF zealots on c.l.c who swear by that.
OTOH some standard text, especially on floating-point, uses math
notation that doesn't come through properly in .txt.
David N. Williams
Here's a portable implementation of Andrew's algorithm for
FNEXTUP based on rawfloat.fs, including tests:
http://www-personal.umich.edu/~williams/archive/forth/utilities/fnextup.fs
I plan to include FNEXTDOWN as well in the next draft of the
proposal, because it's a natural companion. Although it can be
defined simply in terms of FNEXTUP, I think the name needs to be
reserved.
I don't claim to be an expert on the use of these functions, but
here's an example, also requiring rawfloat.fs:
http://www-personal.umich.edu/~williams/archive/forth/utilities/fulps.fs
It defines two words, FULP and F-/ULP, which allow one to
compute relative error directly in ulps.
I'll post links for an application to the recently discussed erf
and erfc accuracy tests as a reply to that thread. Marcel has
some interesting further erf tests -- maybe he'll post them...
-- David
Marcel Hendrix
[..]
> http://www-personal.umich.edu/~williams/archive/forth/utilities/fulps.fs
> It defines two words, FULP and F-/ULP, which allow one to
> compute relative error directly in ulps.
> I'll post links for an application to the recently discussed erf
> and erfc accuracy tests as a reply to that thread. Marcel has
> some interesting further erf tests -- maybe he'll post them...
No big deal ...
CREATE erfdata
0.0e0 F, 0.0e0 F,
0.100000000000000000000000000000e-1 F, 0.112834155558496169159095235481e-1 F,
0.200000000000000000000000000000e-1 F, 0.225645746918449442243658616474e-1 F,
0.300000000000000000000000000000e-1 F, 0.338412223417354333022166542569e-1 F,
0.400000000000000000000000000000e-1 F, 0.451111061451247520897490647292e-1 F,
0.500000000000000000000000000000e-1 F, 0.563719777970166238312711126020e-1 F,
0.600000000000000000000000000000e-1 F, 0.676215943933084420794314523912e-1 F,
0.700000000000000000000000000000e-1 F, 0.788577197708907433569970386680e-1 F,
0.800000000000000000000000000000e-1 F, 0.900781258410181607233921876161e-1 F,
0.900000000000000000000000000000e-1 F, 0.101280593914626883352498163244e0 F,
0.100000000000000000000000000000e0 F, 0.112462916018284892203275071744e0 F,
0.110000000000000000000000000000e0 F, 0.123622896199474309769838562628e0 F,
0.120000000000000000000000000000e0 F, 0.134758351819920069295014126058e0 F,
0.130000000000000000000000000000e0 F, 0.145867114835695749082957608895e0 F,
0.140000000000000000000000000000e0 F, 0.156947033062855800871902246476e0 F,
0.150000000000000000000000000000e0 F, 0.167995971427363494770209522381e0 F,
0.160000000000000000000000000000e0 F, 0.179011813198105671155183728287e0 F,
0.170000000000000000000000000000e0 F, 0.189992461201808793395081985931e0 F,
0.180000000000000000000000000000e0 F, 0.200935839018695792648412710622e0 F,
0.190000000000000000000000000000e0 F, 0.211839892157749742936947248988e0 F,
0.200000000000000000000000000000e0 F, 0.222702589210478454140139006800e0 F,
0.210000000000000000000000000000e0 F, 0.233521922982103572925263000421e0 F,
0.220000000000000000000000000000e0 F, 0.244295911599128684869261694261e0 F,
0.230000000000000000000000000000e0 F, 0.255022599592273160770178527263e0 F,
0.240000000000000000000000000000e0 F, 0.265700058953792029795712090025e0 F,
0.250000000000000000000000000000e0 F, 0.276326390168236932985068267765e0 F,
0.260000000000000000000000000000e0 F, 0.286899723215749152048332055369e0 F,
0.270000000000000000000000000000e0 F, 0.297418218547012757893881942436e0 F,
0.280000000000000000000000000000e0 F, 0.307880068029034016648999599572e0 F,
0.290000000000000000000000000000e0 F, 0.318283495860952262338991948381e0 F,
0.300000000000000000000000000000e0 F, 0.328626759459127427638914047867e0 F,
0.310000000000000000000000000000e0 F, 0.338908150310790248699499364816e0 F,
0.320000000000000000000000000000e0 F, 0.349125994795582757315694633242e0 F,
0.330000000000000000000000000000e0 F, 0.359278654974358972988253953358e0 F,
0.340000000000000000000000000000e0 F, 0.369364529344658637204123996972e0 F,
0.350000000000000000000000000000e0 F, 0.379382053562310320297299983619e0 F,
0.360000000000000000000000000000e0 F, 0.389329701128664204753628450550e0 F,
0.370000000000000000000000000000e0 F, 0.399205984042999234582143710203e0 F,
0.380000000000000000000000000000e0 F, 0.409009453419694044918239128676e0 F,
0.390000000000000000000000000000e0 F, 0.418738700069796075759557782169e0 F,
0.400000000000000000000000000000e0 F, 0.428392355046668455103603845320e0 F,
0.410000000000000000000000000000e0 F, 0.437969090155439536373130097283e0 F,
0.420000000000000000000000000000e0 F, 0.447467618426025319819410593897e0 F,
0.430000000000000000000000000000e0 F, 0.456886694549540304987442003670e0 F,
0.440000000000000000000000000000e0 F, 0.466225115277957539319704643377e0 F,
0.450000000000000000000000000000e0 F, 0.475481719786923675319681525460e0 F,
0.460000000000000000000000000000e0 F, 0.484655390001679655026703122572e0 F,
0.470000000000000000000000000000e0 F, 0.493745050886082136517091596899e0 F,
0.480000000000000000000000000000e0 F, 0.502749670694764894533266228628e0 F,
0.490000000000000000000000000000e0 F, 0.511668261188523099211696708386e0 F,
0.500000000000000000000000000000e0 F, 0.520499877813046537682746653892e0 F,
0.510000000000000000000000000000e0 F, 0.529243619841170429010070312239e0 F,
0.520000000000000000000000000000e0 F, 0.537898630478854431106614413389e0 F,
0.530000000000000000000000000000e0 F, 0.546464096935141688222532672455e0 F,
0.540000000000000000000000000000e0 F, 0.554939250456390260497539494224e0 F,
0.550000000000000000000000000000e0 F, 0.563323366325108955991688363557e0 F,
0.560000000000000000000000000000e0 F, 0.571615763823768395668683848890e0 F,
0.570000000000000000000000000000e0 F, 0.579815806163996030236781957527e0 F,
0.580000000000000000000000000000e0 F, 0.587922900381600743987036943270e0 F,
0.590000000000000000000000000000e0 F, 0.595936497197908576518188514610e0 F,
0.600000000000000000000000000000e0 F, 0.603856090847925922562622436057e0 F,
0.610000000000000000000000000000e0 F, 0.611681218875880289504687484527e0 F,
0.620000000000000000000000000000e0 F, 0.619411461898721260562291908088e0 F,
0.630000000000000000000000000000e0 F, 0.627046443338195690464736260557e0 F,
0.640000000000000000000000000000e0 F, 0.634585829122141313859705694454e0 F,
0.650000000000000000000000000000e0 F, 0.642029327355671841292946466331e0 F,
0.660000000000000000000000000000e0 F, 0.649376687962954222750275011925e0 F,
0.670000000000000000000000000000e0 F, 0.656627702300305046439084463282e0 F,
0.680000000000000000000000000000e0 F, 0.663782202741357985424867698915e0 F,
0.690000000000000000000000000000e0 F, 0.670840062235077784356886812532e0 F,
0.700000000000000000000000000000e0 F, 0.677801193837418472975628809244e0 F,
0.710000000000000000000000000000e0 F, 0.684665550217444285287115347973e0 F,
0.720000000000000000000000000000e0 F, 0.691433123138751138842000477908e0 F,
0.730000000000000000000000000000e0 F, 0.698103942917044475821895794653e0 F,
0.740000000000000000000000000000e0 F, 0.704678077854745777673723006755e0 F,
0.750000000000000000000000000000e0 F, 0.711155633653515131598937834591e0 F,
0.760000000000000000000000000000e0 F, 0.717536752805590846718624840203e0 F,
0.770000000000000000000000000000e0 F, 0.723821613964859289254273344891e0 F,
0.780000000000000000000000000000e0 F, 0.730010431298578831244359445487e0 F,
0.790000000000000000000000000000e0 F, 0.736103453820691090374863741184e0 F,
0.800000000000000000000000000000e0 F, 0.742100964707660486167110586503e0 F,
0.810000000000000000000000000000e0 F, 0.748003280597789559230994456548e0 F,
0.820000000000000000000000000000e0 F, 0.753810750874962507155163885040e0 F,
0.830000000000000000000000000000e0 F, 0.759523756937772996814901189110e0 F,
0.840000000000000000000000000000e0 F, 0.765142711454994534663544476447e0 F,
0.850000000000000000000000000000e0 F, 0.770668057608352532380082742179e0 F,
0.860000000000000000000000000000e0 F, 0.776100268323556715663465156947e0 F,
0.870000000000000000000000000000e0 F, 0.781439845490550711651091968238e0 F,
0.880000000000000000000000000000e0 F, 0.786687319173932540042657464765e0 F,
0.890000000000000000000000000000e0 F, 0.791843246814495351087002484785e0 F,
0.900000000000000000000000000000e0 F, 0.796908212422832128518724785142e0 F,
0.910000000000000000000000000000e0 F, 0.801882825765941237429626001964e0 F,
0.920000000000000000000000000000e0 F, 0.806767721547761677691786836134e0 F,
0.930000000000000000000000000000e0 F, 0.811563558584557736234438237354e0 F,
0.940000000000000000000000000000e0 F, 0.816271018976062450999077401764e0 F,
0.950000000000000000000000000000e0 F, 0.820890807273277941907934455841e0 F,
0.960000000000000000000000000000e0 F, 0.825423649643818267103596022362e0 F,
0.970000000000000000000000000000e0 F, 0.829870293035667064651115188528e0 F,
0.980000000000000000000000000000e0 F, 0.834231504340207880514291707117e0 F,
0.990000000000000000000000000000e0 F, 0.838508069555369803579790230530e0 F,
1.0e0 F, 0.842700792949714869341220635083e0 F,
1.01000000000000000000000000000e0 F, 0.846810496228276697886304168608e0 F,
1.02000000000000000000000000000e0 F, 0.850838017700942042052443797678e0 F,
1.03000000000000000000000000000e0 F, 0.854784211454148381612389223130e0 F,
1.04000000000000000000000000000e0 F, 0.858649946526651453191924862667e0 F,
1.05000000000000000000000000000e0 F, 0.862436106090096697765225802369e0 F,
1.06000000000000000000000000000e0 F, 0.866143586635108082774536825816e0 F,
1.07000000000000000000000000000e0 F, 0.869773297163586659132689545761e0 F,
1.08000000000000000000000000000e0 F, 0.873326158387889589649984139486e0 F,
1.09000000000000000000000000000e0 F, 0.876803101937538279856830082636e0 F,
1.10000000000000000000000000000e0 F, 0.880205069574081699771867766322e0 F,
1.11000000000000000000000000000e0 F, 0.883533012414718050732088822843e0 F,
1.12000000000000000000000000000e0 F, 0.886787890165254649549837295827e0 F,
1.13000000000000000000000000000e0 F, 0.889970670362962317254717716140e0 F,
1.14000000000000000000000000000e0 F, 0.893082327629856715369184913907e0 F,
1.15000000000000000000000000000e0 F, 0.896123842936915012419506280749e0 F,
1.16000000000000000000000000000e0 F, 0.899096202879712030000041279757e0 F,
1.17000000000000000000000000000e0 F, 0.902000398965935653353894127360e0 F,
1.18000000000000000000000000000e0 F, 0.904837426915216837567544967696e0 F,
1.19000000000000000000000000000e0 F, 0.907608285971685037790484905145e0 F,
1.20000000000000000000000000000e0 F, 0.910313978229635380238405775715e0 F,
1.21000000000000000000000000000e0 F, 0.912955507972669409081355862097e0 F,
1.22000000000000000000000000000e0 F, 0.915533881026646830667192778601e0 F,
1.23000000000000000000000000000e0 F, 0.918050104126761367892733003921e0 F,
1.24000000000000000000000000000e0 F, 0.920505184299029669863949485264e0 F,
1.25000000000000000000000000000e0 F, 0.922900128256458230136523481197e0 F,
1.26000000000000000000000000000e0 F, 0.925235941810129484510735701213e0 F,
1.27000000000000000000000000000e0 F, 0.927513629295424719100127785559e0 F,
1.28000000000000000000000000000e0 F, 0.929734193013578152514172319937e0 F,
1.29000000000000000000000000000e0 F, 0.931898632688733592554242962899e0 F,
1.30000000000000000000000000000e0 F, 0.934007944940652436603893327504e0 F,
1.31000000000000000000000000000e0 F, 0.936063122773199513379509123921e0 F,
1.32000000000000000000000000000e0 F, 0.938065155078711378050908607034e0 F,
1.33000000000000000000000000000e0 F, 0.940015026158330197754594721551e0 F,
1.34000000000000000000000000000e0 F, 0.941913715258365323657306592595e0 F,
1.35000000000000000000000000000e0 F, 0.943762196122724061065830197376e0 F,
1.36000000000000000000000000000e0 F, 0.945561436561433041323207220960e0 F,
1.37000000000000000000000000000e0 F, 0.947312398035251987699776423971e0 F,
1.38000000000000000000000000000e0 F, 0.949016035256362570111835075789e0 F,
1.39000000000000000000000000000e0 F, 0.950673295805096476827758302413e0 F,
1.40000000000000000000000000000e0 F, 0.952285119762648810516482691533e0 F,
1.41000000000000000000000000000e0 F, 0.953852439359705454847226275613e0 F,
1.42000000000000000000000000000e0 F, 0.955376178640896168788291303906e0 F,
1.43000000000000000000000000000e0 F, 0.956857253144968859850496783450e0 F,
1.44000000000000000000000000000e0 F, 0.958296569600564774513511293169e0 F,
1.45000000000000000000000000000e0 F, 0.959695025637459232377442841965e0 F,
1.46000000000000000000000000000e0 F, 0.961053509513118027313720864175e0 F,
1.47000000000000000000000000000e0 F, 0.962372899854405729886810572718e0 F,
1.48000000000000000000000000000e0 F, 0.963654065414268855166142007980e0 F,
1.49000000000000000000000000000e0 F, 0.964897864843204212102907237564e0 F,
1.50000000000000000000000000000e0 F, 0.966105146475310727066976261646e0 F,
1.51000000000000000000000000000e0 F, 0.967276748128711635913825991249e0 F,
1.52000000000000000000000000000e0 F, 0.968413496920123165931303771472e0 F,
1.53000000000000000000000000000e0 F, 0.969516209093335679948252100952e0 F,
1.54000000000000000000000000000e0 F, 0.970585689861363727448926556465e0 F,
1.55000000000000000000000000000e0 F, 0.971622733262012538372566218848e0 F,
1.56000000000000000000000000000e0 F, 0.972628122026600200033798766967e0 F,
1.57000000000000000000000000000e0 F, 0.973602627461567070965915933850e0 F,
1.58000000000000000000000000000e0 F, 0.974547009342696901235480736712e0 F,
1.59000000000000000000000000000e0 F, 0.975462015821667639794133856309e0 F,
1.60000000000000000000000000000e0 F, 0.976348383344644007774283447142e0 F,
1.61000000000000000000000000000e0 F, 0.977206836582618593554778863365e0 F,
1.62000000000000000000000000000e0 F, 0.978038088373203471420329149739e0 F,
1.63000000000000000000000000000e0 F, 0.978842839673570150500482990249e0 F,
1.64000000000000000000000000000e0 F, 0.979621779524232013515386782363e0 F,
1.65000000000000000000000000000e0 F, 0.980375585023360294162452056584e0 F,
1.66000000000000000000000000000e0 F, 0.981104921311322055650338501742e0 F,
1.67000000000000000000000000000e0 F, 0.981810441565126558280197511885e0 F,
1.68000000000000000000000000000e0 F, 0.982492787002464827943095818625e0 F,
1.69000000000000000000000000000e0 F, 0.983152586895026146341339142860e0 F,
1.70000000000000000000000000000e0 F, 0.983790458590774563626242588122e0 F,
1.71000000000000000000000000000e0 F, 0.984407007544868370574911769514e0 F,
1.72000000000000000000000000000e0 F, 0.985002827358905745666562259787e0 F,
1.73000000000000000000000000000e0 F, 0.985578499828180497431188071400e0 F,
1.74000000000000000000000000000e0 F, 0.986134594996632938939331448608e0 F,
1.75000000000000000000000000000e0 F, 0.986671671219182443772211100129e0 F,
1.76000000000000000000000000000e0 F, 0.987190275231130125566791385785e0 F,
1.77000000000000000000000000000e0 F, 0.987690942224322340437040798990e0 F,
1.78000000000000000000000000000e0 F, 0.988174195929768317289393751016e0 F,
1.79000000000000000000000000000e0 F, 0.988640548706408159263283944287e0 F,
1.80000000000000000000000000000e0 F, 0.989090501635730714183732810756e0 F,
1.81000000000000000000000000000e0 F, 0.989524544621944366953582063504e0 F,
1.82000000000000000000000000000e0 F, 0.989943156497407646205094917059e0 F,
1.83000000000000000000000000000e0 F, 0.990346805133030645297761617229e0 F,
1.84000000000000000000000000000e0 F, 0.990735947553362618000490567051e0 F,
1.85000000000000000000000000000e0 F, 0.991111030056085706155209611037e0 F,
1.86000000000000000000000000000e0 F, 0.991472488335639574649460683136e0 F,
1.87000000000000000000000000000e0 F, 0.991820747610706752658708197512e0 F,
1.88000000000000000000000000000e0 F, 0.992156222755293694076721035182e0 F,
1.89000000000000000000000000000e0 F, 0.992479318433147959270641806759e0 F,
1.90000000000000000000000000000e0 F, 0.992790429235257469948357539303e0 F,
1.91000000000000000000000000000e0 F, 0.993089939820183484438217176392e0 F,
1.92000000000000000000000000000e0 F, 0.993378225056984767759469410979e0 F,
1.93000000000000000000000000000e0 F, 0.993655650170496375504198387176e0 F,
1.94000000000000000000000000000e0 F, 0.993922570888732519066543703907e0 F,
1.95000000000000000000000000000e0 F, 0.994179333592189118776835823032e0 F,
1.96000000000000000000000000000e0 F, 0.994426275464827868000969423733e0 F,
1.97000000000000000000000000000e0 F, 0.994663724646529912575485409079e0 F,
1.98000000000000000000000000000e0 F, 0.994892000386813583757222594664e0 F,
1.99000000000000000000000000000e0 F, 0.995111413199616997238347240473e0 F,
2.0e0 F, 0.995322265018952734162069256367e0 F,
2.01000000000000000000000000000e0 F, 0.995524849355248241312063095812e0 F,
2.02000000000000000000000000000e0 F, 0.995719451452192015978667962525e0 F,
2.03000000000000000000000000000e0 F, 0.995906348443912066070672206288e0 F,
2.04000000000000000000000000000e0 F, 0.996085809512319547895015186555e0 F,
2.05000000000000000000000000000e0 F, 0.996258096044456873132539403263e0 F,
2.06000000000000000000000000000e0 F, 0.996423461789695933775340133868e0 F,
2.07000000000000000000000000000e0 F, 0.996582153016638410456032750651e0 F,
2.08000000000000000000000000000e0 F, 0.996734408669576397404103730037e0 F,
2.09000000000000000000000000000e0 F, 0.996880460524377788338765083357e0 F,
2.10000000000000000000000000000e0 F, 0.997020533343667014496114983359e0 F,
2.11000000000000000000000000000e0 F, 0.997154845031177801648947935145e0 F,
2.12000000000000000000000000000e0 F, 0.997283606785160610778704424574e0 F,
2.13000000000000000000000000000e0 F, 0.997407023250733340776017430629e0 F,
2.14000000000000000000000000000e0 F, 0.997525292671069695356701128080e0 F,
2.15000000000000000000000000000e0 F, 0.997638607037325344858750510486e0 F,
2.16000000000000000000000000000e0 F, 0.997747152237207641699492752613e0 F,
2.17000000000000000000000000000e0 F, 0.997851108202100171372238604809e0 F,
2.18000000000000000000000000000e0 F, 0.997950649052658834678792689617e0 F,
2.19000000000000000000000000000e0 F, 0.998045943242801457529123598681e0 F,
2.20000000000000000000000000000e0 F, 0.998137153702018108556548243971e0 F,
2.21000000000000000000000000000e0 F, 0.998224437975934368814823456418e0 F,
2.22000000000000000000000000000e0 F, 0.998307948365064739107358899411e0 F,
2.23000000000000000000000000000e0 F, 0.998387832061698186549948381543e0 F,
2.24000000000000000000000000000e0 F, 0.998464231284862520615844968012e0 F,
2.25000000000000000000000000000e0 F, 0.998537283413318848302089203627e0 F,
2.26000000000000000000000000000e0 F, 0.998607121116541786642636179664e0 F,
2.27000000000000000000000000000e0 F, 0.998673872483645407328130841341e0 F,
2.28000000000000000000000000000e0 F, 0.998737661150219051712021176125e0 F,
2.29000000000000000000000000000e0 F, 0.998798606423041184302074966729e0 F,
2.30000000000000000000000000000e0 F, 0.998856823402643348534652540619e0 F,
2.31000000000000000000000000000e0 F, 0.998912423103700050040183186321e0 F,
2.32000000000000000000000000000e0 F, 0.998965512573224019809695702911e0 F,
2.33000000000000000000000000000e0 F, 0.999016195006549802974102506547e0 F,
2.34000000000000000000000000000e0 F, 0.999064569861091978842035580567e0 F,
2.35000000000000000000000000000e0 F, 0.999110732967867545150889457574e0 F,
2.36000000000000000000000000000e0 F, 0.999154776640775095111631875309e0 F,
2.37000000000000000000000000000e0 F, 0.999196789783626380902110236962e0 F,
2.38000000000000000000000000000e0 F, 0.999236857994928693095373816179e0 F,
2.39000000000000000000000000000e0 F, 0.999275063670419193575791325720e0 F,
2.40000000000000000000000000000e0 F, 0.999311486103354921430255067829e0 F,
2.41000000000000000000000000000e0 F, 0.999346201582564648884855889296e0 F,
2.42000000000000000000000000000e0 F, 0.999379283488271099505730786095e0 F,
2.43000000000000000000000000000e0 F, 0.999410802385694255639232279065e0 F,
2.44000000000000000000000000000e0 F, 0.999440826116448578590455613973e0 F,
2.45000000000000000000000000000e0 F, 0.999469419887748945596553592536e0 F,
2.46000000000000000000000000000e0 F, 0.999496646359441974605567857956e0 F,
2.47000000000000000000000000000e0 F, 0.999522565728881163674247158931e0 F,
2.48000000000000000000000000000e0 F, 0.999547235813665918980235727034e0 F,
2.49000000000000000000000000000e0 F, 0.999570712132266086606322840203e0 F,
2.50000000000000000000000000000e0 F, 0.999593047982555041060435784260e0 F,
2.51000000000000000000000000000e0 F, 0.999614294518275720661920525512e0 F,
2.52000000000000000000000000000e0 F, 0.999634500823465239215586949345e0 F,
2.53000000000000000000000000000e0 F, 0.999653713984864847611602181550e0 F,
2.54000000000000000000000000000e0 F, 0.999671979162343070964281585399e0 F,
2.55000000000000000000000000000e0 F, 0.999689339657360809492885128231e0 F,
2.56000000000000000000000000000e0 F, 0.999705836979508067426713308202e0 F,
2.57000000000000000000000000000e0 F, 0.999721510911142766670000652366e0 F,
2.58000000000000000000000000000e0 F, 0.999736399570162813678909351557e0 F,
2.59000000000000000000000000000e0 F, 0.999750539470943221871743838747e0 F,
2.60000000000000000000000000000e0 F, 0.999763965583470650796008996792e0 F,
2.61000000000000000000000000000e0 F, 0.999776711390708210081765550549e0 F,
2.62000000000000000000000000000e0 F, 0.999788808944223793772524050934e0 F,
2.63000000000000000000000000000e0 F, 0.999800288918115561773567671981e0 F,
2.64000000000000000000000000000e0 F, 0.999811180661268472697859895005e0 F,
2.65000000000000000000000000000e0 F, 0.999821512247975999096039245139e0 F,
2.66000000000000000000000000000e0 F, 0.999831310526961324669729427063e0 F,
2.67000000000000000000000000000e0 F, 0.999840601168832436289035175306e0 F,
2.68000000000000000000000000000e0 F, 0.999849408712005584127086981217e0 F,
2.69000000000000000000000000000e0 F, 0.999857756607131593604098286620e0 F,
2.70000000000000000000000000000e0 F, 0.999865667260059475670859881280e0 F,
2.71000000000000000000000000000e0 F, 0.999873162073371699777588724418e0 F,
2.72000000000000000000000000000e0 F, 0.999880261486525369137320906902e0 F,
2.73000000000000000000000000000e0 F, 0.999886985014633373018310896415e0 F,
2.74000000000000000000000000000e0 F, 0.999893351285919388145981320942e0 F,
2.75000000000000000000000000000e0 F, 0.999899378077880363163095608025e0 F,
2.76000000000000000000000000000e0 F, 0.999905082352189848729209758969e0 F,
2.77000000000000000000000000000e0 F, 0.999910480288375233419033074670e0 F,
2.78000000000000000000000000000e0 F, 0.999915587316301614224693273483e0 F,
2.79000000000000000000000000000e0 F, 0.999920418147494672237462812221e0 F,
2.80000000000000000000000000000e0 F, 0.999924986805334540975776754752e0 F,
2.81000000000000000000000000000e0 F, 0.999929306654152248770470340745e0 F,
2.82000000000000000000000000000e0 F, 0.999933390427259889483440616000e0 F,
2.83000000000000000000000000000e0 F, 0.999937250253945229426788983291e0 F,
2.84000000000000000000000000000e0 F, 0.999940897685460994406321233484e0 F,
2.85000000000000000000000000000e0 F, 0.999944343720038601012541387426e0 F,
2.86000000000000000000000000000e0 F, 0.999947598826955602236751136866e0 F,
2.87000000000000000000000000000e0 F, 0.999950672969685610748981162883e0 F,
2.88000000000000000000000000000e0 F, 0.999953575628158945224766772600e0 F,
2.89000000000000000000000000000e0 F, 0.999956315820161717373450749551e0 F,
2.90000000000000000000000000000e0 F, 0.999958902121900541164316132511e0 F,
2.91000000000000000000000000000e0 F, 0.999961342687759502470105416051e0 F,
2.92000000000000000000000000000e0 F, 0.999963645269275478192029159049e0 F,
2.93000000000000000000000000000e0 F, 0.999965817233357340078774458885e0 F,
2.94000000000000000000000000000e0 F, 0.999967865579774021028779347453e0 F,
2.95000000000000000000000000000e0 F, 0.999969796957935861737625829429e0 F,
2.96000000000000000000000000000e0 F, 0.999971617682993094132438878915e0 F,
2.97000000000000000000000000000e0 F, 0.999973333751274756079605293012e0 F,
2.98000000000000000000000000000e0 F, 0.999974950855090770264458608094e0 F,
2.99000000000000000000000000000e0 F, 0.999976474396919359773183773153e0 F,
3.0e0 F, 0.999977909503001414558627223870e0 F,
3.01000000000000000000000000000e0 F, 0.999979261036362867394040426859e0 F,
3.02000000000000000000000000000e0 F, 0.999980533609285585815045302128e0 F,
3.03000000000000000000000000000e0 F, 0.999981731595246738576624671730e0 F,
3.04000000000000000000000000000e0 F, 0.999982859140346051919801542932e0 F,
3.05000000000000000000000000000e0 F, 0.999983920174239833019164107976e0 F,
3.06000000000000000000000000000e0 F, 0.999984918420600105892450039369e0 F,
3.07000000000000000000000000000e0 F, 0.999985857407116679281055505040e0 F,
3.08000000000000000000000000000e0 F, 0.999986740475059447000170883915e0 F,
3.09000000000000000000000000000e0 F, 0.999987570788417709415870356510e0 F,
3.10000000000000000000000000000e0 F, 0.999988351342632800403966293837e0 F,
3.11000000000000000000000000000e0 F, 0.999989084972939807716761399466e0 F,
3.12000000000000000000000000000e0 F, 0.999989774362333686430314532983e0 F,
3.13000000000000000000000000000e0 F, 0.999990422049174585335565365967e0 F,
3.14000000000000000000000000000e0 F, 0.999991030434446735009890260323e0 F,
3.15000000000000000000000000000e0 F, 0.999991601788684784070190825728e0 F,
3.16000000000000000000000000000e0 F, 0.999992138258581016945168598578e0 F,
3.17000000000000000000000000000e0 F, 0.999992641873286442566995483560e0 F,
3.18000000000000000000000000000e0 F, 0.999993114550418308799710348013e0 F,
3.19000000000000000000000000000e0 F, 0.999993558101786172297789983350e0 F,
3.20000000000000000000000000000e0 F, 0.999993974238848237905028257637e0 F,
3.21000000000000000000000000000e0 F, 0.999994364577909275721057139761e0 F,
3.22000000000000000000000000000e0 F, 0.999994730645071027620087574688e0 F,
3.23000000000000000000000000000e0 F, 0.999995073880945628324028274273e0 F,
3.24000000000000000000000000000e0 F, 0.999995395645142189112254826297e0 F,
3.25000000000000000000000000000e0 F, 0.999995697220536324878169524105e0 F,
3.26000000000000000000000000000e0 F, 0.999995979817332047487635384467e0 F,
3.27000000000000000000000000000e0 F, 0.999996244576925100210857162634e0 F,
3.28000000000000000000000000000e0 F, 0.999996492575576469327953008413e0 F,
3.29000000000000000000000000000e0 F, 0.999996724827904479777092287035e0 F,
3.30000000000000000000000000000e0 F, 0.999996942290203561838538196597e0 F,
3.31000000000000000000000000000e0 F, 0.999997145863597465233101167206e0 F,
3.32000000000000000000000000000e0 F, 0.999997336397034395554134612132e0 F,
3.33000000000000000000000000000e0 F, 0.999997514690131255533767551989e0 F,
3.34000000000000000000000000000e0 F, 0.999997681495873890143583639601e0 F,
3.35000000000000000000000000000e0 F, 0.999997837523179959816750326118e0 F,
3.36000000000000000000000000000e0 F, 0.999997983439330800015058664371e0 F,
3.37000000000000000000000000000e0 F, 0.999998119872278367806604007120e0 F,
3.38000000000000000000000000000e0 F, 0.999998247412833126918516822639e0 F,
3.39000000000000000000000000000e0 F, 0.999998366616738481729929719677e0 F,
3.40000000000000000000000000000e0 F, 0.999998478006637137714638242862e0 F,
3.41000000000000000000000000000e0 F, 0.999998582073934540768368551870e0 F,
3.42000000000000000000000000000e0 F, 0.999998679280564330496741742403e0 F,
3.43000000000000000000000000000e0 F, 0.999998770060660532728832182400e0 F,
3.44000000000000000000000000000e0 F, 0.999998854822141014087435215358e0 F,
3.45000000000000000000000000000e0 F, 0.999998933948206526218905503710e0 F,
3.46000000000000000000000000000e0 F, 0.999999007798759479089598360193e0 F,
3.47000000000000000000000000000e0 F, 0.999999076711746401418634470684e0 F,
3.48000000000000000000000000000e0 F, 0.999999141004427871663589193713e0 F,
3.49000000000000000000000000000e0 F, 0.999999200974579534832405212782e0 F,
3.50000000000000000000000000000e0 F, 0.999999256901627658587254476316e0 F,
3.51000000000000000000000000000e0 F, 0.999999309047722526460750589923e0 F,
3.52000000000000000000000000000e0 F, 0.999999357658752816349254991422e0 F,
3.53000000000000000000000000000e0 F, 0.999999402965303968610621548877e0 F,
3.54000000000000000000000000000e0 F, 0.999999445183563409904602391610e0 F,
3.55000000000000000000000000000e0 F, 0.999999484516175366204968200593e0 F,
3.56000000000000000000000000000e0 F, 0.999999521153047871016724875387e0 F,
3.57000000000000000000000000000e0 F, 0.999999555272114452585246711416e0 F,
3.58000000000000000000000000000e0 F, 0.999999587040052866624547428960e0 F,
3.59000000000000000000000000000e0 F, 0.999999616612963128659530455453e0 F,
3.60000000000000000000000000000e0 F, 0.999999644137006992314701184444e0 F,
3.61000000000000000000000000000e0 F, 0.999999669749010916634964353452e0 F,
3.62000000000000000000000000000e0 F, 0.999999693577034466641036818617e0 F,
3.63000000000000000000000000000e0 F, 0.999999715740905996653837826725e0 F,
3.64000000000000000000000000000e0 F, 0.999999736352727375323111211250e0 F,
3.65000000000000000000000000000e0 F, 0.999999755517349424622675552567e0 F,
3.66000000000000000000000000000e0 F, 0.999999773332819662188394430113e0 F,
3.67000000000000000000000000000e0 F, 0.999999789890803857138683454022e0 F,
3.68000000000000000000000000000e0 F, 0.999999805276982833797808638200e0 F,
3.69000000000000000000000000000e0 F, 0.999999819571425885409309834693e0 F,
3.70000000000000000000000000000e0 F, 0.999999832848942090853797625924e0 F,
3.71000000000000000000000000000e0 F, 0.999999845179410761448597354666e0 F,
3.72000000000000000000000000000e0 F, 0.999999856628092181986012402275e0 F,
3.73000000000000000000000000000e0 F, 0.999999867255919750145399560419e0 F,
3.74000000000000000000000000000e0 F, 0.999999877119774561178120406869e0 F,
3.75000000000000000000000000000e0 F, 0.999999886272743430203346740923e0 F,
3.76000000000000000000000000000e0 F, 0.999999894764361292459492770171e0 F,
3.77000000000000000000000000000e0 F, 0.999999902640838872326778435481e0 F,
3.78000000000000000000000000000e0 F, 0.999999909945276464770342238222e0 F,
3.79000000000000000000000000000e0 F, 0.999999916717864627952816999842e0 F,
3.80000000000000000000000000000e0 F, 0.999999922996072543035871301802e0 F,
3.81000000000000000000000000000e0 F, 0.999999928814824756540487444583e0 F,
3.82000000000000000000000000000e0 F, 0.999999934206666981977303113488e0 F,
3.83000000000000000000000000000e0 F, 0.999999939201921600705773558932e0 F,
3.84000000000000000000000000000e0 F, 0.999999943828833467051722383598e0 F,
3.85000000000000000000000000000e0 F, 0.999999948113706589527419217230e0 F,
3.86000000000000000000000000000e0 F, 0.999999952081032228479841109489e0 F,
3.87000000000000000000000000000e0 F, 0.999999955753608920567184625380e0 F,
3.88000000000000000000000000000e0 F, 0.999999959152654912059634443601e0 F,
3.89000000000000000000000000000e0 F, 0.999999962297913456009131347825e0 F,
3.90000000000000000000000000000e0 F, 0.999999965207751402768257721692e0 F,
3.91000000000000000000000000000e0 F, 0.999999967899251489096718973094e0 F,
3.92000000000000000000000000000e0 F, 0.999999970388298708114035029833e0 F,
3.93000000000000000000000000000e0 F, 0.999999972689661120580136629049e0 F,
3.94000000000000000000000000000e0 F, 0.999999974817065447355070745065e0 F,
3.95000000000000000000000000000e0 F, 0.999999976783267763350692539111e0 F,
3.96000000000000000000000000000e0 F, 0.999999978600119594788977794192e0 F,
3.97000000000000000000000000000e0 F, 0.999999980278629704073471565005e0 F,
3.98000000000000000000000000000e0 F, 0.999999981829021830014495012629e0 F,
3.99000000000000000000000000000e0 F, 0.999999983260788635479156623334e0 F,
4.0e0 F, 0.999999984582742099719981147840e0 F,
4.01000000000000000000000000000e0 F, 0.999999985803060578628973811310e0 F,
4.02000000000000000000000000000e0 F, 0.999999986929332742926880683289e0 F,
4.03000000000000000000000000000e0 F, 0.999999987968598591791737966987e0 F,
4.04000000000000000000000000000e0 F, 0.999999988927387727619660545652e0 F,
4.05000000000000000000000000000e0 F, 0.999999989811755066458969764928e0 F,
4.06000000000000000000000000000e0 F, 0.999999990627314148132540129413e0 F,
4.07000000000000000000000000000e0 F, 0.999999991379268200130507554931e0 F,
4.08000000000000000000000000000e0 F, 0.999999992072439099985541182296e0 F,
4.09000000000000000000000000000e0 F, 0.999999992711294372006455502990e0 F,
4.10000000000000000000000000000e0 F, 0.999999993299972345915101627273e0 F,
4.11000000000000000000000000000e0 F, 0.999999993842305597079599254842e0 F,
4.12000000000000000000000000000e0 F, 0.999999994341842780638678614115e0 F,
4.13000000000000000000000000000e0 F, 0.999999994801868964843020600827e0 F,
4.14000000000000000000000000000e0 F, 0.999999995225424562376994050149e0 F,
4.15000000000000000000000000000e0 F, 0.999999995615322952246179748760e0 F,
4.16000000000000000000000000000e0 F, 0.999999995974166879001692277508e0 F,
4.17000000000000000000000000000e0 F, 0.999999996304363710601732706064e0 F,
4.18000000000000000000000000000e0 F, 0.999999996608139631065173452552e0 F,
4.19000000000000000000000000000e0 F, 0.999999996887552839233372260657e0 F,
4.20000000000000000000000000000e0 F, 0.999999997144505820407811384251e0 F,
4.21000000000000000000000000000e0 F, 0.999999997380756753356393392959e0 F,
4.22000000000000000000000000000e0 F, 0.999999997597930111164948394378e0 F,
4.23000000000000000000000000000e0 F, 0.999999997797526510638153870829e0 F,
4.24000000000000000000000000000e0 F, 0.999999997980931861411820942592e0 F,
4.25000000000000000000000000000e0 F, 0.999999998149425862613257479944e0 F,
4.26000000000000000000000000000e0 F, 0.999999998304189891785759031526e0 F,
4.27000000000000000000000000000e0 F, 0.999999998446314327865433847648e0 F,
4.28000000000000000000000000000e0 F, 0.999999998576805347252389970251e0 F,
4.29000000000000000000000000000e0 F, 0.999999998696591229443243860141e0 F,
4.30000000000000000000000000000e0 F, 0.999999998806528206277958696049e0 F,
4.31000000000000000000000000000e0 F, 0.999999998907405886591730810310e0 F,
4.32000000000000000000000000000e0 F, 0.999999998999952285943070625742e0 F,
4.33000000000000000000000000000e0 F, 0.999999999084838489103912519130e0 F,
4.34000000000000000000000000000e0 F, 0.999999999162682971138541822782e0 F,
4.35000000000000000000000000000e0 F, 0.999999999234055601157792481190e0 F,
4.36000000000000000000000000000e0 F, 0.999999999299481351206209929341e0 F,
4.37000000000000000000000000000e0 F, 0.999999999359443731215952324699e0 F,
4.38000000000000000000000000000e0 F, 0.999999999414387969535758671457e0 F,
4.39000000000000000000000000000e0 F, 0.999999999464723957210342368592e0 F,
4.40000000000000000000000000000e0 F, 0.999999999510828972939411158204e0 F,
4.41000000000000000000000000000e0 F, 0.999999999553050204480829840790e0 F,
4.42000000000000000000000000000e0 F, 0.999999999591707081174196866309e0 F,
4.43000000000000000000000000000e0 F, 0.999999999627093431244556363867e0 F,
4.44000000000000000000000000000e0 F, 0.999999999659479476596644416418e0 F,
4.45000000000000000000000000000e0 F, 0.999999999689113676923763682342e0 F,
4.46000000000000000000000000000e0 F, 0.999999999716224434128131400491e0 F,
4.47000000000000000000000000000e0 F, 0.999999999741021667277623092956e0 F,
4.48000000000000000000000000000e0 F, 0.999999999763698267603729136774e0 F,
4.49000000000000000000000000000e0 F, 0.999999999784431442373953599011e0 F,
4.50000000000000000000000000000e0 F, 0.999999999803383955845711252372e0 F,
4.51000000000000000000000000000e0 F, 0.999999999820705274925102741430e0 F,
4.52000000000000000000000000000e0 F, 0.999999999836532626610028664238e0 F,
4.53000000000000000000000000000e0 F, 0.999999999850991973790366257485e0 F,
4.54000000000000000000000000000e0 F, 0.999999999864198915505959672120e0 F,
4.55000000000000000000000000000e0 F, 0.999999999876259517323696760115e0 F,
4.56000000000000000000000000000e0 F, 0.999999999887271077085831737278e0 F,
4.57000000000000000000000000000e0 F, 0.999999999897322830900965564833e0 F,
4.58000000000000000000000000000e0 F, 0.999999999906496603894840018641e0 F,
4.59000000000000000000000000000e0 F, 0.999999999914867409908579731433e0 F,
4.60000000000000000000000000000e0 F, 0.999999999922504004025581681081e0 F,
4.61000000000000000000000000000e0 F, 0.999999999929469391523360046216e0 F,
4.62000000000000000000000000000e0 F, 0.999999999935821296581860068320e0 F,
4.63000000000000000000000000000e0 F, 0.999999999941612593833703413425e0 F,
4.64000000000000000000000000000e0 F, 0.999999999946891705613251788340e0 F,
4.65000000000000000000000000000e0 F, 0.999999999951702967549088757474e0 F,
4.66000000000000000000000000000e0 F, 0.999999999956086964947411723804e0 F,
4.67000000000000000000000000000e0 F, 0.999999999960080842230858548658e0 F,
4.68000000000000000000000000000e0 F, 0.999999999963718587527495376468e0 F,
4.69000000000000000000000000000e0 F, 0.999999999967031294347156301432e0 F,
4.70000000000000000000000000000e0 F, 0.999999999970047402136203396795e0 F,
4.71000000000000000000000000000e0 F, 0.999999999972792917366274934705e0 F,
4.72000000000000000000000000000e0 F, 0.999999999975291616686970292862e0 F,
4.73000000000000000000000000000e0 F, 0.999999999977565233555991084915e0 F,
4.74000000000000000000000000000e0 F, 0.999999999979633629652374485489e0 F,
4.75000000000000000000000000000e0 F, 0.999999999981514952278514689113e0 F,
4.76000000000000000000000000000e0 F, 0.999999999983225778864110488445e0 F,
4.77000000000000000000000000000e0 F, 0.999999999984781249599477465139e0 F,
4.78000000000000000000000000000e0 F, 0.999999999986195189146334073036e0 F,
4.79000000000000000000000000000e0 F, 0.999999999987480218300756940077e0 F,
4.80000000000000000000000000000e0 F, 0.999999999988647856415078039045e0 F,
4.81000000000000000000000000000e0 F, 0.999999999989708615322670025486e0 F,
4.82000000000000000000000000000e0 F, 0.999999999990672085451468209992e0 F,
4.83000000000000000000000000000e0 F, 0.999999999991547014758365897407e0 F,
4.84000000000000000000000000000e0 F, 0.999999999992341381066975482692e0 F,
4.85000000000000000000000000000e0 F, 0.999999999993062458345374197899e0 F,
4.86000000000000000000000000000e0 F, 0.999999999993716877418074910249e0 F,
4.87000000000000000000000000000e0 F, 0.999999999994310681567322356423e0 F,
4.88000000000000000000000000000e0 F, 0.999999999994849377442675180285e0 F,
4.89000000000000000000000000000e0 F, 0.999999999995337981664472470081e0 F,
4.90000000000000000000000000000e0 F, 0.999999999995781063475994218570e0 F,
4.91000000000000000000000000000e0 F, 0.999999999996182783770716950599e0 F,
4.92000000000000000000000000000e0 F, 0.999999999996546930794861023726e0 F,
4.93000000000000000000000000000e0 F, 0.999999999996876952801259878530e0 F,
4.94000000000000000000000000000e0 F, 0.999999999997175987908300720353e0 F,
4.95000000000000000000000000000e0 F, 0.999999999997446891397148719442e0 F,
4.96000000000000000000000000000e0 F, 0.999999999997692260661541069460e0 F,
4.97000000000000000000000000000e0 F, 0.999999999997914458007000963805e0 F,
4.98000000000000000000000000000e0 F, 0.999999999998115631480261460748e0 F,
4.99000000000000000000000000000e0 F, 0.999999999998297733894900324134e0 F,
5.0e0 F, 0.999999999998462540205571965150e0 F,
5.01000000000000000000000000000e0 F, 0.999999999998611663370691457092e0 F,
5.02000000000000000000000000000e0 F, 0.999999999998746568831894794398e0 F,
5.03000000000000000000000000000e0 F, 0.999999999998868587727991854116e0 F,
5.04000000000000000000000000000e0 F, 0.999999999998978928951372376787e0 F,
5.05000000000000000000000000000e0 F, 0.999999999999078690145854538527e0 F,
5.06000000000000000000000000000e0 F, 0.999999999999168867736719118814e0 F,
5.07000000000000000000000000000e0 F, 0.999999999999250366076093262174e0 F,
5.08000000000000000000000000000e0 F, 0.999999999999324005779884039855e0 F,
5.09000000000000000000000000000e0 F, 0.999999999999390531326065053131e0 F,
5.10000000000000000000000000000e0 F, 0.999999999999450617978244470040e0 F,
5.11000000000000000000000000000e0 F, 0.999999999999504878093048843929e0 F,
5.12000000000000000000000000000e0 F, 0.999999999999553866864905673402e0 F,
5.13000000000000000000000000000e0 F, 0.999999999999598087557263701516e0 F,
5.14000000000000000000000000000e0 F, 0.999999999999637996265120898687e0 F,
5.15000000000000000000000000000e0 F, 0.999999999999674006249905918680e0 F,
5.16000000000000000000000000000e0 F, 0.999999999999706491884251872353e0 F,
5.17000000000000000000000000000e0 F, 0.999999999999735792240985991292e0 F,
5.18000000000000000000000000000e0 F, 0.999999999999762214357711603317e0 F,
5.19000000000000000000000000000e0 F, 0.999999999999786036205658104257e0 F,
5.20000000000000000000000000000e0 F, 0.999999999999807509389000276403e0 F,
5.21000000000000000000000000000e0 F, 0.999999999999826861598581937643e0 F,
5.22000000000000000000000000000e0 F, 0.999999999999844298841903523956e0 F,
5.23000000000000000000000000000e0 F, 0.999999999999860007469333172102e0 F,
5.24000000000000000000000000000e0 F, 0.999999999999874156014761780234e0 F,
5.25000000000000000000000000000e0 F, 0.999999999999886896867331128461e0 F,
5.26000000000000000000000000000e0 F, 0.999999999999898367789408243565e0 F,
5.27000000000000000000000000000e0 F, 0.999999999999908693294647572811e0 F,
5.28000000000000000000000000000e0 F, 0.999999999999917985898764872331e0 F,
5.29000000000000000000000000000e0 F, 0.999999999999926347254533528257e0 F,
5.30000000000000000000000000000e0 F, 0.999999999999933869181496592017e0 F,
5.31000000000000000000000000000e0 F, 0.999999999999940634599958109874e0 F,
5.32000000000000000000000000000e0 F, 0.999999999999946718377967996920e0 F,
5.33000000000000000000000000000e0 F, 0.999999999999952188099238983211e0 F,
5.34000000000000000000000000000e0 F, 0.999999999999957104759225833219e0 F,
5.35000000000000000000000000000e0 F, 0.999999999999961523395950406804e0 F,
5.36000000000000000000000000000e0 F, 0.999999999999965493661565956818e0 F,
5.37000000000000000000000000000e0 F, 0.999999999999969060340115542622e0 F,
5.38000000000000000000000000000e0 F, 0.999999999999972263816448174868e0 F,
5.39000000000000000000000000000e0 F, 0.999999999999975140500808254347e0 F,
5.40000000000000000000000000000e0 F, 0.999999999999977723213205322052e0 F,
5.41000000000000000000000000000e0 F, 0.999999999999980041531298703032e0 F,
5.42000000000000000000000000000e0 F, 0.999999999999982122105192191147e0 F,
5.43000000000000000000000000000e0 F, 0.999999999999983988942224634412e0 F,
5.44000000000000000000000000000e0 F, 0.999999999999985663664560529668e0 F,
5.45000000000000000000000000000e0 F, 0.999999999999987165742128129508e0 F,
5.46000000000000000000000000000e0 F, 0.999999999999988512703218914430e0 F,
5.47000000000000000000000000000e0 F, 0.999999999999989720324849585256e0 F,
5.48000000000000000000000000000e0 F, 0.999999999999990802804794151544e0 F,
5.49000000000000000000000000000e0 F, 0.999999999999991772917017553673e0 F,
5.50000000000000000000000000000e0 F, 0.999999999999992642152082025602e0 F,
5.51000000000000000000000000000e0 F, 0.999999999999993420843951679358e0 F,
5.52000000000000000000000000000e0 F, 0.999999999999994118284488289052e0 F,
5.53000000000000000000000000000e0 F, 0.999999999999994742826810799941e0 F,
5.54000000000000000000000000000e0 F, 0.999999999999995301978581616314e0 F,
5.55000000000000000000000000000e0 F, 0.999999999999995802486183252850e0 F,
5.56000000000000000000000000000e0 F, 0.999999999999996250410658574798e0 F,
5.57000000000000000000000000000e0 F, 0.999999999999996651196205787507e0 F,
5.58000000000000000000000000000e0 F, 0.999999999999997009731944821516e0 F,
5.59000000000000000000000000000e0 F, 0.999999999999997330407604116545e0 F,
5.60000000000000000000000000000e0 F, 0.999999999999997617163715416982e0 F,
5.61000000000000000000000000000e0 F, 0.999999999999997873536848487795e0 F,
5.62000000000000000000000000000e0 F, 0.999999999999998102700367128162e0 F,
5.63000000000000000000000000000e0 F, 0.999999999999998307501142030855e0 F,
5.64000000000000000000000000000e0 F, 0.999999999999998490492614480613e0 F,
5.65000000000000000000000000000e0 F, 0.999999999999998653964567214376e0 F,
5.66000000000000000000000000000e0 F, 0.999999999999998799969924625156e0 F,
5.67000000000000000000000000000e0 F, 0.999999999999998930348873555947e0 F,
5.68000000000000000000000000000e0 F, 0.999999999999999046750567905987e0 F,
5.69000000000000000000000000000e0 F, 0.999999999999999150652654890572e0 F,
5.70000000000000000000000000000e0 F, 0.999999999999999243378837813750e0 F,
5.71000000000000000000000000000e0 F, 0.999999999999999326114669408407e0 F,
5.72000000000000000000000000000e0 F, 0.999999999999999399921750968712e0 F,
5.73000000000000000000000000000e0 F, 0.999999999999999465750495461908e0 F,
5.74000000000000000000000000000e0 F, 0.999999999999999524451597393236e0 F,
5.75000000000000000000000000000e0 F, 0.999999999999999576786338257426e0 F,
5.76000000000000000000000000000e0 F, 0.999999999999999623435843805072e0 F,
5.77000000000000000000000000000e0 F, 0.999999999999999665009397956809e0 F,
5.78000000000000000000000000000e0 F, 0.999999999999999702051907899012e0 F,
5.79000000000000000000000000000e0 F, 0.999999999999999735050605588257e0 F,
5.80000000000000000000000000000e0 F, 0.999999999999999764441062484356e0 F,
5.81000000000000000000000000000e0 F, 0.999999999999999790612586738276e0 F,
5.82000000000000000000000000000e0 F, 0.999999999999999813913065204374e0 F,
5.83000000000000000000000000000e0 F, 0.999999999999999834653306456241e0 F,
5.84000000000000000000000000000e0 F, 0.999999999999999853110935398287e0 F,
5.85000000000000000000000000000e0 F, 0.999999999999999869533885023630e0 F,
5.86000000000000000000000000000e0 F, 0.999999999999999884143526320489e0 F,
5.87000000000000000000000000000e0 F, 0.999999999999999897137473226909e0 F,
5.88000000000000000000000000000e0 F, 0.999999999999999908692095834299e0 F,
5.89000000000000000000000000000e0 F, 0.999999999999999918964771705183e0 F,
5.90000000000000000000000000000e0 F, 0.999999999999999928095902164495e0 F,
5.91000000000000000000000000000e0 F, 0.999999999999999936210717714908e0 F,
5.92000000000000000000000000000e0 F, 0.999999999999999943420894286201e0 F,
5.93000000000000000000000000000e0 F, 0.999999999999999949825999830484e0 F,
5.94000000000000000000000000000e0 F, 0.999999999999999955514788795568e0 F,
5.95000000000000000000000000000e0 F, 0.999999999999999960566360226603e0 F,
5.96000000000000000000000000000e0 F, 0.999999999999999965051193641935e0 F,
5.97000000000000000000000000000e0 F, 0.999999999999999969032075385547e0 F,
5.98000000000000000000000000000e0 F, 0.999999999999999972564926859643e0 F,
5.99000000000000000000000000000e0 F, 0.999999999999999975699544872679e0 F,
6.0e0 F, 0.999999999999999978480263287501e0 F,
6.01000000000000000000000000000e0 F, 0.999999999999999980946544209704e0 F,
6.02000000000000000000000000000e0 F, 0.999999999999999983133506107170e0 F,
6.03000000000000000000000000000e0 F, 0.999999999999999985072395488763e0 F,
6.04000000000000000000000000000e0 F, 0.999999999999999986791008084481e0 F,
6.05000000000000000000000000000e0 F, 0.999999999999999988314064853606e0 F,
6.06000000000000000000000000000e0 F, 0.999999999999999989663547594210e0 F,
6.07000000000000000000000000000e0 F, 0.999999999999999990858998430876e0 F,
6.08000000000000000000000000000e0 F, 0.999999999999999991917787011636e0 F,
6.09000000000000000000000000000e0 F, 0.999999999999999992855348845106e0 F,
6.10000000000000000000000000000e0 F, 0.999999999999999993685397849806e0 F,
6.11000000000000000000000000000e0 F, 0.999999999999999994420115865639e0 F,
6.12000000000000000000000000000e0 F, 0.999999999999999995070321588691e0 F,
6.13000000000000000000000000000e0 F, 0.999999999999999995645621131572e0 F,
6.14000000000000000000000000000e0 F, 0.999999999999999996154542179352e0 F,
6.15000000000000000000000000000e0 F, 0.999999999999999996604653503085e0 F,
6.16000000000000000000000000000e0 F, 0.999999999999999997002671406493e0 F,
6.17000000000000000000000000000e0 F, 0.999999999999999997354554514356e0 F,
6.18000000000000000000000000000e0 F, 0.999999999999999997665588161544e0 F,
6.19000000000000000000000000000e0 F, 0.999999999999999997940459507718e0 F,
6.20000000000000000000000000000e0 F, 0.999999999999999998183324382762e0 F,
6.21000000000000000000000000000e0 F, 0.999999999999999998397866760705e0 F,
6.22000000000000000000000000000e0 F, 0.999999999999999998587351663832e0 F,
6.23000000000000000000000000000e0 F, 0.999999999999999998754672212757e0 F,
6.24000000000000000000000000000e0 F, 0.999999999999999998902391461365e0 F,
6.25000000000000000000000000000e0 F, 0.999999999999999999032779586812e0 F,
6.26000000000000000000000000000e0 F, 0.999999999999999999147846943312e0 F,
6.27000000000000000000000000000e0 F, 0.999999999999999999249373433521e0 F,
6.28000000000000000000000000000e0 F, 0.999999999999999999338934602234e0 F,
6.29000000000000000000000000000e0 F, 0.999999999999999999417924813275e0 F,
6.30000000000000000000000000000e0 F, 0.999999999999999999487577831260e0 F,
6.31000000000000000000000000000e0 F, 0.999999999999999999548985094960e0 F,
6.32000000000000000000000000000e0 F, 0.999999999999999999603111937724e0 F,
6.33000000000000000000000000000e0 F, 0.999999999999999999650811982582e0 F,
6.34000000000000000000000000000e0 F, 0.999999999999999999692839914702e0 F,
6.35000000000000000000000000000e0 F, 0.999999999999999999729862811754e0 F,
6.36000000000000000000000000000e0 F, 0.999999999999999999762470192844e0 F,
6.37000000000000000000000000000e0 F, 0.999999999999999999791182929090e0 F,
6.38000000000000000000000000000e0 F, 0.999999999999999999816461143107e0 F,
6.39000000000000000000000000000e0 F, 0.999999999999999999838711210662e0 F,
6.40000000000000000000000000000e0 F, 0.999999999999999999858291965233e0 F,
6.41000000000000000000000000000e0 F, 0.999999999999999999875520195052e0 F,
6.42000000000000000000000000000e0 F, 0.999999999999999999890675512287e0 F,
6.43000000000000000000000000000e0 F, 0.999999999999999999904004665169e0 F,
6.44000000000000000000000000000e0 F, 0.999999999999999999915725355977e0 F,
6.45000000000000000000000000000e0 F, 0.999999999999999999926029620812e0 F,
6.46000000000000000000000000000e0 F, 0.999999999999999999935086820804e0 F,
6.47000000000000000000000000000e0 F, 0.999999999999999999943046288882e0 F,
6.48000000000000000000000000000e0 F, 0.999999999999999999950039671274e0 F,
6.49000000000000000000000000000e0 F, 0.999999999999999999956182998499e0 F,
6.50000000000000000000000000000e0 F, 0.999999999999999999961578516729e0 F,
6.51000000000000000000000000000e0 F, 0.999999999999999999966316306896e0 F,
6.52000000000000000000000000000e0 F, 0.999999999999999999970475715854e0 F,
6.53000000000000000000000000000e0 F, 0.999999999999999999974126621137e0 F,
6.54000000000000000000000000000e0 F, 0.999999999999999999977330548427e0 F,
6.55000000000000000000000000000e0 F, 0.999999999999999999980141658686e0 F,
6.56000000000000000000000000000e0 F, 0.999999999999999999982607619952e0 F,
6.57000000000000000000000000000e0 F, 0.999999999999999999984770377140e0 F,
6.58000000000000000000000000000e0 F, 0.999999999999999999986666831608e0 F,
6.59000000000000000000000000000e0 F, 0.999999999999999999988329440957e0 F,
6.60000000000000000000000000000e0 F, 0.999999999999999999989786748321e0 F,
6.61000000000000000000000000000e0 F, 0.999999999999999999991063849316e0 F,
6.62000000000000000000000000000e0 F, 0.999999999999999999992182803924e0 F,
6.63000000000000000000000000000e0 F, 0.999999999999999999993162999725e0 F,
6.64000000000000000000000000000e0 F, 0.999999999999999999994021472148e0 F,
6.65000000000000000000000000000e0 F, 0.999999999999999999994773186779e0 F,
6.66000000000000000000000000000e0 F, 0.999999999999999999995431288155e0 F,
6.67000000000000000000000000000e0 F, 0.999999999999999999996007318998e0 F,
6.68000000000000000000000000000e0 F, 0.999999999999999999996511413336e0 F,
6.69000000000000000000000000000e0 F, 0.999999999999999999996952466602e0 F,
6.70000000000000000000000000000e0 F, 0.999999999999999999997338285424e0 F,
6.71000000000000000000000000000e0 F, 0.999999999999999999997675719485e0 F,
6.72000000000000000000000000000e0 F, 0.999999999999999999997970777604e0 F,
6.73000000000000000000000000000e0 F, 0.999999999999999999998228729881e0 F,
6.74000000000000000000000000000e0 F, 0.999999999999999999998454197563e0 F,
6.75000000000000000000000000000e0 F, 0.999999999999999999998651232111e0 F,
6.76000000000000000000000000000e0 F, 0.999999999999999999998823384717e0 F,
6.77000000000000000000000000000e0 F, 0.999999999999999999998973767448e0 F,
6.78000000000000000000000000000e0 F, 0.999999999999999999999105106983e0 F,
6.79000000000000000000000000000e0 F, 0.999999999999999999999219791852e0 F,
6.80000000000000000000000000000e0 F, 0.999999999999999999999319913943e0 F,
6.81000000000000000000000000000e0 F, 0.999999999999999999999407304972e0 F,
6.82000000000000000000000000000e0 F, 0.999999999999999999999483568506e0 F,
6.83000000000000000000000000000e0 F, 0.999999999999999999999550108100e0 F,
6.84000000000000000000000000000e0 F, 0.999999999999999999999608151990e0 F,
6.85000000000000000000000000000e0 F, 0.999999999999999999999658774772e0 F,
6.86000000000000000000000000000e0 F, 0.999999999999999999999702916432e0 F,
6.87000000000000000000000000000e0 F, 0.999999999999999999999741399036e0 F,
6.88000000000000000000000000000e0 F, 0.999999999999999999999774941378e0 F,
6.89000000000000000000000000000e0 F, 0.999999999999999999999804171825e0 F,
6.90000000000000000000000000000e0 F, 0.999999999999999999999829639581e0 F,
6.91000000000000000000000000000e0 F, 0.999999999999999999999851824560e0 F,
6.92000000000000000000000000000e0 F, 0.999999999999999999999871146047e0 F,
6.93000000000000000000000000000e0 F, 0.999999999999999999999887970276e0 F,
6.94000000000000000000000000000e0 F, 0.999999999999999999999902617084e0 F,
6.95000000000000000000000000000e0 F, 0.999999999999999999999915365725e0 F,
6.96000000000000000000000000000e0 F, 0.999999999999999999999926459976e0 F,
6.97000000000000000000000000000e0 F, 0.999999999999999999999936112595e0 F,
6.98000000000000000000000000000e0 F, 0.999999999999999999999944509236e0 F,
6.99000000000000000000000000000e0 F, 0.999999999999999999999951811861e0 F,
7.0e0 F, 0.999999999999999999999958161744e0 F,
7.01000000000000000000000000000e0 F, 0.999999999999999999999963682079e0 F,
7.02000000000000000000000000000e0 F, 0.999999999999999999999968480280e0 F,
7.03000000000000000000000000000e0 F, 0.999999999999999999999972649977e0 F,
7.04000000000000000000000000000e0 F, 0.999999999999999999999976272771e0 F,
7.05000000000000000000000000000e0 F, 0.999999999999999999999979419766e0 F,
7.06000000000000000000000000000e0 F, 0.999999999999999999999982152904e0 F,
7.07000000000000000000000000000e0 F, 0.999999999999999999999984526136e0 F,
7.08000000000000000000000000000e0 F, 0.999999999999999999999986586443e0 F,
7.09000000000000000000000000000e0 F, 0.999999999999999999999988374729e0 F,
7.10000000000000000000000000000e0 F, 0.999999999999999999999989926597e0 F,
7.11000000000000000000000000000e0 F, 0.999999999999999999999991273035e0 F,
7.12000000000000000000000000000e0 F, 0.999999999999999999999992441002e0 F,
7.13000000000000000000000000000e0 F, 0.999999999999999999999993453952e0 F,
7.14000000000000000000000000000e0 F, 0.999999999999999999999994332284e0 F,
7.15000000000000000000000000000e0 F, 0.999999999999999999999995093736e0 F,
7.16000000000000000000000000000e0 F, 0.999999999999999999999995753729e0 F,
7.17000000000000000000000000000e0 F, 0.999999999999999999999996325667e0 F,
7.18000000000000000000000000000e0 F, 0.999999999999999999999996821200e0 F,
7.19000000000000000000000000000e0 F, 0.999999999999999999999997250449e0 F,
7.20000000000000000000000000000e0 F, 0.999999999999999999999997622205e0 F,
7.21000000000000000000000000000e0 F, 0.999999999999999999999997944106e0 F,
7.22000000000000000000000000000e0 F, 0.999999999999999999999998222780e0 F,
7.23000000000000000000000000000e0 F, 0.999999999999999999999998463985e0 F,
7.24000000000000000000000000000e0 F, 0.999999999999999999999998672716e0 F,
7.25000000000000000000000000000e0 F, 0.999999999999999999999998853310e0 F,
7.26000000000000000000000000000e0 F, 0.999999999999999999999999009528e0 F,
7.27000000000000000000000000000e0 F, 0.999999999999999999999999144634e0 F,
7.28000000000000000000000000000e0 F, 0.999999999999999999999999261456e0 F,
7.29000000000000000000000000000e0 F, 0.999999999999999999999999362450e0 F,
7.30000000000000000000000000000e0 F, 0.999999999999999999999999449743e0 F,
7.31000000000000000000000000000e0 F, 0.999999999999999999999999525177e0 F,
7.32000000000000000000000000000e0 F, 0.999999999999999999999999590352e0 F,
7.33000000000000000000000000000e0 F, 0.999999999999999999999999646650e0 F,
7.34000000000000000000000000000e0 F, 0.999999999999999999999999695272e0 F,
7.35000000000000000000000000000e0 F, 0.999999999999999999999999737256e0 F,
7.36000000000000000000000000000e0 F, 0.999999999999999999999999773500e0 F,
7.37000000000000000000000000000e0 F, 0.999999999999999999999999804783e0 F,
7.38000000000000000000000000000e0 F, 0.999999999999999999999999831779e0 F,
7.39000000000000000000000000000e0 F, 0.999999999999999999999999855070e0 F,
7.40000000000000000000000000000e0 F, 0.999999999999999999999999875161e0 F,
7.41000000000000000000000000000e0 F, 0.999999999999999999999999892489e0 F,
7.42000000000000000000000000000e0 F, 0.999999999999999999999999907430e0 F,
7.43000000000000000000000000000e0 F, 0.999999999999999999999999920310e0 F,
7.44000000000000000000000000000e0 F, 0.999999999999999999999999931412e0 F,
7.45000000000000000000000000000e0 F, 0.999999999999999999999999940978e0 F,
7.46000000000000000000000000000e0 F, 0.999999999999999999999999949221e0 F,
7.47000000000000000000000000000e0 F, 0.999999999999999999999999956321e0 F,
7.48000000000000000000000000000e0 F, 0.999999999999999999999999962436e0 F,
7.49000000000000000000000000000e0 F, 0.999999999999999999999999967701e0 F,
7.50000000000000000000000000000e0 F, 0.999999999999999999999999972234e0 F,
7.51000000000000000000000000000e0 F, 0.999999999999999999999999976135e0 F,
7.52000000000000000000000000000e0 F, 0.999999999999999999999999979492e0 F,
7.53000000000000000000000000000e0 F, 0.999999999999999999999999982381e0 F,
7.54000000000000000000000000000e0 F, 0.999999999999999999999999984865e0 F,
7.55000000000000000000000000000e0 F, 0.999999999999999999999999987002e0 F,
7.56000000000000000000000000000e0 F, 0.999999999999999999999999988839e0 F,
7.57000000000000000000000000000e0 F, 0.999999999999999999999999990419e0 F,
7.58000000000000000000000000000e0 F, 0.999999999999999999999999991776e0 F,
7.59000000000000000000000000000e0 F, 0.999999999999999999999999992943e0 F,
7.60000000000000000000000000000e0 F, 0.999999999999999999999999993945e0 F,
7.61000000000000000000000000000e0 F, 0.999999999999999999999999994806e0 F,
7.62000000000000000000000000000e0 F, 0.999999999999999999999999995546e0 F,
7.63000000000000000000000000000e0 F, 0.999999999999999999999999996181e0 F,
7.64000000000000000000000000000e0 F, 0.999999999999999999999999996726e0 F,
7.65000000000000000000000000000e0 F, 0.999999999999999999999999997194e0 F,
7.66000000000000000000000000000e0 F, 0.999999999999999999999999997595e0 F,
7.67000000000000000000000000000e0 F, 0.999999999999999999999999997940e0 F,
7.68000000000000000000000000000e0 F, 0.999999999999999999999999998235e0 F,
7.69000000000000000000000000000e0 F, 0.999999999999999999999999998488e0 F,
7.70000000000000000000000000000e0 F, 0.999999999999999999999999998706e0 F,
7.71000000000000000000000000000e0 F, 0.999999999999999999999999998892e0 F,
7.72000000000000000000000000000e0 F, 0.999999999999999999999999999052e0 F,
7.73000000000000000000000000000e0 F, 0.999999999999999999999999999188e0 F,
7.74000000000000000000000000000e0 F, 0.999999999999999999999999999306e0 F,
7.75000000000000000000000000000e0 F, 0.999999999999999999999999999406e0 F,
7.76000000000000000000000000000e0 F, 0.999999999999999999999999999492e0 F,
7.77000000000000000000000000000e0 F, 0.999999999999999999999999999566e0 F,
7.78000000000000000000000000000e0 F, 0.999999999999999999999999999629e0 F,
7.79000000000000000000000000000e0 F, 0.999999999999999999999999999683e0 F,
7.80000000000000000000000000000e0 F, 0.999999999999999999999999999729e0 F,
7.81000000000000000000000000000e0 F, 0.999999999999999999999999999768e0 F,
7.82000000000000000000000000000e0 F, 0.999999999999999999999999999802e0 F,
7.83000000000000000000000000000e0 F, 0.999999999999999999999999999831e0 F,
7.84000000000000000000000000000e0 F, 0.999999999999999999999999999856e0 F,
7.85000000000000000000000000000e0 F, 0.999999999999999999999999999877e0 F,
7.86000000000000000000000000000e0 F, 0.999999999999999999999999999895e0 F,
7.87000000000000000000000000000e0 F, 0.999999999999999999999999999910e0 F,
7.88000000000000000000000000000e0 F, 0.999999999999999999999999999923e0 F,
7.89000000000000000000000000000e0 F, 0.999999999999999999999999999935e0 F,
7.90000000000000000000000000000e0 F, 0.999999999999999999999999999944e0 F,
7.91000000000000000000000000000e0 F, 0.999999999999999999999999999952e0 F,
7.92000000000000000000000000000e0 F, 0.999999999999999999999999999959e0 F,
7.93000000000000000000000000000e0 F, 0.999999999999999999999999999965e0 F,
7.94000000000000000000000000000e0 F, 0.999999999999999999999999999971e0 F,
7.95000000000000000000000000000e0 F, 0.999999999999999999999999999975e0 F,
7.96000000000000000000000000000e0 F, 0.999999999999999999999999999979e0 F,
7.97000000000000000000000000000e0 F, 0.999999999999999999999999999982e0 F,
7.98000000000000000000000000000e0 F, 0.999999999999999999999999999985e0 F,
7.99000000000000000000000000000e0 F, 0.999999999999999999999999999987e0 F,
-- Tests upto 8 for full-range erf
( not using FLOCALS|, F.1, +E.R FLOAT[] )
0e FVALUE emin
0e FVALUE emax
: TESTS ( -- )
-Inf TO emax +Inf TO emin
CR ." ** TESTING erf **"
800 0 DO erfdata I 2* FLOATS + F@ erf
erfdata I 2* 1+ FLOATS + F@ F-/ULP
FDUP emin F< IF FDUP TO emin
CR x F. 2 SPACES FDUP F. ." ulps"
ENDIF
FDUP emax F> IF FDUP TO emax
CR x F. 2 SPACES FDUP F. ." ulps"
ENDIF
FDROP
LOOP
CR ." ** Summary: maximum rel. error = " emax F.
." ulp, minimum = " emin F. ." ulp **" ;
-- output ----------
erf1 -- CGM's code, abs. error < 2e-5
erf -- aph code vsn 1
erf_c -- BSD code port
erff -- BLAS code ok
FORTH> TESTS ( erf1 )
** TESTING erf1 **
0.0e+000 8.4e306 ulps
0.0e+000 8.4e306 ulps
1.0e-002 6.6e12 ulps
2.0e-002 5.1e12 ulps
3.0e-002 3.3e12 ulps
4.0e-002 3.3e12 ulps
5.0e-002 3.2e12 ulps
6.0e-002 2.8e12 ulps
7.0e-002 2.1e12 ulps
1.1e-001 1.6e12 ulps
1.2e-001 7.8e11 ulps
1.3e-001 7.1e11 ulps
1.5e-001 5.0e11 ulps
1.6e-001 4.8e11 ulps
1.7e-001 2.5e11 ulps
1.8e-001 2.3e11 ulps
1.9e-001 1.4e11 ulps
2.0e-001 -472302948.6 ulps
2.1e-001 -1.5e11 ulps
2.2e-001 -2.6e11 ulps
2.4e-001 -2.7e11 ulps
2.5e-001 -3.3e11 ulps
2.7e-001 -4.2e11 ulps
2.9e-001 -4.5e11 ulps
3.0e-001 -5.2e11 ulps
3.8e-001 -6.1e11 ulps
4.4e-001 -6.2e11 ulps
** Summary: maximum rel. error = 8.4e306 ulp, minimum = -6.2e11 ulp ** ok
FORTH> TESTS ( erf )
** TESTING erf **
0.0e+000 0.0 ulps
0.0e+000 0.0 ulps
1.0e-002 0.2 ulps
9.0e-002 0.3 ulps
2.1e-001 0.3 ulps
1.0e+000 0.4 ulps
1.2e+000 0.4 ulps
1.2e+000 0.4 ulps
1.4e+000 0.4 ulps
1.5e+000 0.5 ulps
1.6e+000 0.5 ulps
1.7e+000 0.5 ulps
1.9e+000 0.5 ulps
2.0e+000 0.6 ulps
2.1e+000 0.7 ulps
2.2e+000 0.8 ulps
2.6e+000 0.9 ulps
4.0e+000 0.0 ulps
4.1e+000 0.0 ulps
4.1e+000 0.0 ulps
4.1e+000 0.0 ulps
4.2e+000 0.0 ulps
4.3e+000 0.0 ulps
4.4e+000 0.0 ulps
5.2e+000 0.0 ulps
5.3e+000 0.0 ulps
5.3e+000 0.0 ulps
5.4e+000 0.0 ulps
5.4e+000 0.0 ulps
5.5e+000 0.0 ulps
** Summary: maximum rel. error = 0.9 ulp, minimum = 0.0 ulp ** ok
FORTH> TESTS ( erf_c )
** TESTING erf_c **
0.0e+000 0.0 ulps
0.0e+000 0.0 ulps
1.0e-002 0.1 ulps
3.0e-002 -0.4 ulps
5.0e-002 0.4 ulps
6.0e-002 -0.4 ulps
7.0e-002 0.6 ulps
9.0e-002 -0.6 ulps
1.4e-001 0.8 ulps
2.1e-001 -0.7 ulps
4.7e-001 -0.8 ulps
** Summary: maximum rel. error = 0.8 ulp, minimum = -0.8 ulp ** ok
FORTH> TESTS ( erff )
** TESTING erff **
0.0e+000 0.0 ulps
0.0e+000 0.0 ulps
1.0e-002 -538.8 ulps
9.0e-002 -633.0 ulps
2.1e-001 -691.1 ulps
1.2e+000 91.9 ulps
1.2e+000 131.0 ulps
1.2e+000 146.1 ulps
** Summary: maximum rel. error = 146.1 ulp, minimum = -691.1 ulp ** ok
David N. Williams
> Ed wrote:
> [...]
>> The IEEE strategy of allowing one to accept "overflow" input
>> and then substitute it with various "good results" is disconcerting
>> as it weakens the ability of programs to discriminate. I assume
>> the same 'substitution strategy' applies to your IEEE >FLOAT
>> string input function (?)
>> ...
>
> It may not be as bad as that.
>
> Re-reading the section on IEEE overflows etc again, I note that
> as well as providing default results, one is supposed to "signal"
> the exception. How will an IEEE Forth signal exceptions?
>
> In '94, I/O exceptions are signalled with an "ior" result from the
> relevant function. The text interpreter doesn't have that - it
> simply works on a pass/fail result. So if I enter the following
> in immediate mode:
>
> 1e9000
>
> should it signal the overflow exception - or quietly return an INF
> (or whatever) ? Similarly should calculations entered in immediate
> mode which result in divide-by-zero signal the exception?
>
> IMO before one can define IEEE I/O functions, one must first
> know how IEEE exceptions are to be handled. If one is going to
> use ior values, then that will inform the stack effect of IEEE I/O
> functions. If one plans to use e.g. >FLOAT for IEEE I/O then a
> different method of handling exceptions will be needed.
I certainly agree that how IEEE exceptions are to be handled has
an impact on the rest of the proposal, even if they're to be
deferred in the first go-around.
My assumption has always been that no Forth-like ior's would be
involved, that the kind of signalling that occurs is to set an
fp status flag, produce the default IEEE datum, and keep going,
leaving it to the user to query the status flags or not, as he
chooses. That seems to be the underlying IEEE 754 philosophy,
due to Kahan.
With that approach, there is no impact on stack effects in the
rest of the proposal.
Maybe it would further the discussion of fp exceptions to see an
example of their use in |Z|, copied below from:
http://www-personal.umich.edu/~williams/archive/forth/complex/complex-kahan.fs
The two exception words it uses are copied below from my ancient
fpieee-ext.c pfe module. The actions mentioned in the following
list correspond to FPIEEE [IF] ... [THEN] blocks in the |Z|
code, which is a slavish imitation of a Kahan algorithm.
* At entry, the system fp exception flags are whatever has been
generated from previous operations.
* The FE-INVALID flag is saved at the beginning, and restored at
the end, as part of the NaN-propagation process.
* At the point where irrelevant underflow exceptions might begin
to occur, the FE-UNDERFLOW flag is saved, so if there is no
relevant underflow, the previous underflow state can be
preserved.
* Just before the point where a meaningful underflow exception
might be generated, the original flag is restored. It gets
overwritten if underflow occurs, otherwise it's unchanged.
* As written, the code also works when the fp exception words
are not present. The only effect is that the final exception
state is uncertain.
The fp exception state can then be queried and reacted to by
words that use |Z|.
-- David
------------
/* In the following stack effects, "fexcept" indicates an
* implementation-defined floating-point exception bit mask in a
* C99 fexcept_t type, while "excepts" indicates exception flag
* data corresponding to a combination of the masks in a C int
* type. An ambiguous condition exists if either fexcept_t or
* int does not fit into a cell. The four C99 functions
* fegetexceptflag(), fesetexceptflag(), feraiseexcept(), and
* feclearexcept() return zero if successful and a nonzero value
* if not, which we ignore, because some older systems (BSD?)
* return void instead.
*
* Forth constants FE-DIVBYZERO, FE-INEXACT, FE-INVALID,
* FE-OVERFLOW, FE-UNDERFLOW, FE-ALL-EXCEPT, corresponding to
* the C99 exception masks named with "_" instead of "-". are
* defined under P4_LISTWORDS.
*/
/** FEGETEXCEPTFLAG ( excepts -- fenv )
*/
FCode (p4_fegetexceptflag)
{ fegetexceptflag ( (fexcept_t*) SP, *SP ); }
/** FESETEXCEPTFLAG ( fenv excepts -- )
*/
FCode (p4_fesetexceptflag)
{
fesetexceptflag ( (fexcept_t*) &SP[1], SP[0] );
SP += 2;
}
-------------
-----------------------
: |z| (f: x y -- |z|)
(
This is Kahan's algorithm, which treats infinity, the precision
of 1+sqrt[2], and the roundoff threshhold for large ratios of
the maximum and minimum absolute values of x and y. When the
FPIEEE word list is present, it also treats the IEEE underflow
and invalid flags.
)
fpframe>r
FPIEEE [IF]
FE-INVALID fegetexceptflag ( invalid?)
[THEN]
fabs fswap fabs
zdup f<
IF fswap
THEN \ a >= b >= 0 if not NaN (f: a b )
fdup +Inf 0E f~
IF fnip +Inf (f: a=inf b=inf)
THEN
zdup to-fb to-fa
f- fdup fa (f: t=a-b t a)
fdup +Inf 0E f~ 0E f~ or (f: t)
( t=a.or.a=inf?)
IF \ b is negligible
(f: t) fdrop 0E (f: s=0)
ELSE \ b is not negligible
FPIEEE [IF]
FE-UNDERFLOW fegetexceptflag ( invalid? underflow?)
[THEN]
fa fover (f: t b t)
f< IF \ 2 < a/b
(f: t) fdrop
fa fb f/ (f: u=a/b)
fdup 2/epsilon f<
IF \ a/b < 2/epsilon
fdup f^2 1E f+ fsqrt f+ (f: s=u+sqrt[1+u^2])
ELSE
fdrop 0E (f: s=0)
THEN
ELSE \ 1 <= a/b <= 2
(f: t)
fb f/ (f: v=t/b)
fdup fdup 2E f+ f* (f: v w=[2+v]v)
fdup 2E f+ fsqrt sqrt2 f+ f/ f+ (f: h=w/[sqrt[2+w]+sqrt2]+v)
t2p1 f+ r2p1 f+ (f: s=[t2p1+h]+r2p1)
\ substitute v = u - 1 and w = u^2 - 1 to prove that
\ this s expression is the same as the one above
THEN
(f: s)
fb fswap f/ (f: s=b/s)
FPIEEE [IF] ( invalid? underflow?)
FE-UNDERFLOW fesetexceptflag ( invalid?)
\ any underflow is now deserved
[THEN]
THEN
FPIEEE [IF]
( invalid?) FE-INVALID fesetexceptflag
[THEN]
(f: s) fa f+ (f: s+a)
r>frame ;
--------------------
Andrew Haley
> I certainly agree that how IEEE exceptions are to be handled has an
> impact on the rest of the proposal, even if they're to be deferred
> in the first go-around.
> My assumption has always been that no Forth-like ior's would be
> involved, that the kind of signalling that occurs is to set an fp
> status flag, produce the default IEEE datum, and keep going, leaving
> it to the user to query the status flags or not, as he chooses.
> That seems to be the underlying IEEE 754 philosophy, due to Kahan.
Makes sense. I'm not aware of any great barrier to adopting the
simplest form of floating-point exception handling right now.
> With that approach, there is no impact on stack effects in the
> rest of the proposal.
> Maybe it would further the discussion of fp exceptions to see an
> example of their use in |Z|, copied below from:
> http://www-personal.umich.edu/~williams/archive/forth/complex/complex-kahan.fs
Thanks, that's very useful.
> /** FEGETEXCEPTFLAG ( excepts -- fenv )
> */
> FCode (p4_fegetexceptflag)
> { fegetexceptflag ( (fexcept_t*) SP, *SP ); }
> /** FESETEXCEPTFLAG ( fenv excepts -- )
> */
> FCode (p4_fesetexceptflag)
> {
> fesetexceptflag ( (fexcept_t*) &SP[1], SP[0] );
> SP += 2;
> }
> -------------
With all due respect. :-) Eeew, these names are horrible. They are,
however, the same as those defined in ISO 9899:1999 (C99) Section
7.6.2 so I suppose it makes sense to keep them.
Andrew.
David N. Williams
> David N. Williams <will...@umich.edu> wrote:
>
>> I certainly agree that how IEEE exceptions are to be handled has an
>> impact on the rest of the proposal, even if they're to be deferred
>> in the first go-around.
>
>> My assumption has always been that no Forth-like ior's would be
>> involved, that the kind of signalling that occurs is to set an fp
>> status flag, produce the default IEEE datum, and keep going, leaving
>> it to the user to query the status flags or not, as he chooses.
>> That seems to be the underlying IEEE 754 philosophy, due to Kahan.
>
> Makes sense. I'm not aware of any great barrier to adopting the
> simplest form of floating-point exception handling right now.
There were some objections to the proposal in the first draft. I'll
have to review what they were.
>> With that approach, there is no impact on stack effects in the
>> rest of the proposal.
>> [...]
>> /** FEGETEXCEPTFLAG ( excepts -- fenv )
>> */
>> FCode (p4_fegetexceptflag)
>> { fegetexceptflag ( (fexcept_t*) SP, *SP ); }
>
>> /** FESETEXCEPTFLAG ( fenv excepts -- )
>> */
>> FCode (p4_fesetexceptflag)
>> {
>> fesetexceptflag ( (fexcept_t*) &SP[1], SP[0] );
>> SP += 2;
>> }
>> -------------
>
> With all due respect. :-) Eeew, these names are horrible. They are,
> however, the same as those defined in ISO 9899:1999 (C99) Section
> 7.6.2 so I suppose it makes sense to keep them.
I'm wasn't actually proposing those names, although they were in
the first draft. They're from an implementation nobody used but
myself, chosen, as you mention, to agree with C99.
The naming problem goes to what to call IEEE fp exceptions:
Anton Ertl wrote:
> Andrew Haley <andr...@littlepinkcloud.invalid> writes:
>> I think that you should not use the word "exception" in a Forth
>> standard to refer to anything other than a Forth exception. Maybe
>> call them "FP exceptions"?
>
> I think that's still too easily confused. The first thing I think of
> when I read this is what happens when something throws one of the
> varlues -55, -54, -46..-41.
>
> Either make an acronym of it (e.g., IFENFE for "IEEE FP exception, not
> Forth exception"), or use a different word, e.g, condition.
Maybe what we should take from C99 is that the terminology for
fp exceptions should have "flag" in it, maybe something like
"fe flag".
Or how about just FEF? Then possible names might be:
GET-FEF SET-FEF TEST-FEF RAISE-FEF CLEAR-FEF
Other suggestions?
-- David