1 /*-
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunSoft, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 #include <float.h>
13 #include <stdint.h>
14
15 #include "fpmath.h"
16 #include "math.h"
17 #include "math_private.h"
18
19 #define BIAS (LDBL_MAX_EXP - 1)
20
21 #if LDBL_MANL_SIZE > 32
22 typedef uint64_t manl_t;
23 #else
24 typedef uint32_t manl_t;
25 #endif
26
27 #if LDBL_MANH_SIZE > 32
28 typedef uint64_t manh_t;
29 #else
30 typedef uint32_t manh_t;
31 #endif
32
33 /*
34 * These macros add and remove an explicit integer bit in front of the
35 * fractional mantissa, if the architecture doesn't have such a bit by
36 * default already.
37 */
38 #ifdef LDBL_IMPLICIT_NBIT
39 #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
40 #define HFRAC_BITS LDBL_MANH_SIZE
41 #else
42 #define SET_NBIT(hx) (hx)
43 #define HFRAC_BITS (LDBL_MANH_SIZE - 1)
44 #endif
45
46 #define MANL_SHIFT (LDBL_MANL_SIZE - 1)
47
48 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
49
50 /*
51 * fmodl(x,y)
52 * Return x mod y in exact arithmetic
53 * Method: shift and subtract
54 *
55 * Assumptions:
56 * - The low part of the mantissa fits in a manl_t exactly.
57 * - The high part of the mantissa fits in an int64_t with enough room
58 * for an explicit integer bit in front of the fractional bits.
59 */
60 long double
fmodl(long double x,long double y)61 fmodl(long double x, long double y)
62 {
63 union IEEEl2bits ux, uy;
64 int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
65 manh_t hy;
66 manl_t lx,ly,lz;
67 int ix,iy,n,sx;
68
69 ux.e = x;
70 uy.e = y;
71 sx = ux.bits.sign;
72
73 /* purge off exception values */
74 if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */
75 (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
76 (uy.bits.exp == BIAS + LDBL_MAX_EXP &&
77 ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
78 return nan_mix_op(x, y, *)/nan_mix_op(x, y, *);
79 if(ux.bits.exp<=uy.bits.exp) {
80 if((ux.bits.exp<uy.bits.exp) ||
81 (ux.bits.manh<=uy.bits.manh &&
82 (ux.bits.manh<uy.bits.manh ||
83 ux.bits.manl<uy.bits.manl))) {
84 return x; /* |x|<|y| return x or x-y */
85 }
86 if(ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl) {
87 return Zero[sx]; /* |x|=|y| return x*0*/
88 }
89 }
90
91 /* determine ix = ilogb(x) */
92 if(ux.bits.exp == 0) { /* subnormal x */
93 ux.e *= 0x1.0p512;
94 ix = ux.bits.exp - (BIAS + 512);
95 } else {
96 ix = ux.bits.exp - BIAS;
97 }
98
99 /* determine iy = ilogb(y) */
100 if(uy.bits.exp == 0) { /* subnormal y */
101 uy.e *= 0x1.0p512;
102 iy = uy.bits.exp - (BIAS + 512);
103 } else {
104 iy = uy.bits.exp - BIAS;
105 }
106
107 /* set up {hx,lx}, {hy,ly} and align y to x */
108 hx = SET_NBIT(ux.bits.manh);
109 hy = SET_NBIT(uy.bits.manh);
110 lx = ux.bits.manl;
111 ly = uy.bits.manl;
112
113 /* fix point fmod */
114 n = ix - iy;
115
116 while(n--) {
117 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
118 if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
119 else {
120 if ((hz|lz)==0) /* return sign(x)*0 */
121 return Zero[sx];
122 hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz;
123 }
124 }
125 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
126 if(hz>=0) {hx=hz;lx=lz;}
127
128 /* convert back to floating value and restore the sign */
129 if((hx|lx)==0) /* return sign(x)*0 */
130 return Zero[sx];
131 while(hx<(1ULL<<HFRAC_BITS)) { /* normalize x */
132 hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
133 iy -= 1;
134 }
135 ux.bits.manh = hx; /* The mantissa is truncated here if needed. */
136 ux.bits.manl = lx;
137 if (iy < LDBL_MIN_EXP) {
138 ux.bits.exp = iy + (BIAS + 512);
139 ux.e *= 0x1p-512;
140 } else {
141 ux.bits.exp = iy + BIAS;
142 }
143 x = ux.e * one; /* create necessary signal */
144 return x; /* exact output */
145 }
146