1 /* Software floating-point emulation.
2 Basic one-word fraction declaration and manipulation.
3 Copyright (C) 1997-2023 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
5
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 In addition to the permissions in the GNU Lesser General Public
12 License, the Free Software Foundation gives you unlimited
13 permission to link the compiled version of this file into
14 combinations with other programs, and to distribute those
15 combinations without any restriction coming from the use of this
16 file. (The Lesser General Public License restrictions do apply in
17 other respects; for example, they cover modification of the file,
18 and distribution when not linked into a combine executable.)
19
20 The GNU C Library is distributed in the hope that it will be useful,
21 but WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23 Lesser General Public License for more details.
24
25 You should have received a copy of the GNU Lesser General Public
26 License along with the GNU C Library; if not, see
27 <https://www.gnu.org/licenses/>. */
28
29 #ifndef SOFT_FP_OP_1_H
30 #define SOFT_FP_OP_1_H 1
31
32 #define _FP_FRAC_DECL_1(X) _FP_W_TYPE X##_f _FP_ZERO_INIT
33 #define _FP_FRAC_COPY_1(D, S) (D##_f = S##_f)
34 #define _FP_FRAC_SET_1(X, I) (X##_f = I)
35 #define _FP_FRAC_HIGH_1(X) (X##_f)
36 #define _FP_FRAC_LOW_1(X) (X##_f)
37 #define _FP_FRAC_WORD_1(X, w) (X##_f)
38
39 #define _FP_FRAC_ADDI_1(X, I) (X##_f += I)
40 #define _FP_FRAC_SLL_1(X, N) \
41 do \
42 { \
43 if (__builtin_constant_p (N) && (N) == 1) \
44 X##_f += X##_f; \
45 else \
46 X##_f <<= (N); \
47 } \
48 while (0)
49 #define _FP_FRAC_SRL_1(X, N) (X##_f >>= N)
50
51 /* Right shift with sticky-lsb. */
52 #define _FP_FRAC_SRST_1(X, S, N, sz) __FP_FRAC_SRST_1 (X##_f, S, (N), (sz))
53 #define _FP_FRAC_SRS_1(X, N, sz) __FP_FRAC_SRS_1 (X##_f, (N), (sz))
54
55 #define __FP_FRAC_SRST_1(X, S, N, sz) \
56 do \
57 { \
58 S = (__builtin_constant_p (N) && (N) == 1 \
59 ? X & 1 \
60 : (X << (_FP_W_TYPE_SIZE - (N))) != 0); \
61 X = X >> (N); \
62 } \
63 while (0)
64
65 #define __FP_FRAC_SRS_1(X, N, sz) \
66 (X = (X >> (N) | (__builtin_constant_p (N) && (N) == 1 \
67 ? X & 1 \
68 : (X << (_FP_W_TYPE_SIZE - (N))) != 0)))
69
70 #define _FP_FRAC_ADD_1(R, X, Y) (R##_f = X##_f + Y##_f)
71 #define _FP_FRAC_SUB_1(R, X, Y) (R##_f = X##_f - Y##_f)
72 #define _FP_FRAC_DEC_1(X, Y) (X##_f -= Y##_f)
73 #define _FP_FRAC_CLZ_1(z, X) __FP_CLZ ((z), X##_f)
74
75 /* Predicates. */
76 #define _FP_FRAC_NEGP_1(X) ((_FP_WS_TYPE) X##_f < 0)
77 #define _FP_FRAC_ZEROP_1(X) (X##_f == 0)
78 #define _FP_FRAC_OVERP_1(fs, X) (X##_f & _FP_OVERFLOW_##fs)
79 #define _FP_FRAC_CLEAR_OVERP_1(fs, X) (X##_f &= ~_FP_OVERFLOW_##fs)
80 #define _FP_FRAC_HIGHBIT_DW_1(fs, X) (X##_f & _FP_HIGHBIT_DW_##fs)
81 #define _FP_FRAC_EQ_1(X, Y) (X##_f == Y##_f)
82 #define _FP_FRAC_GE_1(X, Y) (X##_f >= Y##_f)
83 #define _FP_FRAC_GT_1(X, Y) (X##_f > Y##_f)
84
85 #define _FP_ZEROFRAC_1 0
86 #define _FP_MINFRAC_1 1
87 #define _FP_MAXFRAC_1 (~(_FP_WS_TYPE) 0)
88
89 /* Unpack the raw bits of a native fp value. Do not classify or
90 normalize the data. */
91
92 #define _FP_UNPACK_RAW_1(fs, X, val) \
93 do \
94 { \
95 union _FP_UNION_##fs _FP_UNPACK_RAW_1_flo; \
96 _FP_UNPACK_RAW_1_flo.flt = (val); \
97 \
98 X##_f = _FP_UNPACK_RAW_1_flo.bits.frac; \
99 X##_e = _FP_UNPACK_RAW_1_flo.bits.exp; \
100 X##_s = _FP_UNPACK_RAW_1_flo.bits.sign; \
101 } \
102 while (0)
103
104 #define _FP_UNPACK_RAW_1_P(fs, X, val) \
105 do \
106 { \
107 union _FP_UNION_##fs *_FP_UNPACK_RAW_1_P_flo \
108 = (union _FP_UNION_##fs *) (val); \
109 \
110 X##_f = _FP_UNPACK_RAW_1_P_flo->bits.frac; \
111 X##_e = _FP_UNPACK_RAW_1_P_flo->bits.exp; \
112 X##_s = _FP_UNPACK_RAW_1_P_flo->bits.sign; \
113 } \
114 while (0)
115
116 /* Repack the raw bits of a native fp value. */
117
118 #define _FP_PACK_RAW_1(fs, val, X) \
119 do \
120 { \
121 union _FP_UNION_##fs _FP_PACK_RAW_1_flo; \
122 \
123 _FP_PACK_RAW_1_flo.bits.frac = X##_f; \
124 _FP_PACK_RAW_1_flo.bits.exp = X##_e; \
125 _FP_PACK_RAW_1_flo.bits.sign = X##_s; \
126 \
127 (val) = _FP_PACK_RAW_1_flo.flt; \
128 } \
129 while (0)
130
131 #define _FP_PACK_RAW_1_P(fs, val, X) \
132 do \
133 { \
134 union _FP_UNION_##fs *_FP_PACK_RAW_1_P_flo \
135 = (union _FP_UNION_##fs *) (val); \
136 \
137 _FP_PACK_RAW_1_P_flo->bits.frac = X##_f; \
138 _FP_PACK_RAW_1_P_flo->bits.exp = X##_e; \
139 _FP_PACK_RAW_1_P_flo->bits.sign = X##_s; \
140 } \
141 while (0)
142
143
144 /* Multiplication algorithms: */
145
146 /* Basic. Assuming the host word size is >= 2*FRACBITS, we can do the
147 multiplication immediately. */
148
149 #define _FP_MUL_MEAT_DW_1_imm(wfracbits, R, X, Y) \
150 do \
151 { \
152 R##_f = X##_f * Y##_f; \
153 } \
154 while (0)
155
156 #define _FP_MUL_MEAT_1_imm(wfracbits, R, X, Y) \
157 do \
158 { \
159 _FP_MUL_MEAT_DW_1_imm ((wfracbits), R, X, Y); \
160 /* Normalize since we know where the msb of the multiplicands \
161 were (bit B), we know that the msb of the of the product is \
162 at either 2B or 2B-1. */ \
163 _FP_FRAC_SRS_1 (R, (wfracbits)-1, 2*(wfracbits)); \
164 } \
165 while (0)
166
167 /* Given a 1W * 1W => 2W primitive, do the extended multiplication. */
168
169 #define _FP_MUL_MEAT_DW_1_wide(wfracbits, R, X, Y, doit) \
170 do \
171 { \
172 doit (R##_f1, R##_f0, X##_f, Y##_f); \
173 } \
174 while (0)
175
176 #define _FP_MUL_MEAT_1_wide(wfracbits, R, X, Y, doit) \
177 do \
178 { \
179 _FP_FRAC_DECL_2 (_FP_MUL_MEAT_1_wide_Z); \
180 _FP_MUL_MEAT_DW_1_wide ((wfracbits), _FP_MUL_MEAT_1_wide_Z, \
181 X, Y, doit); \
182 /* Normalize since we know where the msb of the multiplicands \
183 were (bit B), we know that the msb of the of the product is \
184 at either 2B or 2B-1. */ \
185 _FP_FRAC_SRS_2 (_FP_MUL_MEAT_1_wide_Z, (wfracbits)-1, \
186 2*(wfracbits)); \
187 R##_f = _FP_MUL_MEAT_1_wide_Z_f0; \
188 } \
189 while (0)
190
191 /* Finally, a simple widening multiply algorithm. What fun! */
192
193 #define _FP_MUL_MEAT_DW_1_hard(wfracbits, R, X, Y) \
194 do \
195 { \
196 _FP_W_TYPE _FP_MUL_MEAT_DW_1_hard_xh, _FP_MUL_MEAT_DW_1_hard_xl; \
197 _FP_W_TYPE _FP_MUL_MEAT_DW_1_hard_yh, _FP_MUL_MEAT_DW_1_hard_yl; \
198 _FP_FRAC_DECL_2 (_FP_MUL_MEAT_DW_1_hard_a); \
199 \
200 /* Split the words in half. */ \
201 _FP_MUL_MEAT_DW_1_hard_xh = X##_f >> (_FP_W_TYPE_SIZE/2); \
202 _FP_MUL_MEAT_DW_1_hard_xl \
203 = X##_f & (((_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE/2)) - 1); \
204 _FP_MUL_MEAT_DW_1_hard_yh = Y##_f >> (_FP_W_TYPE_SIZE/2); \
205 _FP_MUL_MEAT_DW_1_hard_yl \
206 = Y##_f & (((_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE/2)) - 1); \
207 \
208 /* Multiply the pieces. */ \
209 R##_f0 = _FP_MUL_MEAT_DW_1_hard_xl * _FP_MUL_MEAT_DW_1_hard_yl; \
210 _FP_MUL_MEAT_DW_1_hard_a_f0 \
211 = _FP_MUL_MEAT_DW_1_hard_xh * _FP_MUL_MEAT_DW_1_hard_yl; \
212 _FP_MUL_MEAT_DW_1_hard_a_f1 \
213 = _FP_MUL_MEAT_DW_1_hard_xl * _FP_MUL_MEAT_DW_1_hard_yh; \
214 R##_f1 = _FP_MUL_MEAT_DW_1_hard_xh * _FP_MUL_MEAT_DW_1_hard_yh; \
215 \
216 /* Reassemble into two full words. */ \
217 if ((_FP_MUL_MEAT_DW_1_hard_a_f0 += _FP_MUL_MEAT_DW_1_hard_a_f1) \
218 < _FP_MUL_MEAT_DW_1_hard_a_f1) \
219 R##_f1 += (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE/2); \
220 _FP_MUL_MEAT_DW_1_hard_a_f1 \
221 = _FP_MUL_MEAT_DW_1_hard_a_f0 >> (_FP_W_TYPE_SIZE/2); \
222 _FP_MUL_MEAT_DW_1_hard_a_f0 \
223 = _FP_MUL_MEAT_DW_1_hard_a_f0 << (_FP_W_TYPE_SIZE/2); \
224 _FP_FRAC_ADD_2 (R, R, _FP_MUL_MEAT_DW_1_hard_a); \
225 } \
226 while (0)
227
228 #define _FP_MUL_MEAT_1_hard(wfracbits, R, X, Y) \
229 do \
230 { \
231 _FP_FRAC_DECL_2 (_FP_MUL_MEAT_1_hard_z); \
232 _FP_MUL_MEAT_DW_1_hard ((wfracbits), \
233 _FP_MUL_MEAT_1_hard_z, X, Y); \
234 \
235 /* Normalize. */ \
236 _FP_FRAC_SRS_2 (_FP_MUL_MEAT_1_hard_z, \
237 (wfracbits) - 1, 2*(wfracbits)); \
238 R##_f = _FP_MUL_MEAT_1_hard_z_f0; \
239 } \
240 while (0)
241
242
243 /* Division algorithms: */
244
245 /* Basic. Assuming the host word size is >= 2*FRACBITS, we can do the
246 division immediately. Give this macro either _FP_DIV_HELP_imm for
247 C primitives or _FP_DIV_HELP_ldiv for the ISO function. Which you
248 choose will depend on what the compiler does with divrem4. */
249
250 #define _FP_DIV_MEAT_1_imm(fs, R, X, Y, doit) \
251 do \
252 { \
253 _FP_W_TYPE _FP_DIV_MEAT_1_imm_q, _FP_DIV_MEAT_1_imm_r; \
254 X##_f <<= (X##_f < Y##_f \
255 ? R##_e--, _FP_WFRACBITS_##fs \
256 : _FP_WFRACBITS_##fs - 1); \
257 doit (_FP_DIV_MEAT_1_imm_q, _FP_DIV_MEAT_1_imm_r, X##_f, Y##_f); \
258 R##_f = _FP_DIV_MEAT_1_imm_q | (_FP_DIV_MEAT_1_imm_r != 0); \
259 } \
260 while (0)
261
262 /* GCC's longlong.h defines a 2W / 1W => (1W,1W) primitive udiv_qrnnd
263 that may be useful in this situation. This first is for a primitive
264 that requires normalization, the second for one that does not. Look
265 for UDIV_NEEDS_NORMALIZATION to tell which your machine needs. */
266
267 #define _FP_DIV_MEAT_1_udiv_norm(fs, R, X, Y) \
268 do \
269 { \
270 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_norm_nh; \
271 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_norm_nl; \
272 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_norm_q; \
273 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_norm_r; \
274 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_norm_y; \
275 \
276 /* Normalize Y -- i.e. make the most significant bit set. */ \
277 _FP_DIV_MEAT_1_udiv_norm_y = Y##_f << _FP_WFRACXBITS_##fs; \
278 \
279 /* Shift X op correspondingly high, that is, up one full word. */ \
280 if (X##_f < Y##_f) \
281 { \
282 R##_e--; \
283 _FP_DIV_MEAT_1_udiv_norm_nl = 0; \
284 _FP_DIV_MEAT_1_udiv_norm_nh = X##_f; \
285 } \
286 else \
287 { \
288 _FP_DIV_MEAT_1_udiv_norm_nl = X##_f << (_FP_W_TYPE_SIZE - 1); \
289 _FP_DIV_MEAT_1_udiv_norm_nh = X##_f >> 1; \
290 } \
291 \
292 udiv_qrnnd (_FP_DIV_MEAT_1_udiv_norm_q, \
293 _FP_DIV_MEAT_1_udiv_norm_r, \
294 _FP_DIV_MEAT_1_udiv_norm_nh, \
295 _FP_DIV_MEAT_1_udiv_norm_nl, \
296 _FP_DIV_MEAT_1_udiv_norm_y); \
297 R##_f = (_FP_DIV_MEAT_1_udiv_norm_q \
298 | (_FP_DIV_MEAT_1_udiv_norm_r != 0)); \
299 } \
300 while (0)
301
302 #define _FP_DIV_MEAT_1_udiv(fs, R, X, Y) \
303 do \
304 { \
305 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_nh, _FP_DIV_MEAT_1_udiv_nl; \
306 _FP_W_TYPE _FP_DIV_MEAT_1_udiv_q, _FP_DIV_MEAT_1_udiv_r; \
307 if (X##_f < Y##_f) \
308 { \
309 R##_e--; \
310 _FP_DIV_MEAT_1_udiv_nl = X##_f << _FP_WFRACBITS_##fs; \
311 _FP_DIV_MEAT_1_udiv_nh = X##_f >> _FP_WFRACXBITS_##fs; \
312 } \
313 else \
314 { \
315 _FP_DIV_MEAT_1_udiv_nl = X##_f << (_FP_WFRACBITS_##fs - 1); \
316 _FP_DIV_MEAT_1_udiv_nh = X##_f >> (_FP_WFRACXBITS_##fs + 1); \
317 } \
318 udiv_qrnnd (_FP_DIV_MEAT_1_udiv_q, _FP_DIV_MEAT_1_udiv_r, \
319 _FP_DIV_MEAT_1_udiv_nh, _FP_DIV_MEAT_1_udiv_nl, \
320 Y##_f); \
321 R##_f = _FP_DIV_MEAT_1_udiv_q | (_FP_DIV_MEAT_1_udiv_r != 0); \
322 } \
323 while (0)
324
325
326 /* Square root algorithms:
327 We have just one right now, maybe Newton approximation
328 should be added for those machines where division is fast. */
329
330 #define _FP_SQRT_MEAT_1(R, S, T, X, q) \
331 do \
332 { \
333 while ((q) != _FP_WORK_ROUND) \
334 { \
335 T##_f = S##_f + (q); \
336 if (T##_f <= X##_f) \
337 { \
338 S##_f = T##_f + (q); \
339 X##_f -= T##_f; \
340 R##_f += (q); \
341 } \
342 _FP_FRAC_SLL_1 (X, 1); \
343 (q) >>= 1; \
344 } \
345 if (X##_f) \
346 { \
347 if (S##_f < X##_f) \
348 R##_f |= _FP_WORK_ROUND; \
349 R##_f |= _FP_WORK_STICKY; \
350 } \
351 } \
352 while (0)
353
354 /* Assembly/disassembly for converting to/from integral types.
355 No shifting or overflow handled here. */
356
357 #define _FP_FRAC_ASSEMBLE_1(r, X, rsize) ((r) = X##_f)
358 #define _FP_FRAC_DISASSEMBLE_1(X, r, rsize) (X##_f = (r))
359
360
361 /* Convert FP values between word sizes. */
362
363 #define _FP_FRAC_COPY_1_1(D, S) (D##_f = S##_f)
364
365 #endif /* !SOFT_FP_OP_1_H */