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  */
   9
  10 //
  11 // Copyright (c) Microsoft. All rights reserved.
  12 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
  13 //
  14 // Files:
  15 //  - src/classlibnative/bcltype/number.cpp
  16 //
  17 // Ported from C++ to C and adjusted to Mono runtime
  18
  19 #include <glib.h>
  20
  21 #include "number-ms.h"
  22
  23 static const guint64 rgval64Power10[] = {
  24         /* powers of 10 */
  25         0xa000000000000000LL,
  26         0xc800000000000000LL,
  27         0xfa00000000000000LL,
  28         0x9c40000000000000LL,
  29         0xc350000000000000LL,
  30         0xf424000000000000LL,
  31         0x9896800000000000LL,
  32         0xbebc200000000000LL,
  33         0xee6b280000000000LL,
  34         0x9502f90000000000LL,
  35         0xba43b74000000000LL,
  36         0xe8d4a51000000000LL,
  37         0x9184e72a00000000LL,
  38         0xb5e620f480000000LL,
  39         0xe35fa931a0000000LL,
  40
  41         /* powers of 0.1 */
  42         0xcccccccccccccccdLL,
  43         0xa3d70a3d70a3d70bLL,
  44         0x83126e978d4fdf3cLL,
  45         0xd1b71758e219652eLL,
  46         0xa7c5ac471b478425LL,
  47         0x8637bd05af6c69b7LL,
  48         0xd6bf94d5e57a42beLL,
  49         0xabcc77118461ceffLL,
  50         0x89705f4136b4a599LL,
  51         0xdbe6fecebdedd5c2LL,
  52         0xafebff0bcb24ab02LL,
  53         0x8cbccc096f5088cfLL,
  54         0xe12e13424bb40e18LL,
  55         0xb424dc35095cd813LL,
  56         0x901d7cf73ab0acdcLL,
  57 };
  58
  59 static const gint8 rgexp64Power10[] = {
  60         /* exponents for both powers of 10 and 0.1 */
  61         4,
  62         7,
  63         10,
  64         14,
  65         17,
  66         20,
  67         24,
  68         27,
  69         30,
  70         34,
  71         37,
  72         40,
  73         44,
  74         47,
  75         50,
  76 };
  77
  78 static const guint64 rgval64Power10By16[] = {
  79         /* powers of 10^16 */
  80         0x8e1bc9bf04000000LL,
  81         0x9dc5ada82b70b59eLL,
  82         0xaf298d050e4395d6LL,
  83         0xc2781f49ffcfa6d4LL,
  84         0xd7e77a8f87daf7faLL,
  85         0xefb3ab16c59b14a0LL,
  86         0x850fadc09923329cLL,
  87         0x93ba47c980e98cdeLL,
  88         0xa402b9c5a8d3a6e6LL,
  89         0xb616a12b7fe617a8LL,
  90         0xca28a291859bbf90LL,
  91         0xe070f78d39275566LL,
  92         0xf92e0c3537826140LL,
  93         0x8a5296ffe33cc92cLL,
  94         0x9991a6f3d6bf1762LL,
  95         0xaa7eebfb9df9de8aLL,
  96         0xbd49d14aa79dbc7eLL,
  97         0xd226fc195c6a2f88LL,
  98         0xe950df20247c83f8LL,
  99         0x81842f29f2cce373LL,
 100         0x8fcac257558ee4e2LL,
 101
 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,
 124 };
 125
 126 static const gint16 rgexp64Power10By16[] = {
 127         /* exponents for both powers of 10^16 and 0.1^16 */
 128         54,
 129         107,
 130         160,
 131         213,
 132         266,
 133         319,
 134         373,
 135         426,
 136         479,
 137         532,
 138         585,
 139         638,
 140         691,
 141         745,
 142         798,
 143         851,
 144         904,
 145         957,
 146         1010,
 147         1064,
 148         1117,
 149 };
 150
 151 static inline guint64
 152 digits_to_int (guint16 *p, int count)
 153 {
 154         g_assert (1 <= count && count <= 9);
 155         guint8 i = 0;
 156         guint64 res = 0;
 157         switch (count) {
 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');
 167         }
 168         return res;
 169 }
 170
 171 static inline guint64
 172 mul_64_lossy (guint64 a, guint64 b, gint *pexp)
 173 {
 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 */
 177         guint64 val =
 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);
 181
 182         /* normalize */
 183         if ((val & 0x8000000000000000LL) == 0) {
 184                 val <<= 1;
 185                 *pexp -= 1;
 186         }
 187
 188         return val;
 189 }
 190
 191 static inline void
 192 number_to_double (MonoNumber *number, gdouble *value)
 193 {
 194         guint64 val;
 195         guint16 *src;
 196         gint exp, remaining, total, count, scale, absscale, index;
 197
 198         g_assert (number->digits);
 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 }