| 1 | /* |
|---|
| 2 | * Copyright © 2005-2014 Rich Felker, et al. |
|---|
| 3 | * |
|---|
| 4 | * Permission is hereby granted, free of charge, to any person obtaining |
|---|
| 5 | * a copy of this software and associated documentation files (the |
|---|
| 6 | * "Software"), to deal in the Software without restriction, including |
|---|
| 7 | * without limitation the rights to use, copy, modify, merge, publish, |
|---|
| 8 | * distribute, sublicense, and/or sell copies of the Software, and to |
|---|
| 9 | * permit persons to whom the Software is furnished to do so, subject to |
|---|
| 10 | * the following conditions: |
|---|
| 11 | * |
|---|
| 12 | * The above copyright notice and this permission notice shall be |
|---|
| 13 | * included in all copies or substantial portions of the Software. |
|---|
| 14 | * |
|---|
| 15 | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, |
|---|
| 16 | * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF |
|---|
| 17 | * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. |
|---|
| 18 | * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY |
|---|
| 19 | * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, |
|---|
| 20 | * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE |
|---|
| 21 | * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. |
|---|
| 22 | */ |
|---|
| 23 | |
|---|
| 24 | #include <math.h> |
|---|
| 25 | #include <inttypes.h> |
|---|
| 26 | #include "local.h" |
|---|
| 27 | |
|---|
| 28 | union fshape { |
|---|
| 29 | float value; |
|---|
| 30 | uint32_t bits; |
|---|
| 31 | }; |
|---|
| 32 | |
|---|
| 33 | // This is only necessary because the implementation of isnan only works |
|---|
| 34 | // properly when long double == double. |
|---|
| 35 | // See: https://sourceware.org/ml/newlib/2014/msg00684.html |
|---|
| 36 | #ifdef _LDBL_EQ_DOUBLE |
|---|
| 37 | |
|---|
| 38 | float |
|---|
| 39 | nexttowardf (float x, long double y) |
|---|
| 40 | { |
|---|
| 41 | union fshape ux; |
|---|
| 42 | uint32_t e; |
|---|
| 43 | |
|---|
| 44 | if (isnan(x) || isnan(y)) |
|---|
| 45 | return x + y; |
|---|
| 46 | if (x == y) |
|---|
| 47 | return y; |
|---|
| 48 | ux.value = x; |
|---|
| 49 | if (x == 0) { |
|---|
| 50 | ux.bits = 1; |
|---|
| 51 | if (signbit(y)) |
|---|
| 52 | ux.bits |= 0x80000000; |
|---|
| 53 | } else if (x < y) { |
|---|
| 54 | if (signbit(x)) |
|---|
| 55 | ux.bits--; |
|---|
| 56 | else |
|---|
| 57 | ux.bits++; |
|---|
| 58 | } else { |
|---|
| 59 | if (signbit(x)) |
|---|
| 60 | ux.bits++; |
|---|
| 61 | else |
|---|
| 62 | ux.bits--; |
|---|
| 63 | } |
|---|
| 64 | e = ux.bits & 0x7f800000; |
|---|
| 65 | /* raise overflow if ux.value is infinite and x is finite */ |
|---|
| 66 | if (e == 0x7f800000) { |
|---|
| 67 | volatile float force_eval; |
|---|
| 68 | force_eval = x + x; |
|---|
| 69 | } |
|---|
| 70 | /* raise underflow if ux.value is subnormal or zero */ |
|---|
| 71 | if (e == 0) { |
|---|
| 72 | volatile float force_eval; |
|---|
| 73 | force_eval = x*x + ux.value*ux.value; |
|---|
| 74 | } |
|---|
| 75 | return ux.value; |
|---|
| 76 | } |
|---|
| 77 | |
|---|
| 78 | #endif // _LDBL_EQ_DOUBLE |
|---|
