1140609Sdas/*- 2226371Sdas * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> 3140609Sdas * All rights reserved. 4140609Sdas * 5140609Sdas * Redistribution and use in source and binary forms, with or without 6140609Sdas * modification, are permitted provided that the following conditions 7140609Sdas * are met: 8140609Sdas * 1. Redistributions of source code must retain the above copyright 9140609Sdas * notice, this list of conditions and the following disclaimer. 10140609Sdas * 2. Redistributions in binary form must reproduce the above copyright 11140609Sdas * notice, this list of conditions and the following disclaimer in the 12140609Sdas * documentation and/or other materials provided with the distribution. 13140609Sdas * 14140609Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15140609Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16140609Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17140609Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18140609Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19140609Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20140609Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21140609Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22140609Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23140609Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24140609Sdas * SUCH DAMAGE. 25140609Sdas */ 26140609Sdas 27140609Sdas#include <sys/cdefs.h> 28140609Sdas__FBSDID("$FreeBSD: releng/10.2/lib/msun/src/s_fmaf.c 226371 2011-10-15 04:16:58Z das $"); 29140609Sdas 30226371Sdas#include <fenv.h> 31226371Sdas 32177712Sdas#include "math.h" 33226371Sdas#include "math_private.h" 34177712Sdas 35140609Sdas/* 36140609Sdas * Fused multiply-add: Compute x * y + z with a single rounding error. 37140609Sdas * 38140609Sdas * A double has more than twice as much precision than a float, so 39226371Sdas * direct double-precision arithmetic suffices, except where double 40226371Sdas * rounding occurs. 41140609Sdas */ 42140609Sdasfloat 43140609Sdasfmaf(float x, float y, float z) 44140609Sdas{ 45226371Sdas double xy, result; 46226371Sdas uint32_t hr, lr; 47140609Sdas 48226371Sdas xy = (double)x * y; 49226371Sdas result = xy + z; 50226371Sdas EXTRACT_WORDS(hr, lr, result); 51226371Sdas /* Common case: The double precision result is fine. */ 52226371Sdas if ((lr & 0x1fffffff) != 0x10000000 || /* not a halfway case */ 53226371Sdas (hr & 0x7ff00000) == 0x7ff00000 || /* NaN */ 54226371Sdas result - xy == z || /* exact */ 55226371Sdas fegetround() != FE_TONEAREST) /* not round-to-nearest */ 56226371Sdas return (result); 57226371Sdas 58226371Sdas /* 59226371Sdas * If result is inexact, and exactly halfway between two float values, 60226371Sdas * we need to adjust the low-order bit in the direction of the error. 61226371Sdas */ 62226371Sdas fesetround(FE_TOWARDZERO); 63226371Sdas volatile double vxy = xy; /* XXX work around gcc CSE bug */ 64226371Sdas double adjusted_result = vxy + z; 65226371Sdas fesetround(FE_TONEAREST); 66226371Sdas if (result == adjusted_result) 67226371Sdas SET_LOW_WORD(adjusted_result, lr + 1); 68226371Sdas return (adjusted_result); 69140609Sdas} 70