fmod()

Well, I was working on a structure for my own data type, when it came time to do modulus.

I was talking to a friend of mine, who recommended fmod() from math.h, but he told me it was painfully slow.

So I decided to make my own and see if it was faster. Sure enough, the end result operated five times faster than the one in math.h.

Code:

---

double mfmod(double x,double y) { double a; return ((a=x/y)-(int)a)*y; }

---

Also functional, but perhaps slower(?) is the following:

Code:

---

double mfmod(double x,double y) { return x-((int)(x/y))*y; }

---

The timing code was as follows.

Code:

---

#include <stdio.h>
#include <stdlib.h>
#include "fmod.h"
#include <time.h>

int main(int argc, char *argv[])
{
  int x,y;
    double t,u;
    int i;
   
    x=clock();
    for (i=0;i<10000000;i++) u=mfmod(100000.5,3);
    y=clock();
   
    printf("Mine: &#37;f ms\r\n\r\n",(((double)y-(double)x)*1000.0)/(double)CLOCKS_PER_SEC);
   
    x=clock();
    for (i=0;i<10000000;i++) t=fmod(100000.5,3);
    y=clock();
   
    printf("STDC: %f ms\r\n",(((double)y-(double)x)*1000.0)/(double)CLOCKS_PER_SEC);
   
   
   
    printf("t: %f u: %f",t,u);
 
  system("PAUSE");       
  return 0;
}

---

Have I missed something? Why does mine work faster?

I'm guessing that the compiler can optimize yours, which is noticeably error-checking-free, while it is unable to optimize the "hidden" version in the library. Case in point: I compiled (using MinGW's gcc, with optimization set to -O3) and tested your program above: your loop ran in 15 ms, while the library version ran in 735 ms, according to the output.

ETA: I should say that no-error-checking is not necessarily a bad thing here, since the x/y will also set error flags, and I tested and got the same result "1.#IND00" from yours and the library fmod as well.

Ah. I set the optimization as high as it would go.

Mine: 0.000000 ms
STDC: 16.000000 ms

So I guess that is likely a cause.

Either way, as common sense dictates, I'm going to stick with my function for the operator, since it's self contained anyway.

Quote:

---

Originally Posted by Josh@Dreamland

Code:

---

double mfmod(double x,double y) { double a; return ((a=x/y)-(int)a)*y; }

---

---

I'm not sure, but isn't this undefined behavior ? I mean, is there something in the C standard that says the variable a will always have the value of "x / y" before being casted to an int ?

error checking == slow things down ??

Quote:

---

Originally Posted by foxman

I'm not sure, but isn't this undefined behavior ? I mean, is there something in the C standard that says the variable a will always have the value of "x / y" before being casted to an int ?

---

You are quite correct. This is undefined.

A variable cannot be assigned to and used in the same expression, except when it is only use to compute its new value (eg: a=b+a), and when there is a sequence point between the two uses (eg: b < (a=32) && a<7).

Hahaha. My syntax derivatives often generate that kind of reply.

And also, I can tell you and show you what's gonna happen. The memory is already allocated for that double, and thanks to the left to right order of things, by the time it recalls the value it will have already been assigned. It works, and there's a reason behind it working.

If you like, reword it. It just costs about 20 extra milliseconds in the 10mil loop.

audinue--
Yeah. An if is unrealistically slow, which is why I only include error checking in debug mode.

Quote:

---

Originally Posted by audinue

error checking == slow things down ??

---

Yes and no, I guess; here, we're checking that y != 0; I seem to recall that jump on 0 is the easiest jump there is, but it still takes one processor "tick". In the normal scheme of things, who cares; here, we're doing it 10000000 times, so the difference is noticable, although 16 milliseconds in ten million iterations is still not huge.

Either way, gentlemen (!!)

Putting the C Bible down, and looking at what I'm actually adding to the struct (which I'm actually coding in C++ anyway), this is what I'm leaving it as:

Code:

---

double          var::operator% (double y)                          { serr(0,values[0]->type); double a=values[0]->realval/y; a-=(int)a; return a*y; }

---

Just so you know, serr is a macro defined as follows:

Code:

---

#if DEBUGMODE && SHOWERRORS
#define terr(t1,t2) if (t1 != t2) show_error("Performing operations on variables of different trypes.",0);
#define serr(t1,t2) if (t1 == 1 || t2 == 1) show_error("Performing invalid operations on string types.",0); else if (t1 != t2) show_error("Performing operations on variables of different trypes.",0);
#else
#define terr(t1,t2)
#define serr(t1,t2)
#endif

---

Quote:

---

Originally Posted by Josh@Dreamland

And also, I can tell you and show you what's gonna happen. The memory is already allocated for that double, and thanks to the left to right order of things, by the time it recalls the value it will have already been assigned. It works, and there's a reason behind it working.

---

This does not answer to my question. It's not because the - operator is associative from left to right that it necessarily means the left operand will be evaluated before the right operand at the assembly level (except if the standard says it has to). And it's not because it works on your compiler that the result will be the same on every compiler. That's the whole point of undefined behavior.

Anyway, the last version of your code doesn't rely on such thing, so no need to care anymore.

Yep. So I guess in the end all that really matters that it's working, and it's working pretty fast.

