blob: db65acf97b285a3fe5b44a7e935f6577897f125f (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
//===-- lib/extendsfdf2.c - single -> double conversion -----------*- C -*-===//
//
// The LLVM Compiler Infrastructure
//
// This file is dual licensed under the MIT and the University of Illinois Open
// Source Licenses. See LICENSE.TXT for details.
//
//===----------------------------------------------------------------------===//
//
// This file implements a fairly generic conversion from a narrower to a wider
// IEEE-754 floating-point type. The constants and types defined following the
// includes below parameterize the conversion.
//
// This routine can be trivially adapted to support conversions from
// half-precision or to quad-precision. It does not support types that don't
// use the usual IEEE-754 interchange formats; specifically, some work would be
// needed to adapt it to (for example) the Intel 80-bit format or PowerPC
// double-double format.
//
// Note please, however, that this implementation is only intended to support
// *widening* operations; if you need to convert to a *narrower* floating-point
// type (e.g. double -> float), then this routine will not do what you want it
// to.
//
// It also requires that integer types at least as large as both formats
// are available on the target platform; this may pose a problem when trying
// to add support for quad on some 32-bit systems, for example. You also may
// run into trouble finding an appropriate CLZ function for wide source types;
// you will likely need to roll your own on some platforms.
//
// Finally, the following assumptions are made:
//
// 1. floating-point types and integer types have the same endianness on the
// target platform
//
// 2. quiet NaNs, if supported, are indicated by the leading bit of the
// significand field being set
//
//===----------------------------------------------------------------------===//
#include <stdint.h>
#include <limits.h>
typedef float src_t;
typedef uint32_t src_rep_t;
#define SRC_REP_C UINT32_C
static const int srcSigBits = 23;
#define src_rep_t_clz __builtin_clz
typedef double dst_t;
typedef uint64_t dst_rep_t;
#define DST_REP_C UINT64_C
static const int dstSigBits = 52;
// End of specialization parameters. Two helper routines for conversion to and
// from the representation of floating-point data as integer values follow.
static inline src_rep_t srcToRep(src_t x) {
const union { src_t f; src_rep_t i; } rep = {.f = x};
return rep.i;
}
static inline dst_t dstFromRep(dst_rep_t x) {
const union { dst_t f; dst_rep_t i; } rep = {.i = x};
return rep.f;
}
// End helper routines. Conversion implementation follows.
dst_t __extendsfdf2(src_t a) {
// Various constants whose values follow from the type parameters.
// Any reasonable optimizer will fold and propagate all of these.
const int srcBits = sizeof(src_t)*CHAR_BIT;
const int srcExpBits = srcBits - srcSigBits - 1;
const int srcInfExp = (1 << srcExpBits) - 1;
const int srcExpBias = srcInfExp >> 1;
const src_rep_t srcMinNormal = SRC_REP_C(1) << srcSigBits;
const src_rep_t srcInfinity = (src_rep_t)srcInfExp << srcSigBits;
const src_rep_t srcSignMask = SRC_REP_C(1) << (srcSigBits + srcExpBits);
const src_rep_t srcAbsMask = srcSignMask - 1;
const src_rep_t srcQNaN = SRC_REP_C(1) << (srcSigBits - 1);
const src_rep_t srcNaNCode = srcQNaN - 1;
const int dstBits = sizeof(dst_t)*CHAR_BIT;
const int dstExpBits = dstBits - dstSigBits - 1;
const int dstInfExp = (1 << dstExpBits) - 1;
const int dstExpBias = dstInfExp >> 1;
const dst_rep_t dstMinNormal = DST_REP_C(1) << dstSigBits;
// Break a into a sign and representation of the absolute value
const src_rep_t aRep = srcToRep(a);
const src_rep_t aAbs = aRep & srcAbsMask;
const src_rep_t sign = aRep & srcSignMask;
dst_rep_t absResult;
if (aAbs - srcMinNormal < srcInfinity - srcMinNormal) {
// a is a normal number.
// Extend to the destination type by shifting the significand and
// exponent into the proper position and rebiasing the exponent.
absResult = (dst_rep_t)aAbs << (dstSigBits - srcSigBits);
absResult += (dst_rep_t)(dstExpBias - srcExpBias) << dstSigBits;
}
else if (aAbs >= srcInfinity) {
// a is NaN or infinity.
// Conjure the result by beginning with infinity, then setting the qNaN
// bit (if needed) and right-aligning the rest of the trailing NaN
// payload field.
absResult = (dst_rep_t)dstInfExp << dstSigBits;
absResult |= (dst_rep_t)(aAbs & srcQNaN) << (dstSigBits - srcSigBits);
absResult |= aAbs & srcNaNCode;
}
else if (aAbs) {
// a is denormal.
// renormalize the significand and clear the leading bit, then insert
// the correct adjusted exponent in the destination type.
const int scale = src_rep_t_clz(aAbs) - src_rep_t_clz(srcMinNormal);
absResult = (dst_rep_t)aAbs << (dstSigBits - srcSigBits + scale);
absResult ^= dstMinNormal;
const int resultExponent = dstExpBias - srcExpBias - scale + 1;
absResult |= (dst_rep_t)resultExponent << dstSigBits;
}
else {
// a is zero.
absResult = 0;
}
// Apply the signbit to (dst_t)abs(a).
const dst_rep_t result = absResult | (dst_rep_t)sign << (dstBits - srcBits);
return dstFromRep(result);
}
|
