summaryrefslogtreecommitdiff log msg author committer range
path: root/src/math/remainderf.c
blob: c17bb4f434d2f81eacc24f04838a26fc999c727d (plain) (blame)
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 ``` ``````/* origin: FreeBSD /usr/src/lib/msun/src/e_remainderf.c */ /* * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include "libm.h" static const float zero = 0.0; float remainderf(float x, float p) { int32_t hx,hp; uint32_t sx; float p_half; GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hp, p); sx = hx & 0x80000000; hp &= 0x7fffffff; hx &= 0x7fffffff; /* purge off exception values */ if (hp == 0) /* p = 0 */ return (x*p)/(x*p); if (hx >= 0x7f800000 || hp > 0x7f800000) /* x not finite, p is NaN */ // FIXME: why long double? return ((long double)x*p)/((long double)x*p); if (hp <= 0x7effffff) x = fmodf(x, p + p); /* now x < 2p */ if (hx - hp == 0) return zero*x; x = fabsf(x); p = fabsf(p); if (hp < 0x01000000) { if (x + x > p) { x -= p; if (x + x >= p) x -= p; } } else { p_half = (float)0.5*p; if (x > p_half) { x -= p; if (x >= p_half) x -= p; } } GET_FLOAT_WORD(hx, x); if ((hx & 0x7fffffff) == 0) hx = 0; SET_FLOAT_WORD(x, hx ^ sx); return x; } ``````