Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173/* * remainder.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_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 remainder( double x, double y ){
// deal with NaN if( x != x || y != y ) return x + y;
//x = Inf or y = 0, return Invalid per IEEE-754 double fabsx = __builtin_fabs(x); if( fabsx == __builtin_inf() || 0.0 == y ) {#ifdef ARMLIBM_SET_FLAGS return required_add_double( __builtin_inf(), -__builtin_inf() ); //set invalid#else return __builtin_nan("");#endif }
//handle trivial case double fabsy = __builtin_fabs(y); if( fabsy == __builtin_inf() || 0.0 == x ) {#ifdef ARMLIBM_SET_FLAGS required_add_double( fabsx, 0.0 ); // signal underflow (that is, no flag set, // but exception occurs if unmasked) if x is denorm#endif return x; } //we have to work int32_t iquo = 0; int32_t iscx = local_ilogb(fabsx); int32_t iscy = local_ilogb(fabsy); int32_t idiff = iscx - iscy; double x1, y1, z; x1 = fabsx; y1 = fabsy; if( idiff >= 0 ) { int32_t i; if( idiff ) { y1 = local_scalbn( y1, -iscy ); x1 = local_scalbn( x1, -iscx ); for( i = idiff; i != 0; i-- ) { if( x1 >= y1 ) { x1 -= y1; iquo += 1; } iquo += iquo; x1 += x1; } x1 = local_scalbn( x1, iscy ); } if( x1 >= fabsy ) { x1 -= fabsy; iquo += 1; } } if( x1 < 0x1.0p1022 ) { z = x1 + x1; if( (z > fabsy) || ((z == fabsy) & ((iquo & 1) != 0 ))) x1 -= fabsy; } else { // x1 is only this large if y is even bigger, so we can safely divide y by 2 double halfy = 0.5*fabsy; if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 ))) x1 -= fabsy; } if( x < 0 ) x1 = -x1; return x1;}