Conversion, rounding, and truncation
Tom Cunningham
and discussion (and some dross), I thought I'd throw another catalytic
cookie into the ring. Here is a C code fragment:
double d;
long l;
float f;
double pow();
d = pow(2.0,3.0);
l = d;
printf("d = %f, l = %ld\n", d, l);
f = d;
l = f;
printf("f = %f, l = %ld\n", d, l);
The output from code generated by the Sun 3 compiler (4.2BSD), Microsoft C,
and VAX VMS C is:
d = 8.000000, l = 8
d = 8.000000, l = 8
Output from a Unix SYSV compiler is:
d = 8.000000, l = 7
d = 8.000000, l = 8
Now, K&R clearly state on pg. 42 that float to int causes truncation of
any fractional part, and I assume this applies to the double-to-long case
as well. Moreover, they say that double is converted to float by
rounding, which I guess would explain the second SYSV output line.
What I don't understand is why the formatting routines appear to round
the floating point number, since it looks like the pow() in the SYSV
implementation returns something not quite 8.0. Am I missing something
here? If this has already been discussed ad nauseum, someone please
clue me.
Tom Cunningham "Good, fast, cheap -- select two."
USPS: Motorola Inc. 6501 William Cannon Dr. W. Austin, TX 78735-8598
UUCP: {ihnp4,seismo,ctvax,gatech}!ut-sally!oakhill!tomc
Phone: 512-440-2953
Daniel R. Levy
> double d;
> long l;
> float f;
> double pow();
> d = pow(2.0,3.0);
> l = d;
> printf("d = %f, l = %ld\n", d, l);
> f = d;
> l = f;
> printf("f = %f, l = %ld\n", d, l);
>
>The output from code generated by the Sun 3 compiler (4.2BSD), Microsoft C,
>and VAX VMS C is:
> d = 8.000000, l = 8
> d = 8.000000, l = 8
>Output from a Unix SYSV compiler is:
> d = 8.000000, l = 7
> d = 8.000000, l = 8
>implementation returns something not quite 8.0.
You have guessed right. pow(x,y) is basically implemented as exp(y*log(x)),
with domain and range checking too of course. Depending on the exact implemen-
tation of exp() and log(), this can give results which are slightly off even
when an exact solution would be representable (as in your example). As a
guru explained to me, this is done, even if it would be possible to special-
case exactly-integral exponents to pow(), in order to preserve the local
continuity of the function. For an example of why this would be wanted,
consider that a program which iterates to find a desired value of some
complicated, but supposedly continuous, function containing pow() could
conceivably get hung up in an infinite loop if pow(2.0,3.0+delta) < pow(2.0,3.0)
(delta being a small increment greater than zero) due to the latter
being calculated in a different way.
--
------------------------------- Disclaimer: The views contained herein are
| dan levy | yvel nad | my own and are not at all those of my em-
| an engihacker @ | ployer or the administrator of any computer
| at&t computer systems division | upon which I may hack.
| skokie, illinois |
-------------------------------- Path: ..!{akgua,homxb,ihnp4,ltuxa,mvuxa,
go for it! allegra,ulysses,vax135}!ttrdc!levy
David Hough
> with domain and range checking too of course.
> Depending on the exact implemen-
> tation of exp() and log(), this can give results which are slightly off even
> guru explained to me, this is done, even if it would be possible to special-
> case exactly-integral exponents to pow(), in order to preserve the local
> continuity of the function.
It is true that the simplest way to implement pow(x,y) is exp(y*log(x)).
That is not the best way. There is no inherent reason why x**y for
small integer values x and y
can not be computed correctly and without losing monotonicity.
The pow() in recent Sun software releases
and in 4.3 BSD works pretty well. In contrast most Unix implementations
of libm functions are pretty poor - Gene Spafford (gatech!spaf) wrote a book
about them.
The paranoia program distributed by Richard Karpinski
(ucbvax!ucsfcgl!cca.ucsf!dick) tests computer arithmetic, including
pow() at small integer values.
David Hough dho...@sun.com