1 /* Software floating-point emulation.
2    Basic one-word fraction declaration and manipulation.
3    Copyright (C) 1997-2022 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 */
366