| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_fmal.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_fmal.c */ |
| 2 /*- | 2 /*- |
| 3 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> | 3 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> |
| 4 * All rights reserved. | 4 * All rights reserved. |
| 5 * | 5 * |
| 6 * Redistribution and use in source and binary forms, with or without | 6 * Redistribution and use in source and binary forms, with or without |
| 7 * modification, are permitted provided that the following conditions | 7 * modification, are permitted provided that the following conditions |
| 8 * are met: | 8 * are met: |
| 9 * 1. Redistributions of source code must retain the above copyright | 9 * 1. Redistributions of source code must retain the above copyright |
| 10 * notice, this list of conditions and the following disclaimer. | 10 * notice, this list of conditions and the following disclaimer. |
| (...skipping 146 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 157 * Fused multiply-add: Compute x * y + z with a single rounding error. | 157 * Fused multiply-add: Compute x * y + z with a single rounding error. |
| 158 * | 158 * |
| 159 * We use scaling to avoid overflow/underflow, along with the | 159 * We use scaling to avoid overflow/underflow, along with the |
| 160 * canonical precision-doubling technique adapted from: | 160 * canonical precision-doubling technique adapted from: |
| 161 * | 161 * |
| 162 * Dekker, T. A Floating-Point Technique for Extending the | 162 * Dekker, T. A Floating-Point Technique for Extending the |
| 163 * Available Precision. Numer. Math. 18, 224-242 (1971). | 163 * Available Precision. Numer. Math. 18, 224-242 (1971). |
| 164 */ | 164 */ |
| 165 long double fmal(long double x, long double y, long double z) | 165 long double fmal(long double x, long double y, long double z) |
| 166 { | 166 { |
| 167 » #pragma STDC FENV_ACCESS ON | 167 » PRAGMA_STDC_FENV_ACCESS_ON |
| 168 long double xs, ys, zs, adj; | 168 long double xs, ys, zs, adj; |
| 169 struct dd xy, r; | 169 struct dd xy, r; |
| 170 int oround; | 170 int oround; |
| 171 int ex, ey, ez; | 171 int ex, ey, ez; |
| 172 int spread; | 172 int spread; |
| 173 | 173 |
| 174 /* | 174 /* |
| 175 * Handle special cases. The order of operations and the particular | 175 * Handle special cases. The order of operations and the particular |
| 176 * return values here are crucial in handling special cases involving | 176 * return values here are crucial in handling special cases involving |
| 177 * infinities, NaNs, overflows, and signed zeroes correctly. | 177 * infinities, NaNs, overflows, and signed zeroes correctly. |
| (...skipping 106 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 284 return ret; | 284 return ret; |
| 285 } | 285 } |
| 286 | 286 |
| 287 adj = add_adjusted(r.lo, xy.lo); | 287 adj = add_adjusted(r.lo, xy.lo); |
| 288 if (spread + ilogbl(r.hi) > -16383) | 288 if (spread + ilogbl(r.hi) > -16383) |
| 289 return scalbnl(r.hi + adj, spread); | 289 return scalbnl(r.hi + adj, spread); |
| 290 else | 290 else |
| 291 return add_and_denormalize(r.hi, adj, spread); | 291 return add_and_denormalize(r.hi, adj, spread); |
| 292 } | 292 } |
| 293 #endif | 293 #endif |
| OLD | NEW |