- What's the correct fmod()-algorithm - 11 Updates
| Bonita Montero <Bonita.Montero@gmail.com>: Apr 04 04:32PM +0200 If I implement fmod either with ... x / y - trunc( x / y ) ... or with .. x - trunc( x / y ) * y ... this doesn't always give correct values, f.e. if x is very large without a fraction and y is very samll. So what's the correct alogithm for fmod ? |
| Real Troll <real.troll@trolls.com>: Apr 04 03:49PM +0100 On 04/04/2021 15:32, Bonita Montero wrote: > ... this doesn't always give correct values, f.e. if x is very large > without a fraction and y is very samll. So what's the correct alogithm > for fmod ? No problems here using this example: #include <stdio.h> #include <math.h> int main() { float a, b; int c; a = 150.25; b = 27.67; c = 23; printf("Remainder of %f / %d is %lf\n", a, c, fmod(a,c)); printf("Remainder of %f / %f is %lf\n", a, b, fmod(a,b)); return(0); } The result should be: Remainder of 150.250000 / 23 is 12.250000 Remainder of 150.250000 / 27.670000 is 11.900000 |
| Bonita Montero <Bonita.Montero@gmail.com>: Apr 04 05:12PM +0200 > If I implement fmod either with ... > x / y - trunc( x / y ) This won't work for: x = 5290912749423342.0 and y = 10.0 The result is ... 0.18750000000000000 ... although it should be integral. > ... or with .. > x - trunc( x / y ) * y This won't work for: x = 52909127494233416 and y = 10.0 The result is ... -8.0 So ne naive solution doesn't work for such number-combinations. |
| Bonita Montero <Bonita.Montero@gmail.com>: Apr 04 05:14PM +0200 > The result should be: > Remainder of 150.250000 / 23 is 12.250000 > Remainder of 150.250000 / 27.670000 is 11.900000 Sorry, you're naive; look at my latest post. I'm looking for a fmod() -implementation and not for code which shows how to use fmod(), and it must be really tricky. |
| Bart <bc@freeuk.com>: Apr 04 05:23PM +0100 On 04/04/2021 16:12, Bonita Montero wrote: > ... although it should be integral. >> ... or with .. >> x - trunc( x / y ) * y This is the correct formula (but watch out for negative values). > The result is ... > -8.0 > So ne naive solution doesn't work for such number-combinations. I don't have trunc() on Windows. If I emulate it with a simple version the converts to int64 then I get the same result as you. How is trunc() implemented? If it actually returns an int, then it can't represent large intger values anyway (eg. 1e100). And if it simply converts double to long long, then it will not work for values needing more than 52 bits of precision. What do you get by printing the value of trunc(52909127494233416)? If that is wrong, then the fmod result will be wrong. |
| Bonita Montero <Bonita.Montero@gmail.com>: Apr 04 06:46PM +0200 > What do you get by printing the value of trunc(52909127494233416)? > If that is wrong, then the fmod result will be wrong. No, that's the magic of floating point arithemtics. |
| Bart <bc@freeuk.com>: Apr 04 06:25PM +0100 On 04/04/2021 17:46, Bonita Montero wrote: >> What do you get by printing the value of trunc(52909127494233416)? >> If that is wrong, then the fmod result will be wrong. > No, that's the magic of floating point arithemtics. The reason it goes wrong is this: x =52909127494233416.00000000000000000000 x/y =5290912749423342.00000000000000000000 trunc(x/y) =5290912749423342.00000000000000000000 (x/y)*y =52909127494233424.00000000000000000000 x/y (y is 10.0) should be 5290912749423341.6, but it gives ...3342.0. Probably because it's reaching the limits of precision. The difference between first and last figures is 8, rather than 6 (and with the wrong sign). I don't know how the internal fmod works; it might be at the bit level. But it will also start giving problems at a certain point. |
| Bonita Montero <Bonita.Montero@gmail.com>: Apr 04 07:28PM +0200 > x/y (y is 10.0) should be 5290912749423341.6, but it gives ...3342.0. > Probably because it's reaching the limits of precision. ... Right, the difference between x / y and floor( x / y ) becomes larger than the modulo should be or even negative. I'm currently checking the floor()-implementation of glibc, but that's very complex. |
| Paavo Helde <myfirstname@osa.pri.ee>: Apr 04 08:39PM +0300 04.04.2021 19:23 Bart kirjutas: >>> ... or with .. >>> x - trunc( x / y ) * y > This is the correct formula (but watch out for negative values). Arguably the correct definition for intdiv is floor(x/y), not trunc(x/y) (this avoids disruptions at zero), with the matching definition for modulus. Python has got this right, whereas C's fmod() is arguably wrong with negative numbers. >> x = 52909127494233416 and y = 10.0 >> The result is ... >> -8.0 I bet fmod() is carefully coded to return most accurate results for all operand values. If it were trivial to implement it via a simple formula, there would be no need to include it in the standard library. >> So ne naive solution doesn't work for such number-combinations. > I don't have trunc() on Windows. If I emulate it with a simple version > the converts to int64 then I get the same result as you. trunc() is present in MSVC, at least since VS2013. It returns floating-point: https://docs.microsoft.com/en-us/cpp/c-runtime-library/reference/trunc-truncf-truncl?view=msvc-160 |
| "Alf P. Steinbach" <alf.p.steinbach+usenet@gmail.com>: Apr 05 12:37AM +0200 On 04.04.2021 16:32, Bonita Montero wrote: > ... this doesn't always give correct values, f.e. if x is very large > without a fraction and y is very samll. So what's the correct alogithm > for fmod ? cppreference says << The double version of fmod behaves as if implemented as follows double fmod(double x, double y) { #pragma STDC FENV_ACCESS ON double result = std::remainder(std::fabs(x), (y = std::fabs(y))); if (std::signbit(result)) result += y; return std::copysign(result, x); } <url: https://en.cppreference.com/w/cpp/numeric/math/fmod> - Alf |
| Ian Collins <ian-news@hotmail.com>: Apr 05 10:38AM +1200 On 05/04/2021 04:23, Bart wrote: > I don't have trunc() on Windows. If I emulate it with a simple version > the converts to int64 then I get the same result as you. https://code.woboq.org/userspace/glibc/sysdeps/ieee754/dbl-64/s_trunc.c.html It's usually implement with a builtin (__builtin_trunc) in gcc. -- Ian. |
| You received this digest because you're subscribed to updates for this group. You can change your settings on the group membership page. To unsubscribe from this group and stop receiving emails from it send an email to comp.lang.c+++unsubscribe@googlegroups.com. |
No comments:
Post a Comment