# fmodf.c   [plain text]

```/*
*  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.
*
*
*/

#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;

//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

//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.)

//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 >>= 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;
}
*/
```