math_private.h 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
  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. * from: @(#)fdlibm.h 5.1 93/09/24
  13. * $FreeBSD$
  14. */
  15. #ifndef _MATH_PRIVATE_H_
  16. #define _MATH_PRIVATE_H_
  17. #include <sys/types.h>
  18. #if HAVE_MACHINE_ENDIAN_H
  19. #include <machine/endian.h>
  20. #elif HAVE_SYS_ENDIAN_H && HOST_ANDROID
  21. /* Android unified headers don't have machine/endian.h */
  22. #include <sys/endian.h>
  23. #endif
  24. /*
  25. * The original fdlibm code used statements like:
  26. * n0 = ((*(int*)&one)>>29)^1; * index of high word *
  27. * ix0 = *(n0+(int*)&x); * high word of x *
  28. * ix1 = *((1-n0)+(int*)&x); * low word of x *
  29. * to dig two 32 bit words out of the 64 bit IEEE floating point
  30. * value. That is non-ANSI, and, moreover, the gcc instruction
  31. * scheduler gets it wrong. We instead use the following macros.
  32. * Unlike the original code, we determine the endianness at compile
  33. * time, not at run time; I don't see much benefit to selecting
  34. * endianness at run time.
  35. */
  36. /*
  37. * A union which permits us to convert between a double and two 32 bit
  38. * ints.
  39. */
  40. #ifdef __arm__
  41. #if defined(__VFP_FP__) || defined(__ARM_EABI__)
  42. #define IEEE_WORD_ORDER BYTE_ORDER
  43. #else
  44. #define IEEE_WORD_ORDER BIG_ENDIAN
  45. #endif
  46. #else /* __arm__ */
  47. #define IEEE_WORD_ORDER BYTE_ORDER
  48. #endif
  49. #if IEEE_WORD_ORDER == BIG_ENDIAN
  50. typedef union
  51. {
  52. double value;
  53. struct
  54. {
  55. u_int32_t msw;
  56. u_int32_t lsw;
  57. } parts;
  58. struct
  59. {
  60. u_int64_t w;
  61. } xparts;
  62. } ieee_double_shape_type;
  63. #endif
  64. #if IEEE_WORD_ORDER == LITTLE_ENDIAN
  65. typedef union
  66. {
  67. double value;
  68. struct
  69. {
  70. u_int32_t lsw;
  71. u_int32_t msw;
  72. } parts;
  73. struct
  74. {
  75. u_int64_t w;
  76. } xparts;
  77. } ieee_double_shape_type;
  78. #endif
  79. /* Get two 32 bit ints from a double. */
  80. #define EXTRACT_WORDS(ix0,ix1,d) \
  81. do { \
  82. ieee_double_shape_type ew_u; \
  83. ew_u.value = (d); \
  84. (ix0) = ew_u.parts.msw; \
  85. (ix1) = ew_u.parts.lsw; \
  86. } while (0)
  87. /* Get a 64-bit int from a double. */
  88. #define EXTRACT_WORD64(ix,d) \
  89. do { \
  90. ieee_double_shape_type ew_u; \
  91. ew_u.value = (d); \
  92. (ix) = ew_u.xparts.w; \
  93. } while (0)
  94. /* Get the more significant 32 bit int from a double. */
  95. #define GET_HIGH_WORD(i,d) \
  96. do { \
  97. ieee_double_shape_type gh_u; \
  98. gh_u.value = (d); \
  99. (i) = gh_u.parts.msw; \
  100. } while (0)
  101. /* Get the less significant 32 bit int from a double. */
  102. #define GET_LOW_WORD(i,d) \
  103. do { \
  104. ieee_double_shape_type gl_u; \
  105. gl_u.value = (d); \
  106. (i) = gl_u.parts.lsw; \
  107. } while (0)
  108. /* Set a double from two 32 bit ints. */
  109. #define INSERT_WORDS(d,ix0,ix1) \
  110. do { \
  111. ieee_double_shape_type iw_u; \
  112. iw_u.parts.msw = (ix0); \
  113. iw_u.parts.lsw = (ix1); \
  114. (d) = iw_u.value; \
  115. } while (0)
  116. /* Set a double from a 64-bit int. */
  117. #define INSERT_WORD64(d,ix) \
  118. do { \
  119. ieee_double_shape_type iw_u; \
  120. iw_u.xparts.w = (ix); \
  121. (d) = iw_u.value; \
  122. } while (0)
  123. /* Set the more significant 32 bits of a double from an int. */
  124. #define SET_HIGH_WORD(d,v) \
  125. do { \
  126. ieee_double_shape_type sh_u; \
  127. sh_u.value = (d); \
  128. sh_u.parts.msw = (v); \
  129. (d) = sh_u.value; \
  130. } while (0)
  131. /* Set the less significant 32 bits of a double from an int. */
  132. #define SET_LOW_WORD(d,v) \
  133. do { \
  134. ieee_double_shape_type sl_u; \
  135. sl_u.value = (d); \
  136. sl_u.parts.lsw = (v); \
  137. (d) = sl_u.value; \
  138. } while (0)
  139. /*
  140. * A union which permits us to convert between a float and a 32 bit
  141. * int.
  142. */
  143. typedef union
  144. {
  145. float value;
  146. /* FIXME: Assumes 32 bit int. */
  147. unsigned int word;
  148. } ieee_float_shape_type;
  149. /* Get a 32 bit int from a float. */
  150. #define GET_FLOAT_WORD(i,d) \
  151. do { \
  152. ieee_float_shape_type gf_u; \
  153. gf_u.value = (d); \
  154. (i) = gf_u.word; \
  155. } while (0)
  156. /* Set a float from a 32 bit int. */
  157. #define SET_FLOAT_WORD(d,i) \
  158. do { \
  159. ieee_float_shape_type sf_u; \
  160. sf_u.word = (i); \
  161. (d) = sf_u.value; \
  162. } while (0)
  163. /*
  164. * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
  165. * double.
  166. */
  167. #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
  168. do { \
  169. union IEEEl2bits ew_u; \
  170. ew_u.e = (d); \
  171. (ix0) = ew_u.xbits.expsign; \
  172. (ix1) = ew_u.xbits.man; \
  173. } while (0)
  174. /*
  175. * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
  176. * long double.
  177. */
  178. #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
  179. do { \
  180. union IEEEl2bits ew_u; \
  181. ew_u.e = (d); \
  182. (ix0) = ew_u.xbits.expsign; \
  183. (ix1) = ew_u.xbits.manh; \
  184. (ix2) = ew_u.xbits.manl; \
  185. } while (0)
  186. /* Get expsign as a 16 bit int from a long double. */
  187. #define GET_LDBL_EXPSIGN(i,d) \
  188. do { \
  189. union IEEEl2bits ge_u; \
  190. ge_u.e = (d); \
  191. (i) = ge_u.xbits.expsign; \
  192. } while (0)
  193. /*
  194. * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
  195. * mantissa.
  196. */
  197. #define INSERT_LDBL80_WORDS(d,ix0,ix1) \
  198. do { \
  199. union IEEEl2bits iw_u; \
  200. iw_u.xbits.expsign = (ix0); \
  201. iw_u.xbits.man = (ix1); \
  202. (d) = iw_u.e; \
  203. } while (0)
  204. /*
  205. * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
  206. * comprising the mantissa.
  207. */
  208. #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
  209. do { \
  210. union IEEEl2bits iw_u; \
  211. iw_u.xbits.expsign = (ix0); \
  212. iw_u.xbits.manh = (ix1); \
  213. iw_u.xbits.manl = (ix2); \
  214. (d) = iw_u.e; \
  215. } while (0)
  216. /* Set expsign of a long double from a 16 bit int. */
  217. #define SET_LDBL_EXPSIGN(d,v) \
  218. do { \
  219. union IEEEl2bits se_u; \
  220. se_u.e = (d); \
  221. se_u.xbits.expsign = (v); \
  222. (d) = se_u.e; \
  223. } while (0)
  224. #ifdef __i386__
  225. /* Long double constants are broken on i386. */
  226. #define LD80C(m, ex, v) { \
  227. .xbits.man = __CONCAT(m, ULL), \
  228. .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \
  229. }
  230. #else
  231. /* The above works on non-i386 too, but we use this to check v. */
  232. #define LD80C(m, ex, v) { .e = (v), }
  233. #endif
  234. #ifdef FLT_EVAL_METHOD
  235. /*
  236. * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
  237. */
  238. #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
  239. #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
  240. #else
  241. #define STRICT_ASSIGN(type, lval, rval) do { \
  242. volatile type __lval; \
  243. \
  244. if (sizeof(type) >= sizeof(long double)) \
  245. (lval) = (rval); \
  246. else { \
  247. __lval = (rval); \
  248. (lval) = __lval; \
  249. } \
  250. } while (0)
  251. #endif
  252. #endif /* FLT_EVAL_METHOD */
  253. /* Support switching the mode to FP_PE if necessary. */
  254. #if defined(__i386__) && !defined(NO_FPSETPREC)
  255. #define ENTERI() \
  256. long double __retval; \
  257. fp_prec_t __oprec; \
  258. \
  259. if ((__oprec = fpgetprec()) != FP_PE) \
  260. fpsetprec(FP_PE)
  261. #define RETURNI(x) do { \
  262. __retval = (x); \
  263. if (__oprec != FP_PE) \
  264. fpsetprec(__oprec); \
  265. RETURNF(__retval); \
  266. } while (0)
  267. #else
  268. #define ENTERI(x)
  269. #define RETURNI(x) RETURNF(x)
  270. #endif
  271. /* Default return statement if hack*_t() is not used. */
  272. #define RETURNF(v) return (v)
  273. /*
  274. * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
  275. * a == 0, but is slower.
  276. */
  277. #define _2sum(a, b) do { \
  278. __typeof(a) __s, __w; \
  279. \
  280. __w = (a) + (b); \
  281. __s = __w - (a); \
  282. (b) = ((a) - (__w - __s)) + ((b) - __s); \
  283. (a) = __w; \
  284. } while (0)
  285. /*
  286. * 2sumF algorithm.
  287. *
  288. * "Normalize" the terms in the infinite-precision expression a + b for
  289. * the sum of 2 floating point values so that b is as small as possible
  290. * relative to 'a'. (The resulting 'a' is the value of the expression in
  291. * the same precision as 'a' and the resulting b is the rounding error.)
  292. * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
  293. * exponent overflow or underflow must not occur. This uses a Theorem of
  294. * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
  295. * is apparently due to Skewchuk (1997).
  296. *
  297. * For this to always work, assignment of a + b to 'a' must not retain any
  298. * extra precision in a + b. This is required by C standards but broken
  299. * in many compilers. The brokenness cannot be worked around using
  300. * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
  301. * algorithm would be destroyed by non-null strict assignments. (The
  302. * compilers are correct to be broken -- the efficiency of all floating
  303. * point code calculations would be destroyed similarly if they forced the
  304. * conversions.)
  305. *
  306. * Fortunately, a case that works well can usually be arranged by building
  307. * any extra precision into the type of 'a' -- 'a' should have type float_t,
  308. * double_t or long double. b's type should be no larger than 'a's type.
  309. * Callers should use these types with scopes as large as possible, to
  310. * reduce their own extra-precision and efficiciency problems. In
  311. * particular, they shouldn't convert back and forth just to call here.
  312. */
  313. #ifdef DEBUG
  314. #define _2sumF(a, b) do { \
  315. __typeof(a) __w; \
  316. volatile __typeof(a) __ia, __ib, __r, __vw; \
  317. \
  318. __ia = (a); \
  319. __ib = (b); \
  320. assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
  321. \
  322. __w = (a) + (b); \
  323. (b) = ((a) - __w) + (b); \
  324. (a) = __w; \
  325. \
  326. /* The next 2 assertions are weak if (a) is already long double. */ \
  327. assert((long double)__ia + __ib == (long double)(a) + (b)); \
  328. __vw = __ia + __ib; \
  329. __r = __ia - __vw; \
  330. __r += __ib; \
  331. assert(__vw == (a) && __r == (b)); \
  332. } while (0)
  333. #else /* !DEBUG */
  334. #define _2sumF(a, b) do { \
  335. __typeof(a) __w; \
  336. \
  337. __w = (a) + (b); \
  338. (b) = ((a) - __w) + (b); \
  339. (a) = __w; \
  340. } while (0)
  341. #endif /* DEBUG */
  342. /*
  343. * Set x += c, where x is represented in extra precision as a + b.
  344. * x must be sufficiently normalized and sufficiently larger than c,
  345. * and the result is then sufficiently normalized.
  346. *
  347. * The details of ordering are that |a| must be >= |c| (so that (a, c)
  348. * can be normalized without extra work to swap 'a' with c). The details of
  349. * the normalization are that b must be small relative to the normalized 'a'.
  350. * Normalization of (a, c) makes the normalized c tiny relative to the
  351. * normalized a, so b remains small relative to 'a' in the result. However,
  352. * b need not ever be tiny relative to 'a'. For example, b might be about
  353. * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
  354. * That is usually enough, and adding c (which by normalization is about
  355. * 2**53 times smaller than a) cannot change b significantly. However,
  356. * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
  357. * significantly relative to b. The caller must ensure that significant
  358. * cancellation doesn't occur, either by having c of the same sign as 'a',
  359. * or by having |c| a few percent smaller than |a|. Pre-normalization of
  360. * (a, b) may help.
  361. *
  362. * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
  363. * exercise 19). We gain considerable efficiency by requiring the terms to
  364. * be sufficiently normalized and sufficiently increasing.
  365. */
  366. #define _3sumF(a, b, c) do { \
  367. __typeof(a) __tmp; \
  368. \
  369. __tmp = (c); \
  370. _2sumF(__tmp, (a)); \
  371. (b) += (a); \
  372. (a) = __tmp; \
  373. } while (0)
  374. /*
  375. * Common routine to process the arguments to nan(), nanf(), and nanl().
  376. */
  377. void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
  378. #ifdef _COMPLEX_H
  379. /*
  380. * C99 specifies that complex numbers have the same representation as
  381. * an array of two elements, where the first element is the real part
  382. * and the second element is the imaginary part.
  383. */
  384. typedef union {
  385. float complex f;
  386. float a[2];
  387. } float_complex;
  388. typedef union {
  389. double complex f;
  390. double a[2];
  391. } double_complex;
  392. typedef union {
  393. long double complex f;
  394. long double a[2];
  395. } long_double_complex;
  396. #define REALPART(z) ((z).a[0])
  397. #define IMAGPART(z) ((z).a[1])
  398. /*
  399. * Inline functions that can be used to construct complex values.
  400. *
  401. * The C99 standard intends x+I*y to be used for this, but x+I*y is
  402. * currently unusable in general since gcc introduces many overflow,
  403. * underflow, sign and efficiency bugs by rewriting I*y as
  404. * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
  405. * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
  406. * to -0.0+I*0.0.
  407. */
  408. static __inline float complex
  409. cpackf(float x, float y)
  410. {
  411. float_complex z;
  412. REALPART(z) = x;
  413. IMAGPART(z) = y;
  414. return (z.f);
  415. }
  416. static __inline double complex
  417. cpack(double x, double y)
  418. {
  419. double_complex z;
  420. REALPART(z) = x;
  421. IMAGPART(z) = y;
  422. return (z.f);
  423. }
  424. static __inline long double complex
  425. cpackl(long double x, long double y)
  426. {
  427. long_double_complex z;
  428. REALPART(z) = x;
  429. IMAGPART(z) = y;
  430. return (z.f);
  431. }
  432. #endif /* _COMPLEX_H */
  433. #ifdef __GNUCLIKE_ASM
  434. /* Asm versions of some functions. */
  435. #ifdef __amd64__
  436. static __inline int
  437. irint(double x)
  438. {
  439. int n;
  440. asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
  441. return (n);
  442. }
  443. #define HAVE_EFFICIENT_IRINT
  444. #endif
  445. #ifdef __i386__
  446. static __inline int
  447. irint(double x)
  448. {
  449. int n;
  450. asm("fistl %0" : "=m" (n) : "t" (x));
  451. return (n);
  452. }
  453. #define HAVE_EFFICIENT_IRINT
  454. #endif
  455. #if defined(__amd64__) || defined(__i386__)
  456. static __inline int
  457. irintl(long double x)
  458. {
  459. int n;
  460. asm("fistl %0" : "=m" (n) : "t" (x));
  461. return (n);
  462. }
  463. #define HAVE_EFFICIENT_IRINTL
  464. #endif
  465. #endif /* __GNUCLIKE_ASM */
  466. #ifdef DEBUG
  467. #if defined(__amd64__) || defined(__i386__)
  468. #define breakpoint() asm("int $3")
  469. #else
  470. #include <signal.h>
  471. #define breakpoint() raise(SIGTRAP)
  472. #endif
  473. #endif
  474. /* Write a pari script to test things externally. */
  475. #ifdef DOPRINT
  476. #include <stdio.h>
  477. #ifndef DOPRINT_SWIZZLE
  478. #define DOPRINT_SWIZZLE 0
  479. #endif
  480. #ifdef DOPRINT_LD80
  481. #define DOPRINT_START(xp) do { \
  482. uint64_t __lx; \
  483. uint16_t __hx; \
  484. \
  485. /* Hack to give more-problematic args. */ \
  486. EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
  487. __lx ^= DOPRINT_SWIZZLE; \
  488. INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
  489. printf("x = %.21Lg; ", (long double)*xp); \
  490. } while (0)
  491. #define DOPRINT_END1(v) \
  492. printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
  493. #define DOPRINT_END2(hi, lo) \
  494. printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
  495. (long double)(hi), (long double)(lo))
  496. #elif defined(DOPRINT_D64)
  497. #define DOPRINT_START(xp) do { \
  498. uint32_t __hx, __lx; \
  499. \
  500. EXTRACT_WORDS(__hx, __lx, *xp); \
  501. __lx ^= DOPRINT_SWIZZLE; \
  502. INSERT_WORDS(*xp, __hx, __lx); \
  503. printf("x = %.21Lg; ", (long double)*xp); \
  504. } while (0)
  505. #define DOPRINT_END1(v) \
  506. printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
  507. #define DOPRINT_END2(hi, lo) \
  508. printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
  509. (long double)(hi), (long double)(lo))
  510. #elif defined(DOPRINT_F32)
  511. #define DOPRINT_START(xp) do { \
  512. uint32_t __hx; \
  513. \
  514. GET_FLOAT_WORD(__hx, *xp); \
  515. __hx ^= DOPRINT_SWIZZLE; \
  516. SET_FLOAT_WORD(*xp, __hx); \
  517. printf("x = %.21Lg; ", (long double)*xp); \
  518. } while (0)
  519. #define DOPRINT_END1(v) \
  520. printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
  521. #define DOPRINT_END2(hi, lo) \
  522. printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
  523. (long double)(hi), (long double)(lo))
  524. #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
  525. #ifndef DOPRINT_SWIZZLE_HIGH
  526. #define DOPRINT_SWIZZLE_HIGH 0
  527. #endif
  528. #define DOPRINT_START(xp) do { \
  529. uint64_t __lx, __llx; \
  530. uint16_t __hx; \
  531. \
  532. EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
  533. __llx ^= DOPRINT_SWIZZLE; \
  534. __lx ^= DOPRINT_SWIZZLE_HIGH; \
  535. INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
  536. printf("x = %.36Lg; ", (long double)*xp); \
  537. } while (0)
  538. #define DOPRINT_END1(v) \
  539. printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
  540. #define DOPRINT_END2(hi, lo) \
  541. printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
  542. (long double)(hi), (long double)(lo))
  543. #endif /* DOPRINT_LD80 */
  544. #else /* !DOPRINT */
  545. #define DOPRINT_START(xp)
  546. #define DOPRINT_END1(v)
  547. #define DOPRINT_END2(hi, lo)
  548. #endif /* DOPRINT */
  549. #define RETURNP(x) do { \
  550. DOPRINT_END1(x); \
  551. RETURNF(x); \
  552. } while (0)
  553. #define RETURNPI(x) do { \
  554. DOPRINT_END1(x); \
  555. RETURNI(x); \
  556. } while (0)
  557. #define RETURN2P(x, y) do { \
  558. DOPRINT_END2((x), (y)); \
  559. RETURNF((x) + (y)); \
  560. } while (0)
  561. #define RETURN2PI(x, y) do { \
  562. DOPRINT_END2((x), (y)); \
  563. RETURNI((x) + (y)); \
  564. } while (0)
  565. #ifdef STRUCT_RETURN
  566. #define RETURNSP(rp) do { \
  567. if (!(rp)->lo_set) \
  568. RETURNP((rp)->hi); \
  569. RETURN2P((rp)->hi, (rp)->lo); \
  570. } while (0)
  571. #define RETURNSPI(rp) do { \
  572. if (!(rp)->lo_set) \
  573. RETURNPI((rp)->hi); \
  574. RETURN2PI((rp)->hi, (rp)->lo); \
  575. } while (0)
  576. #endif
  577. #define SUM2P(x, y) ({ \
  578. const __typeof (x) __x = (x); \
  579. const __typeof (y) __y = (y); \
  580. \
  581. DOPRINT_END2(__x, __y); \
  582. __x + __y; \
  583. })
  584. /*
  585. * ieee style elementary functions
  586. *
  587. * We rename functions here to improve other sources' diffability
  588. * against fdlibm.
  589. */
  590. #define __ieee754_sqrt sqrt
  591. #define __ieee754_acos acos
  592. #define __ieee754_acosh acosh
  593. #define __ieee754_log log
  594. #define __ieee754_log2 log2
  595. #define __ieee754_atanh atanh
  596. #define __ieee754_asin asin
  597. #define __ieee754_atan2 atan2
  598. #define __ieee754_exp exp
  599. #define __ieee754_cosh cosh
  600. #define __ieee754_fmod fmod
  601. #define __ieee754_pow pow
  602. #define __ieee754_lgamma lgamma
  603. #define __ieee754_gamma gamma
  604. #define __ieee754_lgamma_r lgamma_r
  605. #define __ieee754_gamma_r gamma_r
  606. #define __ieee754_log10 log10
  607. #define __ieee754_sinh sinh
  608. #define __ieee754_hypot hypot
  609. #define __ieee754_j0 j0
  610. #define __ieee754_j1 j1
  611. #define __ieee754_y0 y0
  612. #define __ieee754_y1 y1
  613. #define __ieee754_jn jn
  614. #define __ieee754_yn yn
  615. #define __ieee754_remainder remainder
  616. #define __ieee754_scalb scalb
  617. #define __ieee754_sqrtf sqrtf
  618. #define __ieee754_acosf acosf
  619. #define __ieee754_acoshf acoshf
  620. #define __ieee754_logf logf
  621. #define __ieee754_atanhf atanhf
  622. #define __ieee754_asinf asinf
  623. #define __ieee754_atan2f atan2f
  624. #define __ieee754_expf expf
  625. #define __ieee754_coshf coshf
  626. #define __ieee754_fmodf fmodf
  627. #define __ieee754_powf powf
  628. #define __ieee754_lgammaf lgammaf
  629. #define __ieee754_gammaf gammaf
  630. #define __ieee754_lgammaf_r lgammaf_r
  631. #define __ieee754_gammaf_r gammaf_r
  632. #define __ieee754_log10f log10f
  633. #define __ieee754_log2f log2f
  634. #define __ieee754_sinhf sinhf
  635. #define __ieee754_hypotf hypotf
  636. #define __ieee754_j0f j0f
  637. #define __ieee754_j1f j1f
  638. #define __ieee754_y0f y0f
  639. #define __ieee754_y1f y1f
  640. #define __ieee754_jnf jnf
  641. #define __ieee754_ynf ynf
  642. #define __ieee754_remainderf remainderf
  643. #define __ieee754_scalbf scalbf
  644. /* fdlibm kernel function */
  645. int __kernel_rem_pio2(double*,double*,int,int,int);
  646. /* double precision kernel functions */
  647. #ifndef INLINE_REM_PIO2
  648. int __ieee754_rem_pio2(double,double*);
  649. #endif
  650. double __kernel_sin(double,double,int);
  651. double __kernel_cos(double,double);
  652. double __kernel_tan(double,double,int);
  653. double __ldexp_exp(double,int);
  654. #ifdef _COMPLEX_H
  655. double complex __ldexp_cexp(double complex,int);
  656. #endif
  657. /* float precision kernel functions */
  658. #ifndef INLINE_REM_PIO2F
  659. int __ieee754_rem_pio2f(float,double*);
  660. #endif
  661. #ifndef INLINE_KERNEL_SINDF
  662. float __kernel_sindf(double);
  663. #endif
  664. #ifndef INLINE_KERNEL_COSDF
  665. float __kernel_cosdf(double);
  666. #endif
  667. #ifndef INLINE_KERNEL_TANDF
  668. float __kernel_tandf(double,int);
  669. #endif
  670. float __ldexp_expf(float,int);
  671. #ifdef _COMPLEX_H
  672. float complex __ldexp_cexpf(float complex,int);
  673. #endif
  674. /* long double precision kernel functions */
  675. long double __kernel_sinl(long double, long double, int);
  676. long double __kernel_cosl(long double, long double);
  677. long double __kernel_tanl(long double, long double, int);
  678. #endif /* !_MATH_PRIVATE_H_ */