From: torek@elf.bsdi.com (Chris Torek)
Newsgroups: comp.lang.c
Subject: Re: floor() implementation
Date: 2 Jan 1999 14:59:10 -0800

>"Amund Frislie" <sfris@online.no> writes:
>   Does anyone know where I can find [an] implementation [of] floor() ...

In article <87yanlr303.fsf@pfaffben.user.msu.edu> Ben Pfaff
<pfaffben@pilot.msu.edu> replies:
>floor() can't be implemented in C AFAIK without use of other library
>functions.

Indeed.  Typical implementation, usually written in assembly language,
rely on the floating-point format inside the machine.

Suppose "d" is a double-precision floating point number represented
internally as <sign> <fraction> <exponent>, where <sign> is 0 for
positive, 1 for negative; <fraction> is the bits (binary digits) in
1.bbb...b (where each b is 0 or 1), and <exponent> is a (probably
biased) exponent for 2 (or in at least one case, 16).

Such a value is guaranteed to be a whole integer already if the
exponent exceeds the number of bits.  That is:

	+/-  (1.)111...1  (x 2 sup) 52

can never be "something and a half" if there are only 52 digits
after the binary point, just as:

	1.23456E5

is the same as "123456.0" and has no fractional part (there are
only 5 digits after the decimal point, and the exponent is 5).

So, if your computer's "double" binary format uses 52 bits of
fraction, and the exponent is 52 or greater, floor(x) is just x (x
is itself an integer).  That leaves only those cases where the
exponent is less than 52.  In this case, we can still use the
properties of CPU's arithmetic to separate the pieces.

Suppose you have a value in the range [+0..2^52), and you add to
it a value that is exactly 2^52 (4503599627370496.0, but done in
the computer's internal double format, so that it is 1.0000 x 2^52).
If the result is forced into one of those same "double" formats,
so that it can only hold 52 fractional bits, all of the "irrelevant"
bits will fall out of the sum.  We can then subtract 2^52 again
and get only the relevant bits.  There may be *no* relevant bits,
in which case the result of subtraction should be zero.  This cannot
be written in ordinary C because the compiler might elect to use
a more precise internal representation, in which case irrelevant
fractional bits might remain.

Now, if you have IEEE arithmetic, you also have a "rounding
direction".  The rounding direction must not be set to "towards
negative infinity" because then 2^52 - 2^52 is not +0 but rather
-0.  So again we need to use the properties of the CPU to force
a proper rounding direction (probably trivial in assembly-language,
but again impossible in C).

Note that if the input number is negative, we have a whole different
set of conditions.  We also must worry about inputs being Infinity
and NaN, in IEEE code.

>If you allow use of other library functions then you can implement it
>in terms of modf() I think.

It can, but again you probably need to worry about Inf and NaN on
IEEE systems.  The best bet is usually to write all three in some
mix of machine-specific assembly and/or C.  You may find (and indeed
will find, on the i386) that there is an instruction that does most
of the work. :-)
--
In-Real-Life: Chris Torek, Berkeley Software Design Inc
El Cerrito, CA	Domain:	torek@bsdi.com	+1 510 234 3167
http://claw.bsdi.com/torek/  (not always up)	I report spam to abuse@.