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