1
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
2
#include "numpy/halffloat.h"
3

4
/*
5
 * This chooses between 'ties to even' and 'ties away from zero'.
6
 */
7
#define NPY_HALF_ROUND_TIES_TO_EVEN 1
8
/*
9
 * If these are 1, the conversions try to trigger underflow,
10
 * overflow, and invalid exceptions in the FP system when needed.
11
 */
12
#define NPY_HALF_GENERATE_OVERFLOW 1
13
#define NPY_HALF_GENERATE_UNDERFLOW 1
14
#define NPY_HALF_GENERATE_INVALID 1
15

16
/*
17
 ********************************************************************
18
 *                   HALF-PRECISION ROUTINES                        *
19
 ********************************************************************
20
 */
21

22 1
float npy_half_to_float(npy_half h)
23
{
24
    union { float ret; npy_uint32 retbits; } conv;
25 1
    conv.retbits = npy_halfbits_to_floatbits(h);
26 1
    return conv.ret;
27
}
28

29 1
double npy_half_to_double(npy_half h)
30
{
31
    union { double ret; npy_uint64 retbits; } conv;
32 1
    conv.retbits = npy_halfbits_to_doublebits(h);
33 1
    return conv.ret;
34
}
35

36 1
npy_half npy_float_to_half(float f)
37
{
38
    union { float f; npy_uint32 fbits; } conv;
39 1
    conv.f = f;
40 1
    return npy_floatbits_to_halfbits(conv.fbits);
41
}
42

43 1
npy_half npy_double_to_half(double d)
44
{
45
    union { double d; npy_uint64 dbits; } conv;
46 1
    conv.d = d;
47 1
    return npy_doublebits_to_halfbits(conv.dbits);
48
}
49

50 1
int npy_half_iszero(npy_half h)
51
{
52 1
    return (h&0x7fff) == 0;
53
}
54

55 1
int npy_half_isnan(npy_half h)
56
{
57 1
    return ((h&0x7c00u) == 0x7c00u) && ((h&0x03ffu) != 0x0000u);
58
}
59

60 1
int npy_half_isinf(npy_half h)
61
{
62 1
    return ((h&0x7fffu) == 0x7c00u);
63
}
64

65 1
int npy_half_isfinite(npy_half h)
66
{
67 1
    return ((h&0x7c00u) != 0x7c00u);
68
}
69

70 1
int npy_half_signbit(npy_half h)
71
{
72 1
    return (h&0x8000u) != 0;
73
}
74

75 1
npy_half npy_half_spacing(npy_half h)
76
{
77
    npy_half ret;
78 1
    npy_uint16 h_exp = h&0x7c00u;
79 1
    npy_uint16 h_sig = h&0x03ffu;
80 1
    if (h_exp == 0x7c00u) {
81
#if NPY_HALF_GENERATE_INVALID
82 1
        npy_set_floatstatus_invalid();
83
#endif
84 1
        ret = NPY_HALF_NAN;
85 1
    } else if (h == 0x7bffu) {
86
#if NPY_HALF_GENERATE_OVERFLOW
87 1
        npy_set_floatstatus_overflow();
88
#endif
89 1
        ret = NPY_HALF_PINF;
90 1
    } else if ((h&0x8000u) && h_sig == 0) { /* Negative boundary case */
91 1
        if (h_exp > 0x2c00u) { /* If result is normalized */
92 1
            ret = h_exp - 0x2c00u;
93 1
        } else if(h_exp > 0x0400u) { /* The result is a subnormal, but not the smallest */
94 1
            ret = 1 << ((h_exp >> 10) - 2);
95
        } else {
96
            ret = 0x0001u; /* Smallest subnormal half */
97
        }
98 1
    } else if (h_exp > 0x2800u) { /* If result is still normalized */
99 1
        ret = h_exp - 0x2800u;
100 1
    } else if (h_exp > 0x0400u) { /* The result is a subnormal, but not the smallest */
101 1
        ret = 1 << ((h_exp >> 10) - 1);
102
    } else {
103
        ret = 0x0001u;
104
    }
105

106 1
    return ret;
107
}
108

