Branch data Line data Source code
1 : : // Copyright 2019-2026 David Robillard <d@drobilla.net>
2 : : // SPDX-License-Identifier: ISC
3 : :
4 : : #include "soft_float.h"
5 : :
6 : : #include "attributes.h"
7 : : #include "ieee_float.h"
8 : : #include "int_math.h"
9 : :
10 : : #include <assert.h>
11 : : #include <math.h>
12 : : #include <stdint.h>
13 : :
14 : : /// 10^k for k = min_dec_expt, min_dec_expt + dec_expt_step, ..., max_dec_expt
15 : : static const ExessSoftFloat soft_pow10[] = {
16 : : {0xFA8FD5A0081C0288, -1220}, {0xBAAEE17FA23EBF76, -1193},
17 : : {0x8B16FB203055AC76, -1166}, {0xCF42894A5DCE35EA, -1140},
18 : : {0x9A6BB0AA55653B2D, -1113}, {0xE61ACF033D1A45DF, -1087},
19 : : {0xAB70FE17C79AC6CA, -1060}, {0xFF77B1FCBEBCDC4F, -1034},
20 : : {0xBE5691EF416BD60C, -1007}, {0x8DD01FAD907FFC3C, -980},
21 : : {0xD3515C2831559A83, -954}, {0x9D71AC8FADA6C9B5, -927},
22 : : {0xEA9C227723EE8BCB, -901}, {0xAECC49914078536D, -874},
23 : : {0x823C12795DB6CE57, -847}, {0xC21094364DFB5637, -821},
24 : : {0x9096EA6F3848984F, -794}, {0xD77485CB25823AC7, -768},
25 : : {0xA086CFCD97BF97F4, -741}, {0xEF340A98172AACE5, -715},
26 : : {0xB23867FB2A35B28E, -688}, {0x84C8D4DFD2C63F3B, -661},
27 : : {0xC5DD44271AD3CDBA, -635}, {0x936B9FCEBB25C996, -608},
28 : : {0xDBAC6C247D62A584, -582}, {0xA3AB66580D5FDAF6, -555},
29 : : {0xF3E2F893DEC3F126, -529}, {0xB5B5ADA8AAFF80B8, -502},
30 : : {0x87625F056C7C4A8B, -475}, {0xC9BCFF6034C13053, -449},
31 : : {0x964E858C91BA2655, -422}, {0xDFF9772470297EBD, -396},
32 : : {0xA6DFBD9FB8E5B88F, -369}, {0xF8A95FCF88747D94, -343},
33 : : {0xB94470938FA89BCF, -316}, {0x8A08F0F8BF0F156B, -289},
34 : : {0xCDB02555653131B6, -263}, {0x993FE2C6D07B7FAC, -236},
35 : : {0xE45C10C42A2B3B06, -210}, {0xAA242499697392D3, -183},
36 : : {0xFD87B5F28300CA0E, -157}, {0xBCE5086492111AEB, -130},
37 : : {0x8CBCCC096F5088CC, -103}, {0xD1B71758E219652C, -77},
38 : : {0x9C40000000000000, -50}, {0xE8D4A51000000000, -24},
39 : : {0xAD78EBC5AC620000, 3}, {0x813F3978F8940984, 30},
40 : : {0xC097CE7BC90715B3, 56}, {0x8F7E32CE7BEA5C70, 83},
41 : : {0xD5D238A4ABE98068, 109}, {0x9F4F2726179A2245, 136},
42 : : {0xED63A231D4C4FB27, 162}, {0xB0DE65388CC8ADA8, 189},
43 : : {0x83C7088E1AAB65DB, 216}, {0xC45D1DF942711D9A, 242},
44 : : {0x924D692CA61BE758, 269}, {0xDA01EE641A708DEA, 295},
45 : : {0xA26DA3999AEF774A, 322}, {0xF209787BB47D6B85, 348},
46 : : {0xB454E4A179DD1877, 375}, {0x865B86925B9BC5C2, 402},
47 : : {0xC83553C5C8965D3D, 428}, {0x952AB45CFA97A0B3, 455},
48 : : {0xDE469FBD99A05FE3, 481}, {0xA59BC234DB398C25, 508},
49 : : {0xF6C69A72A3989F5C, 534}, {0xB7DCBF5354E9BECE, 561},
50 : : {0x88FCF317F22241E2, 588}, {0xCC20CE9BD35C78A5, 614},
51 : : {0x98165AF37B2153DF, 641}, {0xE2A0B5DC971F303A, 667},
52 : : {0xA8D9D1535CE3B396, 694}, {0xFB9B7CD9A4A7443C, 720},
53 : : {0xBB764C4CA7A44410, 747}, {0x8BAB8EEFB6409C1A, 774},
54 : : {0xD01FEF10A657842C, 800}, {0x9B10A4E5E9913129, 827},
55 : : {0xE7109BFBA19C0C9D, 853}, {0xAC2820D9623BF429, 880},
56 : : {0x80444B5E7AA7CF85, 907}, {0xBF21E44003ACDD2D, 933},
57 : : {0x8E679C2F5E44FF8F, 960}, {0xD433179D9C8CB841, 986},
58 : : {0x9E19DB92B4E31BA9, 1013}, {0xEB96BF6EBADF77D9, 1039},
59 : : {0xAF87023B9BF0EE6B, 1066}};
60 : :
61 : : EXESS_I_NONBLOCKING ExessSoftFloat
62 : 30667 : soft_float_from_double(const double d)
63 : : {
64 [ - + ]: 30667 : assert(d >= 0.0);
65 : :
66 : 30667 : const uint64_t rep = double_to_rep(d);
67 : 30667 : const uint64_t frac = rep & dbl_mant_mask;
68 : 30667 : const int expt = (int)((rep & dbl_expt_mask) >> dbl_physical_mant_dig);
69 : :
70 [ + + ]: 30667 : if (expt == 0) { // Subnormal
71 : 6 : ExessSoftFloat v = {frac, dbl_subnormal_expt};
72 : 6 : return v;
73 : : }
74 : :
75 : 30661 : const ExessSoftFloat v = {frac + dbl_hidden_bit, expt - dbl_expt_bias};
76 : 30661 : return v;
77 : : }
78 : :
79 : : EXESS_I_NONBLOCKING double
80 : 23630 : soft_float_to_double(const ExessSoftFloat v)
81 : : {
82 : 23630 : return ldexp((double)v.f, v.e);
83 : : }
84 : :
85 : : EXESS_I_NONBLOCKING ExessSoftFloat
86 : 47260 : soft_float_normalize(ExessSoftFloat value)
87 : : {
88 : 47260 : const unsigned amount = exess_clz64(value.f);
89 : :
90 : 47260 : value.f <<= amount;
91 : 47260 : value.e -= (int)amount;
92 : :
93 : 47260 : return value;
94 : : }
95 : :
96 : : EXESS_I_NONBLOCKING ExessSoftFloat
97 : 44183 : soft_float_multiply(const ExessSoftFloat lhs, const ExessSoftFloat rhs)
98 : : {
99 : 44183 : EXESS_I_CONSTEXPR uint64_t mask = 0xFFFFFFFF;
100 : 44183 : EXESS_I_CONSTEXPR uint64_t round = 1UL << 31U;
101 : :
102 : 44183 : const uint64_t l0 = lhs.f >> 32U;
103 : 44183 : const uint64_t l1 = lhs.f & mask;
104 : 44183 : const uint64_t r0 = rhs.f >> 32U;
105 : 44183 : const uint64_t r1 = rhs.f & mask;
106 : 44183 : const uint64_t l0r0 = l0 * r0;
107 : 44183 : const uint64_t l1r0 = l1 * r0;
108 : 44183 : const uint64_t l0r1 = l0 * r1;
109 : 44183 : const uint64_t l1r1 = l1 * r1;
110 : 44183 : const uint64_t mid = (l1r1 >> 32U) + (l0r1 & mask) + (l1r0 & mask) + round;
111 : :
112 : 44183 : const ExessSoftFloat r = {l0r0 + (l0r1 >> 32U) + (l1r0 >> 32U) + (mid >> 32U),
113 : 44183 : lhs.e + rhs.e + 64};
114 : :
115 : 44183 : return r;
116 : : }
117 : :
118 : : EXESS_I_NONBLOCKING ExessSoftFloat
119 : 20553 : soft_float_exact_pow10(const int expt)
120 : : {
121 : 20553 : EXESS_I_CONSTEXPR ExessSoftFloat table[8] = {{0xA000000000000000, -60},
122 : : {0xC800000000000000, -57},
123 : : {0xFA00000000000000, -54},
124 : : {0x9C40000000000000, -50},
125 : : {0xC350000000000000, -47},
126 : : {0xF424000000000000, -44},
127 : : {0x9896800000000000, -40}};
128 : :
129 [ - + ]: 20553 : assert(expt > 0);
130 [ - + ]: 20553 : assert(expt < dec_expt_step);
131 : 20553 : return table[expt - 1];
132 : : }
133 : :
134 : : EXESS_I_NONBLOCKING ExessSoftFloat
135 : 23630 : soft_float_pow10_under(const int max_exponent, int* pow10_exponent)
136 : : {
137 [ - + ]: 23630 : assert(max_exponent >= min_dec_expt);
138 [ - + ]: 23630 : assert(max_exponent < max_dec_expt + dec_expt_step);
139 : :
140 : 23630 : const int cache_offset = -min_dec_expt;
141 : 23630 : const int index = (max_exponent + cache_offset) / dec_expt_step;
142 : :
143 : 23630 : *pow10_exponent = min_dec_expt + (index * dec_expt_step);
144 : :
145 [ - + ]: 23630 : assert(*pow10_exponent <= max_exponent);
146 [ - + ]: 23630 : assert(max_exponent < *pow10_exponent + dec_expt_step);
147 : :
148 : 23630 : return soft_pow10[index];
149 : : }
|
