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
copysign(half x,half y)21 extern 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
frexp(half x,int * eptr)31 extern 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
ilogb(half x)58 extern 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
modf(half x,half * iptr)82 extern 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
nextafter(half x,half y)117 extern half __attribute__((overloadable)) nextafter(half x, half y) {
118   short hx, hy, ix, iy;
119 
120   GET_HALF_WORD(hx, x);
121   GET_HALF_WORD(hy, y);
122   ix = hx & 0x7fff; // |x|
123   iy = hy & 0x7fff; // |y|
124 
125   if ((ix > 0x7c00) || // x is nan
126       (iy > 0x7c00))   // y is nan
127     return x + y;      // return nan
128 
129   if (x == y) return y; // x == y.  return y
130   if (ix == 0) {
131     SET_HALF_WORD(x, (hy & 0x8000) | 1);
132     return x;
133   }
134 
135   if (hx >= 0) {  // x >= 0
136     if (hx > hy)
137       hx -= 1;    // x > y, x -= 1 ulp
138     else
139       hx += 1;    // x < y, x += 1 ulp
140   }
141   else {          // x < 0
142     if (hy>= 0 || hx > hy)
143       hx -= 1;    // x < y, x -= 1 ulp
144     else
145       hx += 1;
146   }
147 
148   SET_HALF_WORD(x, hx);
149   return x;
150 }
151