darwin-ldouble.c 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119
  1. // SPDX-License-Identifier: GPL-2.0+
  2. /*
  3. * Borrowed from GCC 4.2.2 (which still was GPL v2+)
  4. */
  5. /* 128-bit long double support routines for Darwin.
  6. Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
  7. Free Software Foundation, Inc.
  8. This file is part of GCC.
  9. */
  10. /*
  11. * Implementations of floating-point long double basic arithmetic
  12. * functions called by the IBM C compiler when generating code for
  13. * PowerPC platforms. In particular, the following functions are
  14. * implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
  15. * Double-double algorithms are based on the paper "Doubled-Precision
  16. * IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
  17. * 1987. An alternative published reference is "Software for
  18. * Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
  19. * ACM TOMS vol 7 no 3, September 1981, pages 272-283.
  20. */
  21. /*
  22. * Each long double is made up of two IEEE doubles. The value of the
  23. * long double is the sum of the values of the two parts. The most
  24. * significant part is required to be the value of the long double
  25. * rounded to the nearest double, as specified by IEEE. For Inf
  26. * values, the least significant part is required to be one of +0.0 or
  27. * -0.0. No other requirements are made; so, for example, 1.0 may be
  28. * represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
  29. * NaN is don't-care.
  30. *
  31. * This code currently assumes big-endian.
  32. */
  33. #define fabs(x) __builtin_fabs(x)
  34. #define isless(x, y) __builtin_isless(x, y)
  35. #define inf() __builtin_inf()
  36. #define unlikely(x) __builtin_expect((x), 0)
  37. #define nonfinite(a) unlikely(!isless(fabs(a), inf()))
  38. typedef union {
  39. long double ldval;
  40. double dval[2];
  41. } longDblUnion;
  42. /* Add two 'long double' values and return the result. */
  43. long double __gcc_qadd(double a, double aa, double c, double cc)
  44. {
  45. longDblUnion x;
  46. double z, q, zz, xh;
  47. z = a + c;
  48. if (nonfinite(z)) {
  49. z = cc + aa + c + a;
  50. if (nonfinite(z))
  51. return z;
  52. x.dval[0] = z; /* Will always be DBL_MAX. */
  53. zz = aa + cc;
  54. if (fabs(a) > fabs(c))
  55. x.dval[1] = a - z + c + zz;
  56. else
  57. x.dval[1] = c - z + a + zz;
  58. } else {
  59. q = a - z;
  60. zz = q + c + (a - (q + z)) + aa + cc;
  61. /* Keep -0 result. */
  62. if (zz == 0.0)
  63. return z;
  64. xh = z + zz;
  65. if (nonfinite(xh))
  66. return xh;
  67. x.dval[0] = xh;
  68. x.dval[1] = z - xh + zz;
  69. }
  70. return x.ldval;
  71. }
  72. long double __gcc_qsub(double a, double b, double c, double d)
  73. {
  74. return __gcc_qadd(a, b, -c, -d);
  75. }
  76. long double __gcc_qmul(double a, double b, double c, double d)
  77. {
  78. longDblUnion z;
  79. double t, tau, u, v, w;
  80. t = a * c; /* Highest order double term. */
  81. if (unlikely(t == 0) /* Preserve -0. */
  82. || nonfinite(t))
  83. return t;
  84. /* Sum terms of two highest orders. */
  85. /* Use fused multiply-add to get low part of a * c. */
  86. #ifndef __NO_FPRS__
  87. asm("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
  88. #else
  89. tau = fmsub(a, c, t);
  90. #endif
  91. v = a * d;
  92. w = b * c;
  93. tau += v + w; /* Add in other second-order terms. */
  94. u = t + tau;
  95. /* Construct long double result. */
  96. if (nonfinite(u))
  97. return u;
  98. z.dval[0] = u;
  99. z.dval[1] = (t - u) + tau;
  100. return z.ldval;
  101. }