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