109 1
npy_half npy_half_copysign(npy_half x, npy_half y)
110
{
111 1
    return (x&0x7fffu) | (y&0x8000u);
112
}
113

114 1
npy_half npy_half_nextafter(npy_half x, npy_half y)
115
{
116
    npy_half ret;
117

118 1
    if (npy_half_isnan(x) || npy_half_isnan(y)) {
119
        ret = NPY_HALF_NAN;
120 1
    } else if (npy_half_eq_nonan(x, y)) {
121
        ret = x;
122 1
    } else if (npy_half_iszero(x)) {
123 1
        ret = (y&0x8000u) + 1; /* Smallest subnormal half */
124 1
    } else if (!(x&0x8000u)) { /* x > 0 */
125 1
        if ((npy_int16)x > (npy_int16)y) { /* x > y */
126 1
            ret = x-1;
127
        } else {
128 1
            ret = x+1;
129
        }
130
    } else {
131 1
        if (!(y&0x8000u) || (x&0x7fffu) > (y&0x7fffu)) { /* x < y */
132 1
            ret = x-1;
133
        } else {
134 1
            ret = x+1;
135
        }
136
    }
137
#if NPY_HALF_GENERATE_OVERFLOW
138 1
    if (npy_half_isinf(ret) && npy_half_isfinite(x)) {
139 1
        npy_set_floatstatus_overflow();
140
    }
141
#endif
142

143 1
    return ret;
144
}
145

146 1
int npy_half_eq_nonan(npy_half h1, npy_half h2)
147
{
148 1
    return (h1 == h2 || ((h1 | h2) & 0x7fff) == 0);
149
}
150

151 1
int npy_half_eq(npy_half h1, npy_half h2)
152
{
153
    /*
154
     * The equality cases are as follows:
155
     *   - If either value is NaN, never equal.
156
     *   - If the values are equal, equal.
157
     *   - If the values are both signed zeros, equal.
158
     */
159 1
    return (!npy_half_isnan(h1) && !npy_half_isnan(h2)) &&
160 1
           (h1 == h2 || ((h1 | h2) & 0x7fff) == 0);
161
}
162

163 1
int npy_half_ne(npy_half h1, npy_half h2)
164
{
165 1
    return !npy_half_eq(h1, h2);
166
}
167

168 1
int npy_half_lt_nonan(npy_half h1, npy_half h2)
169
{
170 1
    if (h1&0x8000u) {
171 1
        if (h2&0x8000u) {
172 1
            return (h1&0x7fffu) > (h2&0x7fffu);
173
        } else {
174
            /* Signed zeros are equal, have to check for it */
175 1
            return (h1 != 0x8000u) || (h2 != 0x0000u);
176
        }
177
    } else {
178 1
        if (h2&0x8000u) {
179
            return 0;
180
        } else {
181 1
            return (h1&0x7fffu) < (h2&0x7fffu);
182
        }
183
    }
184
}
185

186 1
int npy_half_lt(npy_half h1, npy_half h2)
187
{
188 1
    return (!npy_half_isnan(h1) && !npy_half_isnan(h2)) && npy_half_lt_nonan(h1, h2);
189
}
190

191 1
int npy_half_gt(npy_half h1, npy_half h2)
192
{
193 1
    return npy_half_lt(h2, h1);
194
}
195

196 1
int npy_half_le_nonan(npy_half h1, npy_half h2)
197
{
198 1
    if (h1&0x8000u) {
199 1
        if (h2&0x8000u) {
200 1
            return (h1&0x7fffu) >= (h2&0x7fffu);
201
        } else {
202
            return 1;
203
        }
204
    } else {
205 1
        if (h2&0x8000u) {
206
            /* Signed zeros are equal, have to check for it */
207 1
            return (h1 == 0x0000u) && (h2 == 0x8000u);
208
        } else {
209 1
            return (h1&0x7fffu) <= (h2&0x7fffu);
210
        }
211
    }
212
}
213

