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