1 /* origin: FreeBSD /usr/src/lib/msun/src/s_remquol.c */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
10 * ====================================================
15 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
16 long double remquol(long double x
, long double y
, int *quo
)
18 return remquo(x
, y
, quo
);
20 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
22 #define BIAS (LDBL_MAX_EXP - 1)
24 #if LDBL_MANL_SIZE > 32
25 typedef uint64_t manl_t
;
27 typedef uint32_t manl_t
;
30 #if LDBL_MANH_SIZE > 32
31 typedef uint64_t manh_t
;
33 typedef uint32_t manh_t
;
37 * These macros add and remove an explicit integer bit in front of the
38 * fractional mantissa, if the architecture doesn't have such a bit by
41 #ifdef LDBL_IMPLICIT_NBIT
42 #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
43 #define HFRAC_BITS LDBL_MANH_SIZE
45 #define SET_NBIT(hx) (hx)
46 #define HFRAC_BITS (LDBL_MANH_SIZE - 1)
49 #define MANL_SHIFT (LDBL_MANL_SIZE - 1)
51 static const long double Zero
[] = {0.0, -0.0};
54 * Return the IEEE remainder and set *quo to the last n bits of the
55 * quotient, rounded to the nearest integer. We choose n=31 because
56 * we wind up computing all the integer bits of the quotient anyway as
57 * a side-effect of computing the remainder by the shift and subtract
58 * method. In practice, this is far more bits than are needed to use
59 * remquo in reduction algorithms.
62 * - The low part of the mantissa fits in a manl_t exactly.
63 * - The high part of the mantissa fits in an int64_t with enough room
64 * for an explicit integer bit in front of the fractional bits.
66 long double remquol(long double x
, long double y
, int *quo
)
68 union IEEEl2bits ux
, uy
;
69 int64_t hx
,hz
; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
77 sxy
= sx
^ uy
.bits
.sign
;
78 ux
.bits
.sign
= 0; /* |x| */
79 uy
.bits
.sign
= 0; /* |y| */
82 /* purge off exception values */
83 if ((uy
.bits
.exp
|uy
.bits
.manh
|uy
.bits
.manl
)==0 || /* y=0 */
84 (ux
.bits
.exp
== BIAS
+ LDBL_MAX_EXP
) || /* or x not finite */
85 (uy
.bits
.exp
== BIAS
+ LDBL_MAX_EXP
&&
86 ((uy
.bits
.manh
&~LDBL_NBIT
)|uy
.bits
.manl
)!=0)) /* or y is NaN */
88 if (ux
.bits
.exp
<= uy
.bits
.exp
) {
89 if ((ux
.bits
.exp
< uy
.bits
.exp
) ||
90 (ux
.bits
.manh
<= uy
.bits
.manh
&&
91 (ux
.bits
.manh
< uy
.bits
.manh
||
92 ux
.bits
.manl
< uy
.bits
.manl
))) {
94 goto fixup
; /* |x|<|y| return x or x-y */
96 if (ux
.bits
.manh
== uy
.bits
.manh
&& ux
.bits
.manl
== uy
.bits
.manl
) {
98 return Zero
[sx
]; /* |x|=|y| return x*0*/
102 /* determine ix = ilogb(x) */
103 if (ux
.bits
.exp
== 0) { /* subnormal x */
105 ix
= ux
.bits
.exp
- (BIAS
+ 512);
107 ix
= ux
.bits
.exp
- BIAS
;
110 /* determine iy = ilogb(y) */
111 if (uy
.bits
.exp
== 0) { /* subnormal y */
113 iy
= uy
.bits
.exp
- (BIAS
+ 512);
115 iy
= uy
.bits
.exp
- BIAS
;
118 /* set up {hx,lx}, {hy,ly} and align y to x */
119 hx
= SET_NBIT(ux
.bits
.manh
);
120 hy
= SET_NBIT(uy
.bits
.manh
);
134 hx
= hx
+ hx
+ (lx
>>MANL_SHIFT
);
137 hx
= hz
+ hz
+ (lz
>>MANL_SHIFT
);
153 /* convert back to floating value and restore the sign */
154 if ((hx
|lx
) == 0) { /* return sign(x)*0 */
159 while (hx
< (1ULL<<HFRAC_BITS
)) { /* normalize x */
160 hx
= hx
+ hx
+ (lx
>>MANL_SHIFT
);
164 ux
.bits
.manh
= hx
; /* The integer bit is truncated here if needed. */
166 if (iy
< LDBL_MIN_EXP
) {
167 ux
.bits
.exp
= iy
+ (BIAS
+ 512);
170 ux
.bits
.exp
= iy
+ BIAS
;
176 if (y
< LDBL_MIN
* 2) {
177 if (x
+ x
> y
|| (x
+ x
== y
&& (q
& 1))) {
181 } else if (x
> 0.5*y
|| (x
== 0.5*y
&& (q
& 1))) {