mono/metadata/number-ms.c

   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
 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
 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
 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
 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 }