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