fmod

[joseph.p.doyle]

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.0020000000000 0955
95.021000 - floor(95.021000 / 0.010000) * 0.010000 =
0.0010000000000 0477
95.020000 - floor(95.020000 / 0.010000) * 0.010000 =
0.0000000000000 0000
fmod(95.022000, 0.010000) = 0.0020000000000 0359
fmod(95.021000, 0.010000) = 0.0009999999999 9882
fmod(95.020000, 0.010000) = 0.0099999999999 9404

everything makes sense, except for the last line. why does fmod return
0.01 instead of 0.0?

[pete]

jo************@ comcast.net wrote:

fmod(95.022000, 0.010000) = 0.0020000000000 0359
fmod(95.021000, 0.010000) = 0.0009999999999 9882
fmod(95.020000, 0.010000) = 0.0099999999999 9404

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.0220 00, 0.010000) is 0.002000
fs_fmod(95.0210 00, 0.010000) is 0.001000
fs_fmod(95.0200 00, 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.0220 00, 0.010000));
printf("fs_fmod (95.021000, 0.010000) is %f\n"
, fs_fmod(95.0210 00, 0.010000));
printf("fs_fmod (95.020000, 0.010000) is %f\n"
, fs_fmod(95.0200 00, 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

[Keith Thompson]

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.0020000000000 0955
95.021000 - floor(95.021000 / 0.010000) * 0.010000 =
0.0010000000000 0477
95.020000 - floor(95.020000 / 0.010000) * 0.010000 =
0.0000000000000 0000
fmod(95.022000, 0.010000) = 0.0020000000000 0359
fmod(95.021000, 0.010000) = 0.0009999999999 9882
fmod(95.020000, 0.010000) = 0.0099999999999 9404

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 floating-point.

See section 14 of the comp.lang.c FAQ, <http://www.c-faq.com/>.

--
Keith Thompson (The_Other_Keit h) 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.

[joseph.p.doyle]

Thanks Pete.

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.0000000000045 2 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 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>

[Fred Kleinschmidt]

<jo************ @comcast.net> wrote in message
news:11******** *************@z 14g2000cwz.goog legroups.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.0020000000000 0955
95.021000 - floor(95.021000 / 0.010000) * 0.010000 =
0.0010000000000 0477
95.020000 - floor(95.020000 / 0.010000) * 0.010000 =
0.0000000000000 0000
fmod(95.022000, 0.010000) = 0.0020000000000 0359
fmod(95.021000, 0.010000) = 0.0009999999999 9882
fmod(95.020000, 0.010000) = 0.0099999999999 9404

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.0099999999999 9404 which is "close to" 0.01, but is NOT 0.01.
The roundoff /truncation error encountered is thus 0.0000000000000 0596

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

[joseph.p.doyle]

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
floating-point 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.0099999999999 9404 or
0.0100000000000 0596, 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 floating-point
representation (or implicit type promotions).

In other words, if fmod(95.02, 0.01) returned -0.0000000000000 0404, or
0.0000000000000 0596, 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.

[John Smith]

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_FAILU RE);
}
decfrac = 0.0;
minuend = strtold(argv[1], NULL);
printf("\n\tsta rting: %.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("\tsubtr acting %.12Lf\n", subtrahend);
remainder = minuend - subtrahend;
minuend = remainder;
printf("\tremai nder = %.12Lf\n", remainder);
mantissa[mant_bit] = '1';
decfrac += subtrahend;
continue;
}
else if(minuend - subtrahend > 0.0)
{
printf("\tsubtr acting %.12Lf\n", subtrahend);
remainder = minuend - subtrahend;
minuend = remainder;
printf("\tremai nder = %.12Lf\n", remainder);
mantissa[mant_bit] = '1';
decfrac += subtrahend;
}
else
{
mantissa[mant_bit] = '0';
printf("\n");
}
}
mantissa[mant_bit] = '\0';

printf("\n\tbin ary mantissa: .%s\n", mantissa);
printf("\tdecim al fraction: %.12Lf\n", decfrac);

return 0;
}

[Default User]

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.

[Jordan Abel]

On 2006-02-14, 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
floating-point 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.0099999999999 9404 or
0.0100000000000 0596, 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.7c147ae147a e1p+6 == 95.019999999999 996020960679743 438959121704101 \
5625

0x1.47ae147ae14 7bp-7 == 0.0100000000000 000002081668171 172168513294309 \
377670288085937 5

_those_ are the numbers that fmod sees.

[joseph.p.doyle]

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.