1 /**
2 * \file
3 * System.Double, System.Single and System.Number runtime support
4 *
5 * Author:
6 * Ludovic Henry (ludovic@xamarin.com)
7 *
8 * Copyright 2015 Xamarin, Inc (http://www.xamarin.com)
9 * Licensed under the MIT license. See LICENSE file in the project root for full license information.
10 */
11
12 //
13 // Copyright (c) Microsoft. All rights reserved.
14 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
15 //
16 // Files:
17 // - src/classlibnative/bcltype/number.cpp
18 //
19 // Ported from C++ to C and adjusted to Mono runtime
20
21 #include <glib.h>
22
23 #include "number-ms.h"
24
25 static const guint64 rgval64Power10[] = {
26 /* powers of 10 */
27 0xa000000000000000LL,
28 0xc800000000000000LL,
29 0xfa00000000000000LL,
30 0x9c40000000000000LL,
31 0xc350000000000000LL,
32 0xf424000000000000LL,
33 0x9896800000000000LL,
34 0xbebc200000000000LL,
35 0xee6b280000000000LL,
36 0x9502f90000000000LL,
37 0xba43b74000000000LL,
38 0xe8d4a51000000000LL,
39 0x9184e72a00000000LL,
40 0xb5e620f480000000LL,
41 0xe35fa931a0000000LL,
42
43 /* powers of 0.1 */
44 0xcccccccccccccccdLL,
45 0xa3d70a3d70a3d70bLL,
46 0x83126e978d4fdf3cLL,
47 0xd1b71758e219652eLL,
48 0xa7c5ac471b478425LL,
49 0x8637bd05af6c69b7LL,
50 0xd6bf94d5e57a42beLL,
51 0xabcc77118461ceffLL,
52 0x89705f4136b4a599LL,
53 0xdbe6fecebdedd5c2LL,
54 0xafebff0bcb24ab02LL,
55 0x8cbccc096f5088cfLL,
56 0xe12e13424bb40e18LL,
57 0xb424dc35095cd813LL,
58 0x901d7cf73ab0acdcLL,
59 };
60
61 static const gint8 rgexp64Power10[] = {
62 /* exponents for both powers of 10 and 0.1 */
63 4,
64 7,
65 10,
66 14,
67 17,
68 20,
69 24,
70 27,
71 30,
72 34,
73 37,
74 40,
75 44,
76 47,
77 50,
78 };
79
80 static const guint64 rgval64Power10By16[] = {
81 /* powers of 10^16 */
82 0x8e1bc9bf04000000LL,
83 0x9dc5ada82b70b59eLL,
84 0xaf298d050e4395d6LL,
85 0xc2781f49ffcfa6d4LL,
86 0xd7e77a8f87daf7faLL,
87 0xefb3ab16c59b14a0LL,
88 0x850fadc09923329cLL,
89 0x93ba47c980e98cdeLL,
90 0xa402b9c5a8d3a6e6LL,
91 0xb616a12b7fe617a8LL,
92 0xca28a291859bbf90LL,
93 0xe070f78d39275566LL,
94 0xf92e0c3537826140LL,
95 0x8a5296ffe33cc92cLL,
96 0x9991a6f3d6bf1762LL,
97 0xaa7eebfb9df9de8aLL,
98 0xbd49d14aa79dbc7eLL,
99 0xd226fc195c6a2f88LL,
100 0xe950df20247c83f8LL,
101 0x81842f29f2cce373LL,
102 0x8fcac257558ee4e2LL,
103
104 /* powers of 0.1^16 */
105 0xe69594bec44de160LL,
106 0xcfb11ead453994c3LL,
107 0xbb127c53b17ec165LL,
108 0xa87fea27a539e9b3LL,
109 0x97c560ba6b0919b5LL,
110 0x88b402f7fd7553abLL,
111 0xf64335bcf065d3a0LL,
112 0xddd0467c64bce4c4LL,
113 0xc7caba6e7c5382edLL,
114 0xb3f4e093db73a0b7LL,
115 0xa21727db38cb0053LL,
116 0x91ff83775423cc29LL,
117 0x8380dea93da4bc82LL,
118 0xece53cec4a314f00LL,
119 0xd5605fcdcf32e217LL,
120 0xc0314325637a1978LL,
121 0xad1c8eab5ee43ba2LL,
122 0x9becce62836ac5b0LL,
123 0x8c71dcd9ba0b495cLL,
124 0xfd00b89747823938LL,
125 0xe3e27a444d8d991aLL,
126 };
127
128 static const gint16 rgexp64Power10By16[] = {
129 /* exponents for both powers of 10^16 and 0.1^16 */
130 54,
131 107,
132 160,
133 213,
134 266,
135 319,
136 373,
137 426,
138 479,
139 532,
140 585,
141 638,
142 691,
143 745,
144 798,
145 851,
146 904,
147 957,
148 1010,
149 1064,
150 1117,
151 };
152
153 static inline guint64
digits_to_int(guint16 * p,int count)154 digits_to_int (guint16 *p, int count)
155 {
156 g_assert (1 <= count && count <= 9);
157 guint8 i = 0;
158 guint64 res = 0;
159 switch (count) {
160 case 9: res += 100000000 * (p [i++] - '0');
161 case 8: res += 10000000 * (p [i++] - '0');
162 case 7: res += 1000000 * (p [i++] - '0');
163 case 6: res += 100000 * (p [i++] - '0');
164 case 5: res += 10000 * (p [i++] - '0');
165 case 4: res += 1000 * (p [i++] - '0');
166 case 3: res += 100 * (p [i++] - '0');
167 case 2: res += 10 * (p [i++] - '0');
168 case 1: res += 1 * (p [i++] - '0');
169 }
170 return res;
171 }
172
173 static inline guint64
mul_64_lossy(guint64 a,guint64 b,gint * pexp)174 mul_64_lossy (guint64 a, guint64 b, gint *pexp)
175 {
176 /* it's ok to losse some precision here - it will be called
177 * at most twice during the conversion, so the error won't
178 * propagate to any of the 53 significant bits of the result */
179 guint64 val =
180 ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b >> 32))) )
181 + ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b ))) >> 32)
182 + ((((guint64) (guint32) (a )) * ((guint64) (guint32) (b >> 32))) >> 32);
183
184 /* normalize */
185 if ((val & 0x8000000000000000LL) == 0) {
186 val <<= 1;
187 *pexp -= 1;
188 }
189
190 return val;
191 }
192
193 static inline void
number_to_double(MonoNumber * number,gdouble * value)194 number_to_double (MonoNumber *number, gdouble *value)
195 {
196 guint64 val;
197 guint16 *src;
198 gint exp, remaining, total, count, scale, absscale, index;
199
200 total = 0;
201 src = number->digits;
202 while (*src++) total ++;
203
204 remaining = total;
205
206 src = number->digits;
207 while (*src == '0') {
208 remaining --;
209 src ++;
210 }
211
212 if (remaining == 0) {
213 *value = 0;
214 goto done;
215 }
216
217 count = MIN (remaining, 9);
218 remaining -= count;
219 val = digits_to_int (src, count);
220
221 if (remaining > 0) {
222 count = MIN (remaining, 9);
223 remaining -= count;
224
225 /* get the denormalized power of 10 */
226 guint32 mult = (guint32) (rgval64Power10 [count - 1] >> (64 - rgexp64Power10 [count - 1]));
227 val = ((guint64) (guint32) val) * ((guint64) mult) + digits_to_int (src + 9, count);
228 }
229
230 scale = number->scale - (total - remaining);
231 absscale = abs (scale);
232
233 if (absscale >= 22 * 16) {
234 /* overflow / underflow */
235 *(guint64*) value = (scale > 0) ? 0x7FF0000000000000LL : 0;
236 goto done;
237 }
238
239 exp = 64;
240
241 /* normalize the mantiss */
242 if ((val & 0xFFFFFFFF00000000LL) == 0) { val <<= 32; exp -= 32; }
243 if ((val & 0xFFFF000000000000LL) == 0) { val <<= 16; exp -= 16; }
244 if ((val & 0xFF00000000000000LL) == 0) { val <<= 8; exp -= 8; }
245 if ((val & 0xF000000000000000LL) == 0) { val <<= 4; exp -= 4; }
246 if ((val & 0xC000000000000000LL) == 0) { val <<= 2; exp -= 2; }
247 if ((val & 0x8000000000000000LL) == 0) { val <<= 1; exp -= 1; }
248
249 index = absscale & 15;
250 if (index) {
251 gint multexp = rgexp64Power10 [index - 1];
252 /* the exponents are shared between the inverted and regular table */
253 exp += (scale < 0) ? (-multexp + 1) : multexp;
254
255 guint64 multval = rgval64Power10 [index + ((scale < 0) ? 15 : 0) - 1];
256 val = mul_64_lossy (val, multval, &exp);
257 }
258
259 index = absscale >> 4;
260 if (index) {
261 gint multexp = rgexp64Power10By16 [index - 1];
262 /* the exponents are shared between the inverted and regular table */
263 exp += (scale < 0) ? (-multexp + 1) : multexp;
264
265 guint64 multval = rgval64Power10By16 [index + ((scale < 0) ? 21 : 0) - 1];
266 val = mul_64_lossy (val, multval, &exp);
267 }
268
269 if ((guint32) val & (1 << 10)) {
270 /* IEEE round to even */
271 guint64 tmp = val + ((1 << 10) - 1) + (((guint32) val >> 11) & 1);
272 if (tmp < val) {
273 /* overflow */
274 tmp = (tmp >> 1) | 0x8000000000000000LL;
275 exp += 1;
276 }
277 val = tmp;
278 }
279
280 /* return the exponent to a biased state */
281 exp += 0x3FE;
282
283 /* handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case */
284 if (exp <= 0) {
285 if (exp == -52 && (val >= 0x8000000000000058LL)) {
286 /* round X where {Epsilon > X >= 2.470328229206232730000000E-324} up to Epsilon (instead of down to zero) */
287 val = 0x0000000000000001LL;
288 } else if (exp <= -52) {
289 /* underflow */
290 val = 0;
291 } else {
292 /* denormalized */
293 val >>= (-exp + 11 + 1);
294 }
295 } else if (exp >= 0x7FF) {
296 /* overflow */
297 val = 0x7FF0000000000000LL;
298 } else {
299 /* normal postive exponent case */
300 val = ((guint64) exp << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFFLL);
301 }
302
303 *(guint64*) value = val;
304
305 done:
306 if (number->sign)
307 *(guint64*) value |= 0x8000000000000000LL;
308 }
309
310 gint
mono_double_from_number(gpointer from,MonoDouble * target)311 mono_double_from_number (gpointer from, MonoDouble *target)
312 {
313 MonoDouble_double res;
314 guint e, mant_lo, mant_hi;
315
316 res.d = 0;
317
318 number_to_double ((MonoNumber*) from, &res.d);
319 e = res.s.exp;
320 mant_lo = res.s.mantLo;
321 mant_hi = res.s.mantHi;
322
323 if (e == 0x7ff)
324 return 0;
325
326 if (e == 0 && mant_lo == 0 && mant_hi == 0)
327 res.d = 0;
328
329 *target = res.s;
330 return 1;
331 }
332