rs_f16_math.c revision a53bf64187bcf01ed1c62eacd52f6c6f41b5864b 
1/* Implementations for copysign, ilogb, and nextafter for float16 based on
2 * corresponding float32 implementations in
3 * bionic/libm/upstream-freebsd/lib/msun/src
4 */
5
6/*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 *
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
15 */
16
17#include "rs_core.rsh"
18#include "rs_f16_util.h"
19
20// Based on bionic/libm/upstream-freebsd/lib/msun/src/s_copysignf.c
21extern half __attribute__((overloadable)) copysign(half x, half y) {
22    short hx, hy;
23    GET_HALF_WORD(hx, x);
24    GET_HALF_WORD(hy, y);
25
26    SET_HALF_WORD(x, (hx & 0x7fff) | (hy & 0x8000));
27    return x;
28}
29
30// Based on bionic/libm/upstream-freebsd/lib/msun/src/s_frexpf.c
31extern half __attribute__((overloadable)) frexp(half x, int *eptr) {
32    short hx, ix;
33    static const half two12 = 4096;
34
35    GET_HALF_WORD(hx, x);
36    ix = hx & 0x7fff;
37
38    *eptr = 0;
39    if (ix >= 0x7c00 || ix == 0) return x; // NaN, infinity or zero
40    if (ix <= 0x0400) {
41        // x is subnormal.  Scale it by 2^12 (and adjust eptr accordingly) so
42        // that even the smallest subnormal value becomes normal.
43        x *= two12;
44        GET_HALF_WORD(hx, x);
45        ix = hx & 0x7fff;
46        *eptr = -12;
47    }
48
49    // Adjust eptr by (non-biased exponent of hx + 1).  Set the non-biased
50    // exponent to be equal to -1 so that abs(hx) is between 0.5 and 1.
51    *eptr += (ix >> 10) - 14;
52    hx = (hx & 0x83ff) | 0x3800;
53    SET_HALF_WORD(x, hx);
54    return x;
55}
56
57// Based on bionic/libm/upstream-freebsd/lib/msun/src/s_ilogbf.c
58extern int __attribute__((overloadable)) ilogb(half x) {
59    const int RS_INT_MAX = 0x7fffffff;
60    const int RS_INT_MIN = 0x80000000;
61
62    short hx, ix;
63    GET_HALF_WORD(hx, x);
64    hx &= 0x7fff;
65
66    if (hx < 0x0400) { // subnormal
67        if (hx == 0)
68            return RS_INT_MIN; // for zero
69        for (hx <<= 5, ix = -14; hx > 0; hx <<= 1)
70            ix -= 1;
71        return ix;
72    }
73    else if (hx < 0x7c00) {
74        return (hx >> 10) - 15;
75    }
76    else { // hx >= 0x7c00
77        return RS_INT_MAX; // for NaN and infinity
78    }
79}
80
81// Based on bionic/libm/upstream-freebsd/lib/msun/src/s_modff.c
82extern half __attribute__((overloadable)) modf(half x, half *iptr) {
83    short i0, j0;
84    unsigned short i;
85    GET_HALF_WORD(i0, x);
86    j0 = ((i0 >> 10) & 0x1f) - 15; // exponent of x
87    if (j0 < 10) {
88        if (j0 < 0) { // No integral part
89            SET_HALF_WORD(*iptr, i0 & 0x8000); // *iptr = +/- 0
90            return x;
91        }
92        else {
93            i = 0x03ff >> j0; // mask to check fractional parts of x
94            if ((i0 & i) == 0) { // no bits set in fractional part
95                *iptr = x;
96                SET_HALF_WORD(x, i0 & 0x8000);
97                return x;
98            }
99            else {
100                SET_HALF_WORD(*iptr, i0 & ~i); // zero out fractional parts
101                return x - *iptr;
102            }
103        }
104    }
105    else { // No fractional part
106        unsigned short ix;
107        *iptr = x;
108        if (x != x)
109            return x;
110        GET_HALF_WORD(ix, x);
111        SET_HALF_WORD(x, ix & 0x8000); // x = +/- 0
112        return x;
113    }
114}
115
116// Based on bionic/libm/upstream-freebsd/lib/msun/src/s_nextafterf.c
117extern half __attribute__((overloadable)) nextafter(half x, half y) {
118  volatile half t;
119  short hx, hy, ix, iy;
120
121  GET_HALF_WORD(hx, x);
122  GET_HALF_WORD(hy, y);
123  ix = hx & 0x7fff; // |x|
124  iy = hy & 0x7fff; // |y|
125
126  if ((ix > 0x7c00) || // x is nan
127      (iy > 0x7c00))   // y is nan
128    return x + y;      // return nan
129
130  if (x == y) return y; // x == y.  return y
131  if (ix == 0) {
132    SET_HALF_WORD(x, (hy & 0x8000) | 1);
133    return x;
134  }
135
136  if (hx >= 0) {  // x >= 0
137    if (hx > hy)
138      hx -= 1;    // x > y, x -= 1 ulp
139    else
140      hx += 1;    // x < y, x += 1 ulp
141  }
142  else {          // x < 0
143    if (hy>= 0 || hx > hy)
144      hx -= 1;    // x < y, x -= 1 ulp
145    else
146      hx += 1;
147  }
148
149  SET_HALF_WORD(x, hx);
150  return x;
151}
152