# fmod.c   [plain text]

```/*
*  fmod.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.
*     Modified for armLibm, removed code to set flags as this is a soft-float fallback only.
*
*  Copyright 2007 Apple Inc. All rights reserved.
*
*/

#if defined(__SOFTFP__)

#include <math.h>
#include <float.h>
#include <stdint.h>

static inline int local_ilogb( double x ) __attribute__ ((always_inline));
static inline int local_ilogb( double x )
{
union{ double d; uint64_t u;}u = {x};

u.u &= 0x7fffffffffffffffULL;
int32_t exp = (int32_t) (u.u >> 52);

if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
{ // +-denorm, the other possibilities: +0 +-inf, NaN are screend out by caller
u.u |= 0x3ff0000000000000ULL;
u.d -= 1.0;
exp = (int32_t) (u.u >> 52);

return exp - (1023+1022);
}

return exp - 1023;
}

static inline  double local_scalbn( double x, int n );
static inline  double local_scalbn( double x, int n )
{
union{ uint64_t u; double d;} u;

uint32_t absn = n >> (8*sizeof(n)-1);
absn = (n ^ absn) - absn;

if( absn > 1022 )
{
// step = copysign( 1022, n )
int signMask = n >> (8 * sizeof( int ) - 1);
int step = (1022 ^ signMask) - signMask;

u.u = ( (int64_t) step + 1023ULL) << 52;

if( n < 0 )
{
do
{
x *= u.d;
n -= step;
}while( n < -1022 );
}
else
{
do
{
x *= u.d;
n -= step;
}while( n > 1022 );
}
}

//create 2**n in double
u.u = ( (int64_t) n + 1023ULL) << 52;

//scale by appropriate power of 2
x *= u.d;

//return result
return x;
}

double fmod( double x, double y )
{
union{ double d; uint64_t u;}ux = {x};
union{ double d; uint64_t u;}uy = {y};

uint64_t absx = ux.u & ~0x8000000000000000ULL;
uint64_t absy = uy.u & ~0x8000000000000000ULL;

if( absx-1ULL >= 0x7fefffffffffffffULL || absy-1ULL >= 0x7fefffffffffffffULL )
{
double fabsx = __builtin_fabs(x);
double fabsy = __builtin_fabs(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_inf() || 0.0 == y )
{
return __builtin_nan("");
}

//handle trivial case
if( fabsy == __builtin_inf() || 0.0 == x )
{
return x;
}
}

if( absy >= absx )
{
if( absy == absx )
{
ux.u ^= absx;
return ux.d;
}

return x;
}

int32_t expx = absx >> 52;
int32_t expy = absy >> 52;
int64_t sx = absx & 0x000fffffffffffff;
int64_t sy = absy & 0x000fffffffffffff;

if( 0 == expx )
{
uint32_t shift = __builtin_clzll( absx ) - (64-53);
sx <<= shift;
expx = 1 - shift;
}

if( 0 == expy )
{
uint32_t shift = __builtin_clzll( absy ) - (64-53);
sy <<= shift;
expy = 1 - shift;
}
sx |= 0x0010000000000000ULL;
sy |= 0x0010000000000000ULL;

int32_t idiff = expx - expy;
int32_t shift = 0;
int64_t mask;

do
{
sx <<= shift;
idiff += ~shift;
sx -= sy;
mask = sx >> 63;
sx += sx;
sx += sy & mask;
shift = __builtin_clzll( sx ) - (64-53);
}
while( idiff >= shift && sx != 0LL );

if( idiff < 0 )
{
sx += sy & mask;
sx >>= 1;
idiff = 0;
}

sx <<= idiff;

if( 0 == sx )
{
ux.u &= 0x8000000000000000;
return ux.d;
}

shift = __builtin_clzll( sx ) - (64-53);
sx <<= shift;
expy -= shift;
sx &= 0x000fffffffffffffULL;
sx |= ux.u & 0x8000000000000000ULL;
if( expy > 0 )
{
ux.u = sx | ((int64_t) expy << 52);
return ux.d;
}

expy += 1022;
ux.u = sx | ((int64_t) expy << 52);
return ux.d * 0x1.0p-1022;

}

#endif

```