mono/metadata/decimal.c

   1 /*
   2    decimal.c
   3
   4    conversions and numerical operations for the c# type System.Decimal
   5
   6    Author: Martin Weindel (martin.weindel@t-online.de)
   7
   8    (C) 2001 by Martin Weindel
   9 */
  10
  11 /*
  12  * machine dependent configuration for
  13  * CSharp value type System.Decimal
  14  */
  15
  16 #include <mono/metadata/exception.h>
  17 #include <stdio.h>
  18 #include <memory.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include <math.h>
  22
  23 /* needed for building microsoft dll */
  24 #define DECINLINE __inline
  25
  26 #define LIT_GUINT32(x) x
  27 #define LIT_GUINT64(x) x##L
  28
  29
  30 /* we need a UInt64 type => guint64 */
  31 #include <glib.h>
  32
  33 #include "decimal.h"
  34
  35 /*
  36  * Deal with anon union support.
  37  */
  38 #define ss32 u.ss32
  39 #define signscale u.signscale
  40
  41 /* debugging stuff */
  42 #ifdef _DEBUG
  43 #include <assert.h>
  44 #define PRECONDITION(flag)  assert(flag)
  45 #define POSTCONDITION(flag)  assert(flag)
  46 #define TEST(flag)  assert(flag)
  47 #define INVARIANT_TEST(p) assert(p->signscale.scale >= 0 && p->signscale.scale <= DECIMAL_MAX_SCALE \
  48         && p->signscale.reserved1 == 0 && p->signscale.reserved2 == 0);
  49 #else
  50 #define PRECONDITION(flag)
  51 #define POSTCONDITION(flag)
  52 #define TEST(flag)
  53 #define INVARIANT_TEST(p)
  54 #endif /*#ifdef _DEBUG*/
  55
  56 #define DECIMAL_MAX_SCALE 28
  57 #define DECIMAL_MAX_INTFACTORS 9
  58
  59 #define DECIMAL_SUCCESS 0
  60 #define DECIMAL_FINISHED 1
  61 #define DECIMAL_OVERFLOW 2
  62 #define DECIMAL_INVALID_CHARACTER 2
  63 #define DECIMAL_INTERNAL_ERROR 3
  64 #define DECIMAL_INVALID_BITS 4
  65 #define DECIMAL_DIVIDE_BY_ZERO 5
  66 #define DECIMAL_BUFFER_OVERFLOW 6
  67
  68 /* some MACROS */
  69 #define DECINIT(src) memset(src, 0, sizeof(decimal_repr))
  70
  71 #define DECCOPY(dest, src) memcpy(dest, src, sizeof(decimal_repr))
  72
  73 #define DECSWAP(p1, p2, h) \
  74         h = (p1)->ss32; (p1)->ss32 = (p2)->ss32; (p2)->ss32 = h; \
  75         h = (p1)->hi32; (p1)->hi32 = (p2)->hi32; (p2)->hi32 = h; \
  76         h = (p1)->mid32; (p1)->mid32 = (p2)->mid32; (p2)->mid32 = h; \
  77         h = (p1)->lo32; (p1)->lo32 = (p2)->lo32; (p2)->lo32 = h;
  78
  79 #define DECNEGATE(p1) (p1)->signscale.sign = 1 - (p1)->signscale.sign
  80
  81 #define LIT_DEC128(hi, mid, lo) { (((guint64)mid)<<32 | lo), hi }
  82
  83 #define DECTO128(pd, lo, hi) \
  84         lo = (((guint64)(pd)->mid32) << 32) | (pd)->lo32; \
  85     hi = (pd)->hi32;
  86
  87 /* some constants */
  88 #define LIT_GUINT32_HIGHBIT LIT_GUINT32(0x80000000)
  89 #define LIT_GUINT64_HIGHBIT LIT_GUINT64(0x8000000000000000)
  90
  91 #define DECIMAL_LOG_NEGINF -1000
  92
  93 static guint32 constantsDecadeInt32Factors[DECIMAL_MAX_INTFACTORS+1] = {
  94     LIT_GUINT32(1), LIT_GUINT32(10), LIT_GUINT32(100), LIT_GUINT32(1000),
  95     LIT_GUINT32(10000), LIT_GUINT32(100000), LIT_GUINT32(1000000),
  96     LIT_GUINT32(10000000), LIT_GUINT32(100000000), LIT_GUINT32(1000000000)
  97 };
  98
  99 typedef struct {
 100     guint64 lo;
 101     guint64 hi;
 102 } dec128_repr;
 103
 104 static dec128_repr dec128decadeFactors[DECIMAL_MAX_SCALE+1] = {
 105     LIT_DEC128( 0, 0, 1u), /* == 1 */
 106     LIT_DEC128( 0, 0, 10u), /* == 10 */
 107     LIT_DEC128( 0, 0, 100u), /* == 100 */
 108     LIT_DEC128( 0, 0, 1000u), /* == 1e3m */
 109     LIT_DEC128( 0, 0, 10000u), /* == 1e4m */
 110     LIT_DEC128( 0, 0, 100000u), /* == 1e5m */
 111     LIT_DEC128( 0, 0, 1000000u), /* == 1e6m */
 112     LIT_DEC128( 0, 0, 10000000u), /* == 1e7m */
 113     LIT_DEC128( 0, 0, 100000000u), /* == 1e8m */
 114     LIT_DEC128( 0, 0, 1000000000u), /* == 1e9m */
 115     LIT_DEC128( 0, 2u, 1410065408u), /* == 1e10m */
 116     LIT_DEC128( 0, 23u, 1215752192u), /* == 1e11m */
 117     LIT_DEC128( 0, 232u, 3567587328u), /* == 1e12m */
 118     LIT_DEC128( 0, 2328u, 1316134912u), /* == 1e13m */
 119     LIT_DEC128( 0, 23283u, 276447232u), /* == 1e14m */
 120     LIT_DEC128( 0, 232830u, 2764472320u), /* == 1e15m */
 121     LIT_DEC128( 0, 2328306u, 1874919424u), /* == 1e16m */
 122     LIT_DEC128( 0, 23283064u, 1569325056u), /* == 1e17m */
 123     LIT_DEC128( 0, 232830643u, 2808348672u), /* == 1e18m */
 124     LIT_DEC128( 0, 2328306436u, 2313682944u), /* == 1e19m */
 125     LIT_DEC128( 5u, 1808227885u, 1661992960u), /* == 1e20m */
 126     LIT_DEC128( 54u, 902409669u, 3735027712u), /* == 1e21m */
 127     LIT_DEC128( 542u, 434162106u, 2990538752u), /* == 1e22m */
 128     LIT_DEC128( 5421u, 46653770u, 4135583744u), /* == 1e23m */
 129     LIT_DEC128( 54210u, 466537709u, 2701131776u), /* == 1e24m */
 130     LIT_DEC128( 542101u, 370409800u, 1241513984u), /* == 1e25m */
 131     LIT_DEC128( 5421010u, 3704098002u, 3825205248u), /* == 1e26m */
 132     LIT_DEC128( 54210108u, 2681241660u, 3892314112u), /* == 1e27m */
 133     LIT_DEC128( 542101086u, 1042612833u, 268435456u), /* == 1e28m */
 134 };
 135
 136 /* 192 bit addition: c = a + b
 137    addition is modulo 2**128, any carry is lost */
 138 DECINLINE static void add128(guint64 alo, guint64 ahi,
 139                              guint64 blo, guint64 bhi,
 140                              guint64* pclo, guint64* pchi)
 141 {
 142     alo += blo;
 143     if (alo < blo) ahi++; /* carry */
 144     ahi += bhi;
 145
 146     *pclo = alo;
 147     *pchi = ahi;
 148 }
 149
 150 /* 128 bit subtraction: c = a - b
 151    subtraction is modulo 2**128, any carry is lost */
 152 DECINLINE static void sub128(guint64 alo, guint64 ahi,
 153                              guint64 blo, guint64 bhi,
 154                              guint64* pclo, guint64* pchi)
 155 {
 156     guint64 clo, chi;
 157
 158     clo = alo - blo;
 159     chi = ahi - bhi;
 160     if (alo < blo) chi--; /* borrow */
 161
 162     *pclo = clo;
 163     *pchi = chi;
 164 }
 165
 166 /* 192 bit addition: c = a + b
 167    addition is modulo 2**192, any carry is lost */
 168 DECINLINE static void add192(guint64 alo, guint64 ami, guint64 ahi,
 169                              guint64 blo, guint64 bmi, guint64 bhi,
 170                              guint64* pclo, guint64* pcmi, guint64* pchi)
 171 {
 172     alo += blo;
 173     if (alo < blo) { /* carry low */
 174         ami++;
 175         if (ami == 0) ahi++; /* carry mid */
 176     }
 177     ami += bmi;
 178     if (ami < bmi) ahi++; /* carry mid */
 179     ahi += bhi;
 180     *pclo = alo;
 181     *pcmi = ami;
 182     *pchi = ahi;
 183 }
 184
 185 /* 192 bit subtraction: c = a - b
 186    subtraction is modulo 2**192, any carry is lost */
 187 DECINLINE static void sub192(guint64 alo, guint64 ami, guint64 ahi,
 188                              guint64 blo, guint64 bmi, guint64 bhi,
 189                              guint64* pclo, guint64* pcmi, guint64* pchi)
 190 {
 191     guint64 clo, cmi, chi;
 192
 193     clo = alo - blo;
 194     cmi = ami - bmi;
 195     chi = ahi - bhi;
 196     if (alo < blo) {
 197         if (cmi == 0) chi--; /* borrow mid */
 198         cmi--; /* borrow low */
 199     }
 200     if (ami < bmi) chi--; /* borrow mid */
 201     *pclo = clo;
 202     *pcmi = cmi;
 203     *pchi = chi;
 204 }
 205
 206 /* multiplication c(192bit) = a(96bit) * b(96bit) */
 207 DECINLINE static void mult96by96to192(guint32 alo, guint32 ami, guint32 ahi,
 208                                       guint32 blo, guint32 bmi, guint32 bhi,
 209                                       guint64* pclo, guint64* pcmi, guint64* pchi)
 210 {
 211     guint64 a, b, c, d;
 212     guint32 h0, h1, h2, h3, h4, h5;
 213     int carry0, carry1;
 214
 215     a = ((guint64)alo) * blo;
 216     h0 = (guint32) a;
 217
 218     a >>= 32; carry0 = 0;
 219     b = ((guint64)alo) * bmi;
 220     c = ((guint64)ami) * blo;
 221     a += b; if (a < b) carry0++;
 222     a += c; if (a < c) carry0++;
 223     h1 = (guint32) a;
 224
 225     a >>= 32; carry1 = 0;
 226     b = ((guint64)alo) * bhi;
 227     c = ((guint64)ami) * bmi;
 228     d = ((guint64)ahi) * blo;
 229     a += b; if (a < b) carry1++;
 230     a += c; if (a < c) carry1++;
 231     a += d; if (a < d) carry1++;
 232     h2 = (guint32) a;
 233
 234     a >>= 32; a += carry0; carry0 = 0;
 235     b = ((guint64)ami) * bhi;
 236     c = ((guint64)ahi) * bmi;
 237     a += b; if (a < b) carry0++;
 238     a += c; if (a < c) carry0++;
 239     h3 = (guint32) a;
 240
 241     a >>= 32; a += carry1;
 242     b = ((guint64)ahi) * bhi;
 243     a += b;
 244     h4 = (guint32) a;
 245
 246     a >>= 32; a += carry0;
 247     h5 = (guint32) a;
 248
 249     *pclo = ((guint64)h1) << 32 | h0;
 250     *pcmi = ((guint64)h3) << 32 | h2;
 251     *pchi = ((guint64)h5) << 32 | h4;
 252 }
 253
 254 /* multiplication c(128bit) = a(96bit) * b(32bit) */
 255 DECINLINE static void mult96by32to128(guint32 alo, guint32 ami, guint32 ahi,
 256                                       guint32 factor,
 257                                       guint64* pclo, guint64* pchi)
 258 {
 259     guint64 a;
 260     guint32 h0, h1;
 261
 262     a = ((guint64)alo) * factor;
 263     h0 = (guint32) a;
 264
 265     a >>= 32;
 266     a += ((guint64)ami) * factor;
 267     h1 = (guint32) a;
 268
 269     a >>= 32;
 270     a += ((guint64)ahi) * factor;
 271
 272     *pclo = ((guint64)h1) << 32 | h0;
 273     *pchi = a;
 274 }
 275
 276 /* multiplication c(128bit) *= b(32bit) */
 277 DECINLINE static int mult128by32(guint64* pclo, guint64* pchi, guint32 factor, int roundBit)
 278 {
 279     guint64 a;
 280     guint32 h0, h1;
 281
 282     a = ((guint64)(guint32)(*pclo)) * factor;
 283     if (roundBit) a += factor / 2;
 284     h0 = (guint32) a;
 285
 286     a >>= 32;
 287     a += (*pclo >> 32) * factor;
 288     h1 = (guint32) a;
 289
 290     *pclo = ((guint64)h1) << 32 | h0;
 291
 292     a >>= 32;
 293     a += ((guint64)(guint32)(*pchi)) * factor;
 294     h0 = (guint32) a;
 295
 296     a >>= 32;
 297     a += (*pchi >> 32) * factor;
 298     h1 = (guint32) a;
 299
 300     *pchi = ((guint64)h1) << 32 | h0;
 301
 302     return ((a >> 32) == 0) ? DECIMAL_SUCCESS : DECIMAL_OVERFLOW;
 303 }
 304
 305 DECINLINE static int mult128DecadeFactor(guint64* pclo, guint64* pchi, int powerOfTen)
 306 {
 307     int idx, rc;
 308
 309     while (powerOfTen > 0) {
 310         idx = (powerOfTen >= DECIMAL_MAX_INTFACTORS) ? DECIMAL_MAX_INTFACTORS : powerOfTen;
 311         powerOfTen -= idx;
 312         rc = mult128by32(pclo, pchi, constantsDecadeInt32Factors[idx], 0);
 313         if (rc != DECIMAL_SUCCESS) return rc;
 314     }
 315     return DECIMAL_SUCCESS;
 316 }
 317
 318 /* division: x(128bit) /= factor(32bit)
 319    returns roundBit */
 320 DECINLINE static int div128by32(guint64* plo, guint64* phi, guint32 factor, guint32* pRest)
 321 {
 322     guint64 a, b, c, h;
 323
 324     h = *phi;
 325     a = (guint32)(h >> 32);
 326     b = a / factor;
 327     a -= b * factor;
 328     a <<= 32;
 329     a |= (guint32) h;
 330     c = a / factor;
 331     a -= c * factor;
 332     a <<= 32;
 333     *phi = b << 32 | (guint32)c;
 334
 335     h = *plo;
 336     a |= (guint32)(h >> 32);
 337     b = a / factor;
 338     a -= b * factor;
 339     a <<= 32;
 340     a |= (guint32) h;
 341     c = a / factor;
 342     a -= c * factor;
 343     *plo = b << 32 | (guint32)c;
 344
 345     if (pRest) *pRest = (guint32) a;
 346
 347     a <<= 1;
 348     return (a > factor || (a == factor && (c & 1) == 1)) ? 1 : 0;
 349 }
 350
 351 /* division: x(192bit) /= factor(32bit)
 352    no rest and no rounding*/
 353 DECINLINE static void div192by32(guint64* plo, guint64* pmi, guint64* phi,
 354                                  guint32 factor)
 355 {
 356     guint64 a, b, c, h;
 357
 358     h = *phi;
 359     a = (guint32)(h >> 32);
 360     b = a / factor;
 361     a -= b * factor;
 362     a <<= 32;
 363     a |= (guint32) h;
 364     c = a / factor;
 365     a -= c * factor;
 366     a <<= 32;
 367     *phi = b << 32 | (guint32)c;
 368
 369     h = *pmi;
 370     a |= (guint32)(h >> 32);
 371     b = a / factor;
 372     a -= b * factor;
 373     a <<= 32;
 374     a |= (guint32) h;
 375     c = a / factor;
 376     a -= c * factor;
 377     a <<= 32;
 378     *pmi = b << 32 | (guint32)c;
 379
 380     h = *plo;
 381     a |= (guint32)(h >> 32);
 382     b = a / factor;
 383     a -= b * factor;
 384     a <<= 32;
 385     a |= (guint32) h;
 386     c = a / factor;
 387     a -= c * factor;
 388     a <<= 32;
 389     *plo = b << 32 | (guint32)c;
 390 }
 391
 392 /* returns upper 32bit for a(192bit) /= b(32bit)
 393    a will contain remainder */
 394 static guint32 div192by96to32withRest(guint64* palo, guint64* pami, guint64* pahi,
 395                                       guint32 blo, guint32 bmi, guint32 bhi)
 396 {
 397     guint64 rlo, rmi, rhi; /* remainder */
 398     guint64 tlo, thi; /* term */
 399     guint32 c;
 400
 401     rlo = *palo; rmi = *pami; rhi = *pahi;
 402     if (rhi >= (((guint64)bhi) << 32)) {
 403         c = LIT_GUINT32(0xFFFFFFFF);
 404     } else {
 405         c = (guint32) (rhi / bhi);
 406     }
 407     mult96by32to128(blo, bmi, bhi, c, &tlo, &thi);
 408     sub192(rlo, rmi, rhi, 0, tlo, thi, &rlo, &rmi, &rhi);
 409     while (((gint64)rhi) < 0) {
 410         c--;
 411         add192(rlo, rmi, rhi, 0, (((guint64)bmi)<<32) | blo, bhi, &rlo, &rmi, &rhi);
 412     }
 413     *palo = rlo ; *pami = rmi ; *pahi = rhi;
 414
 415     POSTCONDITION(rhi >> 32 == 0);
 416
 417     return c;
 418 }
 419
 420 /* c(128bit) = a(192bit) / b(96bit)
 421    b must be >= 2^95 */
 422 static void div192by96to128(guint64 alo, guint64 ami, guint64 ahi,
 423                             guint32 blo, guint32 bmi, guint32 bhi,
 424                             guint64* pclo, guint64* pchi)
 425 {
 426     guint64 rlo, rmi, rhi; /* remainder */
 427     guint32 h, c;
 428
 429     PRECONDITION(ahi < (((guint64)bhi) << 32 | bmi)
 430         || (ahi == (((guint64)bhi) << 32 | bmi) && (ami >> 32) > blo));
 431
 432     /* high 32 bit*/
 433     rlo = alo; rmi = ami; rhi = ahi;
 434     h = div192by96to32withRest(&rlo, &rmi, &rhi, blo, bmi, bhi);
 435
 436     /* mid 32 bit*/
 437     rhi = (rhi << 32) | (rmi >> 32); rmi = (rmi << 32) | (rlo >> 32); rlo <<= 32;
 438     *pchi = (((guint64)h) << 32) | div192by96to32withRest(&rlo, &rmi, &rhi, blo, bmi, bhi);
 439
 440     /* low 32 bit */
 441     rhi = (rhi << 32) | (rmi >> 32); rmi = (rmi << 32) | (rlo >> 32); rlo <<= 32;
 442     h = div192by96to32withRest(&rlo, &rmi, &rhi, blo, bmi, bhi);
 443
 444     /* estimate lowest 32 bit (two last bits may be wrong) */
 445     if (rhi >= bhi) {
 446         c = LIT_GUINT32(0xFFFFFFFF);
 447     } else {
 448         rhi <<= 32;
 449         c = (guint32) (rhi / bhi);
 450     }
 451     *pclo = (((guint64)h) << 32) | c;
 452 }
 453
 454 DECINLINE static void roundUp128(guint64* pclo, guint64* pchi) {
 455     if (++(*pclo) == 0) ++(*pchi);
 456 }
 457
 458 static int normalize128(guint64* pclo, guint64* pchi, int* pScale,
 459                         int roundFlag, int roundBit)
 460 {
 461     guint32 overhang = (guint32)(*pchi >> 32);
 462     int scale = *pScale;
 463     int deltaScale;
 464
 465     while (overhang != 0) {
 466         for (deltaScale = 1; deltaScale < DECIMAL_MAX_INTFACTORS; deltaScale++)
 467         {
 468             if (overhang < constantsDecadeInt32Factors[deltaScale]) break;
 469         }
 470
 471         scale -= deltaScale;
 472         if (scale < 0) return DECIMAL_OVERFLOW;
 473
 474         roundBit = div128by32(pclo, pchi, constantsDecadeInt32Factors[deltaScale], 0);
 475
 476         overhang = (guint32)(*pchi >> 32);
 477         if (roundFlag && roundBit && *pclo == (guint64)-1 && (gint32)*pchi == (gint32)-1) {
 478             overhang = 1;
 479         }
 480     }
 481
 482     *pScale = scale;
 483
 484     if (roundFlag && roundBit) {
 485         roundUp128(pclo, pchi);
 486         TEST((*pchi >> 32) == 0);
 487     }
 488
 489     return DECIMAL_SUCCESS;
 490 }
 491
 492 DECINLINE static int maxLeftShift(/*[In, Out]*/decimal_repr* pA)
 493 {
 494     guint64 lo64 = (((guint64)(pA->mid32)) << 32) | pA->lo32;
 495     guint32 hi32 = pA->hi32;
 496     int shift;
 497
 498     for (shift = 0; ((gint32)hi32) >= 0 && shift < 96; shift++) {
 499         hi32 <<= 1;
 500         if (((gint64)lo64) < 0) hi32++;
 501         lo64 <<= 1;
 502     }
 503
 504     pA->lo32 = (guint32) lo64;
 505     pA->mid32 = (guint32)(lo64>>32);
 506     pA->hi32 = hi32;
 507
 508     return shift;
 509 }
 510
 511 DECINLINE static void rshift128(guint64* pclo, guint64* pchi)
 512 {
 513     *pclo >>= 1;
 514     if (*pchi & 1) *pclo |= LIT_GUINT64_HIGHBIT;
 515     *pchi >>= 1;
 516 }
 517
 518 DECINLINE static void lshift96(guint32* pclo, guint32* pcmid, guint32* pchi)
 519 {
 520     *pchi <<= 1;
 521     if (*pcmid & LIT_GUINT32_HIGHBIT) (*pchi)++;
 522     *pcmid <<= 1;
 523     if (*pclo & LIT_GUINT32_HIGHBIT) (*pcmid)++;
 524     *pclo <<= 1;
 525 }
 526
 527 DECINLINE static void lshift128(guint64* pclo, guint64* pchi)
 528 {
 529     *pchi <<= 1;
 530     if (*pclo & LIT_GUINT64_HIGHBIT) (*pchi)++;
 531     *pclo <<= 1;
 532 }
 533
 534 DECINLINE static void rshift192(guint64* pclo, guint64* pcmi, guint64* pchi)
 535 {
 536     *pclo >>= 1;
 537     if (*pcmi & 1) *pclo |= LIT_GUINT64_HIGHBIT;
 538     *pcmi >>= 1;
 539     if (*pchi & 1) *pcmi |= LIT_GUINT64_HIGHBIT;
 540     *pchi >>= 1;
 541 }
 542
 543 /* returns log2(a) or DECIMAL_LOG_NEGINF for a = 0 */
 544 DECINLINE static int log2_32(guint32 a)
 545 {
 546     int tlog2 = 0;
 547
 548     if (a == 0) return DECIMAL_LOG_NEGINF;
 549
 550     if ((a >> 16) != 0) {
 551         a >>= 16;
 552         tlog2 += 16;
 553     }
 554     if ((a >> 8) != 0) {
 555         a >>= 8;
 556         tlog2 += 8;
 557     }
 558     if ((a >> 4) != 0) {
 559         a >>= 4;
 560         tlog2 += 4;
 561     }
 562     if ((a >> 2) != 0) {
 563         a >>= 2;
 564         tlog2 += 2;
 565     }
 566     if ((a >> 1) != 0) {
 567         a >>= 1;
 568         tlog2 += 1;
 569     }
 570     tlog2 += (int) a;
 571
 572     return tlog2;
 573 }
 574
 575 /* returns log2(a) or DECIMAL_LOG_NEGINF for a = 0 */
 576 DECINLINE static int log2_64(guint64 a)
 577 {
 578     int tlog2 = 0;
 579
 580     if (a == 0) return DECIMAL_LOG_NEGINF;
 581
 582     if ((a >> 32) != 0) {
 583         a >>= 32;
 584         tlog2 += 32;
 585     }
 586     if ((a >> 16) != 0) {
 587         a >>= 16;
 588         tlog2 += 16;
 589     }
 590     if ((a >> 8) != 0) {
 591         a >>= 8;
 592         tlog2 += 8;
 593     }
 594     if ((a >> 4) != 0) {
 595         a >>= 4;
 596         tlog2 += 4;
 597     }
 598     if ((a >> 2) != 0) {
 599         a >>= 2;
 600         tlog2 += 2;
 601     }
 602     if ((a >> 1) != 0) {
 603         a >>= 1;
 604         tlog2 += 1;
 605     }
 606     tlog2 += (int) a;
 607
 608     return tlog2;
 609 }
 610
 611 /* returns log2(a) or DECIMAL_LOG_NEGINF for a = 0 */
 612 DECINLINE static int log2_128(guint64 alo, guint64 ahi)
 613 {
 614     if (ahi == 0) return log2_64(alo);
 615     else return log2_64(ahi) + 64;
 616 }
 617
 618 /* returns a upper limit for log2(a) considering scale */
 619 DECINLINE static int log2withScale_128(guint64 alo, guint64 ahi, int scale)
 620 {
 621     int tlog2 = log2_128(alo, ahi);
 622     if (tlog2 < 0) tlog2 = 0;
 623     return tlog2 - (scale * 33219) / 10000;
 624 }
 625
 626 DECINLINE static int pack128toDecimal(/*[Out]*/decimal_repr* pA, guint64 alo, guint64 ahi,
 627                                       int scale, int sign)
 628 {
 629     PRECONDITION((ahi >> 32) == 0);
 630     PRECONDITION(sign == 0 || sign == 1);
 631     PRECONDITION(scale >= 0 && scale <= DECIMAL_MAX_SCALE);
 632
 633     if (scale < 0 || scale > DECIMAL_MAX_SCALE || (ahi >> 32) != 0) {
 634         return DECIMAL_OVERFLOW;
 635     }
 636
 637     pA->lo32 = (guint32) alo;
 638     pA->mid32 = (guint32) (alo >> 32);
 639     pA->hi32 = (guint32) ahi;
 640     pA->signscale.sign = sign;
 641     pA->signscale.scale = scale;
 642
 643     return DECIMAL_SUCCESS;
 644 }
 645
 646 DECINLINE static int adjustScale128(guint64* palo, guint64* pahi, int deltaScale)
 647 {
 648     int idx, rc;
 649
 650     if (deltaScale < 0) {
 651         deltaScale *= -1;
 652         if (deltaScale > DECIMAL_MAX_SCALE) return DECIMAL_INTERNAL_ERROR;
 653         while (deltaScale > 0) {
 654             idx = (deltaScale > DECIMAL_MAX_INTFACTORS) ? DECIMAL_MAX_INTFACTORS : deltaScale;
 655             deltaScale -= idx;
 656             div128by32(palo, pahi, constantsDecadeInt32Factors[idx], 0);
 657         }
 658     } else if (deltaScale > 0) {
 659         if (deltaScale > DECIMAL_MAX_SCALE) return DECIMAL_INTERNAL_ERROR;
 660         while (deltaScale > 0) {
 661             idx = (deltaScale > DECIMAL_MAX_INTFACTORS) ? DECIMAL_MAX_INTFACTORS : deltaScale;
 662             deltaScale -= idx;
 663             rc = mult128by32(palo, pahi, constantsDecadeInt32Factors[idx], 0);
 664             if (rc != DECIMAL_SUCCESS) return rc;
 665         }
 666     }
 667
 668     return DECIMAL_SUCCESS;
 669 }
 670
 671 /* input: c * 10^-(*pScale) * 2^-exp
 672    output: c * 10^-(*pScale) with
 673    minScale <= *pScale <= maxScale and (chi >> 32) == 0 */
 674 DECINLINE static int rescale128(guint64* pclo, guint64* pchi, int* pScale, int texp,
 675                                 int minScale, int maxScale, int roundFlag)
 676 {
 677     guint32 factor, overhang;
 678     int scale, i, rc, roundBit = 0;
 679
 680     PRECONDITION(texp >= 0);
 681
 682     scale = *pScale;
 683
 684     if (texp > 0) {
 685         /* reduce exp */
 686         while (texp > 0 && scale <= maxScale) {
 687             overhang = (guint32)(*pchi >> 32);
 688             while (texp > 0 && ((*pclo & 1) == 0 || overhang > (2<<DECIMAL_MAX_INTFACTORS))) {
 689                 if (--texp == 0) roundBit = (int)(*pclo & 1);
 690                 rshift128(pclo, pchi);
 691                 overhang = (guint32)(*pchi >> 32);
 692             }
 693
 694             if (texp > DECIMAL_MAX_INTFACTORS) i = DECIMAL_MAX_INTFACTORS;
 695             else i = texp;
 696             if (scale + i > maxScale) i = maxScale - scale;
 697             if (i == 0) break;
 698             texp -= i;
 699             scale += i;
 700             factor = constantsDecadeInt32Factors[i] >> i; /* 10^i/2^i=5^i */
 701             mult128by32(pclo, pchi, factor, 0);
 702     /*printf("3: %.17e\n", (((double)chi) * pow(2,64) + clo) * pow(10, -scale) * pow(2, -texp));*/
 703         }
 704
 705         while (texp > 0) {
 706             if (--texp == 0) roundBit = (int)(*pclo & 1);
 707             rshift128(pclo, pchi);
 708         }
 709     }
 710
 711     TEST(texp == 0);
 712
 713     while (scale > maxScale) {
 714         i = scale - maxScale;
 715         if (i > DECIMAL_MAX_INTFACTORS) i = DECIMAL_MAX_INTFACTORS;
 716         scale -= i;
 717         roundBit = div128by32(pclo, pchi, constantsDecadeInt32Factors[i], 0);
 718     }
 719
 720     while (scale < minScale) {
 721         if (!roundFlag) roundBit = 0;
 722         i = minScale - scale;
 723         if (i > DECIMAL_MAX_INTFACTORS) i = DECIMAL_MAX_INTFACTORS;
 724         scale += i;
 725         rc = mult128by32(pclo, pchi, constantsDecadeInt32Factors[i], roundBit);
 726         if (rc != DECIMAL_SUCCESS) return rc;
 727         roundBit = 0;
 728     }
 729
 730     TEST(scale >= 0 && scale <= DECIMAL_MAX_SCALE);
 731
 732     *pScale = scale;
 733
 734     return normalize128(pclo, pchi, pScale, roundFlag, roundBit);
 735 }
 736
 737 /* performs a += b */
 738 gint32 mono_decimalIncr(/*[In, Out]*/decimal_repr* pA, /*[In]*/decimal_repr* pB)
 739 {
 740     MONO_ARCH_SAVE_REGS;
 741
 742     guint64 alo, ahi, blo, bhi;
 743     int log2A, log2B, log2Result, log10Result, rc;
 744     int subFlag, sign, scaleA, scaleB;
 745
 746     DECTO128(pA, alo, ahi);
 747     DECTO128(pB, blo, bhi);
 748
 749     sign = pA->signscale.sign;
 750     subFlag = sign - (int)pB->signscale.sign;
 751     scaleA = pA->signscale.scale;
 752     scaleB = pB->signscale.scale;
 753     if (scaleA == scaleB) {
 754         /* same scale, that's easy */
 755         if (subFlag) {
 756             sub128(alo, ahi, blo, bhi, &alo, &ahi);
 757             if (ahi & LIT_GUINT64_HIGHBIT) {
 758                 alo--;
 759                 alo = ~alo;
 760                 if (alo == 0) ahi--;
 761                 ahi = ~ahi;
 762                 sign = !sign;
 763             }
 764         } else {
 765             add128(alo, ahi, blo, bhi, &alo, &ahi);
 766         }
 767         rc = normalize128(&alo, &ahi, &scaleA, 1, 0);
 768     } else {
 769         /* scales must be adjusted */
 770         /* Estimate log10 and scale of result for adjusting scales */
 771         log2A = log2withScale_128(alo, ahi, scaleA);
 772         log2B = log2withScale_128(blo, bhi, scaleB);
 773         log2Result = (log2A >= log2B) ? log2A : log2B;
 774         if (!subFlag) log2Result++; /* result can have one bit more */
 775         log10Result = (log2Result * 1000) / 3322 + 1;
 776         /* we will calculate in 128bit, so we may need to adjust scale */
 777         if (scaleB > scaleA) scaleA = scaleB;
 778         if (scaleA + log10Result > DECIMAL_MAX_SCALE + 7) {
 779             /* this may not fit in 128bit, so limit it */
 780             scaleA = DECIMAL_MAX_SCALE + 7 - log10Result;
 781         }
 782
 783         rc = adjustScale128(&alo, &ahi, scaleA - (int)pA->signscale.scale);
 784         if (rc != DECIMAL_SUCCESS) return rc;
 785         rc = adjustScale128(&blo, &bhi, scaleA - scaleB);
 786         if (rc != DECIMAL_SUCCESS) return rc;
 787
 788         if (subFlag) {
 789             sub128(alo, ahi, blo, bhi, &alo, &ahi);
 790             if (ahi & LIT_GUINT64_HIGHBIT) {
 791                 alo--;
 792                 alo = ~alo;
 793                 if (alo == 0) ahi--;
 794                 ahi = ~ahi;
 795                 sign = !sign;
 796             }
 797         } else {
 798             add128(alo, ahi, blo, bhi, &alo, &ahi);
 799         }
 800
 801         if (rc != DECIMAL_SUCCESS) return rc;
 802
 803         rc = rescale128(&alo, &ahi,&scaleA, 0, 0, DECIMAL_MAX_SCALE, 1);
 804     }
 805
 806     if (rc != DECIMAL_SUCCESS) return rc;
 807
 808     return pack128toDecimal(pA, alo, ahi, scaleA, sign);
 809 }
 810
 811 /* performs a += factor * constants[idx] */
 812 static int incMultConstant128(guint64* palo, guint64* pahi, int idx, int factor)
 813 {
 814     guint64 blo, bhi, h;
 815
 816     PRECONDITION(idx >= 0 && idx <= DECIMAL_MAX_SCALE);
 817     PRECONDITION(factor > 0 && factor <= 9);
 818
 819     blo = dec128decadeFactors[idx].lo;
 820     h = bhi = dec128decadeFactors[idx].hi;
 821     if (factor != 1) {
 822         mult128by32(&blo, &bhi, factor, 0);
 823         if (h > bhi) return DECIMAL_OVERFLOW;
 824     }
 825     h = *pahi;
 826     add128(*palo, *pahi, blo, bhi, palo, pahi);
 827     if (h > *pahi) return DECIMAL_OVERFLOW;
 828     return DECIMAL_SUCCESS;
 829 }
 830
 831 DECINLINE static void div128DecadeFactor(guint64* palo, guint64* pahi, int powerOfTen)
 832 {
 833     int idx, roundBit = 0;
 834
 835     while (powerOfTen > 0) {
 836         idx = (powerOfTen > DECIMAL_MAX_INTFACTORS) ? DECIMAL_MAX_INTFACTORS : powerOfTen;
 837         powerOfTen -= idx;
 838         roundBit = div128by32(palo, pahi, constantsDecadeInt32Factors[idx], 0);
 839     }
 840
 841     if (roundBit) roundUp128(palo, pahi);
 842 }
 843
 844 /* calc significant digits of mantisse */
 845 DECINLINE static int calcDigits(guint64 alo, guint64 ahi)
 846 {
 847     int tlog2 = 0;
 848     int tlog10;
 849
 850     if (ahi == 0) {
 851         if (alo == 0) {
 852             return 0; /* zero has no signficant digits */
 853         } else {
 854             tlog2 = log2_64(alo);
 855         }
 856     } else {
 857         tlog2 = 64 + log2_64(ahi);
 858     }
 859
 860     tlog10 = (tlog2 * 1000) / 3322;
 861     /* we need an exact floor value of log10(a) */
 862     if (dec128decadeFactors[tlog10].hi > ahi
 863             || (dec128decadeFactors[tlog10].hi == ahi
 864                     && dec128decadeFactors[tlog10].lo > alo)) {
 865         --tlog10;
 866     }
 867     return tlog10+1;
 868 }
 869
 870 gint32 mono_double2decimal(/*[Out]*/decimal_repr* pA, double val, gint32 digits)
 871 {
 872     guint64 alo, ahi;
 873     guint64* p = (guint64*)(&val);
 874     int sigDigits, sign, texp, rc, scale;
 875     guint16 k;
 876
 877     PRECONDITION(digits <= 15);
 878
 879     sign = ((*p & LIT_GUINT64_HIGHBIT) != 0) ? 1 : 0;
 880     k = ((guint16)((*p) >> 52)) & 0x7FF;
 881     alo = (*p & LIT_GUINT64(0xFFFFFFFFFFFFF)) | LIT_GUINT64(0x10000000000000);
 882     ahi = 0;
 883
 884     texp = (k & 0x7FF) - 0x3FF;
 885     if (k == 0x7FF || texp >= 96) return DECIMAL_OVERFLOW; /* NaNs, SNaNs, Infinities or >= 2^96 */
 886     if (k == 0 || texp <= -94) { /* Subnormals, Zeros or < 2^-94 */
 887         DECINIT(pA); /* return zero */
 888         return DECIMAL_SUCCESS;
 889     }
 890
 891     texp -= 52;
 892     if (texp > 0) {
 893         for (; texp > 0; texp--) {
 894             lshift128(&alo, &ahi);
 895         }
 896     }
 897
 898     scale = 0;
 899     rc = rescale128(&alo, &ahi, &scale, -texp, 0, DECIMAL_MAX_SCALE, 0);
 900     if (rc != DECIMAL_SUCCESS) return rc;
 901
 902     sigDigits = calcDigits(alo, ahi);
 903     /* too much digits, then round */
 904     if (sigDigits > digits) {
 905         div128DecadeFactor(&alo, &ahi, sigDigits - digits);
 906         scale -= sigDigits - digits;
 907         /* check value, may be 10^(digits+1) caused by rounding */
 908         if (ahi == dec128decadeFactors[digits].hi
 909             && alo == dec128decadeFactors[digits].lo) {
 910             div128by32(&alo, &ahi, 10, 0);
 911             scale--;
 912         }
 913         if (scale < 0) {
 914             rc = mult128DecadeFactor(&alo, &ahi, -scale);
 915             if (rc != DECIMAL_SUCCESS) return rc;
 916             scale = 0;
 917         }
 918     }
 919
 920     return pack128toDecimal(pA, alo, ahi, scale, sign);
 921 }
 922
 923 /**
 924  * converts a digit string to decimal
 925  * The significant digits must be passed as an integer in buf !
 926  *
 927  * 1. Example:
 928  *   if you want to convert the number 123.456789012345678901234 to decimal
 929  *     buf := "123456789012345678901234"
 930  *     decrDecimal := 3
 931  *     sign := 0
 932  *
 933  * 2. Example:
 934  *   you want to convert -79228162514264337593543950335 to decimal
 935  *     buf := "79228162514264337593543950335"
 936  *     decrDecimal := 29
 937  *     sign := 1
 938  *
 939  * 3. Example:
 940  *   you want to convert -7922816251426433759354395033.250000000000001 to decimal
 941  *     buf := "7922816251426433759354395033250000000000001"
 942  *     decrDecimal := 29
 943  *     sign := 1
 944  *     returns (decimal)-7922816251426433759354395033.3
 945  *
 946  * 4. Example:
 947  *   you want to convert -7922816251426433759354395033.250000000000000 to decimal
 948  *     buf := "7922816251426433759354395033250000000000000"
 949  *     decrDecimal := 29
 950  *     sign := 1
 951  *     returns (decimal)-7922816251426433759354395033.2
 952  *
 953  * 5. Example:
 954  *   you want to convert -7922816251426433759354395033.150000000000000 to decimal
 955  *     buf := "7922816251426433759354395033150000000000000"
 956  *     decrDecimal := 29
 957  *     sign := 1
 958  *     returns (decimal)-7922816251426433759354395033.2
 959  *
 960  * Uses banker's rule for rounding if there are more digits than can be
 961  * represented by the significant
 962  */
 963 gint32 mono_string2decimal(/*[Out]*/decimal_repr* pA, MonoString* str, gint32 decrDecimal, gint32 sign)
 964 {
 965     gushort *buf = mono_string_chars(str);
 966     gushort *p;
 967     guint64 alo, ahi;
 968     int n, rc, i, len, sigLen = -1, firstNonZero;
 969     int scale, roundBit = 0;
 970
 971     alo = ahi = 0;
 972     DECINIT(pA);
 973
 974     for (p = buf, len = 0; *p != 0; len++, p++) { }
 975
 976     for (p = buf, i = 0; *p != 0; i++, p++) {
 977         n = *p - '0';
 978         if (n < 0 || n > 9) {
 979             return DECIMAL_INVALID_CHARACTER;
 980         }
 981         if (n) {
 982             if (sigLen < 0) {
 983                 firstNonZero = i;
 984                 sigLen = (len - firstNonZero > DECIMAL_MAX_SCALE+1)
 985                     ? DECIMAL_MAX_SCALE+1+firstNonZero : len;
 986                 if (decrDecimal > sigLen+1) return DECIMAL_OVERFLOW;
 987             }
 988             if (i >= sigLen) break;
 989             rc = incMultConstant128(&alo, &ahi, sigLen - 1 - i, n);
 990             if (rc != DECIMAL_SUCCESS) {
 991                 return rc;
 992             }
 993         }
 994     }
 995
 996     scale = sigLen - decrDecimal;
 997
 998     if (i < len) { /* too much digits, we must round */
 999         n = buf[i] - '0';
1000         if (n < 0 || n > 9) {
1001             return DECIMAL_INVALID_CHARACTER;
1002         }
1003         if (n > 5) roundBit = 1;
1004         else if (n == 5) { /* we must take a nearer look */
1005             n = buf[i-1] - '0';
1006             for (++i; i < len; ++i) {
1007                 if (buf[i] != '0') break; /* we are greater than .5 */
1008             }
1009             if (i < len /* greater than exactly .5 */
1010                 || n % 2 == 1) { /* exactly .5, use banker's rule for rounding */
1011                 roundBit = 1;
1012             }
1013         }
1014     }
1015
1016     if (ahi != 0) {
1017         rc = normalize128(&alo, &ahi, &scale, 1, roundBit);
1018         if (rc != DECIMAL_SUCCESS) return rc;
1019     }
1020
1021     if (alo == 0 && ahi == 0) {
1022         DECINIT(pA);
1023         return DECIMAL_SUCCESS;
1024     } else {
1025         return pack128toDecimal(pA, alo, ahi, sigLen - decrDecimal, sign);
1026     }
1027 }
1028
1029 /**
1030  * returns minimal number of digit string to represent decimal
1031  * No leading or trailing zeros !
1032  * Examples:
1033  * *pA == 0            =>   buf = "", *pDecPos = 1, *pSign = 0
1034  * *pA == 12.34        =>   buf = "1234", *pDecPos = 2, *pSign = 0
1035  * *pA == -1000.0000   =>   buf = "1", *pDecPos = 4, *pSign = 1
1036  * *pA == -0.00000076  =>   buf = "76", *pDecPos = -6, *pSign = 0
1037  *
1038  * Parameters:
1039  *    pA         decimal instance to convert
1040  *    digits     < 0: use decimals instead
1041  *               = 0: gets mantisse as integer
1042  *               > 0: gets at most <digits> digits, rounded according to banker's rule if necessary
1043  *    decimals   only used if digits < 0
1044  *               >= 0: number of decimal places
1045  *    buf        pointer to result buffer
1046  *    bufSize    size of buffer
1047  *    pDecPos    receives insert position of decimal point relative to start of buffer
1048  *    pSign      receives sign
1049  */
1050
1051 gint32 mono_decimal2string(/*[In]*/decimal_repr* pA, gint32 digits, gint32 decimals,
1052                                    MonoArray* pArray, gint32 bufSize, gint32* pDecPos, gint32* pSign)
1053 {
1054     MONO_ARCH_SAVE_REGS;
1055
1056     guint16 tmp[41];
1057     guint16 *buf = (guint16*) mono_array_addr(pArray, guint16, 0);
1058     guint16 *q, *p = tmp;
1059     decimal_repr aa;
1060     guint64 alo, ahi;
1061     guint32 rest;
1062     gint32 sigDigits, d;
1063     int i, scale, len;
1064
1065     scale = pA->signscale.scale;
1066     DECTO128(pA, alo, ahi);
1067     sigDigits = calcDigits(alo, ahi); /* significant digits */
1068
1069     /* calc needed digits (without leading or trailing zeros) */
1070     d = (digits == 0) ? sigDigits : digits;
1071     if (d < 0) { /* use decimals ? */
1072         if (0 <= decimals && decimals < scale) {
1073             d = sigDigits - scale + decimals;
1074         } else {
1075             d = sigDigits; /* use all you can get */
1076         }
1077     }
1078
1079     if (sigDigits > d) { /* we need to round decimal number */
1080         DECCOPY(&aa, pA);
1081         aa.signscale.scale = DECIMAL_MAX_SCALE;
1082         mono_decimalRound(&aa, DECIMAL_MAX_SCALE - sigDigits + d);
1083         DECTO128(&aa, alo, ahi);
1084         sigDigits += calcDigits(alo, ahi) - d;
1085     }
1086
1087     len = 0;
1088     if (d > 0) {
1089         /* get digits starting from the tail */
1090         for (; (alo != 0 || ahi != 0) && len < 40; len++) {
1091             div128by32(&alo, &ahi, 10, &rest);
1092             *p++ = '0' + (char) rest;
1093         }
1094     }
1095     *p = 0;
1096
1097     if (len >= bufSize) return DECIMAL_BUFFER_OVERFLOW;
1098
1099     /* now we have the minimal count of digits,
1100        extend to wished count of digits or decimals */
1101     q = buf;
1102     if (digits >= 0) { /* count digits */
1103         if (digits >= bufSize) return DECIMAL_BUFFER_OVERFLOW;
1104         if (len == 0) {
1105             /* zero or rounded to zero */
1106             *pDecPos = 1;
1107         } else {
1108             /* copy significant digits */
1109             for (i = 0; i < len; i++) {
1110                 *q++ = *(--p);
1111             }
1112             *pDecPos = sigDigits - scale;
1113         }
1114         /* add trailing zeros */
1115         for (i = len; i < digits; i++) {
1116             *q++ = '0';
1117         }
1118     } else { /* count decimals */
1119         if (scale >= sigDigits) { /* add leading zeros */
1120             if (decimals+2 >= bufSize) return DECIMAL_BUFFER_OVERFLOW;
1121             *pDecPos = 1;
1122             for (i = 0; i <= scale - sigDigits; i++) {
1123                 *q++ = '0';
1124             }
1125         } else {
1126             if (sigDigits - scale + decimals+1 >= bufSize) return DECIMAL_BUFFER_OVERFLOW;
1127             *pDecPos = sigDigits - scale;
1128         }
1129         /* copy significant digits */
1130         for (i = 0; i < len; i++) {
1131             *q++ = *(--p);
1132         }
1133         /* add trailing zeros */
1134         for (i = scale; i < decimals; i++) {
1135             *q++ = '0';
1136         }
1137     }
1138     *q = 0;
1139
1140     *pSign = (sigDigits > 0) ? pA->signscale.sign : 0; /* zero has positive sign */
1141
1142     return DECIMAL_SUCCESS;
1143 }
1144
1145 /**
1146  * converts a decimal to an UInt64 without rounding
1147  */
1148 gint32 mono_decimal2UInt64(/*[In]*/decimal_repr* pA, guint64* pResult)
1149 {
1150     MONO_ARCH_SAVE_REGS;
1151
1152     guint64 alo, ahi;
1153     int scale;
1154
1155     DECTO128(pA, alo, ahi);
1156     scale = pA->signscale.scale;
1157     if (scale > 0) {
1158         div128DecadeFactor(&alo, &ahi, scale);
1159     }
1160
1161     /* overflow if integer too large or < 0 */
1162     if (ahi != 0 || (alo != 0 && pA->signscale.sign)) return DECIMAL_OVERFLOW;
1163
1164     *pResult = alo;
1165     return DECIMAL_SUCCESS;
1166 }
1167
1168 /**
1169  * converts a decimal to an Int64 without rounding
1170  */
1171 gint32 mono_decimal2Int64(/*[In]*/decimal_repr* pA, gint64* pResult)
1172 {
1173     MONO_ARCH_SAVE_REGS;
1174
1175     guint64 alo, ahi;
1176     int sign, scale;
1177
1178     DECTO128(pA, alo, ahi);
1179     scale = pA->signscale.scale;
1180     if (scale > 0) {
1181         div128DecadeFactor(&alo, &ahi, scale);
1182     }
1183
1184     if (ahi != 0) return DECIMAL_OVERFLOW;
1185
1186     sign = pA->signscale.sign;
1187     if (sign && alo != 0) {
1188         if (alo > LIT_GUINT64_HIGHBIT) return DECIMAL_OVERFLOW;
1189         *pResult = (gint64) ~(alo-1);
1190     } else {
1191         if (alo & LIT_GUINT64_HIGHBIT) return DECIMAL_OVERFLOW;
1192         *pResult = (gint64) alo;
1193     }
1194
1195     return DECIMAL_SUCCESS;
1196 }
1197
1198 void mono_decimalFloorAndTrunc(/*[In, Out]*/decimal_repr* pA, gint32 floorFlag)
1199 {
1200     MONO_ARCH_SAVE_REGS;
1201
1202     guint64 alo, ahi;
1203     guint32 factor, rest;
1204     int scale, sign, idx;
1205     int hasRest = 0;
1206
1207     scale = pA->signscale.scale;
1208     if (scale == 0) return; /* nothing to do */
1209
1210     DECTO128(pA, alo, ahi);
1211     sign = pA->signscale.sign;
1212
1213     while (scale > 0) {
1214         idx = (scale > DECIMAL_MAX_INTFACTORS) ? DECIMAL_MAX_INTFACTORS : scale;
1215         factor = constantsDecadeInt32Factors[idx];
1216         scale -= idx;
1217         div128by32(&alo, &ahi, factor, &rest);
1218         hasRest = hasRest || (rest != 0);
1219     }
1220
1221     if (floorFlag && hasRest && sign) { /* floor: if negative, we must round up */
1222         roundUp128(&alo, &ahi);
1223     }
1224
1225     pack128toDecimal(pA, alo, ahi, 0, sign);
1226 }
1227
1228 void mono_decimalRound(/*[In, Out]*/decimal_repr* pA, gint32 decimals)
1229 {
1230     MONO_ARCH_SAVE_REGS;
1231
1232     guint64 alo, ahi;
1233     int scale, sign;
1234
1235     DECTO128(pA, alo, ahi);
1236     scale = pA->signscale.scale;
1237     sign = pA->signscale.sign;
1238     if (scale > decimals) {
1239         div128DecadeFactor(&alo, &ahi, scale - decimals);
1240         scale = decimals;
1241     }
1242
1243     pack128toDecimal(pA, alo, ahi, scale, sign);
1244 }
1245
1246 gint32 mono_decimalMult(/*[In, Out]*/decimal_repr* pA, /*[In]*/decimal_repr* pB)
1247 {
1248     MONO_ARCH_SAVE_REGS;
1249
1250     guint64 low, mid, high;
1251     guint32 factor;
1252     int scale, sign, rc;
1253
1254     mult96by96to192(pA->lo32, pA->mid32, pA->hi32, pB->lo32, pB->mid32, pB->hi32,
1255         &low, &mid, &high);
1256
1257     /* adjust scale and sign */
1258     scale = (int)pA->signscale.scale + (int)pB->signscale.scale;
1259     sign = pA->signscale.sign ^ pB->signscale.sign;
1260
1261     /* first scaling step */
1262     factor = constantsDecadeInt32Factors[DECIMAL_MAX_INTFACTORS];
1263     while (high != 0 || (mid>>32) >= factor) {
1264         if (high < 100) {
1265             factor /= 1000; /* we need some digits for final rounding */
1266             scale -= DECIMAL_MAX_INTFACTORS - 3;
1267         } else {
1268             scale -= DECIMAL_MAX_INTFACTORS;
1269         }
1270
1271         div192by32(&low, &mid, &high, factor);
1272     }
1273
1274     /* second and final scaling */
1275     rc = rescale128(&low, &mid, &scale, 0, 0, DECIMAL_MAX_SCALE, 1);
1276     if (rc != DECIMAL_SUCCESS) return rc;
1277
1278     return pack128toDecimal(pA, low, mid, scale, sign);
1279 }
1280
1281 static int decimalDivSub(/*[In]*/decimal_repr* pA, /*[In]*/decimal_repr* pB,
1282                          guint64* pclo, guint64* pchi, int* pExp)
1283 {
1284     guint64 alo, ami, ahi;
1285     guint64 tlo, tmi, thi;
1286     guint32 blo, bmi, bhi;
1287     int ashift, bshift, extraBit, texp;
1288
1289     ahi = (((guint64)(pA->hi32)) << 32) | pA->mid32;
1290     ami = ((guint64)(pA->lo32)) << 32;
1291     alo = 0;
1292     blo = pB->lo32;
1293     bmi = pB->mid32;
1294     bhi = pB->hi32;
1295
1296     if (blo == 0 && bmi == 0 && bhi == 0) {
1297         return DECIMAL_DIVIDE_BY_ZERO;
1298     }
1299
1300     if (ami == 0 && ahi == 0) {
1301         *pclo = *pchi = 0;
1302         return DECIMAL_FINISHED;
1303     }
1304
1305     /* enlarge dividend to get maximal precision */
1306     for (ashift = 0; (ahi & LIT_GUINT64_HIGHBIT) == 0; ++ashift) {
1307         lshift128(&ami, &ahi);
1308     }
1309
1310     /* ensure that divisor is at least 2^95 */
1311     for (bshift = 0; (bhi & LIT_GUINT32_HIGHBIT) == 0; ++bshift) {
1312         lshift96(&blo, &bmi, &bhi);
1313     }
1314
1315     thi = ((guint64)bhi)<<32 | bmi;
1316     tmi = ((guint64)blo)<<32;
1317     tlo = 0;
1318     if (ahi > thi || (ahi == thi && ami >= tmi)) {
1319         sub192(alo, ami, ahi, tlo, tmi, thi, &alo, &ami, &ahi);
1320         extraBit = 1;
1321     } else {
1322         extraBit = 0;
1323     }
1324
1325     div192by96to128(alo, ami, ahi, blo, bmi, bhi, pclo, pchi);
1326     texp = 128 + ashift - bshift;
1327
1328     if (extraBit) {
1329         rshift128(pclo, pchi);
1330         *pchi += LIT_GUINT64_HIGHBIT;
1331         texp--;
1332     }
1333
1334     /* try loss free right shift */
1335     while (texp > 0 && (*pclo & 1) == 0) {
1336         /* right shift */
1337         rshift128(pclo, pchi);
1338         texp--;
1339     }
1340
1341     *pExp = texp;
1342
1343     return DECIMAL_SUCCESS;
1344 }
1345
1346 gint32 mono_decimalDiv(/*[Out]*/decimal_repr* pC, /*[In]*/decimal_repr* pA, /*[In]*/decimal_repr* pB)
1347 {
1348     MONO_ARCH_SAVE_REGS;
1349
1350     guint64 clo, chi; /* result */
1351     int scale, texp, rc;
1352
1353     rc = decimalDivSub(pA, pB, &clo, &chi, &texp);
1354     if (rc != DECIMAL_SUCCESS) {
1355         if (rc == DECIMAL_FINISHED) rc = DECIMAL_SUCCESS;
1356         return rc;
1357     }
1358
1359     /* adjust scale and sign */
1360     scale = (int)pA->signscale.scale - (int)pB->signscale.scale;
1361
1362     /*test: printf("0: %.17e\n", (((double)chi) * pow(2,64) + clo) * pow(10, -scale) * pow(2, -exp));*/
1363     rc = rescale128(&clo, &chi, &scale, texp, 0, DECIMAL_MAX_SCALE, 1);
1364     if (rc != DECIMAL_SUCCESS) return rc;
1365
1366     return pack128toDecimal(pC, clo, chi, scale, pA->signscale.sign ^ pB->signscale.sign);
1367 }
1368
1369 gint32 mono_decimalIntDiv(/*[Out]*/decimal_repr* pC, /*[In]*/decimal_repr* pA, /*[In]*/decimal_repr* pB)
1370 {
1371     MONO_ARCH_SAVE_REGS;
1372
1373     guint64 clo, chi; /* result */
1374     int scale, texp, rc;
1375
1376     rc = decimalDivSub(pA, pB, &clo, &chi, &texp);
1377     if (rc != DECIMAL_SUCCESS) {
1378         if (rc == DECIMAL_FINISHED) rc = DECIMAL_SUCCESS;
1379         return rc;
1380     }
1381
1382     /* calc scale  */
1383     scale = (int)pA->signscale.scale - (int)pB->signscale.scale;
1384
1385     /* truncate result to integer */
1386     rc = rescale128(&clo, &chi, &scale, texp, 0, 0, 0);
1387     if (rc != DECIMAL_SUCCESS) return rc;
1388
1389     return pack128toDecimal(pC, clo, chi, scale, pA->signscale.sign);
1390 }
1391
1392 /* approximation for log2 of a
1393    If q is the exact value for log2(a), then q <= decimalLog2(a) <= q+1 */
1394 DECINLINE static int decimalLog2(/*[In]*/decimal_repr* pA)
1395 {
1396     int tlog2;
1397     int scale = pA->signscale.scale;
1398
1399     if (pA->hi32 != 0) tlog2 = 64 + log2_32(pA->hi32);
1400     else if (pA->mid32 != 0) tlog2 = 32 + log2_32(pA->mid32);
1401     else tlog2 = log2_32(pA->lo32);
1402
1403     if (tlog2 != DECIMAL_LOG_NEGINF) {
1404         tlog2 -= (scale * 33219) / 10000;
1405     }
1406
1407     return tlog2;
1408 }
1409
1410 DECINLINE static int decimalIsZero(/*[In]*/decimal_repr* pA)
1411 {
1412     return (pA->lo32 == 0 && pA->mid32 == 0 && pA->hi32 == 0);
1413 }
1414
1415 gint32 mono_decimalCompare(/*[In]*/decimal_repr* pA, /*[In]*/decimal_repr* pB)
1416 {
1417     MONO_ARCH_SAVE_REGS;
1418
1419     int log2a, log2b, delta, sign;
1420     decimal_repr aa;
1421
1422     sign = (pA->signscale.sign) ? -1 : 1;
1423
1424     if (pA->signscale.sign ^ pB->signscale.sign) {
1425         return (decimalIsZero(pA) && decimalIsZero(pB)) ? 0 : sign;
1426     }
1427
1428     /* try fast comparison via log2 */
1429     log2a = decimalLog2(pA);
1430     log2b = decimalLog2(pB);
1431     delta = log2a - log2b;
1432      /* decimalLog2 is not exact, so we can say nothing
1433         if abs(delta) <= 1 */
1434     if (delta < -1) return -sign;
1435     if (delta > 1) return sign;
1436
1437     DECCOPY(&aa, pA);
1438     DECNEGATE(&aa);
1439     mono_decimalIncr(&aa, pB);
1440
1441     if (decimalIsZero(&aa)) return 0;
1442
1443     return (aa.signscale.sign) ? 1 : -1;
1444 }
1445
1446 /* d=(-1)^sign * n * 2^(k-52) with sign (1bit), k(11bit), n-2^52(52bit) */
1447 DECINLINE static void buildIEEE754Double(double* pd, int sign, int texp, guint64 mantisse)
1448 {
1449     guint64* p = (guint64*) pd;
1450
1451     PRECONDITION(sign == 0 || sign == 1);
1452     *p = (((guint64)sign) << 63) | (((guint64)((1023+texp)&0x7ff)) << 52) | mantisse;
1453 }
1454
1455 double mono_decimal2double(/*[In]*/decimal_repr* pA)
1456 {
1457     MONO_ARCH_SAVE_REGS;
1458
1459     double d;
1460     guint64 alo, ahi, mantisse;
1461     guint32 overhang, factor, roundBits;
1462     int scale, texp, log5, i;
1463
1464     ahi = (((guint64)(pA->hi32)) << 32) | pA->mid32;
1465     alo = ((guint64)(pA->lo32)) << 32;
1466
1467     /* special case zero */
1468     if (ahi == 0 && alo == 0) return 0.0;
1469
1470     texp = 0;
1471     scale = pA->signscale.scale;
1472
1473     /* transform n * 10^-scale and exp = 0 => m * 2^-exp and scale = 0 */
1474     while (scale > 0) {
1475         while ((ahi & LIT_GUINT64_HIGHBIT) == 0) {
1476             lshift128(&alo, &ahi);
1477             texp++;
1478         }
1479
1480         overhang = (guint32) (ahi >> 32);
1481         if (overhang >= 5) {
1482             /* estimate log5 */
1483             log5 = (log2_32(overhang) * 1000) / 2322; /* ln(5)/ln(2) = 2.3219... */
1484             if (log5 < DECIMAL_MAX_INTFACTORS) {
1485                 /* get maximal factor=5^i, so that overhang / factor >= 1 */
1486                 factor = constantsDecadeInt32Factors[log5] >> log5; /* 5^n = 10^n/2^n */
1487                 i = log5 + overhang / factor;
1488             } else {
1489                 i = DECIMAL_MAX_INTFACTORS; /* we have only constants up to 10^DECIMAL_MAX_INTFACTORS */
1490             }
1491             if (i > scale) i = scale;
1492             factor = constantsDecadeInt32Factors[i] >> i; /* 5^n = 10^n/2^n */
1493             /* n * 10^-scale * 2^-exp => m * 10^-(scale-i) * 2^-(exp+i) with m = n * 5^-i */
1494             div128by32(&alo, &ahi, factor, 0);
1495             scale -= i;
1496             texp += i;
1497         }
1498     }
1499
1500     /* normalize significand (highest bit should be 1) */
1501     while ((ahi & LIT_GUINT64_HIGHBIT) == 0) {
1502         lshift128(&alo, &ahi);
1503         texp++;
1504     }
1505
1506     /* round to nearest even */
1507     roundBits = (guint32)ahi & 0x7ff;
1508     ahi += 0x400;
1509     if ((ahi & LIT_GUINT64_HIGHBIT) == 0) { /* overflow ? */
1510         ahi >>= 1;
1511         texp++;
1512     } else if ((roundBits & 0x400) == 0) ahi &= ~1;
1513
1514     /* 96 bit => 1 implizit bit and 52 explicit bits */
1515     mantisse = (ahi & ~LIT_GUINT64_HIGHBIT) >> 11;
1516
1517     buildIEEE754Double(&d, pA->signscale.sign, -texp+95, mantisse);
1518
1519     return d;
1520 }
1521
1522 /* a *= 10^exp */
1523 gint32 mono_decimalSetExponent(/*[In, Out]*/decimal_repr* pA, gint32 texp)
1524 {
1525     MONO_ARCH_SAVE_REGS;
1526
1527     guint64 alo, ahi;
1528     int rc;
1529     int scale = pA->signscale.scale;
1530
1531     scale -= texp;
1532
1533     if (scale < 0 || scale > DECIMAL_MAX_SCALE) {
1534         DECTO128(pA, alo, ahi);
1535         rc = rescale128(&alo, &ahi, &scale, 0, 0, DECIMAL_MAX_SCALE, 1);
1536         if (rc != DECIMAL_SUCCESS) return rc;
1537         return pack128toDecimal(pA, alo, ahi, scale, pA->signscale.sign);
1538     } else {
1539         pA->signscale.scale = scale;
1540         return DECIMAL_SUCCESS;
1541     }
1542 }