src/libm/e_fmod.c
author Ryan C. Gordon <icculus@icculus.org>
Thu, 21 Apr 2016 03:16:44 -0400
changeset 11729 d1ce8396c356
parent 11685 1788f503254d
permissions -rw-r--r--
Initial shot at a renderer target for Apple's Metal API.

This isn't complete, but is enough to run testsprite2. It's currently
Mac-only; with a little work to figure out how to properly glue in a Metal
layer to a UIView, this will likely work on iOS, too.

This is only wired up to the configure script right now, and disabled by
default. CMake and Xcode still need their bits filled in as appropriate.
     1 /*
     2  * ====================================================
     3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     4  *
     5  * Developed at SunPro, a Sun Microsystems, Inc. business.
     6  * Permission to use, copy, modify, and distribute this
     7  * software is freely granted, provided that this notice
     8  * is preserved.
     9  * ====================================================
    10  */
    11 
    12 /*
    13  * __ieee754_fmod(x,y)
    14  * Return x mod y in exact arithmetic
    15  * Method: shift and subtract
    16  */
    17 
    18 #include "math_libm.h"
    19 #include "math_private.h"
    20 
    21 static const double one = 1.0, Zero[] = {0.0, -0.0,};
    22 
    23 double attribute_hidden __ieee754_fmod(double x, double y)
    24 {
    25 	int32_t n,hx,hy,hz,ix,iy,sx,i;
    26 	u_int32_t lx,ly,lz;
    27 
    28 	EXTRACT_WORDS(hx,lx,x);
    29 	EXTRACT_WORDS(hy,ly,y);
    30 	sx = hx&0x80000000;		/* sign of x */
    31 	hx ^=sx;		/* |x| */
    32 	hy &= 0x7fffffff;	/* |y| */
    33 
    34     /* purge off exception values */
    35 	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
    36 	  ((hy|((ly|-(int32_t)ly)>>31))>0x7ff00000))	/* or y is NaN */
    37 	    return (x*y)/(x*y);
    38 	if(hx<=hy) {
    39 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
    40 	    if(lx==ly)
    41 		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
    42 	}
    43 
    44     /* determine ix = ilogb(x) */
    45 	if(hx<0x00100000) {	/* subnormal x */
    46 	    if(hx==0) {
    47 		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
    48 	    } else {
    49 		for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
    50 	    }
    51 	} else ix = (hx>>20)-1023;
    52 
    53     /* determine iy = ilogb(y) */
    54 	if(hy<0x00100000) {	/* subnormal y */
    55 	    if(hy==0) {
    56 		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
    57 	    } else {
    58 		for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
    59 	    }
    60 	} else iy = (hy>>20)-1023;
    61 
    62     /* set up {hx,lx}, {hy,ly} and align y to x */
    63 	if(ix >= -1022)
    64 	    hx = 0x00100000|(0x000fffff&hx);
    65 	else {		/* subnormal x, shift x to normal */
    66 	    n = -1022-ix;
    67 	    if(n<=31) {
    68 	        hx = (hx<<n)|(lx>>(32-n));
    69 	        lx <<= n;
    70 	    } else {
    71 		hx = lx<<(n-32);
    72 		lx = 0;
    73 	    }
    74 	}
    75 	if(iy >= -1022)
    76 	    hy = 0x00100000|(0x000fffff&hy);
    77 	else {		/* subnormal y, shift y to normal */
    78 	    n = -1022-iy;
    79 	    if(n<=31) {
    80 	        hy = (hy<<n)|(ly>>(32-n));
    81 	        ly <<= n;
    82 	    } else {
    83 		hy = ly<<(n-32);
    84 		ly = 0;
    85 	    }
    86 	}
    87 
    88     /* fix point fmod */
    89 	n = ix - iy;
    90 	while(n--) {
    91 	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
    92 	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
    93 	    else {
    94 	    	if((hz|lz)==0) 		/* return sign(x)*0 */
    95 		    return Zero[(u_int32_t)sx>>31];
    96 	    	hx = hz+hz+(lz>>31); lx = lz+lz;
    97 	    }
    98 	}
    99 	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
   100 	if(hz>=0) {hx=hz;lx=lz;}
   101 
   102     /* convert back to floating value and restore the sign */
   103 	if((hx|lx)==0) 			/* return sign(x)*0 */
   104 	    return Zero[(u_int32_t)sx>>31];
   105 	while(hx<0x00100000) {		/* normalize x */
   106 	    hx = hx+hx+(lx>>31); lx = lx+lx;
   107 	    iy -= 1;
   108 	}
   109 	if(iy>= -1022) {	/* normalize output */
   110 	    hx = ((hx-0x00100000)|((iy+1023)<<20));
   111 	    INSERT_WORDS(x,hx|sx,lx);
   112 	} else {		/* subnormal output */
   113 	    n = -1022 - iy;
   114 	    if(n<=20) {
   115 		lx = (lx>>n)|((u_int32_t)hx<<(32-n));
   116 		hx >>= n;
   117 	    } else if (n<=31) {
   118 		lx = (hx<<(32-n))|(lx>>n); hx = sx;
   119 	    } else {
   120 		lx = hx>>(n-32); hx = sx;
   121 	    }
   122 	    INSERT_WORDS(x,hx|sx,lx);
   123 	    x *= one;		/* create necessary signal */
   124 	}
   125 	return x;		/* exact output */
   126 }
   127 
   128 /*
   129  * wrapper fmod(x,y)
   130  */
   131 #ifndef _IEEE_LIBM
   132 double fmod(double x, double y)
   133 {
   134 	double z = __ieee754_fmod(x, y);
   135 	if (_LIB_VERSION == _IEEE_ || isnan(y) || isnan(x))
   136 		return z;
   137 	if (y == 0.0)
   138 		return __kernel_standard(x, y, 27); /* fmod(x,0) */
   139 	return z;
   140 }
   141 #else
   142 strong_alias(__ieee754_fmod, fmod)
   143 #endif
   144 libm_hidden_def(fmod)