4 #if LDBL_MANT_DIG==64 && LDBL_MAX_EXP==16384
5 /* exact add, assumes exponent_x >= exponent_y */
6 static void add(long double *hi, long double *lo, long double x, long double y)
16 /* exact mul, assumes no over/underflow */
17 static void mul(long double *hi, long double *lo, long double x, long double y)
19 static const long double c = 1.0 + 0x1p32L;
20 long double cx, xh, xl, cy, yh, yl;
29 *lo = (xh*yh - *hi) + xh*yl + xl*yh + xl*yl;
33 assume (long double)(hi+lo) == hi
34 return an adjusted hi so that rounding it to double (or less) precision is correct
36 static long double adjust(long double hi, long double lo)
38 union ldshape uhi, ulo;
46 if ((uhi.i.se & 0x8000) == (ulo.i.se & 0x8000))
49 /* handle underflow and take care of ld80 implicit msb */
50 if (uhi.i.m << 1 == 0) {
59 /* adjusted add so the result is correct when rounded to double (or less) precision */
60 static long double dadd(long double x, long double y)
66 /* adjusted mul so the result is correct when rounded to double (or less) precision */
67 static long double dmul(long double x, long double y)
73 static int getexp(long double x)
77 return u.i.se & 0x7fff;
80 double fma(double x, double y, double z)
82 #pragma STDC FENV_ACCESS ON
83 long double hi, lo1, lo2, xy;
86 /* handle +-inf,nan */
87 if (!isfinite(x) || !isfinite(y))
92 if (x == 0.0 || y == 0.0)
96 if (round == FE_TONEAREST)
101 /* exact mul and add require nearest rounding */
102 /* spurious inexact exceptions may be raised */
103 fesetround(FE_TONEAREST);
104 mul(&xy, &lo1, x, y);
108 add(&hi, &lo2, z, xy);
109 } else if (ez > exy - 12) {
110 add(&hi, &lo2, xy, z);
113 xy + z is 0, but it should be calculated with the
114 original rounding mode so the sign is correct, if the
115 compiler does not support FENV_ACCESS ON it does not
116 know about the changed rounding mode and eliminates
117 the xy + z below without the volatile memory access
122 return (xy + z_) + lo1;
127 the 12 extra bits (1guard, 11round+sticky) are needed so with
129 elo <= ehi - 11, and we use the last 10 bits in adjust so
131 gives correct result when rounded to double
137 the result is stored before return for correct precision and exceptions
139 one corner case is when the underflow flag should be raised because
140 the precise result is an inexact subnormal double, but the calculated
141 long double result is an exact subnormal double
142 (so rounding to double does not raise exceptions)
144 in nearest rounding mode dadd takes care of this: the last bit of the
145 result is adjusted so rounding sees an inexact value when it should
147 in non-nearest rounding mode fenv is used for the workaround
150 if (round == FE_TONEAREST)
151 z = dadd(hi, dadd(lo1, lo2));
153 #if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
154 int e = fetestexcept(FE_INEXACT);
155 feclearexcept(FE_INEXACT);
157 z = hi + (lo1 + lo2);
158 #if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
159 if (getexp(z) < 0x3fff-1022 && fetestexcept(FE_INEXACT))
160 feraiseexcept(FE_UNDERFLOW);
162 feraiseexcept(FE_INEXACT);
168 /* origin: FreeBSD /usr/src/lib/msun/src/s_fma.c */
170 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
171 * All rights reserved.
173 * Redistribution and use in source and binary forms, with or without
174 * modification, are permitted provided that the following conditions
176 * 1. Redistributions of source code must retain the above copyright
177 * notice, this list of conditions and the following disclaimer.
178 * 2. Redistributions in binary form must reproduce the above copyright
179 * notice, this list of conditions and the following disclaimer in the
180 * documentation and/or other materials provided with the distribution.
182 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
183 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
184 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
185 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
186 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
187 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
188 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
189 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
190 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
191 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
196 * A struct dd represents a floating-point number with twice the precision
197 * of a double. We maintain the invariant that "hi" stores the 53 high-order
198 * bits of the result.
206 * Compute a+b exactly, returning the exact result in a struct dd. We assume
207 * that both a and b are finite, but make no assumptions about their relative
210 static inline struct dd dd_add(double a, double b)
217 ret.lo = (a - (ret.hi - s)) + (b - s);
222 * Compute a+b, with a small tweak: The least significant bit of the
223 * result is adjusted into a sticky bit summarizing all the bits that
224 * were lost to rounding. This adjustment negates the effects of double
225 * rounding when the result is added to another number with a higher
226 * exponent. For an explanation of round and sticky bits, see any reference
227 * on FPU design, e.g.,
229 * J. Coonen. An Implementation Guide to a Proposed Standard for
230 * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
232 static inline double add_adjusted(double a, double b)
235 union {double f; uint64_t i;} uhi, ulo;
240 if ((uhi.i & 1) == 0) {
241 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
243 uhi.i += 1 - ((uhi.i ^ ulo.i) >> 62);
251 * Compute ldexp(a+b, scale) with a single rounding error. It is assumed
252 * that the result will be subnormal, and care is taken to ensure that
253 * double rounding does not occur.
255 static inline double add_and_denormalize(double a, double b, int scale)
258 union {double f; uint64_t i;} uhi, ulo;
264 * If we are losing at least two bits of accuracy to denormalization,
265 * then the first lost bit becomes a round bit, and we adjust the
266 * lowest bit of sum.hi to make it a sticky bit summarizing all the
267 * bits in sum.lo. With the sticky bit adjusted, the hardware will
268 * break any ties in the correct direction.
270 * If we are losing only one bit to denormalization, however, we must
271 * break the ties manually.
275 bits_lost = -((int)(uhi.i >> 52) & 0x7ff) - scale + 1;
276 if ((bits_lost != 1) ^ (int)(uhi.i & 1)) {
277 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
279 uhi.i += 1 - (((uhi.i ^ ulo.i) >> 62) & 2);
283 return scalbn(sum.hi, scale);
287 * Compute a*b exactly, returning the exact result in a struct dd. We assume
288 * that both a and b are normalized, so no underflow or overflow will occur.
289 * The current rounding mode must be round-to-nearest.
291 static inline struct dd dd_mul(double a, double b)
293 static const double split = 0x1p27 + 1.0;
295 double ha, hb, la, lb, p, q;
308 q = ha * lb + la * hb;
311 ret.lo = p - ret.hi + q + la * lb;
316 * Fused multiply-add: Compute x * y + z with a single rounding error.
318 * We use scaling to avoid overflow/underflow, along with the
319 * canonical precision-doubling technique adapted from:
321 * Dekker, T. A Floating-Point Technique for Extending the
322 * Available Precision. Numer. Math. 18, 224-242 (1971).
324 * This algorithm is sensitive to the rounding precision. FPUs such
325 * as the i387 must be set in double-precision mode if variables are
326 * to be stored in FP registers in order to avoid incorrect results.
327 * This is the default on FreeBSD, but not on many other systems.
329 * Hardware instructions should be used on architectures that support it,
330 * since this implementation will likely be several times slower.
332 double fma(double x, double y, double z)
334 #pragma STDC FENV_ACCESS ON
335 double xs, ys, zs, adj;
342 * Handle special cases. The order of operations and the particular
343 * return values here are crucial in handling special cases involving
344 * infinities, NaNs, overflows, and signed zeroes correctly.
346 if (!isfinite(x) || !isfinite(y))
350 if (x == 0.0 || y == 0.0)
358 oround = fegetround();
359 spread = ex + ey - ez;
362 * If x * y and z are many orders of magnitude apart, the scaling
363 * will overflow, so we handle these cases specially. Rounding
364 * modes other than FE_TONEAREST are painful.
366 if (spread < -DBL_MANT_DIG) {
368 feraiseexcept(FE_INEXACT);
372 feraiseexcept(FE_UNDERFLOW);
375 default: /* FE_TONEAREST */
379 if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
382 return (nextafter(z, 0));
386 if (x > 0.0 ^ y < 0.0)
389 return (nextafter(z, -INFINITY));
393 if (x > 0.0 ^ y < 0.0)
394 return (nextafter(z, INFINITY));
400 if (spread <= DBL_MANT_DIG * 2)
401 zs = scalbn(zs, -spread);
403 zs = copysign(DBL_MIN, zs);
405 fesetround(FE_TONEAREST);
408 * Basic approach for round-to-nearest:
410 * (xy.hi, xy.lo) = x * y (exact)
411 * (r.hi, r.lo) = xy.hi + z (exact)
412 * adj = xy.lo + r.lo (inexact; low bit is sticky)
413 * result = r.hi + adj (correctly rounded)
416 r = dd_add(xy.hi, zs);
422 * When the addends cancel to 0, ensure that the result has
426 volatile double vzs = zs; /* XXX gcc CSE bug workaround */
427 return xy.hi + vzs + scalbn(xy.lo, spread);
430 if (oround != FE_TONEAREST) {
432 * There is no need to worry about double rounding in directed
434 * But underflow may not be raised properly, example in downward rounding:
435 * fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
438 #if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
439 int e = fetestexcept(FE_INEXACT);
440 feclearexcept(FE_INEXACT);
444 ret = scalbn(r.hi + adj, spread);
445 #if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
446 if (ilogb(ret) < -1022 && fetestexcept(FE_INEXACT))
447 feraiseexcept(FE_UNDERFLOW);
449 feraiseexcept(FE_INEXACT);
454 adj = add_adjusted(r.lo, xy.lo);
455 if (spread + ilogb(r.hi) > -1023)
456 return scalbn(r.hi + adj, spread);
458 return add_and_denormalize(r.hi, adj, spread);