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