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.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?

[pete]

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

[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.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 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.

[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.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 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*********************@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

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

[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_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;
}

[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.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.47ae147ae147bp-7 == 0.0100000000000000002081668171172168513294309\
3776702880859375

_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.

[Martin Ambuhl]

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

[joseph.p.doyle]

anyone care to comment on whether adding epsilon to the numerator would
guarantee that 95.02 would never get represented as 95.0199999999...?

[CBFalconer]

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.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/>

[Keith Thompson]

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.

[Jordan Abel]

On 2006-02-14, 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.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.

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.

[joseph.p.doyle]

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.

[Keith Thompson]

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 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.)

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

[pete]

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

Thanks Pete.

I take it your code is based on an accepted algorithm for computing
floating-point 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