Though I can't see much difference between compilers, the code is actually for a large development environment that will be using GCC as its compiler throughout.

At any rate, all's well. I'm going to finish the other operators and call it a code.

From glibc-2.7 w_fmod.c:

Code:

---

#ifdef __STDC__
        double __fmod(double x, double y)        /* wrapper fmod */
#else
        double __fmod(x,y)                /* wrapper fmod */
        double x,y;
#endif
{
#ifdef _IEEE_LIBM
        return __ieee754_fmod(x,y);
#else
        double z;
        z = __ieee754_fmod(x,y);
        if(_LIB_VERSION == _IEEE_ ||__isnan(y)||__isnan(x)) return z;
        if(y==0.0) {
                return __kernel_standard(x,y,27); /* fmod(x,0) */
        } else
            return z;
#endif
}

---

And __ieee754_fmod() from e_fmod.c in the sysdeps/ieee754/dbl-64 folder:

Code:

---

#ifdef __STDC__
static const double one = 1.0, Zero[] = {0.0, -0.0,};
#else
static double one = 1.0, Zero[] = {0.0, -0.0,};
#endif

#ifdef __STDC__
        double __ieee754_fmod(double x, double y)
#else
        double __ieee754_fmod(x,y)
        double x,y ;
#endif
{
        int32_t n,hx,hy,hz,ix,iy,sx,i;
        u_int32_t lx,ly,lz;

        EXTRACT_WORDS(hx,lx,x);
        EXTRACT_WORDS(hy,ly,y);
        sx = hx&0x80000000;                /* sign of x */
        hx ^=sx;                /* |x| */
        hy &= 0x7fffffff;        /* |y| */

    /* purge off exception values */
        if((hy|ly)==0||(hx>=0x7ff00000)||        /* y=0,or x not finite */
          ((hy|((ly|-ly)>>31))>0x7ff00000))        /* or y is NaN */
            return (x*y)/(x*y);
        if(hx<=hy) {
            if((hx<hy)||(lx<ly)) return x;        /* |x|<|y| return x */
            if(lx==ly)
                return Zero[(u_int32_t)sx>>31];        /* |x|=|y| return x*0*/
        }

    /* determine ix = ilogb(x) */
        if(hx<0x00100000) {        /* subnormal x */
            if(hx==0) {
                for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
            } else {
                for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
            }
        } else ix = (hx>>20)-1023;

    /* determine iy = ilogb(y) */
        if(hy<0x00100000) {        /* subnormal y */
            if(hy==0) {
                for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
            } else {
                for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
            }
        } else iy = (hy>>20)-1023;

    /* set up {hx,lx}, {hy,ly} and align y to x */
        if(ix >= -1022)
            hx = 0x00100000|(0x000fffff&hx);
        else {                /* subnormal x, shift x to normal */
            n = -1022-ix;
            if(n<=31) {
                hx = (hx<<n)|(lx>>(32-n));
                lx <<= n;
            } else {
                hx = lx<<(n-32);
                lx = 0;
            }
        }
        if(iy >= -1022)
            hy = 0x00100000|(0x000fffff&hy);
        else {                /* subnormal y, shift y to normal */
            n = -1022-iy;
            if(n<=31) {
                hy = (hy<<n)|(ly>>(32-n));
                ly <<= n;
            } else {
                hy = ly<<(n-32);
                ly = 0;
            }
        }

    /* fix point fmod */
        n = ix - iy;
        while(n--) {
            hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
            if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
            else {
                    if((hz|lz)==0)                /* return sign(x)*0 */
                    return Zero[(u_int32_t)sx>>31];
                    hx = hz+hz+(lz>>31); lx = lz+lz;
            }
        }
        hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
        if(hz>=0) {hx=hz;lx=lz;}

    /* convert back to floating value and restore the sign */
        if((hx|lx)==0)                        /* return sign(x)*0 */
            return Zero[(u_int32_t)sx>>31];       
        while(hx<0x00100000) {                /* normalize x */
            hx = hx+hx+(lx>>31); lx = lx+lx;
            iy -= 1;
        }
        if(iy>= -1022) {        /* normalize output */
            hx = ((hx-0x00100000)|((iy+1023)<<20));
            INSERT_WORDS(x,hx|sx,lx);
        } else {                /* subnormal output */
            n = -1022 - iy;
            if(n<=20) {
                lx = (lx>>n)|((u_int32_t)hx<<(32-n));
                hx >>= n;
            } else if (n<=31) {
                lx = (hx<<(32-n))|(lx>>n); hx = sx;
            } else {
                lx = hx>>(n-32); hx = sx;
            }
            INSERT_WORDS(x,hx|sx,lx);
            x *= one;                /* create necessary signal */
        }
        return x;                /* exact output */

---

When it comes to the glibc math I assume it to pretty hyperoptimized to no end, though I still get same numbers as you. Could be some rare exceptions being considered in all those if-cases that cause the slowdown, but I can't say much from just the comments alone.

O___o

Those guys certainly put a lot of work into what they do.

Hm. I think I'm comfortable with my one liner, if they're about the same speed. It's the only math function I actually need, since I don't intend to allow bit shifting with floats. Would be a bit silly to include the whole thing for that.

Thanks for your help. ^_^

If that scares you, you should see the trigonometry and exponential functions. I can't even begin to comprehend what goes on there.