math: excess precision fix modf, modff, scalbn, scalbnf
[musl] / src / math / modf.c
1 #include "libm.h"
2
3 double modf(double x, double *iptr)
4 {
5         union {double x; uint64_t n;} u = {x};
6         uint64_t mask;
7         int e;
8
9         e = (int)(u.n>>52 & 0x7ff) - 0x3ff;
10
11         /* no fractional part */
12         if (e >= 52) {
13                 *iptr = x;
14                 if (e == 0x400 && u.n<<12 != 0) /* nan */
15                         return x;
16                 u.n &= (uint64_t)1<<63;
17                 return u.x;
18         }
19
20         /* no integral part*/
21         if (e < 0) {
22                 u.n &= (uint64_t)1<<63;
23                 *iptr = u.x;
24                 return x;
25         }
26
27         mask = (uint64_t)-1>>12 >> e;
28         if ((u.n & mask) == 0) {
29                 *iptr = x;
30                 u.n &= (uint64_t)1<<63;
31                 return u.x;
32         }
33         u.n &= ~mask;
34         *iptr = u.x;
35         STRICT_ASSIGN(double, x, x - *iptr);
36         return x;
37 }