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)
11 // Copyright (c) Microsoft. All rights reserved.
12 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
15 // - src/classlibnative/bcltype/number.cpp
17 // Ported from C++ to C and adjusted to Mono runtime
21 #include "number-ms.h"
23 static const guint64 rgval64Power10[] = {
59 static const gint8 rgexp64Power10[] = {
60 /* exponents for both powers of 10 and 0.1 */
78 static const guint64 rgval64Power10By16[] = {
100 0x8fcac257558ee4e2LL,
102 /* powers of 0.1^16 */
103 0xe69594bec44de160LL,
104 0xcfb11ead453994c3LL,
105 0xbb127c53b17ec165LL,
106 0xa87fea27a539e9b3LL,
107 0x97c560ba6b0919b5LL,
108 0x88b402f7fd7553abLL,
109 0xf64335bcf065d3a0LL,
110 0xddd0467c64bce4c4LL,
111 0xc7caba6e7c5382edLL,
112 0xb3f4e093db73a0b7LL,
113 0xa21727db38cb0053LL,
114 0x91ff83775423cc29LL,
115 0x8380dea93da4bc82LL,
116 0xece53cec4a314f00LL,
117 0xd5605fcdcf32e217LL,
118 0xc0314325637a1978LL,
119 0xad1c8eab5ee43ba2LL,
120 0x9becce62836ac5b0LL,
121 0x8c71dcd9ba0b495cLL,
122 0xfd00b89747823938LL,
123 0xe3e27a444d8d991aLL,
126 static const gint16 rgexp64Power10By16[] = {
127 /* exponents for both powers of 10^16 and 0.1^16 */
151 static inline guint64
152 digits_to_int (guint16 *p, int count)
154 g_assert (1 <= count && count <= 9);
158 case 9: res += 100000000 * (p [i++] - '0');
159 case 8: res += 10000000 * (p [i++] - '0');
160 case 7: res += 1000000 * (p [i++] - '0');
161 case 6: res += 100000 * (p [i++] - '0');
162 case 5: res += 10000 * (p [i++] - '0');
163 case 4: res += 1000 * (p [i++] - '0');
164 case 3: res += 100 * (p [i++] - '0');
165 case 2: res += 10 * (p [i++] - '0');
166 case 1: res += 1 * (p [i++] - '0');
171 static inline guint64
172 mul_64_lossy (guint64 a, guint64 b, gint *pexp)
174 /* it's ok to losse some precision here - it will be called
175 * at most twice during the conversion, so the error won't
176 * propagate to any of the 53 significant bits of the result */
178 ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b >> 32))) )
179 + ((((guint64) (guint32) (a >> 32)) * ((guint64) (guint32) (b ))) >> 32)
180 + ((((guint64) (guint32) (a )) * ((guint64) (guint32) (b >> 32))) >> 32);
183 if ((val & 0x8000000000000000LL) == 0) {
192 number_to_double (MonoNumber *number, gdouble *value)
196 gint exp, remaining, total, count, scale, absscale, index;
198 g_assert (number->digits);
201 src = number->digits;
202 while (*src++) total ++;
206 src = number->digits;
207 while (*src == '0') {
212 if (remaining == 0) {
217 count = MIN (remaining, 9);
219 val = digits_to_int (src, count);
222 count = MIN (remaining, 9);
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);
230 scale = number->scale - (total - remaining);
231 absscale = abs (scale);
233 if (absscale >= 22 * 16) {
234 /* overflow / underflow */
235 *(guint64*) value = (scale > 0) ? 0x7FF0000000000000LL : 0;
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; }
249 index = absscale & 15;
251 gint multexp = rgexp64Power10 [index - 1];
252 /* the exponents are shared between the inverted and regular table */
253 exp += (scale < 0) ? (-multexp + 1) : multexp;
255 guint64 multval = rgval64Power10 [index + ((scale < 0) ? 15 : 0) - 1];
256 val = mul_64_lossy (val, multval, &exp);
259 index = absscale >> 4;
261 gint multexp = rgexp64Power10By16 [index - 1];
262 /* the exponents are shared between the inverted and regular table */
263 exp += (scale < 0) ? (-multexp + 1) : multexp;
265 guint64 multval = rgval64Power10By16 [index + ((scale < 0) ? 21 : 0) - 1];
266 val = mul_64_lossy (val, multval, &exp);
269 if ((guint32) val & (1 << 10)) {
270 /* IEEE round to even */
271 guint64 tmp = val + ((1 << 10) - 1) + (((guint32) val >> 11) & 1);
274 tmp = (tmp >> 1) | 0x8000000000000000LL;
280 /* return the exponent to a biased state */
283 /* handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case */
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) {
293 val >>= (-exp + 11 + 1);
295 } else if (exp >= 0x7FF) {
297 val = 0x7FF0000000000000LL;
299 /* normal postive exponent case */
300 val = ((guint64) exp << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFFLL);
303 *(guint64*) value = val;
307 *(guint64*) value |= 0x8000000000000000LL;
311 mono_double_from_number (gpointer from, MonoDouble *target)
313 MonoDouble_double res;
314 guint e, mant_lo, mant_hi;
318 number_to_double ((MonoNumber*) from, &res.d);
320 mant_lo = res.s.mantLo;
321 mant_hi = res.s.mantHi;
326 if (e == 0 && mant_lo == 0 && mant_hi == 0)