214 1
int npy_half_le(npy_half h1, npy_half h2)
215
{
216 1
    return (!npy_half_isnan(h1) && !npy_half_isnan(h2)) && npy_half_le_nonan(h1, h2);
217
}
218

219 1
int npy_half_ge(npy_half h1, npy_half h2)
220
{
221 1
    return npy_half_le(h2, h1);
222
}
223

224 1
npy_half npy_half_divmod(npy_half h1, npy_half h2, npy_half *modulus)
225
{
226 1
    float fh1 = npy_half_to_float(h1);
227 1
    float fh2 = npy_half_to_float(h2);
228
    float div, mod;
229

230 1
    div = npy_divmodf(fh1, fh2, &mod);
231 1
    *modulus = npy_float_to_half(mod);
232 1
    return npy_float_to_half(div);
233
}
234

235

236

237
/*
238
 ********************************************************************
239
 *                     BIT-LEVEL CONVERSIONS                        *
240
 ********************************************************************
241
 */
242

243 1
npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
244
{
245
    npy_uint32 f_exp, f_sig;
246
    npy_uint16 h_sgn, h_exp, h_sig;
247

248 1
    h_sgn = (npy_uint16) ((f&0x80000000u) >> 16);
249 1
    f_exp = (f&0x7f800000u);
250

251
    /* Exponent overflow/NaN converts to signed inf/NaN */
252 1
    if (f_exp >= 0x47800000u) {
253 1
        if (f_exp == 0x7f800000u) {
254
            /* Inf or NaN */
255 1
            f_sig = (f&0x007fffffu);
256 1
            if (f_sig != 0) {
257
                /* NaN - propagate the flag in the significand... */
258 1
                npy_uint16 ret = (npy_uint16) (0x7c00u + (f_sig >> 13));
259
                /* ...but make sure it stays a NaN */
260 1
                if (ret == 0x7c00u) {
261 0
                    ret++;
262
                }
263 1
                return h_sgn + ret;
264
            } else {
265
                /* signed inf */
266 1
                return (npy_uint16) (h_sgn + 0x7c00u);
267
            }
268
        } else {
269
            /* overflow to signed inf */
270
#if NPY_HALF_GENERATE_OVERFLOW
271 1
            npy_set_floatstatus_overflow();
272
#endif
273 1
            return (npy_uint16) (h_sgn + 0x7c00u);
274
        }
275
    }
276

277
    /* Exponent underflow converts to a subnormal half or signed zero */
278 1
    if (f_exp <= 0x38000000u) {
279
        /*
280
         * Signed zeros, subnormal floats, and floats with small
281
         * exponents all convert to signed zero half-floats.
282
         */
283 1
        if (f_exp < 0x33000000u) {
284
#if NPY_HALF_GENERATE_UNDERFLOW
285
            /* If f != 0, it underflowed to 0 */
286 1
            if ((f&0x7fffffff) != 0) {
287 1
                npy_set_floatstatus_underflow();
288
            }
289
#endif
290
            return h_sgn;
291
        }
292
        /* Make the subnormal significand */
293 1
        f_exp >>= 23;
294 1
        f_sig = (0x00800000u + (f&0x007fffffu));
295
#if NPY_HALF_GENERATE_UNDERFLOW
296
        /* If it's not exactly represented, it underflowed */
297 1
        if ((f_sig&(((npy_uint32)1 << (126 - f_exp)) - 1)) != 0) {
298 1
            npy_set_floatstatus_underflow();
299
        }
300
#endif
301
        /*
302
         * Usually the significand is shifted by 13. For subnormals an
303
         * additional shift needs to occur. This shift is one for the largest
304
         * exponent giving a subnormal `f_exp = 0x38000000 >> 23 = 112`, which
305
         * offsets the new first bit. At most the shift can be 1+10 bits.
306
         */
307 1
        f_sig >>= (113 - f_exp);
308
        /* Handle rounding by adding 1 to the bit beyond half precision */
309
#if NPY_HALF_ROUND_TIES_TO_EVEN
310
        /*
311
         * If the last bit in the half significand is 0 (already even), and
312
         * the remaining bit pattern is 1000...0, then we do not add one
313
         * to the bit after the half significand. However, the (113 - f_exp)
314
         * shift can lose up to 11 bits, so the || checks them in the original.
315
         * In all other cases, we can just add one.
316
         */
317 1
        if (((f_sig&0x00003fffu) != 0x00001000u) || (f&0x000007ffu)) {
318 1
            f_sig += 0x00001000u;
319
        }
320
#else
321
        f_sig += 0x00001000u;
322
#endif
323 1
        h_sig = (npy_uint16) (f_sig >> 13);
324
        /*
325
         * If the rounding causes a bit to spill into h_exp, it will
326
         * increment h_exp from zero to one and h_sig will be zero.
327
         * This is the correct result.
328
         */
329 1
        return (npy_uint16) (h_sgn + h_sig);
330
    }
331

332
    /* Regular case with no overflow or underflow */
333 1
    h_exp = (npy_uint16) ((f_exp - 0x38000000u) >> 13);
334
    /* Handle rounding by adding 1 to the bit beyond half precision */
335 1
    f_sig = (f&0x007fffffu);
336
#if NPY_HALF_ROUND_TIES_TO_EVEN
337
    /*
338
     * If the last bit in the half significand is 0 (already even), and
339
     * the remaining bit pattern is 1000...0, then we do not add one
340
     * to the bit after the half significand.  In all other cases, we do.
341
     */
342 1
    if ((f_sig&0x00003fffu) != 0x00001000u) {
343 1
        f_sig += 0x00001000u;
344
    }
345
#else
346
    f_sig += 0x00001000u;
347
#endif
348 1
    h_sig = (npy_uint16) (f_sig >> 13);
349
    /*
350
     * If the rounding causes a bit to spill into h_exp, it will
351
     * increment h_exp by one and h_sig will be zero.  This is the
352
     * correct result.  h_exp may increment to 15, at greatest, in
353
     * which case the result overflows to a signed inf.
354
     */
355
#if NPY_HALF_GENERATE_OVERFLOW
356 1
    h_sig += h_exp;
357 1
    if (h_sig == 0x7c00u) {
358 1
        npy_set_floatstatus_overflow();
359
    }
360 1
    return h_sgn + h_sig;
361
#else
362
    return h_sgn + h_exp + h_sig;
363
#endif
364
}
365

