root / src / common / bigdecimal.cc @ a3be1d55
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) 20062008 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 < (startscale); ++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 
} 