2 * number-ms.c: System.Double, System.Single and System.Number runtime support
5 * Ludovic Henry (ludovic@xamarin.com)
7 * Copyright 2015 Xamarin, Inc (http://www.xamarin.com)
8 * Licensed under the MIT license. See LICENSE file in the project root for full license information.
12 // Copyright (c) Microsoft. All rights reserved.
13 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
16 // - src/classlibnative/bcltype/number.cpp
18 // Ported from C++ to C and adjusted to Mono runtime
22 #include "number-ms.h"
24 static const guint64 rgval64Power10[] = {
60 static const gint8 rgexp64Power10[] = {
61 /* exponents for both powers of 10 and 0.1 */
79 static const guint64 rgval64Power10By16[] = {
100 0x81842f29f2cce373LL,
101 0x8fcac257558ee4e2LL,
103 /* powers of 0.1^16 */
104 0xe69594bec44de160LL,
105 0xcfb11ead453994c3LL,
106 0xbb127c53b17ec165LL,
107 0xa87fea27a539e9b3LL,
108 0x97c560ba6b0919b5LL,
109 0x88b402f7fd7553abLL,
110 0xf64335bcf065d3a0LL,
111 0xddd0467c64bce4c4LL,
112 0xc7caba6e7c5382edLL,
113 0xb3f4e093db73a0b7LL,
114 0xa21727db38cb0053LL,
115 0x91ff83775423cc29LL,
116 0x8380dea93da4bc82LL,
117 0xece53cec4a314f00LL,
118 0xd5605fcdcf32e217LL,
119 0xc0314325637a1978LL,
120 0xad1c8eab5ee43ba2LL,
121 0x9becce62836ac5b0LL,
122 0x8c71dcd9ba0b495cLL,
123 0xfd00b89747823938LL,
124 0xe3e27a444d8d991aLL,
127 static const gint16 rgexp64Power10By16[] = {
128 /* exponents for both powers of 10^16 and 0.1^16 */
152 static inline guint64
153 digits_to_int (guint16 *p, int count)
155 g_assert (1 <= count && count <= 9);
159 case 9: res += 100000000 * (p [i++] - '0');
160 case 8: res += 10000000 * (p [i++] - '0');
161 case 7: res += 1000000 * (p [i++] - '0');
162 case 6: res += 100000 * (p [i++] - '0');
163 case 5: res += 10000 * (p [i++] - '0');
164 case 4: res += 1000 * (p [i++] - '0');
165 case 3: res += 100 * (p [i++] - '0');
166 case 2: res += 10 * (p [i++] - '0');
167 case 1: res += 1 * (p [i++] - '0');
172 static inline guint64
173 mul_64_lossy (guint64 a, guint64 b, gint *pexp)
175 /* it's ok to losse some precision here - it will be called
176 * at most twice during the conversion, so the error won't
177 * propagate to any of the 53 significant bits of the result */
179 ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b >> 32))) )
180 + ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b ))) >> 32)
181 + ((((guint64) (guint32) (a )) * ((guint64) (guint32) (b >> 32))) >> 32);
184 if ((val & 0x8000000000000000LL) == 0) {
193 number_to_double (MonoNumber *number, gdouble *value)
197 gint exp, remaining, total, count, scale, absscale, index;
200 src = number->digits;
201 while (*src++) total ++;
205 src = number->digits;
206 while (*src == '0') {
211 if (remaining == 0) {
216 count = MIN (remaining, 9);
218 val = digits_to_int (src, count);
221 count = MIN (remaining, 9);
224 /* get the denormalized power of 10 */
225 guint32 mult = (guint32) (rgval64Power10 [count - 1] >> (64 - rgexp64Power10 [count - 1]));
226 val = ((guint64) (guint32) val) * ((guint64) mult) + digits_to_int (src + 9, count);
229 scale = number->scale - (total - remaining);
230 absscale = abs (scale);
232 if (absscale >= 22 * 16) {
233 /* overflow / underflow */
234 *(guint64*) value = (scale > 0) ? 0x7FF0000000000000LL : 0;
240 /* normalize the mantiss */
241 if ((val & 0xFFFFFFFF00000000LL) == 0) { val <<= 32; exp -= 32; }
242 if ((val & 0xFFFF000000000000LL) == 0) { val <<= 16; exp -= 16; }
243 if ((val & 0xFF00000000000000LL) == 0) { val <<= 8; exp -= 8; }
244 if ((val & 0xF000000000000000LL) == 0) { val <<= 4; exp -= 4; }
245 if ((val & 0xC000000000000000LL) == 0) { val <<= 2; exp -= 2; }
246 if ((val & 0x8000000000000000LL) == 0) { val <<= 1; exp -= 1; }
248 index = absscale & 15;
250 gint multexp = rgexp64Power10 [index - 1];
251 /* the exponents are shared between the inverted and regular table */
252 exp += (scale < 0) ? (-multexp + 1) : multexp;
254 guint64 multval = rgval64Power10 [index + ((scale < 0) ? 15 : 0) - 1];
255 val = mul_64_lossy (val, multval, &exp);
258 index = absscale >> 4;
260 gint multexp = rgexp64Power10By16 [index - 1];
261 /* the exponents are shared between the inverted and regular table */
262 exp += (scale < 0) ? (-multexp + 1) : multexp;
264 guint64 multval = rgval64Power10By16 [index + ((scale < 0) ? 21 : 0) - 1];
265 val = mul_64_lossy (val, multval, &exp);
268 if ((guint32) val & (1 << 10)) {
269 /* IEEE round to even */
270 guint64 tmp = val + ((1 << 10) - 1) + (((guint32) val >> 11) & 1);
273 tmp = (tmp >> 1) | 0x8000000000000000LL;
279 /* return the exponent to a biased state */
282 /* handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case */
284 if (exp == -52 && (val >= 0x8000000000000058LL)) {
285 /* round X where {Epsilon > X >= 2.470328229206232730000000E-324} up to Epsilon (instead of down to zero) */
286 val = 0x0000000000000001LL;
287 } else if (exp <= -52) {
292 val >>= (-exp + 11 + 1);
294 } else if (exp >= 0x7FF) {
296 val = 0x7FF0000000000000LL;
298 /* normal postive exponent case */
299 val = ((guint64) exp << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFFLL);
302 *(guint64*) value = val;
306 *(guint64*) value |= 0x8000000000000000LL;
310 mono_double_from_number (gpointer from, MonoDouble *target)
312 MonoDouble_double res;
313 guint e, mant_lo, mant_hi;
317 number_to_double ((MonoNumber*) from, &res.d);
319 mant_lo = res.s.mantLo;
320 mant_hi = res.s.mantHi;
325 if (e == 0 && mant_lo == 0 && mant_hi == 0)