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