1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, 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/*
13 * from: @(#)fdlibm.h 5.1 93/09/24
14 * $FreeBSD: src/lib/msun/src/math_private.h,v 1.20 2005/11/28 04:58:57 bde Exp $
15 */
16
17#ifndef _MATH_PRIVATE_H_
18#define	_MATH_PRIVATE_H_
19
20#include <sys/types.h>
21#include <endian.h>
22
23/*
24 * The original fdlibm code used statements like:
25 *	n0 = ((*(int*)&one)>>29)^1;		* index of high word *
26 *	ix0 = *(n0+(int*)&x);			* high word of x *
27 *	ix1 = *((1-n0)+(int*)&x);		* low word of x *
28 * to dig two 32 bit words out of the 64 bit IEEE floating point
29 * value.  That is non-ANSI, and, moreover, the gcc instruction
30 * scheduler gets it wrong.  We instead use the following macros.
31 * Unlike the original code, we determine the endianness at compile
32 * time, not at run time; I don't see much benefit to selecting
33 * endianness at run time.
34 */
35
36/*
37 * A union which permits us to convert between a double and two 32 bit
38 * ints.
39 */
40
41#if (__BYTE_ORDER == __BIG_ENDIAN) || (defined(__arm__) && !defined(__VFP_FP__))
42
43typedef union
44{
45  double value;
46  struct
47  {
48    u_int32_t msw;
49    u_int32_t lsw;
50  } parts;
51} ieee_double_shape_type;
52
53#endif
54
55#if __BYTE_ORDER == __LITTLE_ENDIAN  && !(defined(__arm__) && !defined(__VFP_FP__))
56
57typedef union
58{
59  double value;
60  struct
61  {
62    u_int32_t lsw;
63    u_int32_t msw;
64  } parts;
65} ieee_double_shape_type;
66
67#endif
68
69/* Get two 32 bit ints from a double.  */
70
71#define EXTRACT_WORDS(ix0,ix1,d)				\
72do {								\
73  ieee_double_shape_type ew_u;					\
74  ew_u.value = (d);						\
75  (ix0) = ew_u.parts.msw;					\
76  (ix1) = ew_u.parts.lsw;					\
77} while (0)
78
79/* Get the more significant 32 bit int from a double.  */
80
81#define GET_HIGH_WORD(i,d)					\
82do {								\
83  ieee_double_shape_type gh_u;					\
84  gh_u.value = (d);						\
85  (i) = gh_u.parts.msw;						\
86} while (0)
87
88/* Get the less significant 32 bit int from a double.  */
89
90#define GET_LOW_WORD(i,d)					\
91do {								\
92  ieee_double_shape_type gl_u;					\
93  gl_u.value = (d);						\
94  (i) = gl_u.parts.lsw;						\
95} while (0)
96
97/* Set a double from two 32 bit ints.  */
98
99#define INSERT_WORDS(d,ix0,ix1)					\
100do {								\
101  ieee_double_shape_type iw_u;					\
102  iw_u.parts.msw = (ix0);					\
103  iw_u.parts.lsw = (ix1);					\
104  (d) = iw_u.value;						\
105} while (0)
106
107/* Set the more significant 32 bits of a double from an int.  */
108
109#define SET_HIGH_WORD(d,v)					\
110do {								\
111  ieee_double_shape_type sh_u;					\
112  sh_u.value = (d);						\
113  sh_u.parts.msw = (v);						\
114  (d) = sh_u.value;						\
115} while (0)
116
117/* Set the less significant 32 bits of a double from an int.  */
118
119#define SET_LOW_WORD(d,v)					\
120do {								\
121  ieee_double_shape_type sl_u;					\
122  sl_u.value = (d);						\
123  sl_u.parts.lsw = (v);						\
124  (d) = sl_u.value;						\
125} while (0)
126
127/*
128 * A union which permits us to convert between a float and a 32 bit
129 * int.
130 */
131
132typedef union
133{
134  float value;
135  /* FIXME: Assumes 32 bit int.  */
136  unsigned int word;
137} ieee_float_shape_type;
138
139/* Get a 32 bit int from a float.  */
140
141#define GET_FLOAT_WORD(i,d)					\
142do {								\
143  ieee_float_shape_type gf_u;					\
144  gf_u.value = (d);						\
145  (i) = gf_u.word;						\
146} while (0)
147
148/* Set a float from a 32 bit int.  */
149
150#define SET_FLOAT_WORD(d,i)					\
151do {								\
152  ieee_float_shape_type sf_u;					\
153  sf_u.word = (i);						\
154  (d) = sf_u.value;						\
155} while (0)
156
157#ifdef _COMPLEX_H
158/*
159 * Inline functions that can be used to construct complex values.
160 *
161 * The C99 standard intends x+I*y to be used for this, but x+I*y is
162 * currently unusable in general since gcc introduces many overflow,
163 * underflow, sign and efficiency bugs by rewriting I*y as
164 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
165 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
166 * to -0.0+I*0.0.
167 */
168static __inline float complex
169cpackf(float x, float y)
170{
171	float complex z;
172
173	__real__ z = x;
174	__imag__ z = y;
175	return (z);
176}
177
178static __inline double complex
179cpack(double x, double y)
180{
181	double complex z;
182
183	__real__ z = x;
184	__imag__ z = y;
185	return (z);
186}
187
188static __inline long double complex
189cpackl(long double x, long double y)
190{
191	long double complex z;
192
193	__real__ z = x;
194	__imag__ z = y;
195	return (z);
196}
197#endif /* _COMPLEX_H */
198
199/*
200 * ieee style elementary functions
201 *
202 * We rename functions here to improve other sources' diffability
203 * against fdlibm.
204 */
205#define	__ieee754_sqrt	sqrt
206#define	__ieee754_acos	acos
207#define	__ieee754_acosh	acosh
208#define	__ieee754_log	log
209#define	__ieee754_atanh	atanh
210#define	__ieee754_asin	asin
211#define	__ieee754_atan2	atan2
212#define	__ieee754_exp	exp
213#define	__ieee754_cosh	cosh
214#define	__ieee754_fmod	fmod
215#define	__ieee754_pow	pow
216#define	__ieee754_lgamma lgamma
217#define	__ieee754_gamma	gamma
218#define	__ieee754_lgamma_r lgamma_r
219#define	__ieee754_gamma_r gamma_r
220#define	__ieee754_log10	log10
221#define	__ieee754_sinh	sinh
222#define	__ieee754_hypot	hypot
223#define	__ieee754_j0	j0
224#define	__ieee754_j1	j1
225#define	__ieee754_y0	y0
226#define	__ieee754_y1	y1
227#define	__ieee754_jn	jn
228#define	__ieee754_yn	yn
229#define	__ieee754_remainder remainder
230#define	__ieee754_scalb	scalb
231#define	__ieee754_sqrtf	sqrtf
232#define	__ieee754_acosf	acosf
233#define	__ieee754_acoshf acoshf
234#define	__ieee754_logf	logf
235#define	__ieee754_atanhf atanhf
236#define	__ieee754_asinf	asinf
237#define	__ieee754_atan2f atan2f
238#define	__ieee754_expf	expf
239#define	__ieee754_coshf	coshf
240#define	__ieee754_fmodf	fmodf
241#define	__ieee754_powf	powf
242#define	__ieee754_lgammaf lgammaf
243#define	__ieee754_gammaf gammaf
244#define	__ieee754_lgammaf_r lgammaf_r
245#define	__ieee754_gammaf_r gammaf_r
246#define	__ieee754_log10f log10f
247#define	__ieee754_sinhf	sinhf
248#define	__ieee754_hypotf hypotf
249#define	__ieee754_j0f	j0f
250#define	__ieee754_j1f	j1f
251#define	__ieee754_y0f	y0f
252#define	__ieee754_y1f	y1f
253#define	__ieee754_jnf	jnf
254#define	__ieee754_ynf	ynf
255#define	__ieee754_remainderf remainderf
256#define	__ieee754_scalbf scalbf
257#define __ieee754_ldexpf ldexpf
258
259/* fdlibm kernel function */
260int	__ieee754_rem_pio2(double,double*);
261double	__kernel_sin(double,double,int);
262double	__kernel_cos(double,double);
263double	__kernel_tan(double,double,int);
264int	__kernel_rem_pio2(double*,double*,int,int,int,const int*);
265
266/* float versions of fdlibm kernel functions */
267int	__ieee754_rem_pio2f(float,float*);
268float	__kernel_sindf(double);
269float	__kernel_cosdf(double);
270float	__kernel_tandf(double,int);
271int	__kernel_rem_pio2f(float*,float*,int,int,int,const int*);
272
273#endif /* !_MATH_PRIVATE_H_ */
274