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         total = 0;
 199         src = number->digits;
 200         while (*src++) total ++;
 201
 202         remaining = total;
 203
 204         src = number->digits;
 205         while (*src == '0') {
 206                 remaining --;
 207                 src ++;
 208         }
 209
 210         if (remaining == 0) {
 211                 *value = 0;
 212                 goto done;
 213         }
 214
 215         count = MIN (remaining, 9);
 216         remaining -= count;
 217         val = digits_to_int (src, count);
 218
 219         if (remaining > 0) {
 220                 count = MIN (remaining, 9);
 221                 remaining -= count;
 222
 223                 /* get the denormalized power of 10 */
 224                 guint32 mult = (guint32) (rgval64Power10 [count - 1] >> (64 - rgexp64Power10 [count - 1]));
 225                 val = ((guint64) (guint32) val) * ((guint64) mult) + digits_to_int (src + 9, count);
 226         }
 227
 228         scale = number->scale - (total - remaining);
 229         absscale = abs (scale);
 230
 231         if (absscale >= 22 * 16) {
 232                 /* overflow / underflow */
 233                 *(guint64*) value = (scale > 0) ? 0x7FF0000000000000LL : 0;
 234                 goto done;
 235         }
 236
 237         exp = 64;
 238
 239         /* normalize the mantiss */
 240         if ((val & 0xFFFFFFFF00000000LL) == 0) { val <<= 32; exp -= 32; }
 241         if ((val & 0xFFFF000000000000LL) == 0) { val <<= 16; exp -= 16; }
 242         if ((val & 0xFF00000000000000LL) == 0) { val <<= 8;  exp -= 8;  }
 243         if ((val & 0xF000000000000000LL) == 0) { val <<= 4;  exp -= 4;  }
 244         if ((val & 0xC000000000000000LL) == 0) { val <<= 2;  exp -= 2;  }
 245         if ((val & 0x8000000000000000LL) == 0) { val <<= 1;  exp -= 1;  }
 246
 247         index = absscale & 15;
 248         if (index) {
 249                 gint multexp = rgexp64Power10 [index - 1];
 250                 /* the exponents are shared between the inverted and regular table */
 251                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 252
 253                 guint64 multval = rgval64Power10 [index + ((scale < 0) ? 15 : 0) - 1];
 254                 val = mul_64_lossy (val, multval, &exp);
 255         }
 256
 257         index = absscale >> 4;
 258         if (index) {
 259                 gint multexp = rgexp64Power10By16 [index - 1];
 260                 /* the exponents are shared between the inverted and regular table */
 261                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 262
 263                 guint64 multval = rgval64Power10By16 [index + ((scale < 0) ? 21 : 0) - 1];
 264                 val = mul_64_lossy (val, multval, &exp);
 265         }
 266
 267         if ((guint32) val & (1 << 10)) {
 268                 /* IEEE round to even */
 269                 guint64 tmp = val + ((1 << 10) - 1) + (((guint32) val >> 11) & 1);
 270                 if (tmp < val) {
 271                         /* overflow */
 272                         tmp = (tmp >> 1) | 0x8000000000000000LL;
 273                         exp += 1;
 274                 }
 275                 val = tmp;
 276         }
 277
 278         /* return the exponent to a biased state */
 279         exp += 0x3FE;
 280
 281         /* handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case */
 282         if (exp <= 0) {
 283                 if (exp == -52 && (val >= 0x8000000000000058LL)) {
 284                         /* round X where {Epsilon > X >= 2.470328229206232730000000E-324} up to Epsilon (instead of down to zero) */
 285                         val = 0x0000000000000001LL;
 286                 } else if (exp <= -52) {
 287                         /* underflow */
 288                         val = 0;
 289                 } else {
 290                         /* denormalized */
 291                         val >>= (-exp + 11 + 1);
 292                 }
 293         } else if (exp >= 0x7FF) {
 294                 /* overflow */
 295                 val = 0x7FF0000000000000LL;
 296         } else {
 297                 /* normal postive exponent case */
 298                 val = ((guint64) exp << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFFLL);
 299         }
 300
 301         *(guint64*) value = val;
 302
 303 done:
 304         if (number->sign)
 305                 *(guint64*) value |= 0x8000000000000000LL;
 306 }
 307
 308 gint
 309 mono_double_from_number (gpointer from, MonoDouble *target)
 310 {
 311         MonoDouble_double res;
 312         guint e, mant_lo, mant_hi;
 313
 314         res.d = 0;
 315
 316         number_to_double ((MonoNumber*) from, &res.d);
 317         e = res.s.exp;
 318         mant_lo = res.s.mantLo;
 319         mant_hi = res.s.mantHi;
 320
 321         if (e == 0x7ff)
 322                 return 0;
 323
 324         if (e == 0 && mant_lo == 0 && mant_hi == 0)
 325                 res.d = 0;
 326
 327         *target = res.s;
 328         return 1;
 329 }