366 1
npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
367
{
368
    npy_uint64 d_exp, d_sig;
369
    npy_uint16 h_sgn, h_exp, h_sig;
370

371 1
    h_sgn = (d&0x8000000000000000ULL) >> 48;
372 1
    d_exp = (d&0x7ff0000000000000ULL);
373

374
    /* Exponent overflow/NaN converts to signed inf/NaN */
375 1
    if (d_exp >= 0x40f0000000000000ULL) {
376 1
        if (d_exp == 0x7ff0000000000000ULL) {
377
            /* Inf or NaN */
378 1
            d_sig = (d&0x000fffffffffffffULL);
379 1
            if (d_sig != 0) {
380
                /* NaN - propagate the flag in the significand... */
381 1
                npy_uint16 ret = (npy_uint16) (0x7c00u + (d_sig >> 42));
382
                /* ...but make sure it stays a NaN */
383 1
                if (ret == 0x7c00u) {
384 0
                    ret++;
385
                }
386 1
                return h_sgn + ret;
387
            } else {
388
                /* signed inf */
389 1
                return h_sgn + 0x7c00u;
390
            }
391
        } else {
392
            /* overflow to signed inf */
393
#if NPY_HALF_GENERATE_OVERFLOW
394 1
            npy_set_floatstatus_overflow();
395
#endif
396 1
            return h_sgn + 0x7c00u;
397
        }
398
    }
399

400
    /* Exponent underflow converts to subnormal half or signed zero */
401 1
    if (d_exp <= 0x3f00000000000000ULL) {
402
        /*
403
         * Signed zeros, subnormal floats, and floats with small
404
         * exponents all convert to signed zero half-floats.
405
         */
406 1
        if (d_exp < 0x3e60000000000000ULL) {
407
#if NPY_HALF_GENERATE_UNDERFLOW
408
            /* If d != 0, it underflowed to 0 */
409 1
            if ((d&0x7fffffffffffffffULL) != 0) {
410 1
                npy_set_floatstatus_underflow();
411
            }
412
#endif
413
            return h_sgn;
414
        }
415
        /* Make the subnormal significand */
416 1
        d_exp >>= 52;
417 1
        d_sig = (0x0010000000000000ULL + (d&0x000fffffffffffffULL));
418
#if NPY_HALF_GENERATE_UNDERFLOW
419
        /* If it's not exactly represented, it underflowed */
420 1
        if ((d_sig&(((npy_uint64)1 << (1051 - d_exp)) - 1)) != 0) {
421 1
            npy_set_floatstatus_underflow();
422
        }
423
#endif
424
        /*
425
         * Unlike floats, doubles have enough room to shift left to align
426
         * the subnormal significand leading to no loss of the last bits.
427
         * The smallest possible exponent giving a subnormal is:
428
         * `d_exp = 0x3e60000000000000 >> 52 = 998`. All larger subnormals are
429
         * shifted with respect to it. This adds a shift of 10+1 bits the final
430
         * right shift when comparing it to the one in the normal branch.
431
         */
432
        assert(d_exp - 998 >= 0);
433 1
        d_sig <<= (d_exp - 998);
434
        /* Handle rounding by adding 1 to the bit beyond half precision */
435
#if NPY_HALF_ROUND_TIES_TO_EVEN
436
        /*
437
         * If the last bit in the half significand is 0 (already even), and
438
         * the remaining bit pattern is 1000...0, then we do not add one
439
         * to the bit after the half significand.  In all other cases, we do.
440
         */
441 1
        if ((d_sig&0x003fffffffffffffULL) != 0x0010000000000000ULL) {
442 1
            d_sig += 0x0010000000000000ULL;
443
        }
444
#else
445
        d_sig += 0x0010000000000000ULL;
446
#endif
447 1
        h_sig = (npy_uint16) (d_sig >> 53);
448
        /*
449
         * If the rounding causes a bit to spill into h_exp, it will
450
         * increment h_exp from zero to one and h_sig will be zero.
451
         * This is the correct result.
452
         */
453 1
        return h_sgn + h_sig;
454
    }
455

456
    /* Regular case with no overflow or underflow */
457 1
    h_exp = (npy_uint16) ((d_exp - 0x3f00000000000000ULL) >> 42);
458
    /* Handle rounding by adding 1 to the bit beyond half precision */
459 1
    d_sig = (d&0x000fffffffffffffULL);
460
#if NPY_HALF_ROUND_TIES_TO_EVEN
461
    /*
462
     * If the last bit in the half significand is 0 (already even), and
463
     * the remaining bit pattern is 1000...0, then we do not add one
464
     * to the bit after the half significand.  In all other cases, we do.
465
     */
466 1
    if ((d_sig&0x000007ffffffffffULL) != 0x0000020000000000ULL) {
467 1
        d_sig += 0x0000020000000000ULL;
468
    }
469
#else
470
    d_sig += 0x0000020000000000ULL;
471
#endif
472 1
    h_sig = (npy_uint16) (d_sig >> 42);
473

474
    /*
475
     * If the rounding causes a bit to spill into h_exp, it will
476
     * increment h_exp by one and h_sig will be zero.  This is the
477
     * correct result.  h_exp may increment to 15, at greatest, in
478
     * which case the result overflows to a signed inf.
479
     */
480
#if NPY_HALF_GENERATE_OVERFLOW
481 1
    h_sig += h_exp;
482 1
    if (h_sig == 0x7c00u) {
483 1
        npy_set_floatstatus_overflow();
484
    }
485 1
    return h_sgn + h_sig;
486
#else
487
    return h_sgn + h_exp + h_sig;
488
#endif
489
}
490

