darling-nix / src / libm / Source / ARM / fmodf.c

/*
 *  fmodf.c
 *  cLibm
 *
 *     Written by Jon Okada, started on December 7th, 1992.                      
 *     Modified by Paul Finlayson (PAF) for MathLib v2.                          
 *     Modified by A. Sazegari (ali) for MathLib v3.                             
 *     Modified and ported by Robert A. Murley (ram) for Mac OS X.               
 *	   Modified for cLibm, fixed a few edge cases, rewrote local_ funcs by Ian Ollmann.
 *
 *  Copyright 2007 Apple Inc. All rights reserved.
 *
 */
 
 #include <math.h>
 #include <float.h>
 #include <stdint.h>

#ifdef ARMLIBM_SET_FLAGS
 #include "required_arithmetic.h"
#endif

static inline int local_ilogbf( float x ) __attribute__ ((always_inline));
static inline int local_ilogbf( float x )
{
	union{ float f; uint32_t u;}u = {x};
	
	u.u &= 0x7fffffff;
	int32_t exp = u.u >> 23;
	
	if( __builtin_expect( (uint32_t) exp - 1U >= 254, 0 ) )
	{ // +-denorm, other possibilities: +-0, +-inf, NaN are screened out by caller
	
		u.u |= 0x3f800000U;
		u.f -= 1.0f;
		exp = u.u >> 23;
		
		return exp - (127+126);
	}

	return exp - 127;
}

static inline float local_scalbnf( float x, int n ) __attribute__ ((always_inline));
static inline float local_scalbnf( float x, int n )
{
	//create 2**n in double
	union{ uint64_t u; double d;} u = {  ( (int64_t) n + 1023) << 52 };
	double d = x;
 
	//scale by appropriate power of 2
	d *= u.d;
 
	//return result
	return (float) d;
}


float fmodf( float x, float y )
{
	union{ float f; uint32_t u; }ux = {x};
	union{ float f; uint32_t u; }uy = {y};

	uint32_t absx = ux.u & 0x7fffffffU;
	uint32_t absy = uy.u & 0x7fffffffU;
	if( absx - 1U >= 0x7f7fffffU || absy - 1U >= 0x7f7fffffU )
	{
		float fabsx = __builtin_fabsf(x);	
		float fabsy = __builtin_fabsf(y);

		// deal with NaN
		if( x != x || y != y )
			return x + y;

		//x = Inf or y = 0, return Invalid per IEEE-754
		if( fabsx == __builtin_inff() || 0.0f == y )
		{
#ifdef ARMLIBM_SET_FLAGS
			return required_add_float( __builtin_inff(), -__builtin_inff() );	//set invalid
#else
			return __builtin_nan("");
#endif
		}

		//handle trivial case
		if( fabsy == __builtin_inff() || 0.0f == x )
		{
#ifdef ARMLIBM_SET_FLAGS
			required_add_float( fabsx, 0.0f );	// signal underflow (that is, no flag set,
												// but exception occurs if unmasked) if x is denorm
#endif
			return x;
		}
	}
 
	if( absy >= absx )
	{
		if( absy == absx )
		{
			ux.u &= 0x80000000;
			return ux.f;
		}
		return x;
	}
 
	//extract exponent, mantissa
	int32_t expx = absx >> 23;
	int32_t expy = absy >> 23;
	int32_t sx = absx & 0x007fffff;			// x significand
	int32_t sy = absy & 0x007fffff;			// y significand
	if( 0 == expx )		//denormal x
	{
		//calculate shift to move leading bit to exponents field
		uint32_t shift = __builtin_clzl( absx ) - (8*sizeof(long) - 24);
		sx <<= shift;	//do the shift
		expx = 1-shift;	//adjust the biased exponent accordingly, -1 here for the implicit shift to make implicit denorm leading bit explicit
	}
	if( 0 == expy )		//denormal y
	{
		//calculate shift to move leading bit to exponents field
		uint32_t shift = __builtin_clzl( absy ) - (8*sizeof(long) - 24);
		sy <<= shift;	//do the shift
		expy = 1-shift;	//adjust the biased exponent accordingly, -1 here for the implicit shift to make implicit denorm leading bit explicit
	}
	//make implicit bits explicit
	sx |= 0x00800000;
	sy |= 0x00800000;
	
		
	// The idea here is to iterate until the exponent of x is the same as y
	// Calculate the number of bits we can safely shave off before we reach parity
	int32_t idiff = expx - expy;	
	int32_t shift = 0;
	int32_t mask;	

	//since idiff is always >= 0 here, it is safe to enter
	//We always shift by shift+1 here, so in the first iteration, the worst that can happen
	do
	{
		// move the leading bit of x to the 23rd bit
		sx <<= shift;

		//Keep track that we did that
		idiff += ~shift;			// idiff -= shift + 1, +1 for the shift below
	
		//The two values both have the 23rd bit set as the leading bit  (24 bit unsigned number)
		//subtract one from the other. This gives us a signed 23 bit number in the range { -0x00ffffff ... 0x00ffffff }
		sx -= sy;		

		//record the sign
		mask = sx >> 31;

		//shift x left 1 to restore a 24 bit quantity
		sx += sx;		//this is potentially 1 shift too far, which we patch up later

		//if negative, we add back sy to restore to postiveness. This is the same as x - y + 0.5y = - 0.5y 
		// instead of x-y. We've effectively hoisted the subtract that would have appeared in the next loop
		// iteration here, and saved a test+branch in exchange for a shift and and. (The add happens either way.)
		sx += sy & mask;	

		//figure out how far we need to shift sx to get the leading bit into the 23rd position
		shift = __builtin_clzl( sx ) - (8*sizeof(long) - 24);
	}
	while( idiff >= shift && sx != 0);
	
	//We may have gone too far 
	if( idiff < 0 )	
	{
		//If so, rewind a little.
		// if sx - sy succeeded, it was the right thing to do, the only thing we did wrong was the shift
		// if sx - sy yielded a negative number, then we shouldn't have done that either
		sx += sy & mask;
		sx >>= 1;
//debug code to make sure we didn't do something dumb here
		idiff = 0;
	}

	//We may have some bits left over to shift in. 
	sx <<= idiff;
	
//convert back to float
	if( 0 == sx )
	{
		ux.u &= 0x80000000;
		return ux.f;
	}
			
	//figure out how far we need to shift in order to move leading bit into exponent field
	shift = __builtin_clzl( sx ) - (8*sizeof(long) - 24);
	sx <<= shift;				// move leading bit into exponent field
	expy -= shift;				// account for the shift in the exponent
	sx &= 0x007fffff;			// remove leading bit
	sx |= ux.u & 0x80000000;	//or in sign
	if( expy > 0 )
	{
		ux.u = sx | (expy << 23);
		return ux.f;
	}
	
	//denormal
	expy += 126;
	ux.u = sx | (expy << 23);
	return ux.f * 0x1.0p-126f;
}

