1 /*	$OpenBSD: math_private.h,v 1.17 2014/06/02 19:31:17 kettenis Exp $	*/
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, 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 /*
14  * from: @(#)fdlibm.h 5.1 93/09/24
15  */
16 
17 #ifndef _MATH_PRIVATE_OPENBSD_H_
18 #define _MATH_PRIVATE_OPENBSD_H_
19 
20 #ifdef _IEEE128_FLOAT
21 
22 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
23 
24 typedef union
25 {
26   long double value;
27   struct {
28     u_int32_t mswhi;
29     u_int32_t mswlo;
30     u_int32_t lswhi;
31     u_int32_t lswlo;
32   } parts32;
33   struct {
34     u_int64_t msw;
35     u_int64_t lsw;
36   } parts64;
37 } ieee_quad_shape_type;
38 
39 #endif
40 
41 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
42 
43 typedef union
44 {
45   long double value;
46   struct {
47     u_int32_t lswlo;
48     u_int32_t lswhi;
49     u_int32_t mswlo;
50     u_int32_t mswhi;
51   } parts32;
52   struct {
53     u_int64_t lsw;
54     u_int64_t msw;
55   } parts64;
56 } ieee_quad_shape_type;
57 
58 #endif
59 
60 /* Get two 64 bit ints from a long double.  */
61 
62 #define GET_LDOUBLE_WORDS64(ix0,ix1,d)				\
63 do {								\
64   ieee_quad_shape_type qw_u;					\
65   qw_u.value = (d);						\
66   (ix0) = qw_u.parts64.msw;					\
67   (ix1) = qw_u.parts64.lsw;					\
68 } while (0)
69 
70 /* Set a long double from two 64 bit ints.  */
71 
72 #define SET_LDOUBLE_WORDS64(d,ix0,ix1)				\
73 do {								\
74   ieee_quad_shape_type qw_u;					\
75   qw_u.parts64.msw = (ix0);					\
76   qw_u.parts64.lsw = (ix1);					\
77   (d) = qw_u.value;						\
78 } while (0)
79 
80 /* Get the more significant 64 bits of a long double mantissa.  */
81 
82 #define GET_LDOUBLE_MSW64(v,d)					\
83 do {								\
84   ieee_quad_shape_type sh_u;					\
85   sh_u.value = (d);						\
86   (v) = sh_u.parts64.msw;					\
87 } while (0)
88 
89 /* Set the more significant 64 bits of a long double mantissa from an int.  */
90 
91 #define SET_LDOUBLE_MSW64(d,v)					\
92 do {								\
93   ieee_quad_shape_type sh_u;					\
94   sh_u.value = (d);						\
95   sh_u.parts64.msw = (v);					\
96   (d) = sh_u.value;						\
97 } while (0)
98 
99 /* Get the least significant 64 bits of a long double mantissa.  */
100 
101 #define GET_LDOUBLE_LSW64(v,d)					\
102 do {								\
103   ieee_quad_shape_type sh_u;					\
104   sh_u.value = (d);						\
105   (v) = sh_u.parts64.lsw;					\
106 } while (0)
107 
108 #define	LDBL_NBIT	0
109 #define LDBL_NBIT_INF   0
110 #define	LDBL_IMPLICIT_NBIT
111 #define	mask_nbit_l(u)	((void)0)
112 #define LDBL_INF_NAN_EXP    32767
113 #define LDBL_EXP_MASK           0x7fff
114 #define LDBL_EXP_SIGN           0x8000
115 
116 #define	LDBL_MANH_SIZE	48
117 #define	LDBL_MANL_SIZE	64
118 
119 static ALWAYS_INLINE int
isnanl_inline(long double x)120 isnanl_inline(long double x)
121 {
122     u_int64_t high, low;
123     GET_LDOUBLE_MSW64(high, x);
124     /* Check exponent for all ones */
125     if (((high >> LDBL_MANH_SIZE) & LDBL_EXP_MASK) != LDBL_EXP_MASK)
126         return 0;
127     GET_LDOUBLE_LSW64(low, x);
128     /* Check for non-zero significand */
129     if (high << (64 - LDBL_MANH_SIZE) == 0 && low == 0)
130         return 0;
131     return 1;
132 }
133 
134 static ALWAYS_INLINE int
issignalingl_inline(long double x)135 issignalingl_inline(long double x)
136 {
137     if (!isnanl_inline(x))
138         return 0;
139     u_int64_t high;
140     GET_LDOUBLE_MSW64(high, x);
141     if (!_IEEE_754_2008_SNAN)
142         return (high & 0x0000800000000000ULL) != 0;
143     else
144         return (high & 0x0000800000000000ULL) == 0;
145 }
146 
147 static ALWAYS_INLINE int
__signbitl(long double x)148 __signbitl(long double x)
149 {
150     int64_t high;
151     GET_LDOUBLE_MSW64(high, x);
152     return high < 0;
153 }
154 
155 /* 128-bit long double */
156 
157 union IEEEl2bits {
158 	long double	e;
159 	struct {
160 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
161 		uint64_t	sign	:1;
162 		uint64_t	exp	:15;
163 		uint64_t	manh	:48;
164 		uint64_t	manl	:64;
165 #endif
166 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
167 		uint64_t	manl	:64;
168 		uint64_t	manh	:48;
169 		uint64_t	exp	:15;
170 		uint64_t	sign	:1;
171 #endif
172 	} bits;
173 	/* TODO andrew: Check the packing here */
174 	struct {
175 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
176 		uint64_t	expsign	:16;
177 		uint64_t	manh	:48;
178 		uint64_t	manl	:64;
179 #endif
180 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
181 		uint64_t	manl	:64;
182 		uint64_t	manh	:48;
183 		uint64_t	expsign	:16;
184 #endif
185 	} xbits;
186 };
187 
188 #define	LDBL_TO_ARRAY32(u, a) do {			\
189 	(a)[0] = (uint32_t)(u).bits.manl;		\
190 	(a)[1] = (uint32_t)((u).bits.manl >> 32);	\
191 	(a)[2] = (uint32_t)(u).bits.manh;		\
192 	(a)[3] = (uint32_t)((u).bits.manh >> 32);	\
193 } while(0)
194 
195 #endif /* _IEEE128_FLOAT */
196 
197 #ifdef _INTEL80_FLOAT
198 
199 /* A union which permits us to convert between a long double and
200    three 32 bit ints.  */
201 
202 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
203 
204 typedef union
205 {
206   long double value;
207   struct {
208 #ifdef __LP64__
209     int padh:32;
210 #endif
211     int exp:16;
212     int padl:16;
213     u_int32_t msw;
214     u_int32_t lsw;
215   } parts;
216 } ieee_extended_shape_type;
217 
218 #endif
219 
220 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
221 
222 typedef union
223 {
224   long double value;
225   struct {
226     u_int32_t lsw;
227     u_int32_t msw;
228     int exp:16;
229     int padl:16;
230 #ifdef __LP64__
231     int padh:32;
232 #endif
233   } parts;
234 } ieee_extended_shape_type;
235 
236 #endif
237 
238 /* Get three 32 bit ints from a double.  */
239 
240 #define GET_LDOUBLE_WORDS(se,ix0,ix1,d)				\
241 do {								\
242   ieee_extended_shape_type ew_u;				\
243   ew_u.value = (d);						\
244   (se) = ew_u.parts.exp;					\
245   (ix0) = ew_u.parts.msw;					\
246   (ix1) = ew_u.parts.lsw;					\
247 } while (0)
248 
249 /* Set a double from two 32 bit ints.  */
250 
251 #define SET_LDOUBLE_WORDS(d,se,ix0,ix1)				\
252 do {								\
253   ieee_extended_shape_type iw_u;				\
254   iw_u.parts.exp = (se);					\
255   iw_u.parts.msw = (ix0);					\
256   iw_u.parts.lsw = (ix1);					\
257   (d) = iw_u.value;						\
258 } while (0)
259 
260 /* Get the more significant 32 bits of a long double mantissa.  */
261 
262 #define GET_LDOUBLE_MSW(v,d)					\
263 do {								\
264   ieee_extended_shape_type sh_u;				\
265   sh_u.value = (d);						\
266   (v) = sh_u.parts.msw;						\
267 } while (0)
268 
269 /* Set the more significant 32 bits of a long double mantissa from an int.  */
270 
271 #define SET_LDOUBLE_MSW(d,v)					\
272 do {								\
273   ieee_extended_shape_type sh_u;				\
274   sh_u.value = (d);						\
275   sh_u.parts.msw = (v);						\
276   (d) = sh_u.value;						\
277 } while (0)
278 
279 /* Get int from the exponent of a long double.  */
280 
281 #define GET_LDOUBLE_EXP(se,d)					\
282 do {								\
283   ieee_extended_shape_type ge_u;				\
284   ge_u.value = (d);						\
285   (se) = ge_u.parts.exp;					\
286 } while (0)
287 
288 /* Set exponent of a long double from an int.  */
289 
290 #define SET_LDOUBLE_EXP(d,se)					\
291 do {								\
292   ieee_extended_shape_type se_u;				\
293   se_u.value = (d);						\
294   se_u.parts.exp = (se);					\
295   (d) = se_u.value;						\
296 } while (0)
297 
298 #define	LDBL_NBIT	0x80000000
299 #ifdef __m68k__
300 #define LDBL_NBIT_INF   0
301 #else
302 #define LDBL_NBIT_INF   LDBL_NBIT
303 #endif
304 #define	mask_nbit_l(u)	((u).bits.manh &= ~LDBL_NBIT)
305 #define LDBL_INF_NAN_EXP    32767
306 #define LDBL_EXP_MASK           0x7fff
307 #define LDBL_EXP_SIGN           0x8000
308 
309 #define	LDBL_MANH_SIZE	32
310 #define	LDBL_MANL_SIZE	32
311 
312 #define	LDBL_TO_ARRAY32(u, a) do {			\
313 	(a)[0] = (uint32_t)(u).bits.manl;		\
314 	(a)[1] = (uint32_t)(u).bits.manh;		\
315 } while (0)
316 
317 static ALWAYS_INLINE int
isnanl_inline(long double x)318 isnanl_inline(long double x)
319 {
320     u_int32_t exp, msw, lsw;
321     GET_LDOUBLE_WORDS(exp, msw, lsw, x);
322     /* Check exponent for all ones */
323     if ((exp & 0x7fff) != 0x7fff)
324         return 0;
325     /* Check for non-zero significand */
326     if (msw == LDBL_NBIT_INF && lsw == 0)
327         return 0;
328     return 1;
329 }
330 
331 static ALWAYS_INLINE int
issignalingl_inline(long double x)332 issignalingl_inline(long double x)
333 {
334     u_int32_t msw;
335     if (!isnanl_inline(x))
336         return 0;
337     GET_LDOUBLE_MSW(msw, x);
338     if (!_IEEE_754_2008_SNAN)
339         return (msw & 0x40000000U) != 0;
340     else
341         return (msw & 0x40000000U) == 0;
342 }
343 
344 static ALWAYS_INLINE int
__signbitl(long double x)345 __signbitl(long double x)
346 {
347     int exp;
348     GET_LDOUBLE_EXP(exp, x);
349     return (exp & 0x8000) != 0;
350 }
351 
352 union IEEEl2bits {
353 	long double	e;
354 	struct {
355 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
356 		uint64_t	sign	:1;
357 		uint64_t	exp	:15;
358 		uint64_t	junk	:16;
359 		uint64_t	manh	:32;
360 		uint64_t	manl	:32;
361 #endif
362 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
363 		uint64_t	manl	:32;
364 		uint64_t	manh	:32;
365 		uint64_t	exp	:15;
366 		uint64_t	sign	:1;
367 		uint64_t	junk	:16;
368 #endif
369 	} bits;
370 	struct {
371 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
372 		uint64_t 	expsign	:16;
373 		uint64_t	junk	:16;
374 		uint64_t        man	:64;
375 #endif
376 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
377 		uint64_t        man	:64;
378 		uint64_t 	expsign	:16;
379 		uint64_t	junk	:16;
380 #endif
381 	} xbits;
382 };
383 
384 #endif /* _INTEL80_FLOAT */
385 
386 #ifdef _DOUBLE_DOUBLE_FLOAT
387 
388 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
389 
390 typedef union
391 {
392   long double value;
393   struct {
394     double      msd;
395     double      lsd;
396   } doubles;
397   struct {
398     u_int32_t mswhi;
399     u_int32_t mswlo;
400     u_int32_t lswhi;
401     u_int32_t lswlo;
402   } parts32;
403   struct {
404     u_int64_t msw;
405     u_int64_t lsw;
406   } parts64;
407   struct {
408     int exp:12;
409     u_int64_t   msw :52;
410     int expl:12;
411     u_int64_t   lsw :52;
412   } parts;
413 } double_double_shape_type;
414 
415 #endif /* __ORDER_BIG_ENDIAN__ */
416 
417 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
418 
419 typedef union
420 {
421   long double value;
422   struct {
423     double      msd;
424     double      lsd;
425   } doubles;
426   struct {
427     u_int32_t mswlo;
428     u_int32_t mswhi;
429     u_int32_t lswlo;
430     u_int32_t lswhi;
431   } parts32;
432   struct {
433     u_int64_t msw;
434     u_int64_t lsw;
435   } parts64;
436   struct {
437     u_int64_t   msw :52;
438     int exp:12;
439     u_int64_t   lsw :52;
440     int expl:12;
441   } parts;
442 } double_double_shape_type;
443 
444 #endif /* __ORDER_LITTLE_ENDIAN__ */
445 
446 /* Get two 64 bit ints from a long double.  */
447 
448 #define GET_LDOUBLE_WORDS64(ix0,ix1,d)				\
449 do {								\
450   double_double_shape_type qw_u;				\
451   qw_u.value = (d);						\
452   (ix0) = qw_u.parts64.msw;					\
453   (ix1) = qw_u.parts64.lsw;					\
454 } while (0)
455 
456 /* Set a long double from two 64 bit ints.  */
457 
458 #define SET_LDOUBLE_WORDS64(d,ix0,ix1)				\
459 do {								\
460   double_double_shape_type qw_u;				\
461   qw_u.parts64.msw = (ix0);					\
462   qw_u.parts64.lsw = (ix1);					\
463   (d) = qw_u.value;						\
464 } while (0)
465 
466 /* Get the more significant 64 bits of a long double mantissa.  */
467 
468 #define GET_LDOUBLE_MSW64(v,d)					\
469 do {								\
470   double_double_shape_type sh_u;				\
471   sh_u.value = (d);						\
472   (v) = sh_u.parts64.msw;					\
473 } while (0)
474 
475 /* Set the more significant 64 bits of a long double mantissa from an int.  */
476 
477 #define SET_LDOUBLE_MSW64(d,v)					\
478 do {								\
479   double_double_shape_type sh_u;				\
480   sh_u.value = (d);						\
481   sh_u.parts64.msw = (v);					\
482   (d) = sh_u.value;						\
483 } while (0)
484 
485 /* Get the least significant 64 bits of a long double mantissa.  */
486 
487 #define GET_LDOUBLE_LSW64(v,d)					\
488 do {								\
489   double_double_shape_type sh_u;				\
490   sh_u.value = (d);						\
491   (v) = sh_u.parts64.lsw;					\
492 } while (0)
493 
494 /* Get int from the exponent of a long double.  */
495 
496 #define GET_LDOUBLE_EXP(se,d)					\
497 do {								\
498   double_double_shape_type ge_u;				\
499   ge_u.value = (d);						\
500   (se) = ge_u.parts.exp;					\
501 } while (0)
502 
503 /* Set exponent of a long double from an int.  */
504 
505 #define SET_LDOUBLE_EXP(d,se)					\
506 do {								\
507   double_double_shape_type se_u;				\
508   se_u.value = (d);						\
509   se_u.parts.exp = (se);					\
510   (d) = se_u.value;						\
511 } while (0)
512 
513 static ALWAYS_INLINE int
isnanl_inline(long double x)514 isnanl_inline(long double x)
515 {
516     uint64_t high;
517     GET_LDOUBLE_MSW64(high, x);
518     /* Check exponent for all ones */
519     if (((high >> 52) & 0x7ff) != 0x7ff)
520         return 0;
521     /* Check for non-zero significand */
522     if (high << 12 == 0)
523         return 0;
524     return 1;
525 }
526 
527 static ALWAYS_INLINE int
issignalingl_inline(long double x)528 issignalingl_inline(long double x)
529 {
530     uint64_t high;
531     if (!isnanl_inline(x))
532         return 0;
533     GET_LDOUBLE_MSW64(high, x);
534     if (!_IEEE_754_2008_SNAN)
535         return (high & 0x0008000000000000ULL) != 0;
536     else
537         return (high & 0x0008000000000000ULL) == 0;
538 }
539 
540 static ALWAYS_INLINE int
__signbitl(long double x)541 __signbitl(long double x)
542 {
543     int exp;
544     GET_LDOUBLE_EXP(exp, x);
545     return exp < 0;
546 }
547 
548 #ifdef __PPC__
549 #ifdef __GNUC__
550 #pragma GCC diagnostic ignored "-Wpragmas"
551 #pragma GCC diagnostic ignored "-Wunknown-warning-option"
552 /* the bitfields confuse the ppc compiler into thinking accessing
553  * manl will be 'out of bounds'
554  */
555 #pragma GCC diagnostic ignored "-Wanalyzer-out-of-bounds"
556 #define yup_all_good
557 #endif
558 #endif
559 
560 union IEEEl2bits {
561 	long double	e;
562 	struct {
563 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
564 		u_int64_t	sign	:1;
565 		u_int64_t	exp	:11;
566 		u_int64_t       manh	:52;
567                 u_int64_t       signl   :1;
568                 u_int64_t       expl    :11;
569                 u_int64_t       manl    :52;
570 #endif
571 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
572 		u_int64_t	manh	:52;
573 		u_int64_t	exp	:11;
574 		u_int64_t	sign	:1;
575                 u_int64_t       manl    :52;
576                 u_int64_t       expl    :11;
577                 u_int64_t       signl   :1;
578 #endif
579 	} bits;
580 	struct {
581 #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
582 		u_int64_t 	expsign	:12;
583 		u_int64_t       manh	:52;
584                 u_int64_t       expsignl:12;
585                 u_int64_t       manl    :52;
586 #endif
587 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
588 		u_int64_t       manh	:52;
589 		u_int64_t 	expsign	:12;
590                 u_int64_t       manl    :52;
591                 u_int64_t       expsignl:12;
592 #endif
593 	} xbits;
594         struct {
595                 double  dh;
596                 double  dl;
597         } dbits;
598 };
599 
600 #define	LDBL_NBIT	        0
601 #define LDBL_NBIT_INF           0
602 #define	LDBL_IMPLICIT_NBIT
603 #define	mask_nbit_l(u)	        ((void)0)
604 #define LDBL_INF_NAN_EXP        2047
605 #define LDBL_EXP_MASK           0x7ff
606 #define LDBL_EXP_SIGN           0x800
607 
608 #define	LDBL_MANH_SIZE	52
609 
610 #define	LDBL_TO_ARRAY32(u, a) do {			\
611 	(a)[0] = (uint32_t)(u).bits.manl;		\
612 	(a)[1] = (uint32_t)((u).bits.manl >> 32);	\
613 	(a)[2] = (uint32_t)(u).bits.manh;		\
614 	(a)[3] = (uint32_t)((u).bits.manh >> 32);	\
615 } while(0)
616 
617 #endif /* _DOUBLE_DOUBLE_FLOAT */
618 
619 /*
620  * Common routine to process the arguments to nan(), nanf(), and nanl().
621  */
622 void __scan_nan(uint32_t *__words, int __num_words, const char *__s);
623 
624 /*
625  * Functions internal to the math package, yet not static.
626  */
627 double __exp__D(double, double);
628 struct Double __log__D(double);
629 long double __p1evll(long double, const long double *, int);
630 long double __polevll(long double, const long double *, int);
631 
632 #endif /* _MATH_PRIVATE_OPENBSD_H_ */
633