summaryrefslogtreecommitdiff log msg author committer range
path: root/src/math/remainder.c
blob: db176c882164ee5ea790c131c3e8a0fe1d5198f0 (plain)
 ```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 65 66 67 68 69 70 ``` ``````/* origin: FreeBSD /usr/src/lib/msun/src/e_remainder.c */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* remainder(x,p) * Return : * returns x REM p = x - [x/p]*p as if in infinite * precise arithmetic, where [x/p] is the (infinite bit) * integer nearest x/p (in half way case choose the even one). * Method : * Based on fmod() return x-[x/p]chopped*p exactlp. */ #include "libm.h" static const double zero = 0.0; double remainder(double x, double p) { int32_t hx,hp; uint32_t sx,lx,lp; double p_half; EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hp, lp, p); sx = hx & 0x80000000; hp &= 0x7fffffff; hx &= 0x7fffffff; /* purge off exception values */ if ((hp|lp) == 0) /* p = 0 */ return (x*p)/(x*p); if (hx >= 0x7ff00000 || /* x not finite */ (hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */ // FIXME: why long double? return ((long double)x*p)/((long double)x*p); if (hp <= 0x7fdfffff) x = fmod(x, p+p); /* now x < 2p */ if (((hx-hp)|(lx-lp)) == 0) return zero*x; x = fabs(x); p = fabs(p); if (hp < 0x00200000) { if (x + x > p) { x -= p; if (x + x >= p) x -= p; } } else { p_half = 0.5*p; if (x > p_half) { x -= p; if (x >= p_half) x -= p; } } GET_HIGH_WORD(hx, x); if ((hx&0x7fffffff) == 0) hx = 0; SET_HIGH_WORD(x, hx^sx); return x; } ``````