Thread: fmod()

  1. #1
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10

    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?
    Last edited by Josh@Dreamland; 07-12-2008 at 05:10 PM.

  2. #2
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    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.
    Last edited by tabstop; 07-12-2008 at 05:24 PM.

  3. #3
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10
    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.

  4. #4
    Chinese pâté foxman's Avatar
    Join Date
    Jul 2007
    Location
    Canada
    Posts
    404
    Quote Originally Posted by Josh@Dreamland View Post
    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 ?
    I hate real numbers.

  5. #5
    Ugly C Lover audinue's Avatar
    Join Date
    Jun 2008
    Location
    Indonesia
    Posts
    489
    error checking == slow things down ??

  6. #6
    Registered User
    Join Date
    Apr 2006
    Posts
    2,149
    Quote Originally Posted by foxman View Post
    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).
    Last edited by King Mir; 07-12-2008 at 07:08 PM.
    It is too clear and so it is hard to see.
    A dunce once searched for fire with a lighted lantern.
    Had he known what fire was,
    He could have cooked his rice much sooner.

  7. #7
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10
    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.
    Last edited by Josh@Dreamland; 07-12-2008 at 07:16 PM.

  8. #8
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    Quote Originally Posted by audinue View Post
    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.

  9. #9
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10
    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&#37; (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

  10. #10
    Chinese pâté foxman's Avatar
    Join Date
    Jul 2007
    Location
    Canada
    Posts
    404
    Quote Originally Posted by Josh@Dreamland View Post
    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.
    I hate real numbers.

  11. #11
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10
    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.

  12. #12
    Registered User OnionKnight's Avatar
    Join Date
    Jan 2005
    Posts
    555
    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.

  13. #13
    Registered User Josh@Dreamland's Avatar
    Join Date
    Dec 2007
    Posts
    10
    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. ^_^

  14. #14
    Registered User OnionKnight's Avatar
    Join Date
    Jan 2005
    Posts
    555
    If that scares you, you should see the trigonometry and exponential functions. I can't even begin to comprehend what goes on there.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. FMOD implimentation, how'd you do it?
    By Jeremy G in forum Game Programming
    Replies: 10
    Last Post: 06-12-2004, 12:04 AM
  2. error with fmod() function
    By minesweeper in forum C++ Programming
    Replies: 3
    Last Post: 08-30-2003, 05:24 AM
  3. Help With Fmod Music Lib!
    By Blizzarddog in forum Game Programming
    Replies: 0
    Last Post: 05-29-2003, 12:29 PM
  4. Fmod and windows programming
    By AtomRiot in forum Windows Programming
    Replies: 1
    Last Post: 01-06-2003, 03:13 PM
  5. fmod problem
    By yusiye in forum C Programming
    Replies: 3
    Last Post: 08-01-2002, 07:28 PM

Tags for this Thread