This code, compiled with visual studio .NET 2003,
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? 17 6215 jo************@comcast.net wrote: 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 jo************@comcast.net writes: This code, compiled with visual studio .NET 2003,
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?
Because none of the literals in your program can be represented
exactly in binary floatingpoint.
See section 14 of the comp.lang.c FAQ, <http://www.cfaq.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.
Thanks Pete.
I take it your code is based on an accepted algorithm for computing
floatingpoint 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/.../xsh/fmod.html
The fmod() function returns the value x  i * y for some integer i such
that, if y is nonzero, 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>
<jo************@comcast.net> wrote in message
news:11*********************@z14g2000cwz.googlegro ups.com... This code, compiled with visual studio .NET 2003,
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?
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
There's a disconnect. I'm fully aware that 0.00999999... is the closest
the machine can get to representing 0.01 after computing the
floatingpoint modulo. However, that's not my question...the question
is, why does fmod(95.02, 0.01) return "something that rounds to 0.01"
(quoting myself), whether that's 0.00999999999999404 or
0.01000000000000596, when the correct decimal modulus of 95.02 and 0.01
is 0.0, and most definitely not 0.01? 95.022 % 0.01 = 0.002. 95.021 %
0.01 = 0.001. 95.020 % 0.01 = 0.000, unless you're using fmod, which
returns something that rounds to 0.01...there's no way that much
difference is due only to limitations on binary floatingpoint
representation (or implicit type promotions).
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. jo************@comcast.net wrote:
<snip code> 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/a...oating_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;
} jo************@comcast.net wrote: 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.
On 20060214, jo************@comcast.net <jo************@comcast.net> wrote: There's a disconnect. I'm fully aware that 0.00999999... is the closest the machine can get to representing 0.01 after computing the floatingpoint modulo. However, that's not my question...the question is, why does fmod(95.02, 0.01) return "something that rounds to 0.01" (quoting myself), whether that's 0.00999999999999404 or 0.01000000000000596, when the correct decimal modulus of 95.02 and 0.01 is 0.0, and most definitely not 0.01? 95.022 % 0.01 = 0.002. 95.021 %
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.47ae147ae147bp7 == 0.0100000000000000002081668171172168513294309\
3776702880859375
_those_ are the numbers that fmod sees.
ah...thanks to all for straightening me out, and sorry for the
confusion. And here I thought I was smart...I was just looking at the
output, and didn't account for the fact that the inputs get adjusted
when represented in binary too.
Jordan Abel wrote: 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.47ae147ae147bp7 == 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.
anyone care to comment on whether adding epsilon to the numerator would
guarantee that 95.02 would never get represented as 95.0199999999...?
Martin Ambuhl wrote: Jordan Abel wrote:
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.47ae147ae147bp7 == 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/> jo************@comcast.net writes: anyone care to comment on whether adding epsilon to the numerator would guarantee that 95.02 would never get represented as 95.0199999999...?
If you expect us to know what you're talking about, please read
<http://cfaj.freeshell.org/google/> before posting another followup.
Thank you.

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.
On 20060214, CBFalconer <cb********@yahoo.com> wrote: Martin Ambuhl wrote: Jordan Abel wrote:
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.47ae147ae147bp7 == 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.
1/2^53 is
..000000000000000111022302462515654042363166809082 03125
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.
Keith Thompson wrote: jo************@comcast.net writes: anyone care to comment on whether adding epsilon to the numerator would guarantee that 95.02 would never get represented as 95.0199999999...?
If you expect us to know what you're talking about, please read <http://cfaj.freeshell.org/google/> before posting another followup. Thank you.
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. jo************@comcast.net writes:
[snip] 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.
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 floatingpoint 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.)

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. jo************@comcast.net wrote: Thanks Pete.
I take it your code is based on an accepted algorithm for computing floatingpoint modulus 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 This discussion thread is closed Replies have been disabled for this discussion. Similar topics
2 posts
views
Thread by Zunbeltz Izaola 
last post: by

8 posts
views
Thread by seia0106 
last post: by

2 posts
views
Thread by Lonnie Princehouse 
last post: by

6 posts
views
Thread by stau 
last post: by

2 posts
views
Thread by Gintautas 
last post: by

14 posts
views
Thread by Aaron Gray 
last post: by

7 posts
views
Thread by bummerland 
last post: by

12 posts
views
Thread by bsabiston 
last post: by
          