491 1
npy_uint32 npy_halfbits_to_floatbits(npy_uint16 h)
492
{
493
    npy_uint16 h_exp, h_sig;
494
    npy_uint32 f_sgn, f_exp, f_sig;
495

496 1
    h_exp = (h&0x7c00u);
497 1
    f_sgn = ((npy_uint32)h&0x8000u) << 16;
498 1
    switch (h_exp) {
499 1
        case 0x0000u: /* 0 or subnormal */
500 1
            h_sig = (h&0x03ffu);
501
            /* Signed zero */
502 1
            if (h_sig == 0) {
503
                return f_sgn;
504
            }
505
            /* Subnormal */
506 1
            h_sig <<= 1;
507 1
            while ((h_sig&0x0400u) == 0) {
508 1
                h_sig <<= 1;
509 1
                h_exp++;
510
            }
511 1
            f_exp = ((npy_uint32)(127 - 15 - h_exp)) << 23;
512 1
            f_sig = ((npy_uint32)(h_sig&0x03ffu)) << 13;
513 1
            return f_sgn + f_exp + f_sig;
514 1
        case 0x7c00u: /* inf or NaN */
515
            /* All-ones exponent and a copy of the significand */
516 1
            return f_sgn + 0x7f800000u + (((npy_uint32)(h&0x03ffu)) << 13);
517 1
        default: /* normalized */
518
            /* Just need to adjust the exponent and shift */
519 1
            return f_sgn + (((npy_uint32)(h&0x7fffu) + 0x1c000u) << 13);
520
    }
521
}
522