//old slower floating point based one.
/*{
	union{ float f; uint32_t u; }ux = {x};
	union{ float f; uint32_t u; }uy = {y};

	float fabsx = __builtin_fabsf(x);	
	float fabsy = __builtin_fabsf(y);
	if( (ux.u & 0x7fffffffU) - 1U >= 0x7f7fffffU || (uy.u & 0x7fffffffU) - 1U >= 0x7f7fffffU )
	{
		// deal with NaN
		if( x != x || y != y )
			return x + y;

		//x = Inf or y = 0, return Invalid per IEEE-754
		if( fabsx == __builtin_inff() || 0.0f == y )
		{
			return required_add_float( __builtin_inff(), -__builtin_inff() );	//set invalid
		}

		//handle trivial case
		if( fabsy == __builtin_inff() || 0.0f == x )
		{
			required_add_float( fabsx, 0.0f );			//signal underflow (that is, no flag set, but exception occurs if unmasked) if x is denorm
			return x;
		}
	}
 
	if( fabsy > fabsx )
		return x;
 
	//we have to work
	int32_t iscx = local_ilogbf(fabsx);
	int32_t iscy = local_ilogbf(fabsy);
	int32_t idiff = iscx - iscy;
	float x1, y1;
	x1 = fabsx;
	y1 = fabsy;
	int32_t i;
		
	if( idiff )
	{
		y1 = local_scalbnf( y1, -iscy );
		x1 = local_scalbnf( x1, -iscx );
		for( i = idiff; i != 0; i-- )
		{
			if( x1 >= y1 )
				x1 -= y1;
			x1 += x1;
		}
		x1 = local_scalbnf( x1, iscy );
	}
	if( x1 >= fabsy )
		x1 -= fabsy;
	 
	if( x < 0 )
		x1 = -x1;
	
	return x1;
}
*/