[comp.lang.c] floating point

firth@sei.cmu.edu.UUCP (04/08/87)

This is not directly relevant to the C language, but
since it continues a topic posted here, I'm sending
it to this group.  It concerns the two issues of
floating-point equality and the reordering of
floating-point code.  It is rather long; I hope you
find it useful.

The background is that I had to decide about these
issues when implementing a set of high-level languages,
of which the most significant in this context were
Algol-60 and Coral-66, both of which are used for
floating calculations.  The language definitions
gave some latitude to the implementor; as implementor
therefore I had to decide what to do with it.

The first step was to look at some user code.  Our
users were mostly mathematicians and engineers, not
computer types.  They used floating point for real
calculations - structures, radar antennae, electron
orbits &c.  Many of these computations were very
complicated, with two and three line expressions
not uncommon. (Good test of your register allocation
strategy, too!)

The Algol users were normally very careful about
tests for convergence and closeness; some of the
Coral code was sloppy.  So, for floating equality,
I had three choices

	(a) fuzzy equality
	(b) bitwise equality
	(c) no equality

Choice (a) seemed wrong.  There was no good definition
of fuzziness, and it would give anomalies such as

	x=y & y=z & x~=z

Choice (b) also seemed wrong - most of the explicit tests
for equality I found were dubious, to say the least.  In
the end, I decided that the compiler would reject a test
for floating equality or nonequality unless one operand
were an explicit zero.

As a second line of defence, the user notes explained this
prohibition, with the comment

	The test 'IF x = y THEN ...' may be rewritten
	'IF x-y = 0 THEN ...', which should be read
	'if the subtraction of y from x causes underflow...'

So the programmer tripped up by the compiler would be alerted
that something subtle was involved here.  In fact, most of
the Algol users were already in the habit of using zero as
the preferred asymptote for a converging series.

The Algol-60 FOR loop is defined to use an ordered comparison
for the termination test, so I decided to leave alone things
like
	FOR x := 0 STEP 0.1 UNTIL 1 DO ...

(which anyway didn't happen often).

Floating point reordering was another hard issue.  I asked the
users about their programs, and the consensus was that they
regarded floating ("real") values not as bit patterns but as
approximations to real numbers, and floating computations as
approximations to mathematical formulae.  They also understood
rounding problems, so for instance would compute known sums by
adding successively larger terms, sum separately positive and
negative terms, and so on.

Moreover, many of their expressions were very complicated
transcriptions of mathematical formulae.  Such code can be
maintained only if the translation is as simple and explicit
as possible.  Forcing users to break up expressions, to insert
parentheses or other symbols, would seriously affect the
maintainability of their code.  So, in the end, I decided
to do NO general reordering, but to honour the conventional
evaluation order, including left associativity of operators.

However, there were some transformations I did feel free to
make.  Given for instance

	(E1) + (E2)

I could evaluate either E1 or E2 first (in the absence of
function calls), so as to use fewer registers.  I could
also make transformations that were exact on the given
target, eg

	(-x)*y => -(x*y)
	(-x)+y => y-x

And, given that I found no user who wanted worse precision, I
could do things like

	x*(y^-1) => x/y

where the result is different but "better".  All this took a
couple of pages of code, and a dozen or so table entries, to
separate actions for floating expressions from actions for
integer expressions.  Not too tedious, and I never had to lie
awake at night worrying that my compiler had totally trashed
the design of the latest phased-array radar antenna.

Thanks for reading all this.  It is obviously only one person's
opinion, but it may help someone.

Robert Firth

henry@utzoo.uucp (Henry Spencer) (03/07/90)

In article <1990Mar5.174628.9141@ncsuvx.ncsu.edu> harish@ecebucolix.ncsu.edu (Harish P. Hiriyannaiah) writes:
>As an aside, what do netters feel about DEC's 64-bit format for floats vs the
>IEEE one ? I contend that the DEC format is superior...

Which DEC 64-bit format?  There are at least two... which brings us back
to how nice it is to have *one* standard, even if DEC doesn't like it.
Even DEC now sells IEEE-format floating-point hardware, in their Mips-based
machines.
-- 
MSDOS, abbrev:  Maybe SomeDay |     Henry Spencer at U of Toronto Zoology
an Operating System.          | uunet!attcan!utzoo!henry henry@zoo.toronto.edu