523 1
npy_uint64 npy_halfbits_to_doublebits(npy_uint16 h)
524
{
525
    npy_uint16 h_exp, h_sig;
526
    npy_uint64 d_sgn, d_exp, d_sig;
527

528 1
    h_exp = (h&0x7c00u);
529 1
    d_sgn = ((npy_uint64)h&0x8000u) << 48;
530 1
    switch (h_exp) {
531 1
        case 0x0000u: /* 0 or subnormal */
532 1
            h_sig = (h&0x03ffu);
533
            /* Signed zero */
534 1
            if (h_sig == 0) {
535
                return d_sgn;
536
            }
537
            /* Subnormal */
538 1
            h_sig <<= 1;
539 1
            while ((h_sig&0x0400u) == 0) {
540 1
                h_sig <<= 1;
541 1
                h_exp++;
542
            }
543 1
            d_exp = ((npy_uint64)(1023 - 15 - h_exp)) << 52;
544 1
            d_sig = ((npy_uint64)(h_sig&0x03ffu)) << 42;
545 1
            return d_sgn + d_exp + d_sig;
546 1
        case 0x7c00u: /* inf or NaN */
547
            /* All-ones exponent and a copy of the significand */
548 1
            return d_sgn + 0x7ff0000000000000ULL +
549 1
                                (((npy_uint64)(h&0x03ffu)) << 42);
550 1
        default: /* normalized */
551
            /* Just need to adjust the exponent and shift */
552 1
            return d_sgn + (((npy_uint64)(h&0x7fffu) + 0xfc000u) << 42);
553
    }
554
}

Read our documentation on viewing source code .

Loading