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