mono/metadata/number-ms.c

   1 /*
   2  * number-ms.c: System.Double, System.Single and System.Number runtime support
   3  *
   4  * Author:
   5  *      Ludovic Henry (ludovic@xamarin.com)
   6  *
   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.
   9  */
  10
  11 //
  12 // Copyright (c) Microsoft. All rights reserved.
  13 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
  14 //
  15 // Files:
  16 //  - src/classlibnative/bcltype/number.cpp
  17 //
  18 // Ported from C++ to C and adjusted to Mono runtime
  19
  20 #include <glib.h>
  21
  22 #include "number-ms.h"
  23
  24 static const guint64 rgval64Power10[] = {
  25         /* powers of 10 */
  26         0xa000000000000000LL,
  27         0xc800000000000000LL,
  28         0xfa00000000000000LL,
  29         0x9c40000000000000LL,
  30         0xc350000000000000LL,
  31         0xf424000000000000LL,
  32         0x9896800000000000LL,
  33         0xbebc200000000000LL,
  34         0xee6b280000000000LL,
  35         0x9502f90000000000LL,
  36         0xba43b74000000000LL,
  37         0xe8d4a51000000000LL,
  38         0x9184e72a00000000LL,
  39         0xb5e620f480000000LL,
  40         0xe35fa931a0000000LL,
  41
  42         /* powers of 0.1 */
  43         0xcccccccccccccccdLL,
  44         0xa3d70a3d70a3d70bLL,
  45         0x83126e978d4fdf3cLL,
  46         0xd1b71758e219652eLL,
  47         0xa7c5ac471b478425LL,
  48         0x8637bd05af6c69b7LL,
  49         0xd6bf94d5e57a42beLL,
  50         0xabcc77118461ceffLL,
  51         0x89705f4136b4a599LL,
  52         0xdbe6fecebdedd5c2LL,
  53         0xafebff0bcb24ab02LL,
  54         0x8cbccc096f5088cfLL,
  55         0xe12e13424bb40e18LL,
  56         0xb424dc35095cd813LL,
  57         0x901d7cf73ab0acdcLL,
  58 };
  59
  60 static const gint8 rgexp64Power10[] = {
  61         /* exponents for both powers of 10 and 0.1 */
  62         4,
  63         7,
  64         10,
  65         14,
  66         17,
  67         20,
  68         24,
  69         27,
  70         30,
  71         34,
  72         37,
  73         40,
  74         44,
  75         47,
  76         50,
  77 };
  78
  79 static const guint64 rgval64Power10By16[] = {
  80         /* powers of 10^16 */
  81         0x8e1bc9bf04000000LL,
  82         0x9dc5ada82b70b59eLL,
  83         0xaf298d050e4395d6LL,
  84         0xc2781f49ffcfa6d4LL,
  85         0xd7e77a8f87daf7faLL,
  86         0xefb3ab16c59b14a0LL,
  87         0x850fadc09923329cLL,
  88         0x93ba47c980e98cdeLL,
  89         0xa402b9c5a8d3a6e6LL,
  90         0xb616a12b7fe617a8LL,
  91         0xca28a291859bbf90LL,
  92         0xe070f78d39275566LL,
  93         0xf92e0c3537826140LL,
  94         0x8a5296ffe33cc92cLL,
  95         0x9991a6f3d6bf1762LL,
  96         0xaa7eebfb9df9de8aLL,
  97         0xbd49d14aa79dbc7eLL,
  98         0xd226fc195c6a2f88LL,
  99         0xe950df20247c83f8LL,
 100         0x81842f29f2cce373LL,
 101         0x8fcac257558ee4e2LL,
 102
 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,
 125 };
 126
 127 static const gint16 rgexp64Power10By16[] = {
 128         /* exponents for both powers of 10^16 and 0.1^16 */
 129         54,
 130         107,
 131         160,
 132         213,
 133         266,
 134         319,
 135         373,
 136         426,
 137         479,
 138         532,
 139         585,
 140         638,
 141         691,
 142         745,
 143         798,
 144         851,
 145         904,
 146         957,
 147         1010,
 148         1064,
 149         1117,
 150 };
 151
 152 static inline guint64
 153 digits_to_int (guint16 *p, int count)
 154 {
 155         g_assert (1 <= count && count <= 9);
 156         guint8 i = 0;
 157         guint64 res = 0;
 158         switch (count) {
 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');
 168         }
 169         return res;
 170 }
 171
 172 static inline guint64
 173 mul_64_lossy (guint64 a, guint64 b, gint *pexp)
 174 {
 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 */
 178         guint64 val =
 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);
 182
 183         /* normalize */
 184         if ((val & 0x8000000000000000LL) == 0) {
 185                 val <<= 1;
 186                 *pexp -= 1;
 187         }
 188
 189         return val;
 190 }
 191
 192 static inline void
 193 number_to_double (MonoNumber *number, gdouble *value)
 194 {
 195         guint64 val;
 196         guint16 *src;
 197         gint exp, remaining, total, count, scale, absscale, index;
 198
 199         total = 0;
 200         src = number->digits;
 201         while (*src++) total ++;
 202
 203         remaining = total;
 204
 205         src = number->digits;
 206         while (*src == '0') {
 207                 remaining --;
 208                 src ++;
 209         }
 210
 211         if (remaining == 0) {
 212                 *value = 0;
 213                 goto done;
 214         }
 215
 216         count = MIN (remaining, 9);
 217         remaining -= count;
 218         val = digits_to_int (src, count);
 219
 220         if (remaining > 0) {
 221                 count = MIN (remaining, 9);
 222                 remaining -= count;
 223
 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);
 227         }
 228
 229         scale = number->scale - (total - remaining);
 230         absscale = abs (scale);
 231
 232         if (absscale >= 22 * 16) {
 233                 /* overflow / underflow */
 234                 *(guint64*) value = (scale > 0) ? 0x7FF0000000000000LL : 0;
 235                 goto done;
 236         }
 237
 238         exp = 64;
 239
 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;  }
 247
 248         index = absscale & 15;
 249         if (index) {
 250                 gint multexp = rgexp64Power10 [index - 1];
 251                 /* the exponents are shared between the inverted and regular table */
 252                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 253
 254                 guint64 multval = rgval64Power10 [index + ((scale < 0) ? 15 : 0) - 1];
 255                 val = mul_64_lossy (val, multval, &exp);
 256         }
 257
 258         index = absscale >> 4;
 259         if (index) {
 260                 gint multexp = rgexp64Power10By16 [index - 1];
 261                 /* the exponents are shared between the inverted and regular table */
 262                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 263
 264                 guint64 multval = rgval64Power10By16 [index + ((scale < 0) ? 21 : 0) - 1];
 265                 val = mul_64_lossy (val, multval, &exp);
 266         }
 267
 268         if ((guint32) val & (1 << 10)) {
 269                 /* IEEE round to even */
 270                 guint64 tmp = val + ((1 << 10) - 1) + (((guint32) val >> 11) & 1);
 271                 if (tmp < val) {
 272                         /* overflow */
 273                         tmp = (tmp >> 1) | 0x8000000000000000LL;
 274                         exp += 1;
 275                 }
 276                 val = tmp;
 277         }
 278
 279         /* return the exponent to a biased state */
 280         exp += 0x3FE;
 281
 282         /* handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case */
 283         if (exp <= 0) {
 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) {
 288                         /* underflow */
 289                         val = 0;
 290                 } else {
 291                         /* denormalized */
 292                         val >>= (-exp + 11 + 1);
 293                 }
 294         } else if (exp >= 0x7FF) {
 295                 /* overflow */
 296                 val = 0x7FF0000000000000LL;
 297         } else {
 298                 /* normal postive exponent case */
 299                 val = ((guint64) exp << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFFLL);
 300         }
 301
 302         *(guint64*) value = val;
 303
 304 done:
 305         if (number->sign)
 306                 *(guint64*) value |= 0x8000000000000000LL;
 307 }
 308
 309 gint
 310 mono_double_from_number (gpointer from, MonoDouble *target)
 311 {
 312         MonoDouble_double res;
 313         guint e, mant_lo, mant_hi;
 314
 315         res.d = 0;
 316
 317         number_to_double ((MonoNumber*) from, &res.d);
 318         e = res.s.exp;
 319         mant_lo = res.s.mantLo;
 320         mant_hi = res.s.mantHi;
 321
 322         if (e == 0x7ff)
 323                 return 0;
 324
 325         if (e == 0 && mant_lo == 0 && mant_hi == 0)
 326                 res.d = 0;
 327
 328         *target = res.s;
 329         return 1;
 330 }