double a = 95.022, b = 0.01;
printf ("%lf - floor(%lf / %lf) * %lf = %.17lf\n", a, a, b, b, a -
floor(a / b) * b);
a = 95.021, b = 0.01;
printf ("%lf - floor(%lf / %lf) * %lf = %.17lf\n", a, a, b, b, a -
floor(a / b) * b);
a = 95.020, b = 0.01;
printf ("%lf - floor(%lf / %lf) * %lf = %.17lf\n", a, a, b, b, a -
floor(a / b) * b);
a = 95.022, b = 0.01;
printf ("fmod(%lf, %lf) = %.17lf\n", a, b, fmod(a, b));
a = 95.021, b = 0.01;
printf ("fmod(%lf, %lf) = %.17lf\n", a, b, fmod(a, b));
a = 95.020, b = 0.01;
printf ("fmod(%lf, %lf) = %.17lf\n", a, b, fmod(a, b));
generates this output:
95.022000 - floor(95.022000 / 0.010000) * 0.010000 =
0.00200000000000955
95.021000 - floor(95.021000 / 0.010000) * 0.010000 =
0.00100000000000477
95.020000 - floor(95.020000 / 0.010000) * 0.010000 =
0.00000000000000000
fmod(95.022000, 0.010000) = 0.00200000000000359
fmod(95.021000, 0.010000) = 0.00099999999999882
fmod(95.020000, 0.010000) = 0.00999999999999404
everything makes sense, except for the last line. why does fmod return
0.01 instead of 0.0?
> fmod(95.022000, 0.010000) = 0.00200000000000359
> fmod(95.021000, 0.010000) = 0.00099999999999882
> fmod(95.020000, 0.010000) = 0.00999999999999404
>
> everything makes sense, except for the last line. why does fmod return
> 0.01 instead of 0.0?
I don't know.
I get similar results with my homemade fmod.
/* BEGIN new.c output */
fs_fmod(95.022000, 0.010000) is 0.002000
fs_fmod(95.021000, 0.010000) is 0.001000
fs_fmod(95.020000, 0.010000) is 0.010000
/* END new.c output */
/* BEGIN new.c */
#include <stdio.h>
#include <float.h>
double fs_fmod(double x, double y);
int main(void)
{
puts("/* BEGIN new.c output */\n");
printf("fs_fmod(95.022000, 0.010000) is %f\n"
, fs_fmod(95.022000, 0.010000));
printf("fs_fmod(95.021000, 0.010000) is %f\n"
, fs_fmod(95.021000, 0.010000));
printf("fs_fmod(95.020000, 0.010000) is %f\n"
, fs_fmod(95.020000, 0.010000));
puts("\n/* END new.c output */");
return 0;
}
double fs_fmod(double x, double y)
{
double a, b;
const double c = x;
if (0 > c) {
x = -x;
}
if (0 > y) {
y = -y;
}
if (y != 0 && DBL_MAX >= y && DBL_MAX >= x) {
while (x >= y) {
a = x / 2;
b = y;
while (a >= b) {
b *= 2;
}
x -= b;
}
} else {
x = 0;
}
return 0 > c ? -x : x;
}
/* END new.c */
--
pete
Because none of the literals in your program can be represented
exactly in binary floating-point.
See section 14 of the comp.lang.c FAQ, <http://www.c-faq.com/>.
--
Keith Thompson (The_Other_Keith) ks...@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.
I take it your code is based on an accepted algorithm for computing
floating-point modulus values. Assuming fmod implementations do the
same thing, or something similar, it seems strange that fmod would
return something that seems that "far off"...if the return value was
0.00000000000452 or something like that, fine. but something that
rounds to 0.01...strange.
Here's a defintion of fmod from
http://www.opengroup.org/onlinepubs/007908799/xsh/fmod.html
The fmod() function returns the value x - i * y for some integer i such
that, if y is non-zero, the result has the same sign as x and magnitude
less than the magnitude of y.
strictly speaking, 0.00999999... is a valid return value since it's
less than 0.01...I haven't cranked through the computation, but maybe
if i went up by one, then the result is less than zero, and it can't
return a negative value for these operands, so somehow you're left with
0.00999999...
calc.exe gets it right <g>
First of all, it doesn't return 0.01. According to you, it returns
0.00999999999999404 which is "close to" 0.01, but is NOT 0.01.
The roundoff /truncation error encountered is thus 0.00000000000000596
Note that the first line is off by a similar amount, the only difference
being that the roundoff/truncation ended up slightly larger than the "exact"
answer rather than slightly smaller.
See comments others have made about the computer's inability to represent
0.01 exactly.
--
Fred L. Kleinschmidt
Boeing Associate Technical Fellow
Technical Architect, Software Reuse Project
In other words, if fmod(95.02, 0.01) returned -0.00000000000000404, or
0.00000000000000596, I wouldn't have posted.
calc.exe returns 0 for 95.02 mod 0.01, not 0.01. I'd guess that's
because calc.exe uses BCD. It still looks like fmod has an unexpected
boundary condition when x % y ~ 0.
> everything makes sense, except for the last line. why does fmod return
> 0.01 instead of 0.0?
>
Keeping Keith Thompson's reply in mind, experiment with the
program below. And while you're at it, take a look at the link.
/* a method for computing the single precision
binary representation of a decimal fraction */
/* see:
http://www.nuvisionmiami.com/books/asm/workbook/floating_tut.htm */
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
int pow2, mant_bit;
long double minuend, subtrahend, remainder, decfrac;
char mantissa[24];
if(argc != 2 || argv[1][0] != '.')
{
printf("usage: >pgm_name <decimal fraction>\n");
exit(EXIT_FAILURE);
}
decfrac = 0.0;
minuend = strtold(argv[1], NULL);
printf("\n\tstarting: %.12Lf\n\n", minuend);
for(pow2=2, mant_bit=0; mant_bit<23; pow2*=2, mant_bit++)
{
printf("%d", mant_bit + 1);
subtrahend = 1.0 / pow2;
if(minuend - subtrahend == 0.0)
{
printf("\tsubtracting %.12Lf\n", subtrahend);
remainder = minuend - subtrahend;
minuend = remainder;
printf("\tremainder = %.12Lf\n", remainder);
mantissa[mant_bit] = '1';
decfrac += subtrahend;
continue;
}
else if(minuend - subtrahend > 0.0)
{
printf("\tsubtracting %.12Lf\n", subtrahend);
remainder = minuend - subtrahend;
minuend = remainder;
printf("\tremainder = %.12Lf\n", remainder);
mantissa[mant_bit] = '1';
decfrac += subtrahend;
}
else
{
mantissa[mant_bit] = '0';
printf("\n");
}
}
mantissa[mant_bit] = '\0';
printf("\n\tbinary mantissa: .%s\n", mantissa);
printf("\tdecimal fraction: %.12Lf\n", decfrac);
return 0;
}
> There's a disconnect.
A disconnect with what? (See below).
Brian
--
Please quote enough of the previous message for context. To do so from
Google, click "show options" and use the Reply shown in the expanded
header.
printf("%a == %.100f\n", 95.02, 95.02);
printf("%a == %.100f\n", .01, .01);
[cutting off excess zeros in the pasted output]
0x1.7c147ae147ae1p+6 == 95.019999999999996020960679743438959121704101\
5625
0x1.47ae147ae147bp-7 == 0.0100000000000000002081668171172168513294309\
3776702880859375
_those_ are the numbers that fmod sees.
> printf("%a == %.100f\n", 95.02, 95.02);
> printf("%a == %.100f\n", .01, .01);
>
> [cutting off excess zeros in the pasted output]
>
> 0x1.7c147ae147ae1p+6 == 95.019999999999996020960679743438959121704101\
> 5625
>
> 0x1.47ae147ae147bp-7 == 0.0100000000000000002081668171172168513294309\
> 3776702880859375
>
> _those_ are the numbers that fmod sees.
Bullshit. Bogus precision. I don't believe a double on your machine
has 160 bits of precision. It appears, from your hex float output that
it has at least 49 and at most 53 bits of precision. Any precision
specifier greater than 17 is, then, utter crap.
But they are accurate representations of the values conveyed. Just
as the fractional value 1/16 is accurately represented by 0.0625.
--
"If you want to post a followup via groups.google.com, don't use
the broken "Reply" link at the bottom of the article. Click on
"show options" at the top of the article, then click on the
"Reply" at the bottom of the article headers." - Keith Thompson
More details at: <http://cfaj.freeshell.org/google/>
Also see <http://www.safalra.com/special/googlegroupsreply/>
If you expect us to know what you're talking about, please read
<http://cfaj.freeshell.org/google/> before posting another followup.
Thank you.
1/2^53 is
.00000000000000011102230246251565404236316680908203125
Just because it can be represented uniquely by 17 digits doesn't mean
that it's not useful to see what the exact value is, down to the last
'5' [any fractional value of a binary floating point number is always
going to end in '5']
> But they are accurate representations of the values conveyed. Just
> as the fractional value 1/16 is accurately represented by 0.0625.
but if you only have 4 bits of precision, any more than 2 digits is
bogus, and it should be 0.06
since 1..16/16 can be uniquely represented as
.06 .12 .18 .25 .31 .37 .43 .50 .56 .62 .68 .75 .81 .87 .93 1.00
in that case.
sorry, I wanted, but wasn't familiar enough with google's groups, to
quote prior messages. Thanks for the pointer.
my original problem was that 95.02 was being represented as
95.019999999 while 0.01 was represented as 0.0100000001 (in general
terms). so of course fmod returns something that rounded to 0.01.
I thought maybe 95.02 + DBL_EPSILON, when represented in binary
floating point, would never be a value less than 95.02, so fmod(95.02 +
DBL_EPSILON, 0.01) (or possibly fmod(95.02 + DBL_EPSILON, 0.01 -
DBL_EPSILON)) would return a value that rounded to 0. I was asking what
anyone thought of that. It's a moot point though, I tried it and it
doesn't work the way I hoped.
Right.
> I thought maybe 95.02 + DBL_EPSILON, when represented in binary
> floating point, would never be a value less than 95.02, so fmod(95.02 +
> DBL_EPSILON, 0.01) (or possibly fmod(95.02 + DBL_EPSILON, 0.01 -
> DBL_EPSILON)) would return a value that rounded to 0. I was asking what
> anyone thought of that. It's a moot point though, I tried it and it
> doesn't work the way I hoped.
DBL_EPSILON is the difference between 1.0 and the next representable
number above 1.0. The absolute resolution of a floating-point number
varies with magnitude; for example, the next representable number
above 2.0 is likely to be 2.0 + 2.0*DBL_EPSILON.
So adding DBL_EPSILON to 95.02 most likely won't do anything useful.
(And similarly, adding DBL_EPSILON to a small value such as 0.1 will
probably skip over several representable values.)
I just made it up.
I guessed that on most systems
that if the original value of y was doubled until
it was between x/2 and x, and then subtracted from x,
and that if that process was repeated until x was less than
the original value of y,
that only a minimum of accuracy would be lost.
--
pete