math: fix tgamma to raise underflow for large negative values
[oweals/musl.git] / src / math / remquof.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_remquof.c */
2 /*-
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 /*
13  * Return the IEEE remainder and set *quo to the last n bits of the
14  * quotient, rounded to the nearest integer.  We choose n=31 because
15  * we wind up computing all the integer bits of the quotient anyway as
16  * a side-effect of computing the remainder by the shift and subtract
17  * method.  In practice, this is far more bits than are needed to use
18  * remquo in reduction algorithms.
19  */
20
21 #include "libm.h"
22
23 static const float Zero[] = {0.0, -0.0,};
24
25 float remquof(float x, float y, int *quo)
26 {
27         int32_t n,hx,hy,hz,ix,iy,sx,i;
28         uint32_t q,sxy;
29
30         GET_FLOAT_WORD(hx, x);
31         GET_FLOAT_WORD(hy, y);
32         sxy = (hx ^ hy) & 0x80000000;
33         sx = hx & 0x80000000;   /* sign of x */
34         hx ^= sx;               /* |x| */
35         hy &= 0x7fffffff;       /* |y| */
36
37         /* purge off exception values */
38         if (hy == 0 || hx >= 0x7f800000 || hy > 0x7f800000) /* y=0,NaN;or x not finite */
39                 return (x*y)/(x*y);
40         if (hx < hy) {       /* |x| < |y| return x or x-y */
41                 q = 0;
42                 goto fixup;
43         } else if(hx==hy) {  /* |x| = |y| return x*0*/
44                 *quo = sxy ? -1 : 1;
45                 return Zero[(uint32_t)sx>>31];
46         }
47
48         /* determine ix = ilogb(x) */
49         if (hx < 0x00800000) {  /* subnormal x */
50                 for (ix = -126, i=hx<<8; i>0; i<<=1) ix--;
51         } else
52                 ix = (hx>>23) - 127;
53
54         /* determine iy = ilogb(y) */
55         if (hy < 0x00800000) {  /* subnormal y */
56                 for (iy = -126, i=hy<<8; i>0; i<<=1) iy--;
57         } else
58                 iy = (hy>>23) - 127;
59
60         /* set up {hx,lx}, {hy,ly} and align y to x */
61         if (ix >= -126)
62                 hx = 0x00800000|(0x007fffff&hx);
63         else {  /* subnormal x, shift x to normal */
64                 n = -126 - ix;
65                 hx <<= n;
66         }
67         if (iy >= -126)
68                 hy = 0x00800000|(0x007fffff&hy);
69         else {  /* subnormal y, shift y to normal */
70                 n = -126 - iy;
71                 hy <<= n;
72         }
73
74         /* fix point fmod */
75         n = ix - iy;
76         q = 0;
77         while (n--) {
78                 hz = hx - hy;
79                 if (hz < 0)
80                         hx = hx << 1;
81                 else {
82                         hx = hz << 1;
83                         q++;
84                 }
85                 q <<= 1;
86         }
87         hz = hx - hy;
88         if (hz >= 0) {
89                 hx = hz;
90                 q++;
91         }
92
93         /* convert back to floating value and restore the sign */
94         if (hx == 0) {                             /* return sign(x)*0 */
95                 q &= 0x7fffffff;
96                 *quo = sxy ? -q : q;
97                 return Zero[(uint32_t)sx>>31];
98         }
99         while (hx < 0x00800000) {  /* normalize x */
100                 hx <<= 1;
101                 iy--;
102         }
103         if (iy >= -126) {          /* normalize output */
104                 hx = (hx-0x00800000)|((iy+127)<<23);
105         } else {                   /* subnormal output */
106                 n = -126 - iy;
107                 hx >>= n;
108         }
109 fixup:
110         SET_FLOAT_WORD(x,hx);
111         y = fabsf(y);
112         if (y < 0x1p-125f) {
113                 if (x + x > y || (x + x == y && (q & 1))) {
114                         q++;
115                         x -= y;
116                 }
117         } else if (x > 0.5f*y || (x == 0.5f*y && (q & 1))) {
118                 q++;
119                 x -= y;
120         }
121         GET_FLOAT_WORD(hx, x);
122         SET_FLOAT_WORD(x, hx ^ sx);
123         q &= 0x7fffffff;
124         *quo = sxy ? -q : q;
125         return x;
126 }