1/**
 2 * This file has no copyright assigned and is placed in the Public Domain.
 3 * This file is part of the mingw-w64 runtime package.
 4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 5 */
 6long double frexpl(long double value, int* exp);
 7
 8#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
 9
10double frexp(double value, int* exp);
11
12/* On ARM `long double` is 64 bits. */
13long double frexpl(long double value, int* exp)
14{
15  return frexp(value, exp);
16}
17
18#elif defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__)
19
20#include <stdint.h>
21
22/* https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format */
23typedef union x87reg_ {
24  struct __attribute__((__packed__)) {
25    uint64_t f64;
26    uint16_t exp : 15;
27    uint16_t sgn :  1;
28  };
29  long double f;
30} x87reg;
31
32long double frexpl(long double value, int* exp)
33{
34  int n;
35  x87reg reg;
36  reg.f = value;
37  if(reg.exp == 0x7FFF) {
38    /* The value is an infinity or NaN.
39     * Store zero in `*exp`. Return the value as is. */
40    *exp = 0;
41    return reg.f;
42  }
43  if(reg.exp != 0) {
44    /* The value is normalized.
45     * Extract and zero out the exponent. */
46    *exp = reg.exp - 0x3FFE;
47    reg.exp = 0x3FFE;
48    return reg.f;
49  }
50  if(reg.f64 == 0) {
51    /* The value is zero.
52     * Store zero in `*exp`. Return the value as is.
53     * Note the signness. */
54    *exp = 0;
55    return reg.f;
56  }
57  /* The value is denormalized.
58   * Extract the exponent, normalize the value, then zero out
59   * the exponent. Note that x87 uses an explicit leading bit. */
60  n = __builtin_clzll(reg.f64);
61  reg.f64 <<= n;
62  *exp = 1 - 0x3FFE - n;
63  reg.exp = 0x3FFE;
64  return reg.f;
65}
66
67#else
68
69#error Please add `frexpl()` implementation for this platform.
70
71#endif