Statistics
| Branch: | Revision:

root / src / common / bigdecimal.cc @ e1750c09

History | View | Annotate | Download (13.4 KB)

1
//==========================================================================
2
//   BIGDECIMAL.CC  - part of
3
//                     OMNeT++/OMNEST
4
//            Discrete System Simulation in C++
5
//
6
//  Author: Tamas Borbely
7
//
8
//==========================================================================
9

    
10
/*--------------------------------------------------------------*
11
  Copyright (C) 2006-2008 OpenSim Ltd.
12

13
  This file is distributed WITHOUT ANY WARRANTY. See the file
14
  `license' for details on this and other legal matters.
15
*--------------------------------------------------------------*/
16

    
17
#include <sstream>
18
#include <assert.h>
19
#include <string.h>
20
#include "opp_ctype.h"
21
#include "platmisc.h"
22
#include "bigdecimal.h"
23
#include "commonutil.h"  //NaN & friends
24

    
25
USING_NAMESPACE
26

    
27
// helpers
28
static inline int64 max(int64 x, int64 y) { return x > y ? x : y; }
29
static inline int64 min(int64 x, int64 y) { return x < y ? x : y; }
30
static inline int64 abs(int64 x) { return x >= 0 ? x : -x; }
31
static inline int sgn(int64 x) { return (x > 0 ? 1 : (x < 0 ? -1 : 0)); }
32

    
33
BigDecimal BigDecimal::Zero(0, 0);
34
BigDecimal BigDecimal::One(1, 0);
35
BigDecimal BigDecimal::MinusOne(-1, 0);
36
BigDecimal BigDecimal::NaN(0, INT_MAX);
37
BigDecimal BigDecimal::PositiveInfinity(1, INT_MAX);
38
BigDecimal BigDecimal::NegativeInfinity(-1, INT_MAX);
39
BigDecimal BigDecimal::Nil;
40

    
41
static int64 powersOfTen[21];
42
static double negativePowersOfTen[21];
43

    
44
class PowersOfTenInitializer
45
{
46
    public:
47
        PowersOfTenInitializer();
48
};
49

    
50
PowersOfTenInitializer initializer;
51

    
52
PowersOfTenInitializer::PowersOfTenInitializer()
53
{
54
    int64 power = 1;
55
    for (unsigned int i = 0; i < sizeof(powersOfTen) / sizeof(*powersOfTen); i++) {
56
        powersOfTen[i] = power;
57
        power *= 10;
58
    }
59

    
60
    double negativePower = 1;
61
    for (unsigned int i = 0; i < sizeof(negativePowersOfTen) / sizeof(*negativePowersOfTen); i++) {
62
        negativePowersOfTen[i] = negativePower;
63
        negativePower /= 10.0;
64
    }
65
}
66

    
67
void BigDecimal::normalize()
68
{
69
    // special values
70
    if (scale == INT_MAX && (intVal == -1 || intVal == 0 || intVal == 1))
71
        return;
72

    
73
    // zero
74
    if (intVal == 0)
75
    {
76
        intVal = 0;
77
        scale = 0;
78
        return;
79
    }
80

    
81
    // underflow, XXX should throw an exception?
82
    if (scale < minScale - INT64_MAX_DIGITS)
83
    {
84
        intVal = 0;
85
        scale = 0;
86
        return;
87
    }
88

    
89
    // overflow
90
    if (scale > maxScale + INT64_MAX_DIGITS)
91
        throw opp_runtime_error("BigDecimal::normalize(): scale too big: %d.", scale); // XXX should be +-Infinity?
92

    
93
    // transform scale between minScale and maxScale
94
    if (scale < minScale)
95
    {
96
        while (scale < minScale)
97
        {
98
            intVal /= 10;
99
            scale++;
100

    
101
            if (intVal == 0)
102
            {
103
                scale = 0;
104
                return;
105
            }
106
        }
107
    }
108
    else if (scale > maxScale)
109
    {
110
        while (scale > maxScale)
111
        {
112
            if (intVal > INT64_MAX/10)
113
                throw opp_runtime_error("BigDecimal::normalize(): arithmetic overflow");
114

    
115
            intVal *= 10;
116
            scale--;
117
        }
118
    }
119

    
120
    // strip trailing zeros
121
    while ((intVal % 10 == 0) && scale < maxScale)
122
    {
123
        intVal /= 10;
124
        scale++;
125
    }
126
}
127

    
128
int64 BigDecimal::getDigits(int scale, int numDigits) const
129
{
130
    assert(!this->isSpecial());
131

    
132
    int start = max(scale, this->scale); // inclusive
133
    int end = scale+numDigits; // exclusive
134

    
135
    if (start >= end)
136
        return 0;
137

    
138
    int64 val = abs(this->intVal);
139
    for (int i = this->scale; i < start; ++i)
140
        val /= 10;
141

    
142
    if (val == 0)
143
        return 0;
144

    
145
    int64 result = 0;
146
    int digit;
147
    int64 multiplier = 1;
148
    for (int i = start; i < end; ++i)
149
    {
150
        digit = val % 10;
151
        val /= 10;
152
        result += multiplier*digit;
153
        multiplier *= 10;
154
    }
155

    
156
    for (int i = 0; i < (start-scale); ++i)
157
        result *= 10;
158

    
159
    return result;
160
}
161

    
162
const BigDecimal& BigDecimal::operator=(double d)
163
{
164
    // check NaN and infinity
165
    if (::isNaN(d))
166
        return *this = NaN;
167
    else if (::isPositiveInfinity(d))
168
        return *this = PositiveInfinity;
169
    else if (::isNegativeInfinity(d))
170
        return *this = NegativeInfinity;
171

    
172
    int sign = 1;
173
    if (d < 0.0)
174
    {
175
        sign = -1;
176
        d = -d;
177
    }
178

    
179
    int exponent;
180
    double mantissa = frexp(d, &exponent); // d = mantissa * 2 ^ exponent
181

    
182
    for (int i=0; i < 52; ++i)
183
    {
184
        mantissa *= 2.0;
185
        exponent--;
186
    }
187

    
188
    int64 intVal = (int64)mantissa;   // d = intVal * 2 ^ exponent
189

    
190

    
191
    // special case for 0.0, next loop would not terminate
192
    if (intVal == 0)
193
    {
194
        this->intVal = 0;
195
        this->scale = 0;
196
        return *this;
197
    }
198

    
199
    // normalize
200
    while ((intVal & 1) == 0)
201
    {
202
        intVal >>= 1;
203
        exponent++;
204
    }
205

    
206
    //printf("intVal=%I64d, exponent=%d\n", intVal, exponent);
207

    
208
    int scale;
209
    if (exponent < 0)
210
    {
211
        scale = exponent;
212
        for (int i = exponent; i < 0; ++i)
213
        {
214
            if (intVal <= INT64_MAX/5)
215
            {
216
                intVal *= 5;
217
            }
218
            else
219
            {
220
                intVal /= 2;
221
                scale++;
222
            }
223
        }
224
    }
225
    else
226
    {
227
        scale = 0;
228
        for (int i = 0; i < exponent; ++i)
229
        {
230
            if (intVal <= INT64_MAX/2)
231
            {
232
                intVal *= 2;
233
            }
234
            else
235
            {
236
                intVal /= 5;
237
                scale++;
238
            }
239
        }
240
    }
241

    
242
    this->intVal = sign * intVal;
243
    this->scale = scale;
244
    this->normalize();
245
    return *this;
246
}
247

    
248
bool BigDecimal::operator<(const BigDecimal &x) const
249
{
250
    if (isSpecial() || x.isSpecial())
251
    {
252
        if (isNil() || x.isNil())
253
            throw opp_runtime_error("BigDecimal::operator<() received Nil.");
254
        else if (isNaN() || x.isNaN())
255
            return false;
256
        else if (x == PositiveInfinity)
257
            return *this != PositiveInfinity;
258
        else if (*this == NegativeInfinity)
259
            return x != NegativeInfinity;
260
        else
261
            return false;
262
    }
263

    
264
    if (scale == x.scale)
265
        return intVal < x.intVal;
266
    if (sgn(intVal) < sgn(x.intVal))
267
        return true;
268
    if (sgn(intVal) > sgn(x.intVal))
269
        return false;
270

    
271
    assert((intVal < 0 && x.intVal < 0) || (intVal > 0 && x.intVal > 0));
272
    bool negatives = intVal < 0;
273

    
274
    // compare absolute values by comparing most significant digits first
275
    bool result = false;
276
    for (int s = max(scale,x.scale); s > min(scale,x.scale)-18; s-=18)
277
    {
278
        int64 digits = this->getDigits(s, 18);
279
        int64 digitsX = x.getDigits(s, 18);
280
        if (digits < digitsX)
281
        {
282
            result = true;
283
            break;
284
        }
285
        else if (digits > digitsX)
286
        {
287
            result = false;
288
            break;
289
        }
290
    }
291

    
292
    return negatives ? !result : result;
293
}
294

    
295
int64 BigDecimal::getMantissaForScale(int reqScale) const
296
{
297
    if (isSpecial())
298
        throw opp_runtime_error("BigDecimal: cannot return mantissa for Nil, NaN or +/-Inf value");
299
    checkScale(reqScale);
300
    int scaleDiff = scale - reqScale;
301
    if (scaleDiff == 0)
302
        return intVal;
303
    else if (scaleDiff < 0)
304
        return intVal / powersOfTen[-scaleDiff];
305
    else
306
        return intVal * powersOfTen[scaleDiff];
307
}
308

    
309
double BigDecimal::dbl() const
310
{
311
    if (isSpecial())
312
    {
313
        if (isNaN())
314
            return ::NaN;
315
        else if (*this == PositiveInfinity)
316
            return POSITIVE_INFINITY;
317
        else if (*this == NegativeInfinity)
318
            return NEGATIVE_INFINITY;
319
        else // Nil
320
            throw opp_runtime_error("BigDecimal::dbl(): received Nil."); // XXX should return NaN?
321
    }
322

    
323
    return (double)intVal * negativePowersOfTen[-scale];
324
}
325

    
326
std::string BigDecimal::str() const
327
{
328
    // delegate to operator<<
329
    std::stringstream out;
330
    out << *this;
331
    return out.str();
332
}
333

    
334
char *BigDecimal::ttoa(char *buf, const BigDecimal &x, char *&endp)
335
{
336
    // special values
337
    if (x.isSpecial())
338
    {
339
        if (x.isNaN())
340
        {
341
            strcpy(buf, "NaN");
342
            endp = buf+3;
343
        }
344
        else if (x == PositiveInfinity)
345
        {
346
            strcpy(buf, "+Inf");
347
            endp = buf+4;
348
        }
349
        else if (x == NegativeInfinity)
350
        {
351
            strcpy(buf, "-Inf");
352
            endp = buf+4;
353
        }
354
        else // Nil
355
            throw opp_runtime_error("BigDecimal::ttoa(): received Nil.");
356
        return buf;
357
    }
358

    
359
    int64 intVal = x.getIntValue();
360
    int scale = x.getScale();
361

    
362
    // prepare for conversion
363
    endp = buf+63;  //19+19+5 should be enough, but play it safe
364
    *endp = '\0';
365
    char *s = endp;
366
    if (intVal==0)
367
        {*--s = '0'; return s;}
368

    
369
    // convert digits
370
    bool negative = intVal<0;
371
    if (negative) intVal = -intVal;
372

    
373
    bool skipzeros = true;
374
    int decimalplace = scale;
375
    do {
376
        int64 res = intVal / 10;
377
        int digit = intVal - (10*res);
378

    
379
        if (skipzeros && (digit!=0 || decimalplace>=0))
380
            skipzeros = false;
381
        if (decimalplace++==0 && s!=endp)
382
            *--s = '.';
383
        if (!skipzeros)
384
            *--s = '0'+digit;
385
        intVal = res;
386
    } while (intVal);
387

    
388
    // add leading zeros, decimal point, etc if needed
389
    if (decimalplace<=0)
390
    {
391
        while (decimalplace++ < 0)
392
            *--s = '0';
393
        *--s = '.';
394
        *--s = '0';
395
    }
396

    
397
    if (negative) *--s = '-';
398
    return s;
399
}
400

    
401
const BigDecimal BigDecimal::parse(const char *s)
402
{
403
    const char *endp;
404
    return parse(s, endp);
405
}
406

    
407
#define OVERFLOW_CHECK(c,s) if (!(c)) throw opp_runtime_error("BigDecimal::parse(\"%s\"): arithmetic overflow", (s));
408

    
409

    
410
const BigDecimal BigDecimal::parse(const char *s, const char *&endp)
411
{
412
    int64 intVal = 0;
413
    int digit;
414
    int digits = 0;
415
    int scale = 0;
416
    int sign = 1;
417
    const char *p = s;
418

    
419
    // skip leading spaces
420
    while (opp_isspace(*p))
421
        ++p;
422

    
423
    // optional signs
424
    if (*p == '-')
425
    {
426
        sign = -1;
427
        ++p;
428
    }
429
    else if (*p == '+')
430
        ++p;
431

    
432
    // parse special numbers
433
    if (opp_isalpha(*p))
434
    {
435
        if (strncasecmp(p, "nan", 3) == 0)
436
        {
437
            endp = p+3;
438
            return NaN;
439
        }
440
        else if (strncasecmp(p, "inf", 3) == 0) // inf and infinity
441
        {
442
            endp = p+3;
443
            if (strncasecmp(endp, "inity", 5) == 0)
444
                endp += 5;
445
            return sign > 0 ? PositiveInfinity : NegativeInfinity;
446
        }
447
        else
448
        {
449
            endp = p;
450
            return Zero;  // XXX should return Nil?
451
        }
452
    }
453
    else if (*p=='1' && *(p+1)=='.' && *(p+2)=='#')
454
    {
455
        if (strncasecmp(p+3, "ind", 3) == 0)
456
        {
457
            endp = p+6;
458
            return NaN;
459
        }
460
        else if (strncasecmp(p+3, "inf", 6) == 0)
461
        {
462
            endp = p+6;
463
            return sign > 0 ? PositiveInfinity : NegativeInfinity;
464
        }
465
    }
466

    
467
    // digits before decimal
468
    while (opp_isdigit(*p))
469
    {
470
        OVERFLOW_CHECK(intVal <= INT64_MAX / 10, s);
471
        intVal *= 10;
472
        digit = ((*p++)-'0');
473
        OVERFLOW_CHECK(intVal <= INT64_MAX - digit, s);
474
        intVal += digit;
475
        digits++;
476
    }
477

    
478
    if (digits == 0 && *p != '.')
479
    {
480
        throw opp_runtime_error("BigDecimal::parse(\"%s\"): missing digits", s);
481
    }
482

    
483
    // digits after decimal
484
    if (*p == '.')
485
    {
486
        p++;
487
        while (opp_isdigit(*p))
488
        {
489
            OVERFLOW_CHECK(intVal <= INT64_MAX / 10, s);
490
            intVal *= 10;
491
            digit = ((*p++)-'0');
492
            OVERFLOW_CHECK(intVal <= INT64_MAX - digit, s);
493
            intVal += digit;
494
            digits++;
495
            scale--;
496
        }
497
    }
498

    
499
    endp = p;
500
    return BigDecimal(sign*intVal, scale);
501
}
502

    
503
const BigDecimal operator+(const BigDecimal& x, const BigDecimal& y)
504
{
505
    // 1. try to add exactly
506
    int scale = std::min(x.scale, y.scale);
507
    int xm = x.scale - scale;
508
    int ym = y.scale - scale;
509

    
510
    const int NUMPOWERS = sizeof(powersOfTen) / sizeof(*powersOfTen);
511

    
512
    if (!x.isSpecial() && !y.isSpecial() && 0 <= xm && xm < NUMPOWERS && 0 <= ym && ym < NUMPOWERS)
513
    {
514
        int64 xmp = powersOfTen[xm];
515
        int64 xv = x.intVal * xmp;
516

    
517
        if (xv / xmp == x.intVal) {
518
            int64 ymp = powersOfTen[ym];
519
            int64 yv = y.intVal * ymp;
520

    
521
            if (yv / ymp == y.intVal) {
522
                bool sameSign = haveSameSign(xv, yv);
523
                int64 intVal = xv + yv;
524

    
525
                if (!sameSign || haveSameSign(intVal, yv))
526
                    return BigDecimal(intVal, scale);
527
            }
528
        }
529
    }
530

    
531
    // 2. add with precision loss
532
    return BigDecimal(x.dbl()+y.dbl());
533
}
534

    
535
const BigDecimal operator-(const BigDecimal& x, const BigDecimal& y)
536
{
537
    // 1. try to subtract exactly
538
    int scale = std::min(x.scale, y.scale);
539
    int xm = x.scale - scale;
540
    int ym = y.scale - scale;
541

    
542
    const int NUMPOWERS = sizeof(powersOfTen) / sizeof(*powersOfTen);
543

    
544
    if (!x.isSpecial() && !y.isSpecial() && 0 <= xm && xm < NUMPOWERS && 0 <= ym && ym < NUMPOWERS)
545
    {
546
        int64 xmp = powersOfTen[xm];
547
        int64 xv = x.intVal * xmp;
548

    
549
        if (xv / xmp == x.intVal) {
550
            int64 ymp = powersOfTen[ym];
551
            int64 yv = y.intVal * ymp;
552

    
553
            if (yv / ymp == y.intVal) {
554
                bool differentSign = !haveSameSign(xv, yv);
555
                int64 intVal = xv - yv;
556

    
557
                if (!differentSign || !haveSameSign(intVal, yv))
558
                    return BigDecimal(intVal, scale);
559
            }
560
        }
561
    }
562

    
563
    // 2. subtract with precision loss
564
    return BigDecimal(x.dbl()-y.dbl());
565
}