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