third-party-mirror / mingw-w64 / 9b207d4da4806fdaea3bf9f4a5d66dc6056670b6 / . / mingw-w64-crt / math / fma.c

/** | |

* This file has no copyright assigned and is placed in the Public Domain. | |

* This file is part of the mingw-w64 runtime package. | |

* No warranty is given; refer to the file DISCLAIMER.PD within this package. | |

*/ | |

double fma(double x, double y, double z); | |

#if defined(_ARM_) || defined(__arm__) | |

/* Use hardware FMA on ARM. */ | |

double fma(double x, double y, double z){ | |

__asm__ ( | |

"fmacd %0, %1, %2 \n" | |

: "+w"(z) | |

: "w"(x), "w"(y) | |

); | |

return z; | |

} | |

#elif defined(_ARM64_) || defined(__aarch64__) | |

/* Use hardware FMA on ARM64. */ | |

double fma(double x, double y, double z){ | |

__asm__ ( | |

"fmadd %d0, %d1, %d2, %d0 \n" | |

: "+w"(z) | |

: "w"(x), "w"(y) | |

); | |

return z; | |

} | |

#elif defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__) | |

#include <math.h> | |

#include <stdint.h> | |

/* This is in accordance with the IEC 559 double-precision format. | |

* Be advised that due to the hidden bit, the higher half actually has 26 bits. | |

* Multiplying two 27-bit numbers will cause a 1-ULP error, which we cannot | |

* avoid. It is kept in the very last position. | |

*/ | |

typedef union iec559_double_ { | |

struct __attribute__((__packed__)) { | |

uint64_t mlo : 27; | |

uint64_t mhi : 25; | |

uint64_t exp : 11; | |

uint64_t sgn : 1; | |

}; | |

double f; | |

} iec559_double; | |

static inline void break_down(iec559_double *restrict lo, iec559_double *restrict hi, double x) { | |

hi->f = x; | |

/* Erase low-order significant bits. `hi->f` now has only 26 significant bits. */ | |

hi->mlo = 0; | |

/* Store the low-order half. It will be normalized by the hardware. */ | |

lo->f = x - hi->f; | |

/* Preserve signness in case of zero. */ | |

lo->sgn = hi->sgn; | |

} | |

double fma(double x, double y, double z) { | |

/* | |

POSIX-2013: | |

1. If x or y are NaN, a NaN shall be returned. | |

2. If x multiplied by y is an exact infinity and z is also an infinity | |

but with the opposite sign, a domain error shall occur, and either a NaN | |

(if supported), or an implementation-defined value shall be returned. | |

3. If one of x and y is infinite, the other is zero, and z is not a NaN, | |

a domain error shall occur, and either a NaN (if supported), or an | |

implementation-defined value shall be returned. | |

4. If one of x and y is infinite, the other is zero, and z is a NaN, a NaN | |

shall be returned and a domain error may occur. | |

5. If x* y is not 0*Inf nor Inf*0 and z is a NaN, a NaN shall be returned. | |

*/ | |

/* Check whether the result is finite. */ | |

double ret = x * y + z; | |

if(!isfinite(ret)) { | |

return ret; /* If this naive check doesn't yield a finite value, the FMA isn't | |

likely to return one either. Forward the value as is. */ | |

} | |

iec559_double xlo, xhi, ylo, yhi; | |

break_down(&xlo, &xhi, x); | |

break_down(&ylo, &yhi, y); | |

/* The order of these four statements is essential. Don't move them around. */ | |

ret = z; | |

ret += xhi.f * yhi.f; /* The most significant item comes first. */ | |

ret += xhi.f * ylo.f + xlo.f * yhi.f; /* They are equally significant. */ | |

ret += xlo.f * ylo.f; /* The least significant item comes last. */ | |

return ret; | |

} | |

#else | |

#error Please add FMA implementation for this platform. | |

#endif |