1 /*
2 *
3 * Ruby BigDecimal(Variable decimal precision) extension library.
4 *
5 * Copyright(C) 2002 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp)
6 *
7 */
8
9 /* #define BIGDECIMAL_DEBUG 1 */
10 #ifdef BIGDECIMAL_DEBUG
11 # define BIGDECIMAL_ENABLE_VPRINT 1
12 #endif
13 #include "bigdecimal.h"
14 #include "ruby/util.h"
15
16 #ifndef BIGDECIMAL_DEBUG
17 # define NDEBUG
18 #endif
19 #include <assert.h>
20
21 #include <ctype.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <errno.h>
26 #include <math.h>
27 #include "math.h"
28
29 #ifdef HAVE_IEEEFP_H
30 #include <ieeefp.h>
31 #endif
32
33 /* #define ENABLE_NUMERIC_STRING */
34
35 #define MUL_OVERFLOW_SIGNED_INTEGER_P(a, b, min, max) ( \
36 (a) == 0 ? 0 : \
37 (a) == -1 ? (b) < -(max) : \
38 (a) > 0 ? \
39 ((b) > 0 ? (max) / (a) < (b) : (min) / (a) > (b)) : \
40 ((b) > 0 ? (min) / (a) < (b) : (max) / (a) > (b)))
41 #define SIGNED_VALUE_MAX INTPTR_MAX
42 #define SIGNED_VALUE_MIN INTPTR_MIN
43 #define MUL_OVERFLOW_SIGNED_VALUE_P(a, b) MUL_OVERFLOW_SIGNED_INTEGER_P(a, b, SIGNED_VALUE_MIN, SIGNED_VALUE_MAX)
44
45 VALUE rb_cBigDecimal;
46 VALUE rb_mBigMath;
47
48 static ID id_BigDecimal_exception_mode;
49 static ID id_BigDecimal_rounding_mode;
50 static ID id_BigDecimal_precision_limit;
51
52 static ID id_up;
53 static ID id_down;
54 static ID id_truncate;
55 static ID id_half_up;
56 static ID id_default;
57 static ID id_half_down;
58 static ID id_half_even;
59 static ID id_banker;
60 static ID id_ceiling;
61 static ID id_ceil;
62 static ID id_floor;
63 static ID id_to_r;
64 static ID id_eq;
65 static ID id_half;
66
67 /* MACRO's to guard objects from GC by keeping them in stack */
68 #define ENTER(n) volatile VALUE RB_UNUSED_VAR(vStack[n]);int iStack=0
69 #define PUSH(x) (vStack[iStack++] = (VALUE)(x))
70 #define SAVE(p) PUSH((p)->obj)
71 #define GUARD_OBJ(p,y) ((p)=(y), SAVE(p))
72
73 #define BASE_FIG RMPD_COMPONENT_FIGURES
74 #define BASE RMPD_BASE
75
76 #define HALF_BASE (BASE/2)
77 #define BASE1 (BASE/10)
78
79 #ifndef DBLE_FIG
80 #define DBLE_FIG (DBL_DIG+1) /* figure of double */
81 #endif
82
83 #ifndef RRATIONAL_ZERO_P
84 # define RRATIONAL_ZERO_P(x) (FIXNUM_P(rb_rational_num(x)) && \
85 FIX2LONG(rb_rational_num(x)) == 0)
86 #endif
87
88 #ifndef RRATIONAL_NEGATIVE_P
89 # define RRATIONAL_NEGATIVE_P(x) RTEST(rb_funcall((x), '<', 1, INT2FIX(0)))
90 #endif
91
92 #ifndef DECIMAL_SIZE_OF_BITS
93 #define DECIMAL_SIZE_OF_BITS(n) (((n) * 3010 + 9998) / 9999)
94 /* an approximation of ceil(n * log10(2)), upto 65536 at least */
95 #endif
96
97 #ifdef PRIsVALUE
98 # define RB_OBJ_CLASSNAME(obj) rb_obj_class(obj)
99 # define RB_OBJ_STRING(obj) (obj)
100 #else
101 # define PRIsVALUE "s"
102 # define RB_OBJ_CLASSNAME(obj) rb_obj_classname(obj)
103 # define RB_OBJ_STRING(obj) StringValueCStr(obj)
104 #endif
105
106 #ifndef HAVE_RB_RATIONAL_NUM
107 static inline VALUE
rb_rational_num(VALUE rat)108 rb_rational_num(VALUE rat)
109 {
110 #ifdef HAVE_TYPE_STRUCT_RRATIONAL
111 return RRATIONAL(rat)->num;
112 #else
113 return rb_funcall(rat, rb_intern("numerator"), 0);
114 #endif
115 }
116 #endif
117
118 #ifndef HAVE_RB_RATIONAL_DEN
119 static inline VALUE
rb_rational_den(VALUE rat)120 rb_rational_den(VALUE rat)
121 {
122 #ifdef HAVE_TYPE_STRUCT_RRATIONAL
123 return RRATIONAL(rat)->den;
124 #else
125 return rb_funcall(rat, rb_intern("denominator"), 0);
126 #endif
127 }
128 #endif
129
130 #define BIGDECIMAL_POSITIVE_P(bd) ((bd)->sign > 0)
131 #define BIGDECIMAL_NEGATIVE_P(bd) ((bd)->sign < 0)
132
133 /*
134 * ================== Ruby Interface part ==========================
135 */
136 #define DoSomeOne(x,y,f) rb_num_coerce_bin(x,y,f)
137
138 /*
139 * VP routines used in BigDecimal part
140 */
141 static unsigned short VpGetException(void);
142 static void VpSetException(unsigned short f);
143 static void VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v);
144 static int VpLimitRound(Real *c, size_t ixDigit);
145 static Real *VpCopy(Real *pv, Real const* const x);
146
147 #ifdef BIGDECIMAL_ENABLE_VPRINT
148 static int VPrint(FILE *fp,const char *cntl_chr,Real *a);
149 #endif
150
151 /*
152 * **** BigDecimal part ****
153 */
154
155 static void
BigDecimal_delete(void * pv)156 BigDecimal_delete(void *pv)
157 {
158 VpFree(pv);
159 }
160
161 static size_t
BigDecimal_memsize(const void * ptr)162 BigDecimal_memsize(const void *ptr)
163 {
164 const Real *pv = ptr;
165 return (sizeof(*pv) + pv->MaxPrec * sizeof(BDIGIT));
166 }
167
168 static const rb_data_type_t BigDecimal_data_type = {
169 "BigDecimal",
170 { 0, BigDecimal_delete, BigDecimal_memsize, },
171 #ifdef RUBY_TYPED_FREE_IMMEDIATELY
172 0, 0, RUBY_TYPED_FREE_IMMEDIATELY
173 #endif
174 };
175
176 static inline int
is_kind_of_BigDecimal(VALUE const v)177 is_kind_of_BigDecimal(VALUE const v)
178 {
179 return rb_typeddata_is_kind_of(v, &BigDecimal_data_type);
180 }
181
182 static VALUE
ToValue(Real * p)183 ToValue(Real *p)
184 {
185 if (VpIsNaN(p)) {
186 VpException(VP_EXCEPTION_NaN, "Computation results to 'NaN'(Not a Number)", 0);
187 }
188 else if (VpIsPosInf(p)) {
189 VpException(VP_EXCEPTION_INFINITY, "Computation results to 'Infinity'", 0);
190 }
191 else if (VpIsNegInf(p)) {
192 VpException(VP_EXCEPTION_INFINITY, "Computation results to '-Infinity'", 0);
193 }
194 return p->obj;
195 }
196
197 NORETURN(static void cannot_be_coerced_into_BigDecimal(VALUE, VALUE));
198
199 static void
cannot_be_coerced_into_BigDecimal(VALUE exc_class,VALUE v)200 cannot_be_coerced_into_BigDecimal(VALUE exc_class, VALUE v)
201 {
202 VALUE str;
203
204 if (rb_special_const_p(v)) {
205 str = rb_inspect(v);
206 }
207 else {
208 str = rb_class_name(rb_obj_class(v));
209 }
210
211 str = rb_str_cat2(rb_str_dup(str), " can't be coerced into BigDecimal");
212 rb_exc_raise(rb_exc_new3(exc_class, str));
213 }
214
215 static inline VALUE BigDecimal_div2(VALUE, VALUE, VALUE);
216
217 static Real*
GetVpValueWithPrec(VALUE v,long prec,int must)218 GetVpValueWithPrec(VALUE v, long prec, int must)
219 {
220 ENTER(1);
221 Real *pv;
222 VALUE num, bg;
223 char szD[128];
224 VALUE orig = Qundef;
225 double d;
226
227 again:
228 switch(TYPE(v)) {
229 case T_FLOAT:
230 if (prec < 0) goto unable_to_coerce_without_prec;
231 if (prec > DBL_DIG+1) goto SomeOneMayDoIt;
232 d = RFLOAT_VALUE(v);
233 if (!isfinite(d)) {
234 pv = VpCreateRbObject(1, NULL);
235 VpDtoV(pv, d);
236 return pv;
237 }
238 if (d != 0.0) {
239 v = rb_funcall(v, id_to_r, 0);
240 goto again;
241 }
242 if (1/d < 0.0) {
243 return VpCreateRbObject(prec, "-0");
244 }
245 return VpCreateRbObject(prec, "0");
246
247 case T_RATIONAL:
248 if (prec < 0) goto unable_to_coerce_without_prec;
249
250 if (orig == Qundef ? (orig = v, 1) : orig != v) {
251 num = rb_rational_num(v);
252 pv = GetVpValueWithPrec(num, -1, must);
253 if (pv == NULL) goto SomeOneMayDoIt;
254
255 v = BigDecimal_div2(ToValue(pv), rb_rational_den(v), LONG2NUM(prec));
256 goto again;
257 }
258
259 v = orig;
260 goto SomeOneMayDoIt;
261
262 case T_DATA:
263 if (is_kind_of_BigDecimal(v)) {
264 pv = DATA_PTR(v);
265 return pv;
266 }
267 else {
268 goto SomeOneMayDoIt;
269 }
270 break;
271
272 case T_FIXNUM:
273 sprintf(szD, "%ld", FIX2LONG(v));
274 return VpCreateRbObject(VpBaseFig() * 2 + 1, szD);
275
276 #ifdef ENABLE_NUMERIC_STRING
277 case T_STRING:
278 StringValueCStr(v);
279 rb_check_safe_obj(v);
280 return VpCreateRbObject(RSTRING_LEN(v) + VpBaseFig() + 1,
281 RSTRING_PTR(v));
282 #endif /* ENABLE_NUMERIC_STRING */
283
284 case T_BIGNUM:
285 bg = rb_big2str(v, 10);
286 PUSH(bg);
287 return VpCreateRbObject(strlen(RSTRING_PTR(bg)) + VpBaseFig() + 1,
288 RSTRING_PTR(bg));
289 default:
290 goto SomeOneMayDoIt;
291 }
292
293 SomeOneMayDoIt:
294 if (must) {
295 cannot_be_coerced_into_BigDecimal(rb_eTypeError, v);
296 }
297 return NULL; /* NULL means to coerce */
298
299 unable_to_coerce_without_prec:
300 if (must) {
301 rb_raise(rb_eArgError,
302 "%"PRIsVALUE" can't be coerced into BigDecimal without a precision",
303 RB_OBJ_CLASSNAME(v));
304 }
305 return NULL;
306 }
307
308 static Real*
GetVpValue(VALUE v,int must)309 GetVpValue(VALUE v, int must)
310 {
311 return GetVpValueWithPrec(v, -1, must);
312 }
313
314 /* call-seq:
315 * BigDecimal.double_fig
316 *
317 * The BigDecimal.double_fig class method returns the number of digits a
318 * Float number is allowed to have. The result depends upon the CPU and OS
319 * in use.
320 */
321 static VALUE
BigDecimal_double_fig(VALUE self)322 BigDecimal_double_fig(VALUE self)
323 {
324 return INT2FIX(VpDblFig());
325 }
326
327 /* call-seq:
328 * big_decimal.precs -> array
329 *
330 * Returns an Array of two Integer values.
331 *
332 * The first value is the current number of significant digits in the
333 * BigDecimal. The second value is the maximum number of significant digits
334 * for the BigDecimal.
335 *
336 * BigDecimal('5').precs #=> [9, 18]
337 */
338
339 static VALUE
BigDecimal_prec(VALUE self)340 BigDecimal_prec(VALUE self)
341 {
342 ENTER(1);
343 Real *p;
344 VALUE obj;
345
346 GUARD_OBJ(p, GetVpValue(self, 1));
347 obj = rb_assoc_new(INT2NUM(p->Prec*VpBaseFig()),
348 INT2NUM(p->MaxPrec*VpBaseFig()));
349 return obj;
350 }
351
352 /*
353 * call-seq: hash
354 *
355 * Creates a hash for this BigDecimal.
356 *
357 * Two BigDecimals with equal sign,
358 * fractional part and exponent have the same hash.
359 */
360 static VALUE
BigDecimal_hash(VALUE self)361 BigDecimal_hash(VALUE self)
362 {
363 ENTER(1);
364 Real *p;
365 st_index_t hash;
366
367 GUARD_OBJ(p, GetVpValue(self, 1));
368 hash = (st_index_t)p->sign;
369 /* hash!=2: the case for 0(1),NaN(0) or +-Infinity(3) is sign itself */
370 if(hash == 2 || hash == (st_index_t)-2) {
371 hash ^= rb_memhash(p->frac, sizeof(BDIGIT)*p->Prec);
372 hash += p->exponent;
373 }
374 return ST2FIX(hash);
375 }
376
377 /*
378 * call-seq: _dump
379 *
380 * Method used to provide marshalling support.
381 *
382 * inf = BigDecimal('Infinity')
383 * #=> Infinity
384 * BigDecimal._load(inf._dump)
385 * #=> Infinity
386 *
387 * See the Marshal module.
388 */
389 static VALUE
BigDecimal_dump(int argc,VALUE * argv,VALUE self)390 BigDecimal_dump(int argc, VALUE *argv, VALUE self)
391 {
392 ENTER(5);
393 Real *vp;
394 char *psz;
395 VALUE dummy;
396 volatile VALUE dump;
397
398 rb_scan_args(argc, argv, "01", &dummy);
399 GUARD_OBJ(vp,GetVpValue(self, 1));
400 dump = rb_str_new(0, VpNumOfChars(vp, "E")+50);
401 psz = RSTRING_PTR(dump);
402 sprintf(psz, "%"PRIuSIZE":", VpMaxPrec(vp)*VpBaseFig());
403 VpToString(vp, psz+strlen(psz), 0, 0);
404 rb_str_resize(dump, strlen(psz));
405 return dump;
406 }
407
408 /*
409 * Internal method used to provide marshalling support. See the Marshal module.
410 */
411 static VALUE
BigDecimal_load(VALUE self,VALUE str)412 BigDecimal_load(VALUE self, VALUE str)
413 {
414 ENTER(2);
415 Real *pv;
416 unsigned char *pch;
417 unsigned char ch;
418 unsigned long m=0;
419
420 pch = (unsigned char *)StringValueCStr(str);
421 rb_check_safe_obj(str);
422 /* First get max prec */
423 while((*pch) != (unsigned char)'\0' && (ch = *pch++) != (unsigned char)':') {
424 if(!ISDIGIT(ch)) {
425 rb_raise(rb_eTypeError, "load failed: invalid character in the marshaled string");
426 }
427 m = m*10 + (unsigned long)(ch-'0');
428 }
429 if (m > VpBaseFig()) m -= VpBaseFig();
430 GUARD_OBJ(pv, VpNewRbClass(m, (char *)pch, self));
431 m /= VpBaseFig();
432 if (m && pv->MaxPrec > m) {
433 pv->MaxPrec = m+1;
434 }
435 return ToValue(pv);
436 }
437
438 static unsigned short
check_rounding_mode_option(VALUE const opts)439 check_rounding_mode_option(VALUE const opts)
440 {
441 VALUE mode;
442 char const *s;
443 long l;
444
445 assert(RB_TYPE_P(opts, T_HASH));
446
447 if (NIL_P(opts))
448 goto noopt;
449
450 mode = rb_hash_lookup2(opts, ID2SYM(id_half), Qundef);
451 if (mode == Qundef || NIL_P(mode))
452 goto noopt;
453
454 if (SYMBOL_P(mode))
455 mode = rb_sym2str(mode);
456 else if (!RB_TYPE_P(mode, T_STRING)) {
457 VALUE str_mode = rb_check_string_type(mode);
458 if (NIL_P(str_mode)) goto invalid;
459 mode = str_mode;
460 }
461 s = RSTRING_PTR(mode);
462 l = RSTRING_LEN(mode);
463 switch (l) {
464 case 2:
465 if (strncasecmp(s, "up", 2) == 0)
466 return VP_ROUND_HALF_UP;
467 break;
468 case 4:
469 if (strncasecmp(s, "even", 4) == 0)
470 return VP_ROUND_HALF_EVEN;
471 else if (strncasecmp(s, "down", 4) == 0)
472 return VP_ROUND_HALF_DOWN;
473 break;
474 default:
475 break;
476 }
477 invalid:
478 if (NIL_P(mode))
479 rb_raise(rb_eArgError, "invalid rounding mode: nil");
480 else
481 rb_raise(rb_eArgError, "invalid rounding mode: %"PRIsVALUE, mode);
482
483 noopt:
484 return VpGetRoundMode();
485 }
486
487 static unsigned short
check_rounding_mode(VALUE const v)488 check_rounding_mode(VALUE const v)
489 {
490 unsigned short sw;
491 ID id;
492 switch (TYPE(v)) {
493 case T_SYMBOL:
494 id = SYM2ID(v);
495 if (id == id_up)
496 return VP_ROUND_UP;
497 if (id == id_down || id == id_truncate)
498 return VP_ROUND_DOWN;
499 if (id == id_half_up || id == id_default)
500 return VP_ROUND_HALF_UP;
501 if (id == id_half_down)
502 return VP_ROUND_HALF_DOWN;
503 if (id == id_half_even || id == id_banker)
504 return VP_ROUND_HALF_EVEN;
505 if (id == id_ceiling || id == id_ceil)
506 return VP_ROUND_CEIL;
507 if (id == id_floor)
508 return VP_ROUND_FLOOR;
509 rb_raise(rb_eArgError, "invalid rounding mode");
510
511 default:
512 break;
513 }
514
515 sw = NUM2USHORT(v);
516 if (!VpIsRoundMode(sw)) {
517 rb_raise(rb_eArgError, "invalid rounding mode");
518 }
519 return sw;
520 }
521
522 /* call-seq:
523 * BigDecimal.mode(mode, value)
524 *
525 * Controls handling of arithmetic exceptions and rounding. If no value
526 * is supplied, the current value is returned.
527 *
528 * Six values of the mode parameter control the handling of arithmetic
529 * exceptions:
530 *
531 * BigDecimal::EXCEPTION_NaN
532 * BigDecimal::EXCEPTION_INFINITY
533 * BigDecimal::EXCEPTION_UNDERFLOW
534 * BigDecimal::EXCEPTION_OVERFLOW
535 * BigDecimal::EXCEPTION_ZERODIVIDE
536 * BigDecimal::EXCEPTION_ALL
537 *
538 * For each mode parameter above, if the value set is false, computation
539 * continues after an arithmetic exception of the appropriate type.
540 * When computation continues, results are as follows:
541 *
542 * EXCEPTION_NaN:: NaN
543 * EXCEPTION_INFINITY:: +Infinity or -Infinity
544 * EXCEPTION_UNDERFLOW:: 0
545 * EXCEPTION_OVERFLOW:: +Infinity or -Infinity
546 * EXCEPTION_ZERODIVIDE:: +Infinity or -Infinity
547 *
548 * One value of the mode parameter controls the rounding of numeric values:
549 * BigDecimal::ROUND_MODE. The values it can take are:
550 *
551 * ROUND_UP, :up:: round away from zero
552 * ROUND_DOWN, :down, :truncate:: round towards zero (truncate)
553 * ROUND_HALF_UP, :half_up, :default:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round away from zero. (default)
554 * ROUND_HALF_DOWN, :half_down:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards zero.
555 * ROUND_HALF_EVEN, :half_even, :banker:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards the even neighbor (Banker's rounding)
556 * ROUND_CEILING, :ceiling, :ceil:: round towards positive infinity (ceil)
557 * ROUND_FLOOR, :floor:: round towards negative infinity (floor)
558 *
559 */
560 static VALUE
BigDecimal_mode(int argc,VALUE * argv,VALUE self)561 BigDecimal_mode(int argc, VALUE *argv, VALUE self)
562 {
563 VALUE which;
564 VALUE val;
565 unsigned long f,fo;
566
567 rb_scan_args(argc, argv, "11", &which, &val);
568 f = (unsigned long)NUM2INT(which);
569
570 if (f & VP_EXCEPTION_ALL) {
571 /* Exception mode setting */
572 fo = VpGetException();
573 if (val == Qnil) return INT2FIX(fo);
574 if (val != Qfalse && val!=Qtrue) {
575 rb_raise(rb_eArgError, "second argument must be true or false");
576 return Qnil; /* Not reached */
577 }
578 if (f & VP_EXCEPTION_INFINITY) {
579 VpSetException((unsigned short)((val == Qtrue) ? (fo | VP_EXCEPTION_INFINITY) :
580 (fo & (~VP_EXCEPTION_INFINITY))));
581 }
582 fo = VpGetException();
583 if (f & VP_EXCEPTION_NaN) {
584 VpSetException((unsigned short)((val == Qtrue) ? (fo | VP_EXCEPTION_NaN) :
585 (fo & (~VP_EXCEPTION_NaN))));
586 }
587 fo = VpGetException();
588 if (f & VP_EXCEPTION_UNDERFLOW) {
589 VpSetException((unsigned short)((val == Qtrue) ? (fo | VP_EXCEPTION_UNDERFLOW) :
590 (fo & (~VP_EXCEPTION_UNDERFLOW))));
591 }
592 fo = VpGetException();
593 if(f & VP_EXCEPTION_ZERODIVIDE) {
594 VpSetException((unsigned short)((val == Qtrue) ? (fo | VP_EXCEPTION_ZERODIVIDE) :
595 (fo & (~VP_EXCEPTION_ZERODIVIDE))));
596 }
597 fo = VpGetException();
598 return INT2FIX(fo);
599 }
600 if (VP_ROUND_MODE == f) {
601 /* Rounding mode setting */
602 unsigned short sw;
603 fo = VpGetRoundMode();
604 if (NIL_P(val)) return INT2FIX(fo);
605 sw = check_rounding_mode(val);
606 fo = VpSetRoundMode(sw);
607 return INT2FIX(fo);
608 }
609 rb_raise(rb_eTypeError, "first argument for BigDecimal.mode invalid");
610 return Qnil;
611 }
612
613 static size_t
GetAddSubPrec(Real * a,Real * b)614 GetAddSubPrec(Real *a, Real *b)
615 {
616 size_t mxs;
617 size_t mx = a->Prec;
618 SIGNED_VALUE d;
619
620 if (!VpIsDef(a) || !VpIsDef(b)) return (size_t)-1L;
621 if (mx < b->Prec) mx = b->Prec;
622 if (a->exponent != b->exponent) {
623 mxs = mx;
624 d = a->exponent - b->exponent;
625 if (d < 0) d = -d;
626 mx = mx + (size_t)d;
627 if (mx < mxs) {
628 return VpException(VP_EXCEPTION_INFINITY, "Exponent overflow", 0);
629 }
630 }
631 return mx;
632 }
633
634 static SIGNED_VALUE
GetPrecisionInt(VALUE v)635 GetPrecisionInt(VALUE v)
636 {
637 SIGNED_VALUE n;
638 n = NUM2INT(v);
639 if (n < 0) {
640 rb_raise(rb_eArgError, "negative precision");
641 }
642 return n;
643 }
644
645 VP_EXPORT Real *
VpNewRbClass(size_t mx,const char * str,VALUE klass)646 VpNewRbClass(size_t mx, const char *str, VALUE klass)
647 {
648 VALUE obj = TypedData_Wrap_Struct(klass, &BigDecimal_data_type, 0);
649 Real *pv = VpAlloc(mx, str, 1, 1);
650 RTYPEDDATA_DATA(obj) = pv;
651 pv->obj = obj;
652 RB_OBJ_FREEZE(obj);
653 return pv;
654 }
655
656 VP_EXPORT Real *
VpCreateRbObject(size_t mx,const char * str)657 VpCreateRbObject(size_t mx, const char *str)
658 {
659 return VpNewRbClass(mx, str, rb_cBigDecimal);
660 }
661
662 #define VpAllocReal(prec) (Real *)VpMemAlloc(offsetof(Real, frac) + (prec) * sizeof(BDIGIT))
663 #define VpReallocReal(ptr, prec) (Real *)VpMemRealloc((ptr), offsetof(Real, frac) + (prec) * sizeof(BDIGIT))
664
665 static Real *
VpCopy(Real * pv,Real const * const x)666 VpCopy(Real *pv, Real const* const x)
667 {
668 assert(x != NULL);
669
670 pv = VpReallocReal(pv, x->MaxPrec);
671 pv->MaxPrec = x->MaxPrec;
672 pv->Prec = x->Prec;
673 pv->exponent = x->exponent;
674 pv->sign = x->sign;
675 pv->flag = x->flag;
676 MEMCPY(pv->frac, x->frac, BDIGIT, pv->MaxPrec);
677
678 return pv;
679 }
680
681 /* Returns True if the value is Not a Number. */
682 static VALUE
BigDecimal_IsNaN(VALUE self)683 BigDecimal_IsNaN(VALUE self)
684 {
685 Real *p = GetVpValue(self, 1);
686 if (VpIsNaN(p)) return Qtrue;
687 return Qfalse;
688 }
689
690 /* Returns nil, -1, or +1 depending on whether the value is finite,
691 * -Infinity, or +Infinity.
692 */
693 static VALUE
BigDecimal_IsInfinite(VALUE self)694 BigDecimal_IsInfinite(VALUE self)
695 {
696 Real *p = GetVpValue(self, 1);
697 if (VpIsPosInf(p)) return INT2FIX(1);
698 if (VpIsNegInf(p)) return INT2FIX(-1);
699 return Qnil;
700 }
701
702 /* Returns True if the value is finite (not NaN or infinite). */
703 static VALUE
BigDecimal_IsFinite(VALUE self)704 BigDecimal_IsFinite(VALUE self)
705 {
706 Real *p = GetVpValue(self, 1);
707 if (VpIsNaN(p)) return Qfalse;
708 if (VpIsInf(p)) return Qfalse;
709 return Qtrue;
710 }
711
712 static void
BigDecimal_check_num(Real * p)713 BigDecimal_check_num(Real *p)
714 {
715 if (VpIsNaN(p)) {
716 VpException(VP_EXCEPTION_NaN, "Computation results to 'NaN'(Not a Number)", 1);
717 }
718 else if (VpIsPosInf(p)) {
719 VpException(VP_EXCEPTION_INFINITY, "Computation results to 'Infinity'", 1);
720 }
721 else if (VpIsNegInf(p)) {
722 VpException(VP_EXCEPTION_INFINITY, "Computation results to '-Infinity'", 1);
723 }
724 }
725
726 static VALUE BigDecimal_split(VALUE self);
727
728 /* Returns the value as an Integer.
729 *
730 * If the BigDecimal is infinity or NaN, raises FloatDomainError.
731 */
732 static VALUE
BigDecimal_to_i(VALUE self)733 BigDecimal_to_i(VALUE self)
734 {
735 ENTER(5);
736 ssize_t e, nf;
737 Real *p;
738
739 GUARD_OBJ(p, GetVpValue(self, 1));
740 BigDecimal_check_num(p);
741
742 e = VpExponent10(p);
743 if (e <= 0) return INT2FIX(0);
744 nf = VpBaseFig();
745 if (e <= nf) {
746 return LONG2NUM((long)(VpGetSign(p) * (BDIGIT_DBL_SIGNED)p->frac[0]));
747 }
748 else {
749 VALUE a = BigDecimal_split(self);
750 VALUE digits = RARRAY_AREF(a, 1);
751 VALUE numerator = rb_funcall(digits, rb_intern("to_i"), 0);
752 VALUE ret;
753 ssize_t dpower = e - (ssize_t)RSTRING_LEN(digits);
754
755 if (BIGDECIMAL_NEGATIVE_P(p)) {
756 numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1));
757 }
758 if (dpower < 0) {
759 ret = rb_funcall(numerator, rb_intern("div"), 1,
760 rb_funcall(INT2FIX(10), rb_intern("**"), 1,
761 INT2FIX(-dpower)));
762 }
763 else {
764 ret = rb_funcall(numerator, '*', 1,
765 rb_funcall(INT2FIX(10), rb_intern("**"), 1,
766 INT2FIX(dpower)));
767 }
768 if (RB_TYPE_P(ret, T_FLOAT)) {
769 rb_raise(rb_eFloatDomainError, "Infinity");
770 }
771 return ret;
772 }
773 }
774
775 /* Returns a new Float object having approximately the same value as the
776 * BigDecimal number. Normal accuracy limits and built-in errors of binary
777 * Float arithmetic apply.
778 */
779 static VALUE
BigDecimal_to_f(VALUE self)780 BigDecimal_to_f(VALUE self)
781 {
782 ENTER(1);
783 Real *p;
784 double d;
785 SIGNED_VALUE e;
786 char *buf;
787 volatile VALUE str;
788
789 GUARD_OBJ(p, GetVpValue(self, 1));
790 if (VpVtoD(&d, &e, p) != 1)
791 return rb_float_new(d);
792 if (e > (SIGNED_VALUE)(DBL_MAX_10_EXP+BASE_FIG))
793 goto overflow;
794 if (e < (SIGNED_VALUE)(DBL_MIN_10_EXP-BASE_FIG))
795 goto underflow;
796
797 str = rb_str_new(0, VpNumOfChars(p, "E"));
798 buf = RSTRING_PTR(str);
799 VpToString(p, buf, 0, 0);
800 errno = 0;
801 d = strtod(buf, 0);
802 if (errno == ERANGE) {
803 if (d == 0.0) goto underflow;
804 if (fabs(d) >= HUGE_VAL) goto overflow;
805 }
806 return rb_float_new(d);
807
808 overflow:
809 VpException(VP_EXCEPTION_OVERFLOW, "BigDecimal to Float conversion", 0);
810 if (BIGDECIMAL_NEGATIVE_P(p))
811 return rb_float_new(VpGetDoubleNegInf());
812 else
813 return rb_float_new(VpGetDoublePosInf());
814
815 underflow:
816 VpException(VP_EXCEPTION_UNDERFLOW, "BigDecimal to Float conversion", 0);
817 if (BIGDECIMAL_NEGATIVE_P(p))
818 return rb_float_new(-0.0);
819 else
820 return rb_float_new(0.0);
821 }
822
823
824 /* Converts a BigDecimal to a Rational.
825 */
826 static VALUE
BigDecimal_to_r(VALUE self)827 BigDecimal_to_r(VALUE self)
828 {
829 Real *p;
830 ssize_t sign, power, denomi_power;
831 VALUE a, digits, numerator;
832
833 p = GetVpValue(self, 1);
834 BigDecimal_check_num(p);
835
836 sign = VpGetSign(p);
837 power = VpExponent10(p);
838 a = BigDecimal_split(self);
839 digits = RARRAY_AREF(a, 1);
840 denomi_power = power - RSTRING_LEN(digits);
841 numerator = rb_funcall(digits, rb_intern("to_i"), 0);
842
843 if (sign < 0) {
844 numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1));
845 }
846 if (denomi_power < 0) {
847 return rb_Rational(numerator,
848 rb_funcall(INT2FIX(10), rb_intern("**"), 1,
849 INT2FIX(-denomi_power)));
850 }
851 else {
852 return rb_Rational1(rb_funcall(numerator, '*', 1,
853 rb_funcall(INT2FIX(10), rb_intern("**"), 1,
854 INT2FIX(denomi_power))));
855 }
856 }
857
858 /* The coerce method provides support for Ruby type coercion. It is not
859 * enabled by default.
860 *
861 * This means that binary operations like + * / or - can often be performed
862 * on a BigDecimal and an object of another type, if the other object can
863 * be coerced into a BigDecimal value.
864 *
865 * e.g.
866 * a = BigDecimal("1.0")
867 * b = a / 2.0 #=> 0.5
868 *
869 * Note that coercing a String to a BigDecimal is not supported by default;
870 * it requires a special compile-time option when building Ruby.
871 */
872 static VALUE
BigDecimal_coerce(VALUE self,VALUE other)873 BigDecimal_coerce(VALUE self, VALUE other)
874 {
875 ENTER(2);
876 VALUE obj;
877 Real *b;
878
879 if (RB_TYPE_P(other, T_FLOAT)) {
880 GUARD_OBJ(b, GetVpValueWithPrec(other, DBL_DIG+1, 1));
881 obj = rb_assoc_new(ToValue(b), self);
882 }
883 else {
884 if (RB_TYPE_P(other, T_RATIONAL)) {
885 Real* pv = DATA_PTR(self);
886 GUARD_OBJ(b, GetVpValueWithPrec(other, pv->Prec*VpBaseFig(), 1));
887 }
888 else {
889 GUARD_OBJ(b, GetVpValue(other, 1));
890 }
891 obj = rb_assoc_new(b->obj, self);
892 }
893
894 return obj;
895 }
896
897 /*
898 * call-seq:
899 * +big_decimal -> big_decimal
900 *
901 * Return self.
902 *
903 * +BigDecimal('5') #=> 0.5e1
904 */
905
906 static VALUE
BigDecimal_uplus(VALUE self)907 BigDecimal_uplus(VALUE self)
908 {
909 return self;
910 }
911
912 /*
913 * Document-method: BigDecimal#add
914 * Document-method: BigDecimal#+
915 *
916 * call-seq:
917 * add(value, digits)
918 *
919 * Add the specified value.
920 *
921 * e.g.
922 * c = a.add(b,n)
923 * c = a + b
924 *
925 * digits:: If specified and less than the number of significant digits of the
926 * result, the result is rounded to that number of digits, according
927 * to BigDecimal.mode.
928 */
929 static VALUE
BigDecimal_add(VALUE self,VALUE r)930 BigDecimal_add(VALUE self, VALUE r)
931 {
932 ENTER(5);
933 Real *c, *a, *b;
934 size_t mx;
935
936 GUARD_OBJ(a, GetVpValue(self, 1));
937 if (RB_TYPE_P(r, T_FLOAT)) {
938 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
939 }
940 else if (RB_TYPE_P(r, T_RATIONAL)) {
941 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
942 }
943 else {
944 b = GetVpValue(r, 0);
945 }
946
947 if (!b) return DoSomeOne(self,r,'+');
948 SAVE(b);
949
950 if (VpIsNaN(b)) return b->obj;
951 if (VpIsNaN(a)) return a->obj;
952
953 mx = GetAddSubPrec(a, b);
954 if (mx == (size_t)-1L) {
955 GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
956 VpAddSub(c, a, b, 1);
957 }
958 else {
959 GUARD_OBJ(c, VpCreateRbObject(mx * (VpBaseFig() + 1), "0"));
960 if(!mx) {
961 VpSetInf(c, VpGetSign(a));
962 }
963 else {
964 VpAddSub(c, a, b, 1);
965 }
966 }
967 return ToValue(c);
968 }
969
970 /* call-seq:
971 * a - b -> bigdecimal
972 *
973 * Subtract the specified value.
974 *
975 * e.g.
976 * c = a - b
977 *
978 * The precision of the result value depends on the type of +b+.
979 *
980 * If +b+ is a Float, the precision of the result is Float::DIG+1.
981 *
982 * If +b+ is a BigDecimal, the precision of the result is +b+'s precision of
983 * internal representation from platform. So, it's return value is platform
984 * dependent.
985 *
986 */
987 static VALUE
BigDecimal_sub(VALUE self,VALUE r)988 BigDecimal_sub(VALUE self, VALUE r)
989 {
990 ENTER(5);
991 Real *c, *a, *b;
992 size_t mx;
993
994 GUARD_OBJ(a, GetVpValue(self,1));
995 if (RB_TYPE_P(r, T_FLOAT)) {
996 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
997 }
998 else if (RB_TYPE_P(r, T_RATIONAL)) {
999 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
1000 }
1001 else {
1002 b = GetVpValue(r,0);
1003 }
1004
1005 if (!b) return DoSomeOne(self,r,'-');
1006 SAVE(b);
1007
1008 if (VpIsNaN(b)) return b->obj;
1009 if (VpIsNaN(a)) return a->obj;
1010
1011 mx = GetAddSubPrec(a,b);
1012 if (mx == (size_t)-1L) {
1013 GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
1014 VpAddSub(c, a, b, -1);
1015 }
1016 else {
1017 GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
1018 if (!mx) {
1019 VpSetInf(c,VpGetSign(a));
1020 }
1021 else {
1022 VpAddSub(c, a, b, -1);
1023 }
1024 }
1025 return ToValue(c);
1026 }
1027
1028 static VALUE
BigDecimalCmp(VALUE self,VALUE r,char op)1029 BigDecimalCmp(VALUE self, VALUE r,char op)
1030 {
1031 ENTER(5);
1032 SIGNED_VALUE e;
1033 Real *a, *b=0;
1034 GUARD_OBJ(a, GetVpValue(self, 1));
1035 switch (TYPE(r)) {
1036 case T_DATA:
1037 if (!is_kind_of_BigDecimal(r)) break;
1038 /* fall through */
1039 case T_FIXNUM:
1040 /* fall through */
1041 case T_BIGNUM:
1042 GUARD_OBJ(b, GetVpValue(r, 0));
1043 break;
1044
1045 case T_FLOAT:
1046 GUARD_OBJ(b, GetVpValueWithPrec(r, DBL_DIG+1, 0));
1047 break;
1048
1049 case T_RATIONAL:
1050 GUARD_OBJ(b, GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 0));
1051 break;
1052
1053 default:
1054 break;
1055 }
1056 if (b == NULL) {
1057 ID f = 0;
1058
1059 switch (op) {
1060 case '*':
1061 return rb_num_coerce_cmp(self, r, rb_intern("<=>"));
1062
1063 case '=':
1064 return RTEST(rb_num_coerce_cmp(self, r, rb_intern("=="))) ? Qtrue : Qfalse;
1065
1066 case 'G':
1067 f = rb_intern(">=");
1068 break;
1069
1070 case 'L':
1071 f = rb_intern("<=");
1072 break;
1073
1074 case '>':
1075 /* fall through */
1076 case '<':
1077 f = (ID)op;
1078 break;
1079
1080 default:
1081 break;
1082 }
1083 return rb_num_coerce_relop(self, r, f);
1084 }
1085 SAVE(b);
1086 e = VpComp(a, b);
1087 if (e == 999)
1088 return (op == '*') ? Qnil : Qfalse;
1089 switch (op) {
1090 case '*':
1091 return INT2FIX(e); /* any op */
1092
1093 case '=':
1094 if (e == 0) return Qtrue;
1095 return Qfalse;
1096
1097 case 'G':
1098 if (e >= 0) return Qtrue;
1099 return Qfalse;
1100
1101 case '>':
1102 if (e > 0) return Qtrue;
1103 return Qfalse;
1104
1105 case 'L':
1106 if (e <= 0) return Qtrue;
1107 return Qfalse;
1108
1109 case '<':
1110 if (e < 0) return Qtrue;
1111 return Qfalse;
1112
1113 default:
1114 break;
1115 }
1116
1117 rb_bug("Undefined operation in BigDecimalCmp()");
1118
1119 UNREACHABLE;
1120 }
1121
1122 /* Returns True if the value is zero. */
1123 static VALUE
BigDecimal_zero(VALUE self)1124 BigDecimal_zero(VALUE self)
1125 {
1126 Real *a = GetVpValue(self, 1);
1127 return VpIsZero(a) ? Qtrue : Qfalse;
1128 }
1129
1130 /* Returns self if the value is non-zero, nil otherwise. */
1131 static VALUE
BigDecimal_nonzero(VALUE self)1132 BigDecimal_nonzero(VALUE self)
1133 {
1134 Real *a = GetVpValue(self, 1);
1135 return VpIsZero(a) ? Qnil : self;
1136 }
1137
1138 /* The comparison operator.
1139 * a <=> b is 0 if a == b, 1 if a > b, -1 if a < b.
1140 */
1141 static VALUE
BigDecimal_comp(VALUE self,VALUE r)1142 BigDecimal_comp(VALUE self, VALUE r)
1143 {
1144 return BigDecimalCmp(self, r, '*');
1145 }
1146
1147 /*
1148 * Tests for value equality; returns true if the values are equal.
1149 *
1150 * The == and === operators and the eql? method have the same implementation
1151 * for BigDecimal.
1152 *
1153 * Values may be coerced to perform the comparison:
1154 *
1155 * BigDecimal('1.0') == 1.0 #=> true
1156 */
1157 static VALUE
BigDecimal_eq(VALUE self,VALUE r)1158 BigDecimal_eq(VALUE self, VALUE r)
1159 {
1160 return BigDecimalCmp(self, r, '=');
1161 }
1162
1163 /* call-seq:
1164 * a < b
1165 *
1166 * Returns true if a is less than b.
1167 *
1168 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
1169 */
1170 static VALUE
BigDecimal_lt(VALUE self,VALUE r)1171 BigDecimal_lt(VALUE self, VALUE r)
1172 {
1173 return BigDecimalCmp(self, r, '<');
1174 }
1175
1176 /* call-seq:
1177 * a <= b
1178 *
1179 * Returns true if a is less than or equal to b.
1180 *
1181 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
1182 */
1183 static VALUE
BigDecimal_le(VALUE self,VALUE r)1184 BigDecimal_le(VALUE self, VALUE r)
1185 {
1186 return BigDecimalCmp(self, r, 'L');
1187 }
1188
1189 /* call-seq:
1190 * a > b
1191 *
1192 * Returns true if a is greater than b.
1193 *
1194 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
1195 */
1196 static VALUE
BigDecimal_gt(VALUE self,VALUE r)1197 BigDecimal_gt(VALUE self, VALUE r)
1198 {
1199 return BigDecimalCmp(self, r, '>');
1200 }
1201
1202 /* call-seq:
1203 * a >= b
1204 *
1205 * Returns true if a is greater than or equal to b.
1206 *
1207 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce)
1208 */
1209 static VALUE
BigDecimal_ge(VALUE self,VALUE r)1210 BigDecimal_ge(VALUE self, VALUE r)
1211 {
1212 return BigDecimalCmp(self, r, 'G');
1213 }
1214
1215 /*
1216 * call-seq:
1217 * -big_decimal -> big_decimal
1218 *
1219 * Return the negation of self.
1220 *
1221 * -BigDecimal('5') #=> -0.5e1
1222 */
1223
1224 static VALUE
BigDecimal_neg(VALUE self)1225 BigDecimal_neg(VALUE self)
1226 {
1227 ENTER(5);
1228 Real *c, *a;
1229 GUARD_OBJ(a, GetVpValue(self, 1));
1230 GUARD_OBJ(c, VpCreateRbObject(a->Prec *(VpBaseFig() + 1), "0"));
1231 VpAsgn(c, a, -1);
1232 return ToValue(c);
1233 }
1234
1235 /*
1236 * Document-method: BigDecimal#mult
1237 *
1238 * call-seq: mult(value, digits)
1239 *
1240 * Multiply by the specified value.
1241 *
1242 * e.g.
1243 * c = a.mult(b,n)
1244 * c = a * b
1245 *
1246 * digits:: If specified and less than the number of significant digits of the
1247 * result, the result is rounded to that number of digits, according
1248 * to BigDecimal.mode.
1249 */
1250 static VALUE
BigDecimal_mult(VALUE self,VALUE r)1251 BigDecimal_mult(VALUE self, VALUE r)
1252 {
1253 ENTER(5);
1254 Real *c, *a, *b;
1255 size_t mx;
1256
1257 GUARD_OBJ(a, GetVpValue(self, 1));
1258 if (RB_TYPE_P(r, T_FLOAT)) {
1259 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
1260 }
1261 else if (RB_TYPE_P(r, T_RATIONAL)) {
1262 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
1263 }
1264 else {
1265 b = GetVpValue(r,0);
1266 }
1267
1268 if (!b) return DoSomeOne(self, r, '*');
1269 SAVE(b);
1270
1271 mx = a->Prec + b->Prec;
1272 GUARD_OBJ(c, VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
1273 VpMult(c, a, b);
1274 return ToValue(c);
1275 }
1276
1277 static VALUE
BigDecimal_divide(Real ** c,Real ** res,Real ** div,VALUE self,VALUE r)1278 BigDecimal_divide(Real **c, Real **res, Real **div, VALUE self, VALUE r)
1279 /* For c = self.div(r): with round operation */
1280 {
1281 ENTER(5);
1282 Real *a, *b;
1283 size_t mx;
1284
1285 GUARD_OBJ(a, GetVpValue(self, 1));
1286 if (RB_TYPE_P(r, T_FLOAT)) {
1287 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
1288 }
1289 else if (RB_TYPE_P(r, T_RATIONAL)) {
1290 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
1291 }
1292 else {
1293 b = GetVpValue(r, 0);
1294 }
1295
1296 if (!b) return DoSomeOne(self, r, '/');
1297 SAVE(b);
1298
1299 *div = b;
1300 mx = a->Prec + vabs(a->exponent);
1301 if (mx < b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent);
1302 mx++; /* NOTE: An additional digit is needed for the compatibility to
1303 the version 1.2.1 and the former. */
1304 mx = (mx + 1) * VpBaseFig();
1305 GUARD_OBJ((*c), VpCreateRbObject(mx, "#0"));
1306 GUARD_OBJ((*res), VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
1307 VpDivd(*c, *res, a, b);
1308 return Qnil;
1309 }
1310
1311 /* call-seq:
1312 * a / b -> bigdecimal
1313 * quo(value) -> bigdecimal
1314 *
1315 * Divide by the specified value.
1316 *
1317 * See BigDecimal#div.
1318 */
1319 static VALUE
BigDecimal_div(VALUE self,VALUE r)1320 BigDecimal_div(VALUE self, VALUE r)
1321 /* For c = self/r: with round operation */
1322 {
1323 ENTER(5);
1324 Real *c=NULL, *res=NULL, *div = NULL;
1325 r = BigDecimal_divide(&c, &res, &div, self, r);
1326 if (!NIL_P(r)) return r; /* coerced by other */
1327 SAVE(c); SAVE(res); SAVE(div);
1328 /* a/b = c + r/b */
1329 /* c xxxxx
1330 r 00000yyyyy ==> (y/b)*BASE >= HALF_BASE
1331 */
1332 /* Round */
1333 if (VpHasVal(div)) { /* frac[0] must be zero for NaN,INF,Zero */
1334 VpInternalRound(c, 0, c->frac[c->Prec-1], (BDIGIT)(VpBaseVal() * (BDIGIT_DBL)res->frac[0] / div->frac[0]));
1335 }
1336 return ToValue(c);
1337 }
1338
1339 /*
1340 * %: mod = a%b = a - (a.to_f/b).floor * b
1341 * div = (a.to_f/b).floor
1342 */
1343 static VALUE
BigDecimal_DoDivmod(VALUE self,VALUE r,Real ** div,Real ** mod)1344 BigDecimal_DoDivmod(VALUE self, VALUE r, Real **div, Real **mod)
1345 {
1346 ENTER(8);
1347 Real *c=NULL, *d=NULL, *res=NULL;
1348 Real *a, *b;
1349 size_t mx;
1350
1351 GUARD_OBJ(a, GetVpValue(self, 1));
1352 if (RB_TYPE_P(r, T_FLOAT)) {
1353 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
1354 }
1355 else if (RB_TYPE_P(r, T_RATIONAL)) {
1356 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
1357 }
1358 else {
1359 b = GetVpValue(r, 0);
1360 }
1361
1362 if (!b) return Qfalse;
1363 SAVE(b);
1364
1365 if (VpIsNaN(a) || VpIsNaN(b)) goto NaN;
1366 if (VpIsInf(a) && VpIsInf(b)) goto NaN;
1367 if (VpIsZero(b)) {
1368 rb_raise(rb_eZeroDivError, "divided by 0");
1369 }
1370 if (VpIsInf(a)) {
1371 GUARD_OBJ(d, VpCreateRbObject(1, "0"));
1372 VpSetInf(d, (SIGNED_VALUE)(VpGetSign(a) == VpGetSign(b) ? 1 : -1));
1373 GUARD_OBJ(c, VpCreateRbObject(1, "NaN"));
1374 *div = d;
1375 *mod = c;
1376 return Qtrue;
1377 }
1378 if (VpIsInf(b)) {
1379 GUARD_OBJ(d, VpCreateRbObject(1, "0"));
1380 *div = d;
1381 *mod = a;
1382 return Qtrue;
1383 }
1384 if (VpIsZero(a)) {
1385 GUARD_OBJ(c, VpCreateRbObject(1, "0"));
1386 GUARD_OBJ(d, VpCreateRbObject(1, "0"));
1387 *div = d;
1388 *mod = c;
1389 return Qtrue;
1390 }
1391
1392 mx = a->Prec + vabs(a->exponent);
1393 if (mx<b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent);
1394 mx = (mx + 1) * VpBaseFig();
1395 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1396 GUARD_OBJ(res, VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
1397 VpDivd(c, res, a, b);
1398 mx = c->Prec * (VpBaseFig() + 1);
1399 GUARD_OBJ(d, VpCreateRbObject(mx, "0"));
1400 VpActiveRound(d, c, VP_ROUND_DOWN, 0);
1401 VpMult(res, d, b);
1402 VpAddSub(c, a, res, -1);
1403 if (!VpIsZero(c) && (VpGetSign(a) * VpGetSign(b) < 0)) {
1404 VpAddSub(res, d, VpOne(), -1);
1405 GUARD_OBJ(d, VpCreateRbObject(GetAddSubPrec(c, b)*(VpBaseFig() + 1), "0"));
1406 VpAddSub(d, c, b, 1);
1407 *div = res;
1408 *mod = d;
1409 } else {
1410 *div = d;
1411 *mod = c;
1412 }
1413 return Qtrue;
1414
1415 NaN:
1416 GUARD_OBJ(c, VpCreateRbObject(1, "NaN"));
1417 GUARD_OBJ(d, VpCreateRbObject(1, "NaN"));
1418 *div = d;
1419 *mod = c;
1420 return Qtrue;
1421 }
1422
1423 /* call-seq:
1424 * a % b
1425 * a.modulo(b)
1426 *
1427 * Returns the modulus from dividing by b.
1428 *
1429 * See BigDecimal#divmod.
1430 */
1431 static VALUE
BigDecimal_mod(VALUE self,VALUE r)1432 BigDecimal_mod(VALUE self, VALUE r) /* %: a%b = a - (a.to_f/b).floor * b */
1433 {
1434 ENTER(3);
1435 Real *div = NULL, *mod = NULL;
1436
1437 if (BigDecimal_DoDivmod(self, r, &div, &mod)) {
1438 SAVE(div); SAVE(mod);
1439 return ToValue(mod);
1440 }
1441 return DoSomeOne(self, r, '%');
1442 }
1443
1444 static VALUE
BigDecimal_divremain(VALUE self,VALUE r,Real ** dv,Real ** rv)1445 BigDecimal_divremain(VALUE self, VALUE r, Real **dv, Real **rv)
1446 {
1447 ENTER(10);
1448 size_t mx;
1449 Real *a = NULL, *b = NULL, *c = NULL, *res = NULL, *d = NULL, *rr = NULL, *ff = NULL;
1450 Real *f = NULL;
1451
1452 GUARD_OBJ(a, GetVpValue(self, 1));
1453 if (RB_TYPE_P(r, T_FLOAT)) {
1454 b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
1455 }
1456 else if (RB_TYPE_P(r, T_RATIONAL)) {
1457 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
1458 }
1459 else {
1460 b = GetVpValue(r, 0);
1461 }
1462
1463 if (!b) return DoSomeOne(self, r, rb_intern("remainder"));
1464 SAVE(b);
1465
1466 mx = (a->MaxPrec + b->MaxPrec) *VpBaseFig();
1467 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1468 GUARD_OBJ(res, VpCreateRbObject((mx+1) * 2 + (VpBaseFig() + 1), "#0"));
1469 GUARD_OBJ(rr, VpCreateRbObject((mx+1) * 2 + (VpBaseFig() + 1), "#0"));
1470 GUARD_OBJ(ff, VpCreateRbObject((mx+1) * 2 + (VpBaseFig() + 1), "#0"));
1471
1472 VpDivd(c, res, a, b);
1473
1474 mx = c->Prec *(VpBaseFig() + 1);
1475
1476 GUARD_OBJ(d, VpCreateRbObject(mx, "0"));
1477 GUARD_OBJ(f, VpCreateRbObject(mx, "0"));
1478
1479 VpActiveRound(d, c, VP_ROUND_DOWN, 0); /* 0: round off */
1480
1481 VpFrac(f, c);
1482 VpMult(rr, f, b);
1483 VpAddSub(ff, res, rr, 1);
1484
1485 *dv = d;
1486 *rv = ff;
1487 return Qnil;
1488 }
1489
1490 /* call-seq:
1491 * remainder(value)
1492 *
1493 * Returns the remainder from dividing by the value.
1494 *
1495 * x.remainder(y) means x-y*(x/y).truncate
1496 */
1497 static VALUE
BigDecimal_remainder(VALUE self,VALUE r)1498 BigDecimal_remainder(VALUE self, VALUE r) /* remainder */
1499 {
1500 VALUE f;
1501 Real *d, *rv = 0;
1502 f = BigDecimal_divremain(self, r, &d, &rv);
1503 if (!NIL_P(f)) return f;
1504 return ToValue(rv);
1505 }
1506
1507 /* call-seq:
1508 * divmod(value)
1509 *
1510 * Divides by the specified value, and returns the quotient and modulus
1511 * as BigDecimal numbers. The quotient is rounded towards negative infinity.
1512 *
1513 * For example:
1514 *
1515 * require 'bigdecimal'
1516 *
1517 * a = BigDecimal("42")
1518 * b = BigDecimal("9")
1519 *
1520 * q, m = a.divmod(b)
1521 *
1522 * c = q * b + m
1523 *
1524 * a == c #=> true
1525 *
1526 * The quotient q is (a/b).floor, and the modulus is the amount that must be
1527 * added to q * b to get a.
1528 */
1529 static VALUE
BigDecimal_divmod(VALUE self,VALUE r)1530 BigDecimal_divmod(VALUE self, VALUE r)
1531 {
1532 ENTER(5);
1533 Real *div = NULL, *mod = NULL;
1534
1535 if (BigDecimal_DoDivmod(self, r, &div, &mod)) {
1536 SAVE(div); SAVE(mod);
1537 return rb_assoc_new(ToValue(div), ToValue(mod));
1538 }
1539 return DoSomeOne(self,r,rb_intern("divmod"));
1540 }
1541
1542 /*
1543 * See BigDecimal#quo
1544 */
1545 static inline VALUE
BigDecimal_div2(VALUE self,VALUE b,VALUE n)1546 BigDecimal_div2(VALUE self, VALUE b, VALUE n)
1547 {
1548 ENTER(5);
1549 SIGNED_VALUE ix;
1550
1551 if (NIL_P(n)) { /* div in Float sense */
1552 Real *div = NULL;
1553 Real *mod;
1554 if (BigDecimal_DoDivmod(self, b, &div, &mod)) {
1555 return BigDecimal_to_i(ToValue(div));
1556 }
1557 return DoSomeOne(self, b, rb_intern("div"));
1558 }
1559
1560 /* div in BigDecimal sense */
1561 ix = GetPrecisionInt(n);
1562 if (ix == 0) {
1563 return BigDecimal_div(self, b);
1564 }
1565 else {
1566 Real *res = NULL;
1567 Real *av = NULL, *bv = NULL, *cv = NULL;
1568 size_t mx = ix + VpBaseFig()*2;
1569 size_t pl = VpSetPrecLimit(0);
1570
1571 GUARD_OBJ(cv, VpCreateRbObject(mx + VpBaseFig(), "0"));
1572 GUARD_OBJ(av, GetVpValue(self, 1));
1573 GUARD_OBJ(bv, GetVpValue(b, 1));
1574 mx = av->Prec + bv->Prec + 2;
1575 if (mx <= cv->MaxPrec) mx = cv->MaxPrec + 1;
1576 GUARD_OBJ(res, VpCreateRbObject((mx * 2 + 2)*VpBaseFig(), "#0"));
1577 VpDivd(cv, res, av, bv);
1578 VpSetPrecLimit(pl);
1579 VpLeftRound(cv, VpGetRoundMode(), ix);
1580 return ToValue(cv);
1581 }
1582 }
1583
1584 /*
1585 * Document-method: BigDecimal#div
1586 *
1587 * call-seq:
1588 * div(value, digits) -> bigdecimal or integer
1589 *
1590 * Divide by the specified value.
1591 *
1592 * digits:: If specified and less than the number of significant digits of the
1593 * result, the result is rounded to that number of digits, according
1594 * to BigDecimal.mode.
1595 *
1596 * If digits is 0, the result is the same as for the / operator
1597 * or #quo.
1598 *
1599 * If digits is not specified, the result is an integer,
1600 * by analogy with Float#div; see also BigDecimal#divmod.
1601 *
1602 * Examples:
1603 *
1604 * a = BigDecimal("4")
1605 * b = BigDecimal("3")
1606 *
1607 * a.div(b, 3) # => 0.133e1
1608 *
1609 * a.div(b, 0) # => 0.1333333333333333333e1
1610 * a / b # => 0.1333333333333333333e1
1611 * a.quo(b) # => 0.1333333333333333333e1
1612 *
1613 * a.div(b) # => 1
1614 */
1615 static VALUE
BigDecimal_div3(int argc,VALUE * argv,VALUE self)1616 BigDecimal_div3(int argc, VALUE *argv, VALUE self)
1617 {
1618 VALUE b,n;
1619
1620 rb_scan_args(argc, argv, "11", &b, &n);
1621
1622 return BigDecimal_div2(self, b, n);
1623 }
1624
1625 static VALUE
BigDecimal_add2(VALUE self,VALUE b,VALUE n)1626 BigDecimal_add2(VALUE self, VALUE b, VALUE n)
1627 {
1628 ENTER(2);
1629 Real *cv;
1630 SIGNED_VALUE mx = GetPrecisionInt(n);
1631 if (mx == 0) return BigDecimal_add(self, b);
1632 else {
1633 size_t pl = VpSetPrecLimit(0);
1634 VALUE c = BigDecimal_add(self, b);
1635 VpSetPrecLimit(pl);
1636 GUARD_OBJ(cv, GetVpValue(c, 1));
1637 VpLeftRound(cv, VpGetRoundMode(), mx);
1638 return ToValue(cv);
1639 }
1640 }
1641
1642 /* call-seq:
1643 * sub(value, digits) -> bigdecimal
1644 *
1645 * Subtract the specified value.
1646 *
1647 * e.g.
1648 * c = a.sub(b,n)
1649 *
1650 * digits:: If specified and less than the number of significant digits of the
1651 * result, the result is rounded to that number of digits, according
1652 * to BigDecimal.mode.
1653 *
1654 */
1655 static VALUE
BigDecimal_sub2(VALUE self,VALUE b,VALUE n)1656 BigDecimal_sub2(VALUE self, VALUE b, VALUE n)
1657 {
1658 ENTER(2);
1659 Real *cv;
1660 SIGNED_VALUE mx = GetPrecisionInt(n);
1661 if (mx == 0) return BigDecimal_sub(self, b);
1662 else {
1663 size_t pl = VpSetPrecLimit(0);
1664 VALUE c = BigDecimal_sub(self, b);
1665 VpSetPrecLimit(pl);
1666 GUARD_OBJ(cv, GetVpValue(c, 1));
1667 VpLeftRound(cv, VpGetRoundMode(), mx);
1668 return ToValue(cv);
1669 }
1670 }
1671
1672
1673 static VALUE
BigDecimal_mult2(VALUE self,VALUE b,VALUE n)1674 BigDecimal_mult2(VALUE self, VALUE b, VALUE n)
1675 {
1676 ENTER(2);
1677 Real *cv;
1678 SIGNED_VALUE mx = GetPrecisionInt(n);
1679 if (mx == 0) return BigDecimal_mult(self, b);
1680 else {
1681 size_t pl = VpSetPrecLimit(0);
1682 VALUE c = BigDecimal_mult(self, b);
1683 VpSetPrecLimit(pl);
1684 GUARD_OBJ(cv, GetVpValue(c, 1));
1685 VpLeftRound(cv, VpGetRoundMode(), mx);
1686 return ToValue(cv);
1687 }
1688 }
1689
1690 /*
1691 * call-seq:
1692 * big_decimal.abs -> big_decimal
1693 *
1694 * Returns the absolute value, as a BigDecimal.
1695 *
1696 * BigDecimal('5').abs #=> 0.5e1
1697 * BigDecimal('-3').abs #=> 0.3e1
1698 */
1699
1700 static VALUE
BigDecimal_abs(VALUE self)1701 BigDecimal_abs(VALUE self)
1702 {
1703 ENTER(5);
1704 Real *c, *a;
1705 size_t mx;
1706
1707 GUARD_OBJ(a, GetVpValue(self, 1));
1708 mx = a->Prec *(VpBaseFig() + 1);
1709 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1710 VpAsgn(c, a, 1);
1711 VpChangeSign(c, 1);
1712 return ToValue(c);
1713 }
1714
1715 /* call-seq:
1716 * sqrt(n)
1717 *
1718 * Returns the square root of the value.
1719 *
1720 * Result has at least n significant digits.
1721 */
1722 static VALUE
BigDecimal_sqrt(VALUE self,VALUE nFig)1723 BigDecimal_sqrt(VALUE self, VALUE nFig)
1724 {
1725 ENTER(5);
1726 Real *c, *a;
1727 size_t mx, n;
1728
1729 GUARD_OBJ(a, GetVpValue(self, 1));
1730 mx = a->Prec * (VpBaseFig() + 1);
1731
1732 n = GetPrecisionInt(nFig) + VpDblFig() + BASE_FIG;
1733 if (mx <= n) mx = n;
1734 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1735 VpSqrt(c, a);
1736 return ToValue(c);
1737 }
1738
1739 /* Return the integer part of the number, as a BigDecimal.
1740 */
1741 static VALUE
BigDecimal_fix(VALUE self)1742 BigDecimal_fix(VALUE self)
1743 {
1744 ENTER(5);
1745 Real *c, *a;
1746 size_t mx;
1747
1748 GUARD_OBJ(a, GetVpValue(self, 1));
1749 mx = a->Prec *(VpBaseFig() + 1);
1750 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1751 VpActiveRound(c, a, VP_ROUND_DOWN, 0); /* 0: round off */
1752 return ToValue(c);
1753 }
1754
1755 /* call-seq:
1756 * round(n, mode)
1757 *
1758 * Round to the nearest integer (by default), returning the result as a
1759 * BigDecimal.
1760 *
1761 * BigDecimal('3.14159').round #=> 3
1762 * BigDecimal('8.7').round #=> 9
1763 * BigDecimal('-9.9').round #=> -10
1764 *
1765 * If n is specified and positive, the fractional part of the result has no
1766 * more than that many digits.
1767 *
1768 * If n is specified and negative, at least that many digits to the left of the
1769 * decimal point will be 0 in the result.
1770 *
1771 * BigDecimal('3.14159').round(3) #=> 3.142
1772 * BigDecimal('13345.234').round(-2) #=> 13300.0
1773 *
1774 * The value of the optional mode argument can be used to determine how
1775 * rounding is performed; see BigDecimal.mode.
1776 */
1777 static VALUE
BigDecimal_round(int argc,VALUE * argv,VALUE self)1778 BigDecimal_round(int argc, VALUE *argv, VALUE self)
1779 {
1780 ENTER(5);
1781 Real *c, *a;
1782 int iLoc = 0;
1783 VALUE vLoc;
1784 VALUE vRound;
1785 size_t mx, pl;
1786
1787 unsigned short sw = VpGetRoundMode();
1788
1789 switch (rb_scan_args(argc, argv, "02", &vLoc, &vRound)) {
1790 case 0:
1791 iLoc = 0;
1792 break;
1793 case 1:
1794 if (RB_TYPE_P(vLoc, T_HASH)) {
1795 sw = check_rounding_mode_option(vLoc);
1796 }
1797 else {
1798 iLoc = NUM2INT(vLoc);
1799 }
1800 break;
1801 case 2:
1802 iLoc = NUM2INT(vLoc);
1803 if (RB_TYPE_P(vRound, T_HASH)) {
1804 sw = check_rounding_mode_option(vRound);
1805 }
1806 else {
1807 sw = check_rounding_mode(vRound);
1808 }
1809 break;
1810 default:
1811 break;
1812 }
1813
1814 pl = VpSetPrecLimit(0);
1815 GUARD_OBJ(a, GetVpValue(self, 1));
1816 mx = a->Prec * (VpBaseFig() + 1);
1817 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1818 VpSetPrecLimit(pl);
1819 VpActiveRound(c, a, sw, iLoc);
1820 if (argc == 0) {
1821 return BigDecimal_to_i(ToValue(c));
1822 }
1823 return ToValue(c);
1824 }
1825
1826 /* call-seq:
1827 * truncate(n)
1828 *
1829 * Truncate to the nearest integer (by default), returning the result as a
1830 * BigDecimal.
1831 *
1832 * BigDecimal('3.14159').truncate #=> 3
1833 * BigDecimal('8.7').truncate #=> 8
1834 * BigDecimal('-9.9').truncate #=> -9
1835 *
1836 * If n is specified and positive, the fractional part of the result has no
1837 * more than that many digits.
1838 *
1839 * If n is specified and negative, at least that many digits to the left of the
1840 * decimal point will be 0 in the result.
1841 *
1842 * BigDecimal('3.14159').truncate(3) #=> 3.141
1843 * BigDecimal('13345.234').truncate(-2) #=> 13300.0
1844 */
1845 static VALUE
BigDecimal_truncate(int argc,VALUE * argv,VALUE self)1846 BigDecimal_truncate(int argc, VALUE *argv, VALUE self)
1847 {
1848 ENTER(5);
1849 Real *c, *a;
1850 int iLoc;
1851 VALUE vLoc;
1852 size_t mx, pl = VpSetPrecLimit(0);
1853
1854 if (rb_scan_args(argc, argv, "01", &vLoc) == 0) {
1855 iLoc = 0;
1856 }
1857 else {
1858 iLoc = NUM2INT(vLoc);
1859 }
1860
1861 GUARD_OBJ(a, GetVpValue(self, 1));
1862 mx = a->Prec * (VpBaseFig() + 1);
1863 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1864 VpSetPrecLimit(pl);
1865 VpActiveRound(c, a, VP_ROUND_DOWN, iLoc); /* 0: truncate */
1866 if (argc == 0) {
1867 return BigDecimal_to_i(ToValue(c));
1868 }
1869 return ToValue(c);
1870 }
1871
1872 /* Return the fractional part of the number, as a BigDecimal.
1873 */
1874 static VALUE
BigDecimal_frac(VALUE self)1875 BigDecimal_frac(VALUE self)
1876 {
1877 ENTER(5);
1878 Real *c, *a;
1879 size_t mx;
1880
1881 GUARD_OBJ(a, GetVpValue(self, 1));
1882 mx = a->Prec * (VpBaseFig() + 1);
1883 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1884 VpFrac(c, a);
1885 return ToValue(c);
1886 }
1887
1888 /* call-seq:
1889 * floor(n)
1890 *
1891 * Return the largest integer less than or equal to the value, as a BigDecimal.
1892 *
1893 * BigDecimal('3.14159').floor #=> 3
1894 * BigDecimal('-9.1').floor #=> -10
1895 *
1896 * If n is specified and positive, the fractional part of the result has no
1897 * more than that many digits.
1898 *
1899 * If n is specified and negative, at least that
1900 * many digits to the left of the decimal point will be 0 in the result.
1901 *
1902 * BigDecimal('3.14159').floor(3) #=> 3.141
1903 * BigDecimal('13345.234').floor(-2) #=> 13300.0
1904 */
1905 static VALUE
BigDecimal_floor(int argc,VALUE * argv,VALUE self)1906 BigDecimal_floor(int argc, VALUE *argv, VALUE self)
1907 {
1908 ENTER(5);
1909 Real *c, *a;
1910 int iLoc;
1911 VALUE vLoc;
1912 size_t mx, pl = VpSetPrecLimit(0);
1913
1914 if (rb_scan_args(argc, argv, "01", &vLoc)==0) {
1915 iLoc = 0;
1916 }
1917 else {
1918 iLoc = NUM2INT(vLoc);
1919 }
1920
1921 GUARD_OBJ(a, GetVpValue(self, 1));
1922 mx = a->Prec * (VpBaseFig() + 1);
1923 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1924 VpSetPrecLimit(pl);
1925 VpActiveRound(c, a, VP_ROUND_FLOOR, iLoc);
1926 #ifdef BIGDECIMAL_DEBUG
1927 VPrint(stderr, "floor: c=%\n", c);
1928 #endif
1929 if (argc == 0) {
1930 return BigDecimal_to_i(ToValue(c));
1931 }
1932 return ToValue(c);
1933 }
1934
1935 /* call-seq:
1936 * ceil(n)
1937 *
1938 * Return the smallest integer greater than or equal to the value, as a BigDecimal.
1939 *
1940 * BigDecimal('3.14159').ceil #=> 4
1941 * BigDecimal('-9.1').ceil #=> -9
1942 *
1943 * If n is specified and positive, the fractional part of the result has no
1944 * more than that many digits.
1945 *
1946 * If n is specified and negative, at least that
1947 * many digits to the left of the decimal point will be 0 in the result.
1948 *
1949 * BigDecimal('3.14159').ceil(3) #=> 3.142
1950 * BigDecimal('13345.234').ceil(-2) #=> 13400.0
1951 */
1952 static VALUE
BigDecimal_ceil(int argc,VALUE * argv,VALUE self)1953 BigDecimal_ceil(int argc, VALUE *argv, VALUE self)
1954 {
1955 ENTER(5);
1956 Real *c, *a;
1957 int iLoc;
1958 VALUE vLoc;
1959 size_t mx, pl = VpSetPrecLimit(0);
1960
1961 if (rb_scan_args(argc, argv, "01", &vLoc) == 0) {
1962 iLoc = 0;
1963 } else {
1964 iLoc = NUM2INT(vLoc);
1965 }
1966
1967 GUARD_OBJ(a, GetVpValue(self, 1));
1968 mx = a->Prec * (VpBaseFig() + 1);
1969 GUARD_OBJ(c, VpCreateRbObject(mx, "0"));
1970 VpSetPrecLimit(pl);
1971 VpActiveRound(c, a, VP_ROUND_CEIL, iLoc);
1972 if (argc == 0) {
1973 return BigDecimal_to_i(ToValue(c));
1974 }
1975 return ToValue(c);
1976 }
1977
1978 /* call-seq:
1979 * to_s(s)
1980 *
1981 * Converts the value to a string.
1982 *
1983 * The default format looks like 0.xxxxEnn.
1984 *
1985 * The optional parameter s consists of either an integer; or an optional '+'
1986 * or ' ', followed by an optional number, followed by an optional 'E' or 'F'.
1987 *
1988 * If there is a '+' at the start of s, positive values are returned with
1989 * a leading '+'.
1990 *
1991 * A space at the start of s returns positive values with a leading space.
1992 *
1993 * If s contains a number, a space is inserted after each group of that many
1994 * fractional digits.
1995 *
1996 * If s ends with an 'E', engineering notation (0.xxxxEnn) is used.
1997 *
1998 * If s ends with an 'F', conventional floating point notation is used.
1999 *
2000 * Examples:
2001 *
2002 * BigDecimal('-123.45678901234567890').to_s('5F')
2003 * #=> '-123.45678 90123 45678 9'
2004 *
2005 * BigDecimal('123.45678901234567890').to_s('+8F')
2006 * #=> '+123.45678901 23456789'
2007 *
2008 * BigDecimal('123.45678901234567890').to_s(' F')
2009 * #=> ' 123.4567890123456789'
2010 */
2011 static VALUE
BigDecimal_to_s(int argc,VALUE * argv,VALUE self)2012 BigDecimal_to_s(int argc, VALUE *argv, VALUE self)
2013 {
2014 ENTER(5);
2015 int fmt = 0; /* 0: E format, 1: F format */
2016 int fPlus = 0; /* 0: default, 1: set ' ' before digits, 2: set '+' before digits. */
2017 Real *vp;
2018 volatile VALUE str;
2019 char *psz;
2020 char ch;
2021 size_t nc, mc = 0;
2022 SIGNED_VALUE m;
2023 VALUE f;
2024
2025 GUARD_OBJ(vp, GetVpValue(self, 1));
2026
2027 if (rb_scan_args(argc, argv, "01", &f) == 1) {
2028 if (RB_TYPE_P(f, T_STRING)) {
2029 psz = StringValueCStr(f);
2030 rb_check_safe_obj(f);
2031 if (*psz == ' ') {
2032 fPlus = 1;
2033 psz++;
2034 }
2035 else if (*psz == '+') {
2036 fPlus = 2;
2037 psz++;
2038 }
2039 while ((ch = *psz++) != 0) {
2040 if (ISSPACE(ch)) {
2041 continue;
2042 }
2043 if (!ISDIGIT(ch)) {
2044 if (ch == 'F' || ch == 'f') {
2045 fmt = 1; /* F format */
2046 }
2047 break;
2048 }
2049 mc = mc*10 + ch - '0';
2050 }
2051 }
2052 else {
2053 m = NUM2INT(f);
2054 if (m <= 0) {
2055 rb_raise(rb_eArgError, "argument must be positive");
2056 }
2057 mc = (size_t)m;
2058 }
2059 }
2060 if (fmt) {
2061 nc = VpNumOfChars(vp, "F");
2062 }
2063 else {
2064 nc = VpNumOfChars(vp, "E");
2065 }
2066 if (mc > 0) {
2067 nc += (nc + mc - 1) / mc + 1;
2068 }
2069
2070 str = rb_str_new(0, nc);
2071 psz = RSTRING_PTR(str);
2072
2073 if (fmt) {
2074 VpToFString(vp, psz, mc, fPlus);
2075 }
2076 else {
2077 VpToString (vp, psz, mc, fPlus);
2078 }
2079 rb_str_resize(str, strlen(psz));
2080 return str;
2081 }
2082
2083 /* Splits a BigDecimal number into four parts, returned as an array of values.
2084 *
2085 * The first value represents the sign of the BigDecimal, and is -1 or 1, or 0
2086 * if the BigDecimal is Not a Number.
2087 *
2088 * The second value is a string representing the significant digits of the
2089 * BigDecimal, with no leading zeros.
2090 *
2091 * The third value is the base used for arithmetic (currently always 10) as an
2092 * Integer.
2093 *
2094 * The fourth value is an Integer exponent.
2095 *
2096 * If the BigDecimal can be represented as 0.xxxxxx*10**n, then xxxxxx is the
2097 * string of significant digits with no leading zeros, and n is the exponent.
2098 *
2099 * From these values, you can translate a BigDecimal to a float as follows:
2100 *
2101 * sign, significant_digits, base, exponent = a.split
2102 * f = sign * "0.#{significant_digits}".to_f * (base ** exponent)
2103 *
2104 * (Note that the to_f method is provided as a more convenient way to translate
2105 * a BigDecimal to a Float.)
2106 */
2107 static VALUE
BigDecimal_split(VALUE self)2108 BigDecimal_split(VALUE self)
2109 {
2110 ENTER(5);
2111 Real *vp;
2112 VALUE obj,str;
2113 ssize_t e, s;
2114 char *psz1;
2115
2116 GUARD_OBJ(vp, GetVpValue(self, 1));
2117 str = rb_str_new(0, VpNumOfChars(vp, "E"));
2118 psz1 = RSTRING_PTR(str);
2119 VpSzMantissa(vp, psz1);
2120 s = 1;
2121 if(psz1[0] == '-') {
2122 size_t len = strlen(psz1 + 1);
2123
2124 memmove(psz1, psz1 + 1, len);
2125 psz1[len] = '\0';
2126 s = -1;
2127 }
2128 if (psz1[0] == 'N') s = 0; /* NaN */
2129 e = VpExponent10(vp);
2130 obj = rb_ary_new2(4);
2131 rb_ary_push(obj, INT2FIX(s));
2132 rb_ary_push(obj, str);
2133 rb_str_resize(str, strlen(psz1));
2134 rb_ary_push(obj, INT2FIX(10));
2135 rb_ary_push(obj, INT2NUM(e));
2136 return obj;
2137 }
2138
2139 /* Returns the exponent of the BigDecimal number, as an Integer.
2140 *
2141 * If the number can be represented as 0.xxxxxx*10**n where xxxxxx is a string
2142 * of digits with no leading zeros, then n is the exponent.
2143 */
2144 static VALUE
BigDecimal_exponent(VALUE self)2145 BigDecimal_exponent(VALUE self)
2146 {
2147 ssize_t e = VpExponent10(GetVpValue(self, 1));
2148 return INT2NUM(e);
2149 }
2150
2151 /* Returns a string representation of self.
2152 *
2153 * BigDecimal("1234.5678").inspect
2154 * #=> "0.12345678e4"
2155 */
2156 static VALUE
BigDecimal_inspect(VALUE self)2157 BigDecimal_inspect(VALUE self)
2158 {
2159 ENTER(5);
2160 Real *vp;
2161 volatile VALUE str;
2162 size_t nc;
2163
2164 GUARD_OBJ(vp, GetVpValue(self, 1));
2165 nc = VpNumOfChars(vp, "E");
2166
2167 str = rb_str_new(0, nc);
2168 VpToString(vp, RSTRING_PTR(str), 0, 0);
2169 rb_str_resize(str, strlen(RSTRING_PTR(str)));
2170 return str;
2171 }
2172
2173 static VALUE BigMath_s_exp(VALUE, VALUE, VALUE);
2174 static VALUE BigMath_s_log(VALUE, VALUE, VALUE);
2175
2176 #define BigMath_exp(x, n) BigMath_s_exp(rb_mBigMath, (x), (n))
2177 #define BigMath_log(x, n) BigMath_s_log(rb_mBigMath, (x), (n))
2178
2179 inline static int
is_integer(VALUE x)2180 is_integer(VALUE x)
2181 {
2182 return (RB_TYPE_P(x, T_FIXNUM) || RB_TYPE_P(x, T_BIGNUM));
2183 }
2184
2185 inline static int
is_negative(VALUE x)2186 is_negative(VALUE x)
2187 {
2188 if (FIXNUM_P(x)) {
2189 return FIX2LONG(x) < 0;
2190 }
2191 else if (RB_TYPE_P(x, T_BIGNUM)) {
2192 return FIX2INT(rb_big_cmp(x, INT2FIX(0))) < 0;
2193 }
2194 else if (RB_TYPE_P(x, T_FLOAT)) {
2195 return RFLOAT_VALUE(x) < 0.0;
2196 }
2197 return RTEST(rb_funcall(x, '<', 1, INT2FIX(0)));
2198 }
2199
2200 #define is_positive(x) (!is_negative(x))
2201
2202 inline static int
is_zero(VALUE x)2203 is_zero(VALUE x)
2204 {
2205 VALUE num;
2206
2207 switch (TYPE(x)) {
2208 case T_FIXNUM:
2209 return FIX2LONG(x) == 0;
2210
2211 case T_BIGNUM:
2212 return Qfalse;
2213
2214 case T_RATIONAL:
2215 num = rb_rational_num(x);
2216 return FIXNUM_P(num) && FIX2LONG(num) == 0;
2217
2218 default:
2219 break;
2220 }
2221
2222 return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(0)));
2223 }
2224
2225 inline static int
is_one(VALUE x)2226 is_one(VALUE x)
2227 {
2228 VALUE num, den;
2229
2230 switch (TYPE(x)) {
2231 case T_FIXNUM:
2232 return FIX2LONG(x) == 1;
2233
2234 case T_BIGNUM:
2235 return Qfalse;
2236
2237 case T_RATIONAL:
2238 num = rb_rational_num(x);
2239 den = rb_rational_den(x);
2240 return FIXNUM_P(den) && FIX2LONG(den) == 1 &&
2241 FIXNUM_P(num) && FIX2LONG(num) == 1;
2242
2243 default:
2244 break;
2245 }
2246
2247 return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(1)));
2248 }
2249
2250 inline static int
is_even(VALUE x)2251 is_even(VALUE x)
2252 {
2253 switch (TYPE(x)) {
2254 case T_FIXNUM:
2255 return (FIX2LONG(x) % 2) == 0;
2256
2257 case T_BIGNUM:
2258 {
2259 unsigned long l;
2260 rb_big_pack(x, &l, 1);
2261 return l % 2 == 0;
2262 }
2263
2264 default:
2265 break;
2266 }
2267
2268 return 0;
2269 }
2270
2271 static VALUE
rmpd_power_by_big_decimal(Real const * x,Real const * exp,ssize_t const n)2272 rmpd_power_by_big_decimal(Real const* x, Real const* exp, ssize_t const n)
2273 {
2274 VALUE log_x, multiplied, y;
2275 volatile VALUE obj = exp->obj;
2276
2277 if (VpIsZero(exp)) {
2278 return ToValue(VpCreateRbObject(n, "1"));
2279 }
2280
2281 log_x = BigMath_log(x->obj, SSIZET2NUM(n+1));
2282 multiplied = BigDecimal_mult2(exp->obj, log_x, SSIZET2NUM(n+1));
2283 y = BigMath_exp(multiplied, SSIZET2NUM(n));
2284 RB_GC_GUARD(obj);
2285
2286 return y;
2287 }
2288
2289 /* call-seq:
2290 * power(n)
2291 * power(n, prec)
2292 *
2293 * Returns the value raised to the power of n.
2294 *
2295 * Note that n must be an Integer.
2296 *
2297 * Also available as the operator **.
2298 */
2299 static VALUE
BigDecimal_power(int argc,VALUE * argv,VALUE self)2300 BigDecimal_power(int argc, VALUE*argv, VALUE self)
2301 {
2302 ENTER(5);
2303 VALUE vexp, prec;
2304 Real* exp = NULL;
2305 Real *x, *y;
2306 ssize_t mp, ma, n;
2307 SIGNED_VALUE int_exp;
2308 double d;
2309
2310 rb_scan_args(argc, argv, "11", &vexp, &prec);
2311
2312 GUARD_OBJ(x, GetVpValue(self, 1));
2313 n = NIL_P(prec) ? (ssize_t)(x->Prec*VpBaseFig()) : NUM2SSIZET(prec);
2314
2315 if (VpIsNaN(x)) {
2316 y = VpCreateRbObject(n, "0");
2317 RB_GC_GUARD(y->obj);
2318 VpSetNaN(y);
2319 return ToValue(y);
2320 }
2321
2322 retry:
2323 switch (TYPE(vexp)) {
2324 case T_FIXNUM:
2325 break;
2326
2327 case T_BIGNUM:
2328 break;
2329
2330 case T_FLOAT:
2331 d = RFLOAT_VALUE(vexp);
2332 if (d == round(d)) {
2333 if (FIXABLE(d)) {
2334 vexp = LONG2FIX((long)d);
2335 }
2336 else {
2337 vexp = rb_dbl2big(d);
2338 }
2339 goto retry;
2340 }
2341 exp = GetVpValueWithPrec(vexp, DBL_DIG+1, 1);
2342 break;
2343
2344 case T_RATIONAL:
2345 if (is_zero(rb_rational_num(vexp))) {
2346 if (is_positive(vexp)) {
2347 vexp = INT2FIX(0);
2348 goto retry;
2349 }
2350 }
2351 else if (is_one(rb_rational_den(vexp))) {
2352 vexp = rb_rational_num(vexp);
2353 goto retry;
2354 }
2355 exp = GetVpValueWithPrec(vexp, n, 1);
2356 break;
2357
2358 case T_DATA:
2359 if (is_kind_of_BigDecimal(vexp)) {
2360 VALUE zero = INT2FIX(0);
2361 VALUE rounded = BigDecimal_round(1, &zero, vexp);
2362 if (RTEST(BigDecimal_eq(vexp, rounded))) {
2363 vexp = BigDecimal_to_i(vexp);
2364 goto retry;
2365 }
2366 exp = DATA_PTR(vexp);
2367 break;
2368 }
2369 /* fall through */
2370 default:
2371 rb_raise(rb_eTypeError,
2372 "wrong argument type %"PRIsVALUE" (expected scalar Numeric)",
2373 RB_OBJ_CLASSNAME(vexp));
2374 }
2375
2376 if (VpIsZero(x)) {
2377 if (is_negative(vexp)) {
2378 y = VpCreateRbObject(n, "#0");
2379 RB_GC_GUARD(y->obj);
2380 if (BIGDECIMAL_NEGATIVE_P(x)) {
2381 if (is_integer(vexp)) {
2382 if (is_even(vexp)) {
2383 /* (-0) ** (-even_integer) -> Infinity */
2384 VpSetPosInf(y);
2385 }
2386 else {
2387 /* (-0) ** (-odd_integer) -> -Infinity */
2388 VpSetNegInf(y);
2389 }
2390 }
2391 else {
2392 /* (-0) ** (-non_integer) -> Infinity */
2393 VpSetPosInf(y);
2394 }
2395 }
2396 else {
2397 /* (+0) ** (-num) -> Infinity */
2398 VpSetPosInf(y);
2399 }
2400 return ToValue(y);
2401 }
2402 else if (is_zero(vexp)) {
2403 return ToValue(VpCreateRbObject(n, "1"));
2404 }
2405 else {
2406 return ToValue(VpCreateRbObject(n, "0"));
2407 }
2408 }
2409
2410 if (is_zero(vexp)) {
2411 return ToValue(VpCreateRbObject(n, "1"));
2412 }
2413 else if (is_one(vexp)) {
2414 return self;
2415 }
2416
2417 if (VpIsInf(x)) {
2418 if (is_negative(vexp)) {
2419 if (BIGDECIMAL_NEGATIVE_P(x)) {
2420 if (is_integer(vexp)) {
2421 if (is_even(vexp)) {
2422 /* (-Infinity) ** (-even_integer) -> +0 */
2423 return ToValue(VpCreateRbObject(n, "0"));
2424 }
2425 else {
2426 /* (-Infinity) ** (-odd_integer) -> -0 */
2427 return ToValue(VpCreateRbObject(n, "-0"));
2428 }
2429 }
2430 else {
2431 /* (-Infinity) ** (-non_integer) -> -0 */
2432 return ToValue(VpCreateRbObject(n, "-0"));
2433 }
2434 }
2435 else {
2436 return ToValue(VpCreateRbObject(n, "0"));
2437 }
2438 }
2439 else {
2440 y = VpCreateRbObject(n, "0");
2441 if (BIGDECIMAL_NEGATIVE_P(x)) {
2442 if (is_integer(vexp)) {
2443 if (is_even(vexp)) {
2444 VpSetPosInf(y);
2445 }
2446 else {
2447 VpSetNegInf(y);
2448 }
2449 }
2450 else {
2451 /* TODO: support complex */
2452 rb_raise(rb_eMathDomainError,
2453 "a non-integral exponent for a negative base");
2454 }
2455 }
2456 else {
2457 VpSetPosInf(y);
2458 }
2459 return ToValue(y);
2460 }
2461 }
2462
2463 if (exp != NULL) {
2464 return rmpd_power_by_big_decimal(x, exp, n);
2465 }
2466 else if (RB_TYPE_P(vexp, T_BIGNUM)) {
2467 VALUE abs_value = BigDecimal_abs(self);
2468 if (is_one(abs_value)) {
2469 return ToValue(VpCreateRbObject(n, "1"));
2470 }
2471 else if (RTEST(rb_funcall(abs_value, '<', 1, INT2FIX(1)))) {
2472 if (is_negative(vexp)) {
2473 y = VpCreateRbObject(n, "0");
2474 if (is_even(vexp)) {
2475 VpSetInf(y, VpGetSign(x));
2476 }
2477 else {
2478 VpSetInf(y, -VpGetSign(x));
2479 }
2480 return ToValue(y);
2481 }
2482 else if (BIGDECIMAL_NEGATIVE_P(x) && is_even(vexp)) {
2483 return ToValue(VpCreateRbObject(n, "-0"));
2484 }
2485 else {
2486 return ToValue(VpCreateRbObject(n, "0"));
2487 }
2488 }
2489 else {
2490 if (is_positive(vexp)) {
2491 y = VpCreateRbObject(n, "0");
2492 if (is_even(vexp)) {
2493 VpSetInf(y, VpGetSign(x));
2494 }
2495 else {
2496 VpSetInf(y, -VpGetSign(x));
2497 }
2498 return ToValue(y);
2499 }
2500 else if (BIGDECIMAL_NEGATIVE_P(x) && is_even(vexp)) {
2501 return ToValue(VpCreateRbObject(n, "-0"));
2502 }
2503 else {
2504 return ToValue(VpCreateRbObject(n, "0"));
2505 }
2506 }
2507 }
2508
2509 int_exp = FIX2LONG(vexp);
2510 ma = int_exp;
2511 if (ma < 0) ma = -ma;
2512 if (ma == 0) ma = 1;
2513
2514 if (VpIsDef(x)) {
2515 mp = x->Prec * (VpBaseFig() + 1);
2516 GUARD_OBJ(y, VpCreateRbObject(mp * (ma + 1), "0"));
2517 }
2518 else {
2519 GUARD_OBJ(y, VpCreateRbObject(1, "0"));
2520 }
2521 VpPower(y, x, int_exp);
2522 if (!NIL_P(prec) && VpIsDef(y)) {
2523 VpMidRound(y, VpGetRoundMode(), n);
2524 }
2525 return ToValue(y);
2526 }
2527
2528 /* call-seq:
2529 * a ** n -> bigdecimal
2530 *
2531 * Returns the value raised to the power of n.
2532 *
2533 * See BigDecimal#power.
2534 */
2535 static VALUE
BigDecimal_power_op(VALUE self,VALUE exp)2536 BigDecimal_power_op(VALUE self, VALUE exp)
2537 {
2538 return BigDecimal_power(1, &exp, self);
2539 }
2540
2541 /* :nodoc:
2542 *
2543 * private method for dup and clone the provided BigDecimal +other+
2544 */
2545 static VALUE
BigDecimal_initialize_copy(VALUE self,VALUE other)2546 BigDecimal_initialize_copy(VALUE self, VALUE other)
2547 {
2548 Real *pv = rb_check_typeddata(self, &BigDecimal_data_type);
2549 Real *x = rb_check_typeddata(other, &BigDecimal_data_type);
2550
2551 if (self != other) {
2552 DATA_PTR(self) = VpCopy(pv, x);
2553 }
2554 return self;
2555 }
2556
2557 static VALUE
BigDecimal_clone(VALUE self)2558 BigDecimal_clone(VALUE self)
2559 {
2560 return self;
2561 }
2562
2563 static int
opts_exception_p(VALUE opts)2564 opts_exception_p(VALUE opts)
2565 {
2566 static ID kwds[1];
2567 VALUE exception;
2568 if (!kwds[0]) {
2569 kwds[0] = rb_intern_const("exception");
2570 }
2571 rb_get_kwargs(opts, kwds, 0, 1, &exception);
2572 return exception != Qfalse;
2573 }
2574
2575 static Real *
BigDecimal_new(int argc,VALUE * argv)2576 BigDecimal_new(int argc, VALUE *argv)
2577 {
2578 size_t mf;
2579 VALUE opts = Qnil;
2580 VALUE nFig;
2581 VALUE iniValue;
2582 double d;
2583 int exc;
2584
2585 argc = rb_scan_args(argc, argv, "11:", &iniValue, &nFig, &opts);
2586 exc = opts_exception_p(opts);
2587
2588 if (argc == 1) {
2589 mf = 0;
2590 }
2591 else {
2592 /* expand GetPrecisionInt for exception suppression */
2593 ssize_t n = NUM2INT(nFig);
2594 if (n < 0) {
2595 if (!exc) {
2596 return NULL;
2597 }
2598 rb_raise(rb_eArgError, "negative precision");
2599 }
2600 mf = (size_t)n;
2601 }
2602
2603 if (SPECIAL_CONST_P(iniValue)) {
2604 switch (iniValue) {
2605 case Qnil:
2606 if (!exc) return NULL;
2607 rb_raise(rb_eTypeError, "can't convert nil into BigDecimal");
2608 case Qtrue:
2609 if (!exc) return NULL;
2610 rb_raise(rb_eTypeError, "can't convert true into BigDecimal");
2611 case Qfalse:
2612 if (!exc) return NULL;
2613 rb_raise(rb_eTypeError, "can't convert false into BigDecimal");
2614 default:
2615 break;
2616 }
2617 }
2618
2619 switch (TYPE(iniValue)) {
2620 case T_DATA:
2621 if (is_kind_of_BigDecimal(iniValue)) {
2622 return DATA_PTR(iniValue);
2623 }
2624 break;
2625
2626 case T_FIXNUM:
2627 /* fall through */
2628 case T_BIGNUM:
2629 return GetVpValue(iniValue, 1);
2630
2631 case T_FLOAT:
2632 d = RFLOAT_VALUE(iniValue);
2633 if (!isfinite(d)) {
2634 Real *pv = VpCreateRbObject(1, NULL);
2635 VpDtoV(pv, d);
2636 return pv;
2637 }
2638 if (mf > DBL_DIG+1) {
2639 if (!exc) {
2640 return NULL;
2641 }
2642 rb_raise(rb_eArgError, "precision too large.");
2643 }
2644 /* fall through */
2645 case T_RATIONAL:
2646 if (NIL_P(nFig)) {
2647 if (!exc) {
2648 return NULL;
2649 }
2650 rb_raise(rb_eArgError,
2651 "can't omit precision for a %"PRIsVALUE".",
2652 RB_OBJ_CLASSNAME(iniValue));
2653 }
2654 return GetVpValueWithPrec(iniValue, mf, 1);
2655
2656 case T_STRING:
2657 /* fall through */
2658 default:
2659 break;
2660 }
2661 /* TODO: support to_d */
2662 if (!exc) {
2663 iniValue = rb_check_convert_type(iniValue, T_STRING, "String", "to_str");
2664 if (NIL_P(iniValue)) return NULL;
2665 }
2666 StringValueCStr(iniValue);
2667 return VpAlloc(mf, RSTRING_PTR(iniValue), 1, exc);
2668 }
2669
2670 /* call-seq:
2671 * BigDecimal(initial, digits, exception: true)
2672 *
2673 * Create a new BigDecimal object.
2674 *
2675 * initial:: The initial value, as an Integer, a Float, a Rational,
2676 * a BigDecimal, or a String.
2677 *
2678 * If it is a String, spaces are ignored and unrecognized characters
2679 * terminate the value.
2680 *
2681 * digits:: The number of significant digits, as an Integer. If omitted or 0,
2682 * the number of significant digits is determined from the initial
2683 * value.
2684 *
2685 * The actual number of significant digits used in computation is
2686 * usually larger than the specified number.
2687 *
2688 * exception:: Whether an exception should be raised on invalid arguments.
2689 * +true+ by default, if passed +false+, just returns +nil+
2690 * for invalid.
2691 *
2692 *
2693 * ==== Exceptions
2694 *
2695 * TypeError:: If the +initial+ type is neither Integer, Float,
2696 * Rational, nor BigDecimal, this exception is raised.
2697 *
2698 * TypeError:: If the +digits+ is not an Integer, this exception is raised.
2699 *
2700 * ArgumentError:: If +initial+ is a Float, and the +digits+ is larger than
2701 * Float::DIG + 1, this exception is raised.
2702 *
2703 * ArgumentError:: If the +initial+ is a Float or Rational, and the +digits+
2704 * value is omitted, this exception is raised.
2705 */
2706 static VALUE
f_BigDecimal(int argc,VALUE * argv,VALUE self)2707 f_BigDecimal(int argc, VALUE *argv, VALUE self)
2708 {
2709 ENTER(1);
2710 Real *pv;
2711 VALUE obj;
2712
2713 obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, 0);
2714 pv = BigDecimal_new(argc, argv);
2715 if (pv == NULL) return Qnil;
2716 SAVE(pv);
2717 if (ToValue(pv)) pv = VpCopy(NULL, pv);
2718 RTYPEDDATA_DATA(obj) = pv;
2719 RB_OBJ_FREEZE(obj);
2720 return pv->obj = obj;
2721 }
2722
2723 /* call-seq:
2724 * BigDecimal.limit(digits)
2725 *
2726 * Limit the number of significant digits in newly created BigDecimal
2727 * numbers to the specified value. Rounding is performed as necessary,
2728 * as specified by BigDecimal.mode.
2729 *
2730 * A limit of 0, the default, means no upper limit.
2731 *
2732 * The limit specified by this method takes less priority over any limit
2733 * specified to instance methods such as ceil, floor, truncate, or round.
2734 */
2735 static VALUE
BigDecimal_limit(int argc,VALUE * argv,VALUE self)2736 BigDecimal_limit(int argc, VALUE *argv, VALUE self)
2737 {
2738 VALUE nFig;
2739 VALUE nCur = INT2NUM(VpGetPrecLimit());
2740
2741 if (rb_scan_args(argc, argv, "01", &nFig) == 1) {
2742 int nf;
2743 if (NIL_P(nFig)) return nCur;
2744 nf = NUM2INT(nFig);
2745 if (nf < 0) {
2746 rb_raise(rb_eArgError, "argument must be positive");
2747 }
2748 VpSetPrecLimit(nf);
2749 }
2750 return nCur;
2751 }
2752
2753 /* Returns the sign of the value.
2754 *
2755 * Returns a positive value if > 0, a negative value if < 0, and a
2756 * zero if == 0.
2757 *
2758 * The specific value returned indicates the type and sign of the BigDecimal,
2759 * as follows:
2760 *
2761 * BigDecimal::SIGN_NaN:: value is Not a Number
2762 * BigDecimal::SIGN_POSITIVE_ZERO:: value is +0
2763 * BigDecimal::SIGN_NEGATIVE_ZERO:: value is -0
2764 * BigDecimal::SIGN_POSITIVE_INFINITE:: value is +Infinity
2765 * BigDecimal::SIGN_NEGATIVE_INFINITE:: value is -Infinity
2766 * BigDecimal::SIGN_POSITIVE_FINITE:: value is positive
2767 * BigDecimal::SIGN_NEGATIVE_FINITE:: value is negative
2768 */
2769 static VALUE
BigDecimal_sign(VALUE self)2770 BigDecimal_sign(VALUE self)
2771 { /* sign */
2772 int s = GetVpValue(self, 1)->sign;
2773 return INT2FIX(s);
2774 }
2775
2776 /*
2777 * call-seq: BigDecimal.save_exception_mode { ... }
2778 *
2779 * Execute the provided block, but preserve the exception mode
2780 *
2781 * BigDecimal.save_exception_mode do
2782 * BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false)
2783 * BigDecimal.mode(BigDecimal::EXCEPTION_NaN, false)
2784 *
2785 * BigDecimal(BigDecimal('Infinity'))
2786 * BigDecimal(BigDecimal('-Infinity'))
2787 * BigDecimal(BigDecimal('NaN'))
2788 * end
2789 *
2790 * For use with the BigDecimal::EXCEPTION_*
2791 *
2792 * See BigDecimal.mode
2793 */
2794 static VALUE
BigDecimal_save_exception_mode(VALUE self)2795 BigDecimal_save_exception_mode(VALUE self)
2796 {
2797 unsigned short const exception_mode = VpGetException();
2798 int state;
2799 VALUE ret = rb_protect(rb_yield, Qnil, &state);
2800 VpSetException(exception_mode);
2801 if (state) rb_jump_tag(state);
2802 return ret;
2803 }
2804
2805 /*
2806 * call-seq: BigDecimal.save_rounding_mode { ... }
2807 *
2808 * Execute the provided block, but preserve the rounding mode
2809 *
2810 * BigDecimal.save_rounding_mode do
2811 * BigDecimal.mode(BigDecimal::ROUND_MODE, :up)
2812 * puts BigDecimal.mode(BigDecimal::ROUND_MODE)
2813 * end
2814 *
2815 * For use with the BigDecimal::ROUND_*
2816 *
2817 * See BigDecimal.mode
2818 */
2819 static VALUE
BigDecimal_save_rounding_mode(VALUE self)2820 BigDecimal_save_rounding_mode(VALUE self)
2821 {
2822 unsigned short const round_mode = VpGetRoundMode();
2823 int state;
2824 VALUE ret = rb_protect(rb_yield, Qnil, &state);
2825 VpSetRoundMode(round_mode);
2826 if (state) rb_jump_tag(state);
2827 return ret;
2828 }
2829
2830 /*
2831 * call-seq: BigDecimal.save_limit { ... }
2832 *
2833 * Execute the provided block, but preserve the precision limit
2834 *
2835 * BigDecimal.limit(100)
2836 * puts BigDecimal.limit
2837 * BigDecimal.save_limit do
2838 * BigDecimal.limit(200)
2839 * puts BigDecimal.limit
2840 * end
2841 * puts BigDecimal.limit
2842 *
2843 */
2844 static VALUE
BigDecimal_save_limit(VALUE self)2845 BigDecimal_save_limit(VALUE self)
2846 {
2847 size_t const limit = VpGetPrecLimit();
2848 int state;
2849 VALUE ret = rb_protect(rb_yield, Qnil, &state);
2850 VpSetPrecLimit(limit);
2851 if (state) rb_jump_tag(state);
2852 return ret;
2853 }
2854
2855 /* call-seq:
2856 * BigMath.exp(decimal, numeric) -> BigDecimal
2857 *
2858 * Computes the value of e (the base of natural logarithms) raised to the
2859 * power of +decimal+, to the specified number of digits of precision.
2860 *
2861 * If +decimal+ is infinity, returns Infinity.
2862 *
2863 * If +decimal+ is NaN, returns NaN.
2864 */
2865 static VALUE
BigMath_s_exp(VALUE klass,VALUE x,VALUE vprec)2866 BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
2867 {
2868 ssize_t prec, n, i;
2869 Real* vx = NULL;
2870 VALUE one, d, y;
2871 int negative = 0;
2872 int infinite = 0;
2873 int nan = 0;
2874 double flo;
2875
2876 prec = NUM2SSIZET(vprec);
2877 if (prec <= 0) {
2878 rb_raise(rb_eArgError, "Zero or negative precision for exp");
2879 }
2880
2881 /* TODO: the following switch statement is almost same as one in the
2882 * BigDecimalCmp function. */
2883 switch (TYPE(x)) {
2884 case T_DATA:
2885 if (!is_kind_of_BigDecimal(x)) break;
2886 vx = DATA_PTR(x);
2887 negative = BIGDECIMAL_NEGATIVE_P(vx);
2888 infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
2889 nan = VpIsNaN(vx);
2890 break;
2891
2892 case T_FIXNUM:
2893 /* fall through */
2894 case T_BIGNUM:
2895 vx = GetVpValue(x, 0);
2896 break;
2897
2898 case T_FLOAT:
2899 flo = RFLOAT_VALUE(x);
2900 negative = flo < 0;
2901 infinite = isinf(flo);
2902 nan = isnan(flo);
2903 if (!infinite && !nan) {
2904 vx = GetVpValueWithPrec(x, DBL_DIG+1, 0);
2905 }
2906 break;
2907
2908 case T_RATIONAL:
2909 vx = GetVpValueWithPrec(x, prec, 0);
2910 break;
2911
2912 default:
2913 break;
2914 }
2915 if (infinite) {
2916 if (negative) {
2917 return ToValue(GetVpValueWithPrec(INT2FIX(0), prec, 1));
2918 }
2919 else {
2920 Real* vy;
2921 vy = VpCreateRbObject(prec, "#0");
2922 VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
2923 RB_GC_GUARD(vy->obj);
2924 return ToValue(vy);
2925 }
2926 }
2927 else if (nan) {
2928 Real* vy;
2929 vy = VpCreateRbObject(prec, "#0");
2930 VpSetNaN(vy);
2931 RB_GC_GUARD(vy->obj);
2932 return ToValue(vy);
2933 }
2934 else if (vx == NULL) {
2935 cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
2936 }
2937 x = vx->obj;
2938
2939 n = prec + rmpd_double_figures();
2940 negative = BIGDECIMAL_NEGATIVE_P(vx);
2941 if (negative) {
2942 VpSetSign(vx, 1);
2943 }
2944
2945 one = ToValue(VpCreateRbObject(1, "1"));
2946 y = one;
2947 d = y;
2948 i = 1;
2949
2950 while (!VpIsZero((Real*)DATA_PTR(d))) {
2951 SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
2952 SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
2953 ssize_t m = n - vabs(ey - ed);
2954
2955 rb_thread_check_ints();
2956
2957 if (m <= 0) {
2958 break;
2959 }
2960 else if ((size_t)m < rmpd_double_figures()) {
2961 m = rmpd_double_figures();
2962 }
2963
2964 d = BigDecimal_mult(d, x); /* d <- d * x */
2965 d = BigDecimal_div2(d, SSIZET2NUM(i), SSIZET2NUM(m)); /* d <- d / i */
2966 y = BigDecimal_add(y, d); /* y <- y + d */
2967 ++i; /* i <- i + 1 */
2968 }
2969
2970 if (negative) {
2971 return BigDecimal_div2(one, y, vprec);
2972 }
2973 else {
2974 vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
2975 return BigDecimal_round(1, &vprec, y);
2976 }
2977
2978 RB_GC_GUARD(one);
2979 RB_GC_GUARD(x);
2980 RB_GC_GUARD(y);
2981 RB_GC_GUARD(d);
2982 }
2983
2984 /* call-seq:
2985 * BigMath.log(decimal, numeric) -> BigDecimal
2986 *
2987 * Computes the natural logarithm of +decimal+ to the specified number of
2988 * digits of precision, +numeric+.
2989 *
2990 * If +decimal+ is zero or negative, raises Math::DomainError.
2991 *
2992 * If +decimal+ is positive infinity, returns Infinity.
2993 *
2994 * If +decimal+ is NaN, returns NaN.
2995 */
2996 static VALUE
BigMath_s_log(VALUE klass,VALUE x,VALUE vprec)2997 BigMath_s_log(VALUE klass, VALUE x, VALUE vprec)
2998 {
2999 ssize_t prec, n, i;
3000 SIGNED_VALUE expo;
3001 Real* vx = NULL;
3002 VALUE vn, one, two, w, x2, y, d;
3003 int zero = 0;
3004 int negative = 0;
3005 int infinite = 0;
3006 int nan = 0;
3007 double flo;
3008 long fix;
3009
3010 if (!is_integer(vprec)) {
3011 rb_raise(rb_eArgError, "precision must be an Integer");
3012 }
3013
3014 prec = NUM2SSIZET(vprec);
3015 if (prec <= 0) {
3016 rb_raise(rb_eArgError, "Zero or negative precision for exp");
3017 }
3018
3019 /* TODO: the following switch statement is almost same as one in the
3020 * BigDecimalCmp function. */
3021 switch (TYPE(x)) {
3022 case T_DATA:
3023 if (!is_kind_of_BigDecimal(x)) break;
3024 vx = DATA_PTR(x);
3025 zero = VpIsZero(vx);
3026 negative = BIGDECIMAL_NEGATIVE_P(vx);
3027 infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
3028 nan = VpIsNaN(vx);
3029 break;
3030
3031 case T_FIXNUM:
3032 fix = FIX2LONG(x);
3033 zero = fix == 0;
3034 negative = fix < 0;
3035 goto get_vp_value;
3036
3037 case T_BIGNUM:
3038 i = FIX2INT(rb_big_cmp(x, INT2FIX(0)));
3039 zero = i == 0;
3040 negative = i < 0;
3041 get_vp_value:
3042 if (zero || negative) break;
3043 vx = GetVpValue(x, 0);
3044 break;
3045
3046 case T_FLOAT:
3047 flo = RFLOAT_VALUE(x);
3048 zero = flo == 0;
3049 negative = flo < 0;
3050 infinite = isinf(flo);
3051 nan = isnan(flo);
3052 if (!zero && !negative && !infinite && !nan) {
3053 vx = GetVpValueWithPrec(x, DBL_DIG+1, 1);
3054 }
3055 break;
3056
3057 case T_RATIONAL:
3058 zero = RRATIONAL_ZERO_P(x);
3059 negative = RRATIONAL_NEGATIVE_P(x);
3060 if (zero || negative) break;
3061 vx = GetVpValueWithPrec(x, prec, 1);
3062 break;
3063
3064 case T_COMPLEX:
3065 rb_raise(rb_eMathDomainError,
3066 "Complex argument for BigMath.log");
3067
3068 default:
3069 break;
3070 }
3071 if (infinite && !negative) {
3072 Real* vy;
3073 vy = VpCreateRbObject(prec, "#0");
3074 RB_GC_GUARD(vy->obj);
3075 VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
3076 return ToValue(vy);
3077 }
3078 else if (nan) {
3079 Real* vy;
3080 vy = VpCreateRbObject(prec, "#0");
3081 RB_GC_GUARD(vy->obj);
3082 VpSetNaN(vy);
3083 return ToValue(vy);
3084 }
3085 else if (zero || negative) {
3086 rb_raise(rb_eMathDomainError,
3087 "Zero or negative argument for log");
3088 }
3089 else if (vx == NULL) {
3090 cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
3091 }
3092 x = ToValue(vx);
3093
3094 RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
3095 RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2"));
3096
3097 n = prec + rmpd_double_figures();
3098 RB_GC_GUARD(vn) = SSIZET2NUM(n);
3099 expo = VpExponent10(vx);
3100 if (expo < 0 || expo >= 3) {
3101 char buf[DECIMAL_SIZE_OF_BITS(SIZEOF_VALUE * CHAR_BIT) + 4];
3102 snprintf(buf, sizeof(buf), "1E%"PRIdVALUE, -expo);
3103 x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
3104 }
3105 else {
3106 expo = 0;
3107 }
3108 w = BigDecimal_sub(x, one);
3109 x = BigDecimal_div2(w, BigDecimal_add(x, one), vn);
3110 RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn);
3111 RB_GC_GUARD(y) = x;
3112 RB_GC_GUARD(d) = y;
3113 i = 1;
3114 while (!VpIsZero((Real*)DATA_PTR(d))) {
3115 SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
3116 SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
3117 ssize_t m = n - vabs(ey - ed);
3118 if (m <= 0) {
3119 break;
3120 }
3121 else if ((size_t)m < rmpd_double_figures()) {
3122 m = rmpd_double_figures();
3123 }
3124
3125 x = BigDecimal_mult2(x2, x, vn);
3126 i += 2;
3127 d = BigDecimal_div2(x, SSIZET2NUM(i), SSIZET2NUM(m));
3128 y = BigDecimal_add(y, d);
3129 }
3130
3131 y = BigDecimal_mult(y, two);
3132 if (expo != 0) {
3133 VALUE log10, vexpo, dy;
3134 log10 = BigMath_s_log(klass, INT2FIX(10), vprec);
3135 vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1));
3136 dy = BigDecimal_mult(log10, vexpo);
3137 y = BigDecimal_add(y, dy);
3138 }
3139
3140 return y;
3141 }
3142
3143 VALUE
rmpd_util_str_to_d(VALUE str)3144 rmpd_util_str_to_d(VALUE str)
3145 {
3146 ENTER(1);
3147 char const *c_str;
3148 Real *pv;
3149
3150 c_str = StringValueCStr(str);
3151 GUARD_OBJ(pv, VpAlloc(0, c_str, 0, 1));
3152 pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv);
3153 RB_OBJ_FREEZE(pv->obj);
3154 return pv->obj;
3155 }
3156
3157 /* Document-class: BigDecimal
3158 * BigDecimal provides arbitrary-precision floating point decimal arithmetic.
3159 *
3160 * == Introduction
3161 *
3162 * Ruby provides built-in support for arbitrary precision integer arithmetic.
3163 *
3164 * For example:
3165 *
3166 * 42**13 #=> 1265437718438866624512
3167 *
3168 * BigDecimal provides similar support for very large or very accurate floating
3169 * point numbers.
3170 *
3171 * Decimal arithmetic is also useful for general calculation, because it
3172 * provides the correct answers people expect--whereas normal binary floating
3173 * point arithmetic often introduces subtle errors because of the conversion
3174 * between base 10 and base 2.
3175 *
3176 * For example, try:
3177 *
3178 * sum = 0
3179 * 10_000.times do
3180 * sum = sum + 0.0001
3181 * end
3182 * print sum #=> 0.9999999999999062
3183 *
3184 * and contrast with the output from:
3185 *
3186 * require 'bigdecimal'
3187 *
3188 * sum = BigDecimal("0")
3189 * 10_000.times do
3190 * sum = sum + BigDecimal("0.0001")
3191 * end
3192 * print sum #=> 0.1E1
3193 *
3194 * Similarly:
3195 *
3196 * (BigDecimal("1.2") - BigDecimal("1.0")) == BigDecimal("0.2") #=> true
3197 *
3198 * (1.2 - 1.0) == 0.2 #=> false
3199 *
3200 * == Special features of accurate decimal arithmetic
3201 *
3202 * Because BigDecimal is more accurate than normal binary floating point
3203 * arithmetic, it requires some special values.
3204 *
3205 * === Infinity
3206 *
3207 * BigDecimal sometimes needs to return infinity, for example if you divide
3208 * a value by zero.
3209 *
3210 * BigDecimal("1.0") / BigDecimal("0.0") #=> Infinity
3211 * BigDecimal("-1.0") / BigDecimal("0.0") #=> -Infinity
3212 *
3213 * You can represent infinite numbers to BigDecimal using the strings
3214 * <code>'Infinity'</code>, <code>'+Infinity'</code> and
3215 * <code>'-Infinity'</code> (case-sensitive)
3216 *
3217 * === Not a Number
3218 *
3219 * When a computation results in an undefined value, the special value +NaN+
3220 * (for 'not a number') is returned.
3221 *
3222 * Example:
3223 *
3224 * BigDecimal("0.0") / BigDecimal("0.0") #=> NaN
3225 *
3226 * You can also create undefined values.
3227 *
3228 * NaN is never considered to be the same as any other value, even NaN itself:
3229 *
3230 * n = BigDecimal('NaN')
3231 * n == 0.0 #=> false
3232 * n == n #=> false
3233 *
3234 * === Positive and negative zero
3235 *
3236 * If a computation results in a value which is too small to be represented as
3237 * a BigDecimal within the currently specified limits of precision, zero must
3238 * be returned.
3239 *
3240 * If the value which is too small to be represented is negative, a BigDecimal
3241 * value of negative zero is returned.
3242 *
3243 * BigDecimal("1.0") / BigDecimal("-Infinity") #=> -0.0
3244 *
3245 * If the value is positive, a value of positive zero is returned.
3246 *
3247 * BigDecimal("1.0") / BigDecimal("Infinity") #=> 0.0
3248 *
3249 * (See BigDecimal.mode for how to specify limits of precision.)
3250 *
3251 * Note that +-0.0+ and +0.0+ are considered to be the same for the purposes of
3252 * comparison.
3253 *
3254 * Note also that in mathematics, there is no particular concept of negative
3255 * or positive zero; true mathematical zero has no sign.
3256 *
3257 * == bigdecimal/util
3258 *
3259 * When you require +bigdecimal/util+, the #to_d method will be
3260 * available on BigDecimal and the native Integer, Float, Rational,
3261 * and String classes:
3262 *
3263 * require 'bigdecimal/util'
3264 *
3265 * 42.to_d # => 0.42e2
3266 * 0.5.to_d # => 0.5e0
3267 * (2/3r).to_d(3) # => 0.667e0
3268 * "0.5".to_d # => 0.5e0
3269 *
3270 * == License
3271 *
3272 * Copyright (C) 2002 by Shigeo Kobayashi <shigeo@tinyforest.gr.jp>.
3273 *
3274 * BigDecimal is released under the Ruby and 2-clause BSD licenses.
3275 * See LICENSE.txt for details.
3276 *
3277 * Maintained by mrkn <mrkn@mrkn.jp> and ruby-core members.
3278 *
3279 * Documented by zzak <zachary@zacharyscott.net>, mathew <meta@pobox.com>, and
3280 * many other contributors.
3281 */
3282 void
Init_bigdecimal(void)3283 Init_bigdecimal(void)
3284 {
3285 VALUE arg;
3286
3287 id_BigDecimal_exception_mode = rb_intern_const("BigDecimal.exception_mode");
3288 id_BigDecimal_rounding_mode = rb_intern_const("BigDecimal.rounding_mode");
3289 id_BigDecimal_precision_limit = rb_intern_const("BigDecimal.precision_limit");
3290
3291 /* Initialize VP routines */
3292 VpInit(0UL);
3293
3294 /* Class and method registration */
3295 rb_cBigDecimal = rb_define_class("BigDecimal", rb_cNumeric);
3296
3297 /* Global function */
3298 rb_define_global_function("BigDecimal", f_BigDecimal, -1);
3299
3300 /* Class methods */
3301 rb_undef_method(CLASS_OF(rb_cBigDecimal), "allocate");
3302 rb_undef_method(CLASS_OF(rb_cBigDecimal), "new");
3303 rb_define_singleton_method(rb_cBigDecimal, "mode", BigDecimal_mode, -1);
3304 rb_define_singleton_method(rb_cBigDecimal, "limit", BigDecimal_limit, -1);
3305 rb_define_singleton_method(rb_cBigDecimal, "double_fig", BigDecimal_double_fig, 0);
3306 rb_define_singleton_method(rb_cBigDecimal, "_load", BigDecimal_load, 1);
3307
3308 rb_define_singleton_method(rb_cBigDecimal, "save_exception_mode", BigDecimal_save_exception_mode, 0);
3309 rb_define_singleton_method(rb_cBigDecimal, "save_rounding_mode", BigDecimal_save_rounding_mode, 0);
3310 rb_define_singleton_method(rb_cBigDecimal, "save_limit", BigDecimal_save_limit, 0);
3311
3312 /* Constants definition */
3313
3314 #ifndef RUBY_BIGDECIMAL_VERSION
3315 # error RUBY_BIGDECIMAL_VERSION is not defined
3316 #endif
3317 /*
3318 * The version of bigdecimal library
3319 */
3320 rb_define_const(rb_cBigDecimal, "VERSION", rb_str_new2(RUBY_BIGDECIMAL_VERSION));
3321
3322 /*
3323 * Base value used in internal calculations. On a 32 bit system, BASE
3324 * is 10000, indicating that calculation is done in groups of 4 digits.
3325 * (If it were larger, BASE**2 wouldn't fit in 32 bits, so you couldn't
3326 * guarantee that two groups could always be multiplied together without
3327 * overflow.)
3328 */
3329 rb_define_const(rb_cBigDecimal, "BASE", INT2FIX((SIGNED_VALUE)VpBaseVal()));
3330
3331 /* Exceptions */
3332
3333 /*
3334 * 0xff: Determines whether overflow, underflow or zero divide result in
3335 * an exception being thrown. See BigDecimal.mode.
3336 */
3337 rb_define_const(rb_cBigDecimal, "EXCEPTION_ALL", INT2FIX(VP_EXCEPTION_ALL));
3338
3339 /*
3340 * 0x02: Determines what happens when the result of a computation is not a
3341 * number (NaN). See BigDecimal.mode.
3342 */
3343 rb_define_const(rb_cBigDecimal, "EXCEPTION_NaN", INT2FIX(VP_EXCEPTION_NaN));
3344
3345 /*
3346 * 0x01: Determines what happens when the result of a computation is
3347 * infinity. See BigDecimal.mode.
3348 */
3349 rb_define_const(rb_cBigDecimal, "EXCEPTION_INFINITY", INT2FIX(VP_EXCEPTION_INFINITY));
3350
3351 /*
3352 * 0x04: Determines what happens when the result of a computation is an
3353 * underflow (a result too small to be represented). See BigDecimal.mode.
3354 */
3355 rb_define_const(rb_cBigDecimal, "EXCEPTION_UNDERFLOW", INT2FIX(VP_EXCEPTION_UNDERFLOW));
3356
3357 /*
3358 * 0x01: Determines what happens when the result of a computation is an
3359 * overflow (a result too large to be represented). See BigDecimal.mode.
3360 */
3361 rb_define_const(rb_cBigDecimal, "EXCEPTION_OVERFLOW", INT2FIX(VP_EXCEPTION_OVERFLOW));
3362
3363 /*
3364 * 0x10: Determines what happens when a division by zero is performed.
3365 * See BigDecimal.mode.
3366 */
3367 rb_define_const(rb_cBigDecimal, "EXCEPTION_ZERODIVIDE", INT2FIX(VP_EXCEPTION_ZERODIVIDE));
3368
3369 /*
3370 * 0x100: Determines what happens when a result must be rounded in order to
3371 * fit in the appropriate number of significant digits. See
3372 * BigDecimal.mode.
3373 */
3374 rb_define_const(rb_cBigDecimal, "ROUND_MODE", INT2FIX(VP_ROUND_MODE));
3375
3376 /* 1: Indicates that values should be rounded away from zero. See
3377 * BigDecimal.mode.
3378 */
3379 rb_define_const(rb_cBigDecimal, "ROUND_UP", INT2FIX(VP_ROUND_UP));
3380
3381 /* 2: Indicates that values should be rounded towards zero. See
3382 * BigDecimal.mode.
3383 */
3384 rb_define_const(rb_cBigDecimal, "ROUND_DOWN", INT2FIX(VP_ROUND_DOWN));
3385
3386 /* 3: Indicates that digits >= 5 should be rounded up, others rounded down.
3387 * See BigDecimal.mode. */
3388 rb_define_const(rb_cBigDecimal, "ROUND_HALF_UP", INT2FIX(VP_ROUND_HALF_UP));
3389
3390 /* 4: Indicates that digits >= 6 should be rounded up, others rounded down.
3391 * See BigDecimal.mode.
3392 */
3393 rb_define_const(rb_cBigDecimal, "ROUND_HALF_DOWN", INT2FIX(VP_ROUND_HALF_DOWN));
3394 /* 5: Round towards +Infinity. See BigDecimal.mode. */
3395 rb_define_const(rb_cBigDecimal, "ROUND_CEILING", INT2FIX(VP_ROUND_CEIL));
3396
3397 /* 6: Round towards -Infinity. See BigDecimal.mode. */
3398 rb_define_const(rb_cBigDecimal, "ROUND_FLOOR", INT2FIX(VP_ROUND_FLOOR));
3399
3400 /* 7: Round towards the even neighbor. See BigDecimal.mode. */
3401 rb_define_const(rb_cBigDecimal, "ROUND_HALF_EVEN", INT2FIX(VP_ROUND_HALF_EVEN));
3402
3403 /* 0: Indicates that a value is not a number. See BigDecimal.sign. */
3404 rb_define_const(rb_cBigDecimal, "SIGN_NaN", INT2FIX(VP_SIGN_NaN));
3405
3406 /* 1: Indicates that a value is +0. See BigDecimal.sign. */
3407 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_ZERO", INT2FIX(VP_SIGN_POSITIVE_ZERO));
3408
3409 /* -1: Indicates that a value is -0. See BigDecimal.sign. */
3410 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_ZERO", INT2FIX(VP_SIGN_NEGATIVE_ZERO));
3411
3412 /* 2: Indicates that a value is positive and finite. See BigDecimal.sign. */
3413 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_FINITE", INT2FIX(VP_SIGN_POSITIVE_FINITE));
3414
3415 /* -2: Indicates that a value is negative and finite. See BigDecimal.sign. */
3416 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_FINITE", INT2FIX(VP_SIGN_NEGATIVE_FINITE));
3417
3418 /* 3: Indicates that a value is positive and infinite. See BigDecimal.sign. */
3419 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_INFINITE", INT2FIX(VP_SIGN_POSITIVE_INFINITE));
3420
3421 /* -3: Indicates that a value is negative and infinite. See BigDecimal.sign. */
3422 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_INFINITE", INT2FIX(VP_SIGN_NEGATIVE_INFINITE));
3423
3424 arg = rb_str_new2("+Infinity");
3425 /* Positive infinity value. */
3426 rb_define_const(rb_cBigDecimal, "INFINITY", f_BigDecimal(1, &arg, rb_cBigDecimal));
3427 arg = rb_str_new2("NaN");
3428 /* 'Not a Number' value. */
3429 rb_define_const(rb_cBigDecimal, "NAN", f_BigDecimal(1, &arg, rb_cBigDecimal));
3430
3431
3432 /* instance methods */
3433 rb_define_method(rb_cBigDecimal, "initialize_copy", BigDecimal_initialize_copy, 1);
3434 rb_define_method(rb_cBigDecimal, "precs", BigDecimal_prec, 0);
3435
3436 rb_define_method(rb_cBigDecimal, "add", BigDecimal_add2, 2);
3437 rb_define_method(rb_cBigDecimal, "sub", BigDecimal_sub2, 2);
3438 rb_define_method(rb_cBigDecimal, "mult", BigDecimal_mult2, 2);
3439 rb_define_method(rb_cBigDecimal, "div", BigDecimal_div3, -1);
3440 rb_define_method(rb_cBigDecimal, "hash", BigDecimal_hash, 0);
3441 rb_define_method(rb_cBigDecimal, "to_s", BigDecimal_to_s, -1);
3442 rb_define_method(rb_cBigDecimal, "to_i", BigDecimal_to_i, 0);
3443 rb_define_method(rb_cBigDecimal, "to_int", BigDecimal_to_i, 0);
3444 rb_define_method(rb_cBigDecimal, "to_r", BigDecimal_to_r, 0);
3445 rb_define_method(rb_cBigDecimal, "split", BigDecimal_split, 0);
3446 rb_define_method(rb_cBigDecimal, "+", BigDecimal_add, 1);
3447 rb_define_method(rb_cBigDecimal, "-", BigDecimal_sub, 1);
3448 rb_define_method(rb_cBigDecimal, "+@", BigDecimal_uplus, 0);
3449 rb_define_method(rb_cBigDecimal, "-@", BigDecimal_neg, 0);
3450 rb_define_method(rb_cBigDecimal, "*", BigDecimal_mult, 1);
3451 rb_define_method(rb_cBigDecimal, "/", BigDecimal_div, 1);
3452 rb_define_method(rb_cBigDecimal, "quo", BigDecimal_div, 1);
3453 rb_define_method(rb_cBigDecimal, "%", BigDecimal_mod, 1);
3454 rb_define_method(rb_cBigDecimal, "modulo", BigDecimal_mod, 1);
3455 rb_define_method(rb_cBigDecimal, "remainder", BigDecimal_remainder, 1);
3456 rb_define_method(rb_cBigDecimal, "divmod", BigDecimal_divmod, 1);
3457 rb_define_method(rb_cBigDecimal, "clone", BigDecimal_clone, 0);
3458 rb_define_method(rb_cBigDecimal, "dup", BigDecimal_clone, 0);
3459 rb_define_method(rb_cBigDecimal, "to_f", BigDecimal_to_f, 0);
3460 rb_define_method(rb_cBigDecimal, "abs", BigDecimal_abs, 0);
3461 rb_define_method(rb_cBigDecimal, "sqrt", BigDecimal_sqrt, 1);
3462 rb_define_method(rb_cBigDecimal, "fix", BigDecimal_fix, 0);
3463 rb_define_method(rb_cBigDecimal, "round", BigDecimal_round, -1);
3464 rb_define_method(rb_cBigDecimal, "frac", BigDecimal_frac, 0);
3465 rb_define_method(rb_cBigDecimal, "floor", BigDecimal_floor, -1);
3466 rb_define_method(rb_cBigDecimal, "ceil", BigDecimal_ceil, -1);
3467 rb_define_method(rb_cBigDecimal, "power", BigDecimal_power, -1);
3468 rb_define_method(rb_cBigDecimal, "**", BigDecimal_power_op, 1);
3469 rb_define_method(rb_cBigDecimal, "<=>", BigDecimal_comp, 1);
3470 rb_define_method(rb_cBigDecimal, "==", BigDecimal_eq, 1);
3471 rb_define_method(rb_cBigDecimal, "===", BigDecimal_eq, 1);
3472 rb_define_method(rb_cBigDecimal, "eql?", BigDecimal_eq, 1);
3473 rb_define_method(rb_cBigDecimal, "<", BigDecimal_lt, 1);
3474 rb_define_method(rb_cBigDecimal, "<=", BigDecimal_le, 1);
3475 rb_define_method(rb_cBigDecimal, ">", BigDecimal_gt, 1);
3476 rb_define_method(rb_cBigDecimal, ">=", BigDecimal_ge, 1);
3477 rb_define_method(rb_cBigDecimal, "zero?", BigDecimal_zero, 0);
3478 rb_define_method(rb_cBigDecimal, "nonzero?", BigDecimal_nonzero, 0);
3479 rb_define_method(rb_cBigDecimal, "coerce", BigDecimal_coerce, 1);
3480 rb_define_method(rb_cBigDecimal, "inspect", BigDecimal_inspect, 0);
3481 rb_define_method(rb_cBigDecimal, "exponent", BigDecimal_exponent, 0);
3482 rb_define_method(rb_cBigDecimal, "sign", BigDecimal_sign, 0);
3483 rb_define_method(rb_cBigDecimal, "nan?", BigDecimal_IsNaN, 0);
3484 rb_define_method(rb_cBigDecimal, "infinite?", BigDecimal_IsInfinite, 0);
3485 rb_define_method(rb_cBigDecimal, "finite?", BigDecimal_IsFinite, 0);
3486 rb_define_method(rb_cBigDecimal, "truncate", BigDecimal_truncate, -1);
3487 rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1);
3488
3489 rb_mBigMath = rb_define_module("BigMath");
3490 rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2);
3491 rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2);
3492
3493 id_up = rb_intern_const("up");
3494 id_down = rb_intern_const("down");
3495 id_truncate = rb_intern_const("truncate");
3496 id_half_up = rb_intern_const("half_up");
3497 id_default = rb_intern_const("default");
3498 id_half_down = rb_intern_const("half_down");
3499 id_half_even = rb_intern_const("half_even");
3500 id_banker = rb_intern_const("banker");
3501 id_ceiling = rb_intern_const("ceiling");
3502 id_ceil = rb_intern_const("ceil");
3503 id_floor = rb_intern_const("floor");
3504 id_to_r = rb_intern_const("to_r");
3505 id_eq = rb_intern_const("==");
3506 id_half = rb_intern_const("half");
3507 }
3508
3509 /*
3510 *
3511 * ============================================================================
3512 *
3513 * vp_ routines begin from here.
3514 *
3515 * ============================================================================
3516 *
3517 */
3518 #ifdef BIGDECIMAL_DEBUG
3519 static int gfDebug = 1; /* Debug switch */
3520 #if 0
3521 static int gfCheckVal = 1; /* Value checking flag in VpNmlz() */
3522 #endif
3523 #endif /* BIGDECIMAL_DEBUG */
3524
3525 static Real *VpConstOne; /* constant 1.0 */
3526 static Real *VpPt5; /* constant 0.5 */
3527 #define maxnr 100UL /* Maximum iterations for calculating sqrt. */
3528 /* used in VpSqrt() */
3529
3530 /* ETC */
3531 #define MemCmp(x,y,z) memcmp(x,y,z)
3532 #define StrCmp(x,y) strcmp(x,y)
3533
3534 enum op_sw {
3535 OP_SW_ADD = 1, /* + */
3536 OP_SW_SUB, /* - */
3537 OP_SW_MULT, /* * */
3538 OP_SW_DIV /* / */
3539 };
3540
3541 static int VpIsDefOP(Real *c, Real *a, Real *b, enum op_sw sw);
3542 static int AddExponent(Real *a, SIGNED_VALUE n);
3543 static BDIGIT VpAddAbs(Real *a,Real *b,Real *c);
3544 static BDIGIT VpSubAbs(Real *a,Real *b,Real *c);
3545 static size_t VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, BDIGIT *av, BDIGIT *bv);
3546 static int VpNmlz(Real *a);
3547 static void VpFormatSt(char *psz, size_t fFmt);
3548 static int VpRdup(Real *m, size_t ind_m);
3549
3550 #ifdef BIGDECIMAL_DEBUG
3551 static int gnAlloc = 0; /* Memory allocation counter */
3552 #endif /* BIGDECIMAL_DEBUG */
3553
3554 VP_EXPORT void *
VpMemAlloc(size_t mb)3555 VpMemAlloc(size_t mb)
3556 {
3557 void *p = xmalloc(mb);
3558 if (!p) {
3559 VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
3560 }
3561 memset(p, 0, mb);
3562 #ifdef BIGDECIMAL_DEBUG
3563 gnAlloc++; /* Count allocation call */
3564 #endif /* BIGDECIMAL_DEBUG */
3565 return p;
3566 }
3567
3568 VP_EXPORT void *
VpMemRealloc(void * ptr,size_t mb)3569 VpMemRealloc(void *ptr, size_t mb)
3570 {
3571 void *p = xrealloc(ptr, mb);
3572 if (!p) {
3573 VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
3574 }
3575 return p;
3576 }
3577
3578 VP_EXPORT void
VpFree(Real * pv)3579 VpFree(Real *pv)
3580 {
3581 if (pv != NULL) {
3582 xfree(pv);
3583 #ifdef BIGDECIMAL_DEBUG
3584 gnAlloc--; /* Decrement allocation count */
3585 if (gnAlloc == 0) {
3586 printf(" *************** All memories allocated freed ****************\n");
3587 /*getchar();*/
3588 }
3589 if (gnAlloc < 0) {
3590 printf(" ??????????? Too many memory free calls(%d) ?????????????\n", gnAlloc);
3591 /*getchar();*/
3592 }
3593 #endif /* BIGDECIMAL_DEBUG */
3594 }
3595 }
3596
3597 /*
3598 * EXCEPTION Handling.
3599 */
3600
3601 #define rmpd_set_thread_local_exception_mode(mode) \
3602 rb_thread_local_aset( \
3603 rb_thread_current(), \
3604 id_BigDecimal_exception_mode, \
3605 INT2FIX((int)(mode)) \
3606 )
3607
3608 static unsigned short
VpGetException(void)3609 VpGetException (void)
3610 {
3611 VALUE const vmode = rb_thread_local_aref(
3612 rb_thread_current(),
3613 id_BigDecimal_exception_mode
3614 );
3615
3616 if (NIL_P(vmode)) {
3617 rmpd_set_thread_local_exception_mode(RMPD_EXCEPTION_MODE_DEFAULT);
3618 return RMPD_EXCEPTION_MODE_DEFAULT;
3619 }
3620
3621 return NUM2USHORT(vmode);
3622 }
3623
3624 static void
VpSetException(unsigned short f)3625 VpSetException(unsigned short f)
3626 {
3627 rmpd_set_thread_local_exception_mode(f);
3628 }
3629
3630 /*
3631 * Precision limit.
3632 */
3633
3634 #define rmpd_set_thread_local_precision_limit(limit) \
3635 rb_thread_local_aset( \
3636 rb_thread_current(), \
3637 id_BigDecimal_precision_limit, \
3638 SIZET2NUM(limit) \
3639 )
3640 #define RMPD_PRECISION_LIMIT_DEFAULT ((size_t)0)
3641
3642 /* These 2 functions added at v1.1.7 */
3643 VP_EXPORT size_t
VpGetPrecLimit(void)3644 VpGetPrecLimit(void)
3645 {
3646 VALUE const vlimit = rb_thread_local_aref(
3647 rb_thread_current(),
3648 id_BigDecimal_precision_limit
3649 );
3650
3651 if (NIL_P(vlimit)) {
3652 rmpd_set_thread_local_precision_limit(RMPD_PRECISION_LIMIT_DEFAULT);
3653 return RMPD_PRECISION_LIMIT_DEFAULT;
3654 }
3655
3656 return NUM2SIZET(vlimit);
3657 }
3658
3659 VP_EXPORT size_t
VpSetPrecLimit(size_t n)3660 VpSetPrecLimit(size_t n)
3661 {
3662 size_t const s = VpGetPrecLimit();
3663 rmpd_set_thread_local_precision_limit(n);
3664 return s;
3665 }
3666
3667 /*
3668 * Rounding mode.
3669 */
3670
3671 #define rmpd_set_thread_local_rounding_mode(mode) \
3672 rb_thread_local_aset( \
3673 rb_thread_current(), \
3674 id_BigDecimal_rounding_mode, \
3675 INT2FIX((int)(mode)) \
3676 )
3677
3678 VP_EXPORT unsigned short
VpGetRoundMode(void)3679 VpGetRoundMode(void)
3680 {
3681 VALUE const vmode = rb_thread_local_aref(
3682 rb_thread_current(),
3683 id_BigDecimal_rounding_mode
3684 );
3685
3686 if (NIL_P(vmode)) {
3687 rmpd_set_thread_local_rounding_mode(RMPD_ROUNDING_MODE_DEFAULT);
3688 return RMPD_ROUNDING_MODE_DEFAULT;
3689 }
3690
3691 return NUM2USHORT(vmode);
3692 }
3693
3694 VP_EXPORT int
VpIsRoundMode(unsigned short n)3695 VpIsRoundMode(unsigned short n)
3696 {
3697 switch (n) {
3698 case VP_ROUND_UP:
3699 case VP_ROUND_DOWN:
3700 case VP_ROUND_HALF_UP:
3701 case VP_ROUND_HALF_DOWN:
3702 case VP_ROUND_CEIL:
3703 case VP_ROUND_FLOOR:
3704 case VP_ROUND_HALF_EVEN:
3705 return 1;
3706
3707 default:
3708 return 0;
3709 }
3710 }
3711
3712 VP_EXPORT unsigned short
VpSetRoundMode(unsigned short n)3713 VpSetRoundMode(unsigned short n)
3714 {
3715 if (VpIsRoundMode(n)) {
3716 rmpd_set_thread_local_rounding_mode(n);
3717 return n;
3718 }
3719
3720 return VpGetRoundMode();
3721 }
3722
3723 /*
3724 * 0.0 & 1.0 generator
3725 * These gZero_..... and gOne_..... can be any name
3726 * referenced from nowhere except Zero() and One().
3727 * gZero_..... and gOne_..... must have global scope
3728 * (to let the compiler know they may be changed in outside
3729 * (... but not actually..)).
3730 */
3731 volatile const double gOne_ABCED9B4_CE73__00400511F31D = 1.0;
3732
3733 static double
One(void)3734 One(void)
3735 {
3736 return gOne_ABCED9B4_CE73__00400511F31D;
3737 }
3738
3739 /*
3740 ----------------------------------------------------------------
3741 Value of sign in Real structure is reserved for future use.
3742 short sign;
3743 ==0 : NaN
3744 1 : Positive zero
3745 -1 : Negative zero
3746 2 : Positive number
3747 -2 : Negative number
3748 3 : Positive infinite number
3749 -3 : Negative infinite number
3750 ----------------------------------------------------------------
3751 */
3752
3753 VP_EXPORT double
VpGetDoubleNaN(void)3754 VpGetDoubleNaN(void) /* Returns the value of NaN */
3755 {
3756 return nan("");
3757 }
3758
3759 VP_EXPORT double
VpGetDoublePosInf(void)3760 VpGetDoublePosInf(void) /* Returns the value of +Infinity */
3761 {
3762 return HUGE_VAL;
3763 }
3764
3765 VP_EXPORT double
VpGetDoubleNegInf(void)3766 VpGetDoubleNegInf(void) /* Returns the value of -Infinity */
3767 {
3768 return -HUGE_VAL;
3769 }
3770
3771 VP_EXPORT double
VpGetDoubleNegZero(void)3772 VpGetDoubleNegZero(void) /* Returns the value of -0 */
3773 {
3774 static double nzero = 1000.0;
3775 if (nzero != 0.0) nzero = (One()/VpGetDoubleNegInf());
3776 return nzero;
3777 }
3778
3779 #if 0 /* unused */
3780 VP_EXPORT int
3781 VpIsNegDoubleZero(double v)
3782 {
3783 double z = VpGetDoubleNegZero();
3784 return MemCmp(&v,&z,sizeof(v))==0;
3785 }
3786 #endif
3787
3788 VP_EXPORT int
VpException(unsigned short f,const char * str,int always)3789 VpException(unsigned short f, const char *str,int always)
3790 {
3791 unsigned short const exception_mode = VpGetException();
3792
3793 if (f == VP_EXCEPTION_OP || f == VP_EXCEPTION_MEMORY) always = 1;
3794
3795 if (always || (exception_mode & f)) {
3796 switch(f) {
3797 /* case VP_EXCEPTION_OVERFLOW: */
3798 case VP_EXCEPTION_ZERODIVIDE:
3799 case VP_EXCEPTION_INFINITY:
3800 case VP_EXCEPTION_NaN:
3801 case VP_EXCEPTION_UNDERFLOW:
3802 case VP_EXCEPTION_OP:
3803 rb_raise(rb_eFloatDomainError, "%s", str);
3804 break;
3805 case VP_EXCEPTION_MEMORY:
3806 default:
3807 rb_fatal("%s", str);
3808 }
3809 }
3810 return 0; /* 0 Means VpException() raised no exception */
3811 }
3812
3813 /* Throw exception or returns 0,when resulting c is Inf or NaN */
3814 /* sw=1:+ 2:- 3:* 4:/ */
3815 static int
VpIsDefOP(Real * c,Real * a,Real * b,enum op_sw sw)3816 VpIsDefOP(Real *c, Real *a, Real *b, enum op_sw sw)
3817 {
3818 if (VpIsNaN(a) || VpIsNaN(b)) {
3819 /* at least a or b is NaN */
3820 VpSetNaN(c);
3821 goto NaN;
3822 }
3823
3824 if (VpIsInf(a)) {
3825 if (VpIsInf(b)) {
3826 switch(sw) {
3827 case OP_SW_ADD: /* + */
3828 if (VpGetSign(a) == VpGetSign(b)) {
3829 VpSetInf(c, VpGetSign(a));
3830 goto Inf;
3831 }
3832 else {
3833 VpSetNaN(c);
3834 goto NaN;
3835 }
3836 case OP_SW_SUB: /* - */
3837 if (VpGetSign(a) != VpGetSign(b)) {
3838 VpSetInf(c, VpGetSign(a));
3839 goto Inf;
3840 }
3841 else {
3842 VpSetNaN(c);
3843 goto NaN;
3844 }
3845 case OP_SW_MULT: /* * */
3846 VpSetInf(c, VpGetSign(a)*VpGetSign(b));
3847 goto Inf;
3848 case OP_SW_DIV: /* / */
3849 VpSetNaN(c);
3850 goto NaN;
3851 }
3852 VpSetNaN(c);
3853 goto NaN;
3854 }
3855 /* Inf op Finite */
3856 switch(sw) {
3857 case OP_SW_ADD: /* + */
3858 case OP_SW_SUB: /* - */
3859 VpSetInf(c, VpGetSign(a));
3860 break;
3861 case OP_SW_MULT: /* * */
3862 if (VpIsZero(b)) {
3863 VpSetNaN(c);
3864 goto NaN;
3865 }
3866 VpSetInf(c, VpGetSign(a)*VpGetSign(b));
3867 break;
3868 case OP_SW_DIV: /* / */
3869 VpSetInf(c, VpGetSign(a)*VpGetSign(b));
3870 }
3871 goto Inf;
3872 }
3873
3874 if (VpIsInf(b)) {
3875 switch(sw) {
3876 case OP_SW_ADD: /* + */
3877 VpSetInf(c, VpGetSign(b));
3878 break;
3879 case OP_SW_SUB: /* - */
3880 VpSetInf(c, -VpGetSign(b));
3881 break;
3882 case OP_SW_MULT: /* * */
3883 if (VpIsZero(a)) {
3884 VpSetNaN(c);
3885 goto NaN;
3886 }
3887 VpSetInf(c, VpGetSign(a)*VpGetSign(b));
3888 break;
3889 case OP_SW_DIV: /* / */
3890 VpSetZero(c, VpGetSign(a)*VpGetSign(b));
3891 }
3892 goto Inf;
3893 }
3894 return 1; /* Results OK */
3895
3896 Inf:
3897 if (VpIsPosInf(c)) {
3898 return VpException(VP_EXCEPTION_INFINITY, "Computation results to 'Infinity'", 0);
3899 }
3900 else {
3901 return VpException(VP_EXCEPTION_INFINITY, "Computation results to '-Infinity'", 0);
3902 }
3903
3904 NaN:
3905 return VpException(VP_EXCEPTION_NaN, "Computation results to 'NaN'", 0);
3906 }
3907
3908 /*
3909 ----------------------------------------------------------------
3910 */
3911
3912 /*
3913 * returns number of chars needed to represent vp in specified format.
3914 */
3915 VP_EXPORT size_t
VpNumOfChars(Real * vp,const char * pszFmt)3916 VpNumOfChars(Real *vp,const char *pszFmt)
3917 {
3918 SIGNED_VALUE ex;
3919 size_t nc;
3920
3921 if (vp == NULL) return BASE_FIG*2+6;
3922 if (!VpIsDef(vp)) return 32; /* not sure,may be OK */
3923
3924 switch(*pszFmt) {
3925 case 'F':
3926 nc = BASE_FIG*(vp->Prec + 1)+2;
3927 ex = vp->exponent;
3928 if (ex < 0) {
3929 nc += BASE_FIG*(size_t)(-ex);
3930 }
3931 else {
3932 if ((size_t)ex > vp->Prec) {
3933 nc += BASE_FIG*((size_t)ex - vp->Prec);
3934 }
3935 }
3936 break;
3937 case 'E':
3938 /* fall through */
3939 default:
3940 nc = BASE_FIG*(vp->Prec + 2)+6; /* 3: sign + exponent chars */
3941 }
3942 return nc;
3943 }
3944
3945 /*
3946 * Initializer for Vp routines and constants used.
3947 * [Input]
3948 * BaseVal: Base value(assigned to BASE) for Vp calculation.
3949 * It must be the form BaseVal=10**n.(n=1,2,3,...)
3950 * If Base <= 0L,then the BASE will be calculated so
3951 * that BASE is as large as possible satisfying the
3952 * relation MaxVal <= BASE*(BASE+1). Where the value
3953 * MaxVal is the largest value which can be represented
3954 * by one BDIGIT word in the computer used.
3955 *
3956 * [Returns]
3957 * 1+DBL_DIG ... OK
3958 */
3959 VP_EXPORT size_t
VpInit(BDIGIT BaseVal)3960 VpInit(BDIGIT BaseVal)
3961 {
3962 /* Setup +/- Inf NaN -0 */
3963 VpGetDoubleNegZero();
3964
3965 /* Allocates Vp constants. */
3966 VpConstOne = VpAlloc(1UL, "1", 1, 1);
3967 VpPt5 = VpAlloc(1UL, ".5", 1, 1);
3968
3969 #ifdef BIGDECIMAL_DEBUG
3970 gnAlloc = 0;
3971 #endif /* BIGDECIMAL_DEBUG */
3972
3973 #ifdef BIGDECIMAL_DEBUG
3974 if (gfDebug) {
3975 printf("VpInit: BaseVal = %"PRIuBDIGIT"\n", BaseVal);
3976 printf("\tBASE = %"PRIuBDIGIT"\n", BASE);
3977 printf("\tHALF_BASE = %"PRIuBDIGIT"\n", HALF_BASE);
3978 printf("\tBASE1 = %"PRIuBDIGIT"\n", BASE1);
3979 printf("\tBASE_FIG = %u\n", BASE_FIG);
3980 printf("\tDBLE_FIG = %d\n", DBLE_FIG);
3981 }
3982 #endif /* BIGDECIMAL_DEBUG */
3983
3984 return rmpd_double_figures();
3985 }
3986
3987 VP_EXPORT Real *
VpOne(void)3988 VpOne(void)
3989 {
3990 return VpConstOne;
3991 }
3992
3993 /* If exponent overflows,then raise exception or returns 0 */
3994 static int
AddExponent(Real * a,SIGNED_VALUE n)3995 AddExponent(Real *a, SIGNED_VALUE n)
3996 {
3997 SIGNED_VALUE e = a->exponent;
3998 SIGNED_VALUE m = e+n;
3999 SIGNED_VALUE eb, mb;
4000 if (e > 0) {
4001 if (n > 0) {
4002 if (MUL_OVERFLOW_SIGNED_VALUE_P(m, (SIGNED_VALUE)BASE_FIG) ||
4003 MUL_OVERFLOW_SIGNED_VALUE_P(e, (SIGNED_VALUE)BASE_FIG))
4004 goto overflow;
4005 mb = m*(SIGNED_VALUE)BASE_FIG;
4006 eb = e*(SIGNED_VALUE)BASE_FIG;
4007 if (eb - mb > 0) goto overflow;
4008 }
4009 }
4010 else if (n < 0) {
4011 if (MUL_OVERFLOW_SIGNED_VALUE_P(m, (SIGNED_VALUE)BASE_FIG) ||
4012 MUL_OVERFLOW_SIGNED_VALUE_P(e, (SIGNED_VALUE)BASE_FIG))
4013 goto underflow;
4014 mb = m*(SIGNED_VALUE)BASE_FIG;
4015 eb = e*(SIGNED_VALUE)BASE_FIG;
4016 if (mb - eb > 0) goto underflow;
4017 }
4018 a->exponent = m;
4019 return 1;
4020
4021 /* Overflow/Underflow ==> Raise exception or returns 0 */
4022 underflow:
4023 VpSetZero(a, VpGetSign(a));
4024 return VpException(VP_EXCEPTION_UNDERFLOW, "Exponent underflow", 0);
4025
4026 overflow:
4027 VpSetInf(a, VpGetSign(a));
4028 return VpException(VP_EXCEPTION_OVERFLOW, "Exponent overflow", 0);
4029 }
4030
4031 Real *
rmpd_parse_special_string(const char * str)4032 rmpd_parse_special_string(const char *str)
4033 {
4034 static const struct {
4035 const char *str;
4036 size_t len;
4037 int sign;
4038 } table[] = {
4039 { SZ_INF, sizeof(SZ_INF) - 1, VP_SIGN_POSITIVE_INFINITE },
4040 { SZ_PINF, sizeof(SZ_PINF) - 1, VP_SIGN_POSITIVE_INFINITE },
4041 { SZ_NINF, sizeof(SZ_NINF) - 1, VP_SIGN_NEGATIVE_INFINITE },
4042 { SZ_NaN, sizeof(SZ_NaN) - 1, VP_SIGN_NaN }
4043 };
4044 static const size_t table_length = sizeof(table) / sizeof(table[0]);
4045 size_t i;
4046
4047 for (i = 0; i < table_length; ++i) {
4048 const char *p;
4049 if (strncmp(str, table[i].str, table[i].len) != 0) {
4050 continue;
4051 }
4052
4053 p = str + table[i].len;
4054 while (*p && ISSPACE(*p)) ++p;
4055 if (*p == '\0') {
4056 Real *vp = VpAllocReal(1);
4057 vp->MaxPrec = 1;
4058 switch (table[i].sign) {
4059 default:
4060 UNREACHABLE; break;
4061 case VP_SIGN_POSITIVE_INFINITE:
4062 VpSetPosInf(vp);
4063 return vp;
4064 case VP_SIGN_NEGATIVE_INFINITE:
4065 VpSetNegInf(vp);
4066 return vp;
4067 case VP_SIGN_NaN:
4068 VpSetNaN(vp);
4069 return vp;
4070 }
4071 }
4072 }
4073
4074 return NULL;
4075 }
4076
4077 /*
4078 * Allocates variable.
4079 * [Input]
4080 * mx ... allocation unit, if zero then mx is determined by szVal.
4081 * The mx is the number of effective digits can to be stored.
4082 * szVal ... value assigned(char). If szVal==NULL,then zero is assumed.
4083 * If szVal[0]=='#' then Max. Prec. will not be considered(1.1.7),
4084 * full precision specified by szVal is allocated.
4085 *
4086 * [Returns]
4087 * Pointer to the newly allocated variable, or
4088 * NULL be returned if memory allocation is failed,or any error.
4089 */
4090 VP_EXPORT Real *
VpAlloc(size_t mx,const char * szVal,int strict_p,int exc)4091 VpAlloc(size_t mx, const char *szVal, int strict_p, int exc)
4092 {
4093 const char *orig_szVal = szVal;
4094 size_t i, j, ni, ipf, nf, ipe, ne, dot_seen, exp_seen, nalloc;
4095 char v, *psz;
4096 int sign=1;
4097 Real *vp = NULL;
4098 size_t mf = VpGetPrecLimit();
4099 VALUE buf;
4100
4101 mx = (mx + BASE_FIG - 1) / BASE_FIG; /* Determine allocation unit. */
4102 if (mx == 0) ++mx;
4103
4104 if (szVal) {
4105 /* Skipping leading spaces */
4106 while (ISSPACE(*szVal)) szVal++;
4107
4108 /* Processing the leading one `#` */
4109 if (*szVal != '#') {
4110 if (mf) {
4111 mf = (mf + BASE_FIG - 1) / BASE_FIG + 2; /* Needs 1 more for div */
4112 if (mx > mf) {
4113 mx = mf;
4114 }
4115 }
4116 }
4117 else {
4118 ++szVal;
4119 }
4120 }
4121 else {
4122 return_zero:
4123 /* necessary to be able to store */
4124 /* at least mx digits. */
4125 /* szVal==NULL ==> allocate zero value. */
4126 vp = VpAllocReal(mx);
4127 /* xmalloc() alway returns(or throw interruption) */
4128 vp->MaxPrec = mx; /* set max precision */
4129 VpSetZero(vp, 1); /* initialize vp to zero. */
4130 return vp;
4131 }
4132
4133 /* Check on Inf & NaN */
4134 if ((vp = rmpd_parse_special_string(szVal)) != NULL) {
4135 return vp;
4136 }
4137
4138 /* Scanning digits */
4139
4140 /* A buffer for keeping scanned digits */
4141 buf = rb_str_tmp_new(strlen(szVal) + 1);
4142 psz = RSTRING_PTR(buf);
4143
4144 /* cursor: i for psz, and j for szVal */
4145 i = j = 0;
4146
4147 /* Scanning: sign part */
4148 v = psz[i] = szVal[j];
4149 if ((v == '-') || (v == '+')) {
4150 sign = -(v == '-');
4151 ++i;
4152 ++j;
4153 }
4154
4155 /* Scanning: integer part */
4156 ni = 0; /* number of digits in the integer part */
4157 while ((v = psz[i] = szVal[j]) != '\0') {
4158 if (!strict_p && ISSPACE(v)) {
4159 v = psz[i] = '\0';
4160 break;
4161 }
4162 if (v == '_') {
4163 if (ni > 0) {
4164 v = szVal[j+1];
4165 if (v == '\0' || ISSPACE(v) || ISDIGIT(v)) {
4166 ++j;
4167 continue;
4168 }
4169 if (!strict_p) {
4170 v = psz[i] = '\0';
4171 break;
4172 }
4173 }
4174 goto invalid_value;
4175 }
4176 if (!ISDIGIT(v)) {
4177 break;
4178 }
4179 ++ni;
4180 ++i;
4181 ++j;
4182 }
4183
4184 /* Scanning: fractional part */
4185 nf = 0; /* number of digits in the fractional part */
4186 ne = 0; /* number of digits in the exponential part */
4187 ipf = 0; /* index of the beginning of the fractional part */
4188 ipe = 0; /* index of the beginning of the exponential part */
4189 dot_seen = 0;
4190 exp_seen = 0;
4191
4192 if (v != '\0') {
4193 /* Scanning fractional part */
4194 if ((psz[i] = szVal[j]) == '.') {
4195 dot_seen = 1;
4196 ++i;
4197 ++j;
4198 ipf = i;
4199 while ((v = psz[i] = szVal[j]) != '\0') {
4200 if (!strict_p && ISSPACE(v)) {
4201 v = psz[i] = '\0';
4202 break;
4203 }
4204 if (v == '_') {
4205 if (nf > 0 && ISDIGIT(szVal[j+1])) {
4206 ++j;
4207 continue;
4208 }
4209 if (!strict_p) {
4210 v = psz[i] = '\0';
4211 if (nf == 0) {
4212 dot_seen = 0;
4213 }
4214 break;
4215 }
4216 goto invalid_value;
4217 }
4218 if (!ISDIGIT(v)) break;
4219 ++i;
4220 ++j;
4221 ++nf;
4222 }
4223 }
4224
4225 /* Scanning exponential part */
4226 if (v != '\0') {
4227 switch ((psz[i] = szVal[j])) {
4228 case '\0':
4229 break;
4230 case 'e': case 'E':
4231 case 'd': case 'D':
4232 exp_seen = 1;
4233 ++i;
4234 ++j;
4235 ipe = i;
4236 v = psz[i] = szVal[j];
4237 if ((v == '-') || (v == '+')) {
4238 ++i;
4239 ++j;
4240 }
4241 while ((v = psz[i] = szVal[j]) != '\0') {
4242 if (!strict_p && ISSPACE(v)) {
4243 v = psz[i] = '\0';
4244 break;
4245 }
4246 if (v == '_') {
4247 if (ne > 0 && ISDIGIT(szVal[j+1])) {
4248 ++j;
4249 continue;
4250 }
4251 if (!strict_p) {
4252 v = psz[i] = '\0';
4253 if (ne == 0) {
4254 exp_seen = 0;
4255 }
4256 break;
4257 }
4258 goto invalid_value;
4259 }
4260 if (!ISDIGIT(v)) break;
4261 ++i;
4262 ++j;
4263 ++ne;
4264 }
4265 break;
4266 default:
4267 break;
4268 }
4269 }
4270
4271 if (v != '\0') {
4272 /* Scanning trailing spaces */
4273 while (ISSPACE(szVal[j])) ++j;
4274
4275 /* Invalid character */
4276 if (szVal[j] && strict_p) {
4277 goto invalid_value;
4278 }
4279 }
4280 }
4281
4282 psz[i] = '\0';
4283
4284 if (((ni == 0 || dot_seen) && nf == 0) || (exp_seen && ne == 0)) {
4285 VALUE str;
4286 invalid_value:
4287 if (!strict_p) {
4288 goto return_zero;
4289 }
4290 if (!exc) {
4291 return NULL;
4292 }
4293 str = rb_str_new2(orig_szVal);
4294 rb_raise(rb_eArgError, "invalid value for BigDecimal(): \"%"PRIsVALUE"\"", str);
4295 }
4296
4297 nalloc = (ni + nf + BASE_FIG - 1) / BASE_FIG + 1; /* set effective allocation */
4298 /* units for szVal[] */
4299 if (mx == 0) mx = 1;
4300 nalloc = Max(nalloc, mx);
4301 mx = nalloc;
4302 vp = VpAllocReal(mx);
4303 /* xmalloc() alway returns(or throw interruption) */
4304 vp->MaxPrec = mx; /* set max precision */
4305 VpSetZero(vp, sign);
4306 VpCtoV(vp, psz, ni, psz + ipf, nf, psz + ipe, ne);
4307 rb_str_resize(buf, 0);
4308 return vp;
4309 }
4310
4311 /*
4312 * Assignment(c=a).
4313 * [Input]
4314 * a ... RHSV
4315 * isw ... switch for assignment.
4316 * c = a when isw > 0
4317 * c = -a when isw < 0
4318 * if c->MaxPrec < a->Prec,then round operation
4319 * will be performed.
4320 * [Output]
4321 * c ... LHSV
4322 */
4323 VP_EXPORT size_t
VpAsgn(Real * c,Real * a,int isw)4324 VpAsgn(Real *c, Real *a, int isw)
4325 {
4326 size_t n;
4327 if (VpIsNaN(a)) {
4328 VpSetNaN(c);
4329 return 0;
4330 }
4331 if (VpIsInf(a)) {
4332 VpSetInf(c, isw * VpGetSign(a));
4333 return 0;
4334 }
4335
4336 /* check if the RHS is zero */
4337 if (!VpIsZero(a)) {
4338 c->exponent = a->exponent; /* store exponent */
4339 VpSetSign(c, isw * VpGetSign(a)); /* set sign */
4340 n = (a->Prec < c->MaxPrec) ? (a->Prec) : (c->MaxPrec);
4341 c->Prec = n;
4342 memcpy(c->frac, a->frac, n * sizeof(BDIGIT));
4343 /* Needs round ? */
4344 if (isw != 10) {
4345 /* Not in ActiveRound */
4346 if(c->Prec < a->Prec) {
4347 VpInternalRound(c, n, (n>0) ? a->frac[n-1] : 0, a->frac[n]);
4348 }
4349 else {
4350 VpLimitRound(c,0);
4351 }
4352 }
4353 }
4354 else {
4355 /* The value of 'a' is zero. */
4356 VpSetZero(c, isw * VpGetSign(a));
4357 return 1;
4358 }
4359 return c->Prec * BASE_FIG;
4360 }
4361
4362 /*
4363 * c = a + b when operation = 1 or 2
4364 * c = a - b when operation = -1 or -2.
4365 * Returns number of significant digits of c
4366 */
4367 VP_EXPORT size_t
VpAddSub(Real * c,Real * a,Real * b,int operation)4368 VpAddSub(Real *c, Real *a, Real *b, int operation)
4369 {
4370 short sw, isw;
4371 Real *a_ptr, *b_ptr;
4372 size_t n, na, nb, i;
4373 BDIGIT mrv;
4374
4375 #ifdef BIGDECIMAL_DEBUG
4376 if (gfDebug) {
4377 VPrint(stdout, "VpAddSub(enter) a=% \n", a);
4378 VPrint(stdout, " b=% \n", b);
4379 printf(" operation=%d\n", operation);
4380 }
4381 #endif /* BIGDECIMAL_DEBUG */
4382
4383 if (!VpIsDefOP(c, a, b, (operation > 0) ? OP_SW_ADD : OP_SW_SUB)) return 0; /* No significant digits */
4384
4385 /* check if a or b is zero */
4386 if (VpIsZero(a)) {
4387 /* a is zero,then assign b to c */
4388 if (!VpIsZero(b)) {
4389 VpAsgn(c, b, operation);
4390 }
4391 else {
4392 /* Both a and b are zero. */
4393 if (VpGetSign(a) < 0 && operation * VpGetSign(b) < 0) {
4394 /* -0 -0 */
4395 VpSetZero(c, -1);
4396 }
4397 else {
4398 VpSetZero(c, 1);
4399 }
4400 return 1; /* 0: 1 significant digits */
4401 }
4402 return c->Prec * BASE_FIG;
4403 }
4404 if (VpIsZero(b)) {
4405 /* b is zero,then assign a to c. */
4406 VpAsgn(c, a, 1);
4407 return c->Prec*BASE_FIG;
4408 }
4409
4410 if (operation < 0) sw = -1;
4411 else sw = 1;
4412
4413 /* compare absolute value. As a result,|a_ptr|>=|b_ptr| */
4414 if (a->exponent > b->exponent) {
4415 a_ptr = a;
4416 b_ptr = b;
4417 } /* |a|>|b| */
4418 else if (a->exponent < b->exponent) {
4419 a_ptr = b;
4420 b_ptr = a;
4421 } /* |a|<|b| */
4422 else {
4423 /* Exponent part of a and b is the same,then compare fraction */
4424 /* part */
4425 na = a->Prec;
4426 nb = b->Prec;
4427 n = Min(na, nb);
4428 for (i=0; i < n; ++i) {
4429 if (a->frac[i] > b->frac[i]) {
4430 a_ptr = a;
4431 b_ptr = b;
4432 goto end_if;
4433 }
4434 else if (a->frac[i] < b->frac[i]) {
4435 a_ptr = b;
4436 b_ptr = a;
4437 goto end_if;
4438 }
4439 }
4440 if (na > nb) {
4441 a_ptr = a;
4442 b_ptr = b;
4443 goto end_if;
4444 }
4445 else if (na < nb) {
4446 a_ptr = b;
4447 b_ptr = a;
4448 goto end_if;
4449 }
4450 /* |a| == |b| */
4451 if (VpGetSign(a) + sw *VpGetSign(b) == 0) {
4452 VpSetZero(c, 1); /* abs(a)=abs(b) and operation = '-' */
4453 return c->Prec * BASE_FIG;
4454 }
4455 a_ptr = a;
4456 b_ptr = b;
4457 }
4458
4459 end_if:
4460 isw = VpGetSign(a) + sw *VpGetSign(b);
4461 /*
4462 * isw = 0 ...( 1)+(-1),( 1)-( 1),(-1)+(1),(-1)-(-1)
4463 * = 2 ...( 1)+( 1),( 1)-(-1)
4464 * =-2 ...(-1)+(-1),(-1)-( 1)
4465 * If isw==0, then c =(Sign a_ptr)(|a_ptr|-|b_ptr|)
4466 * else c =(Sign ofisw)(|a_ptr|+|b_ptr|)
4467 */
4468 if (isw) { /* addition */
4469 VpSetSign(c, 1);
4470 mrv = VpAddAbs(a_ptr, b_ptr, c);
4471 VpSetSign(c, isw / 2);
4472 }
4473 else { /* subtraction */
4474 VpSetSign(c, 1);
4475 mrv = VpSubAbs(a_ptr, b_ptr, c);
4476 if (a_ptr == a) {
4477 VpSetSign(c,VpGetSign(a));
4478 }
4479 else {
4480 VpSetSign(c, VpGetSign(a_ptr) * sw);
4481 }
4482 }
4483 VpInternalRound(c, 0, (c->Prec > 0) ? c->frac[c->Prec-1] : 0, mrv);
4484
4485 #ifdef BIGDECIMAL_DEBUG
4486 if (gfDebug) {
4487 VPrint(stdout, "VpAddSub(result) c=% \n", c);
4488 VPrint(stdout, " a=% \n", a);
4489 VPrint(stdout, " b=% \n", b);
4490 printf(" operation=%d\n", operation);
4491 }
4492 #endif /* BIGDECIMAL_DEBUG */
4493 return c->Prec * BASE_FIG;
4494 }
4495
4496 /*
4497 * Addition of two values with variable precision
4498 * a and b assuming abs(a)>abs(b).
4499 * c = abs(a) + abs(b) ; where |a|>=|b|
4500 */
4501 static BDIGIT
VpAddAbs(Real * a,Real * b,Real * c)4502 VpAddAbs(Real *a, Real *b, Real *c)
4503 {
4504 size_t word_shift;
4505 size_t ap;
4506 size_t bp;
4507 size_t cp;
4508 size_t a_pos;
4509 size_t b_pos, b_pos_with_word_shift;
4510 size_t c_pos;
4511 BDIGIT av, bv, carry, mrv;
4512
4513 #ifdef BIGDECIMAL_DEBUG
4514 if (gfDebug) {
4515 VPrint(stdout, "VpAddAbs called: a = %\n", a);
4516 VPrint(stdout, " b = %\n", b);
4517 }
4518 #endif /* BIGDECIMAL_DEBUG */
4519
4520 word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
4521 a_pos = ap;
4522 b_pos = bp;
4523 c_pos = cp;
4524
4525 if (word_shift == (size_t)-1L) return 0; /* Overflow */
4526 if (b_pos == (size_t)-1L) goto Assign_a;
4527
4528 mrv = av + bv; /* Most right val. Used for round. */
4529
4530 /* Just assign the last few digits of b to c because a has no */
4531 /* corresponding digits to be added. */
4532 if (b_pos > 0) {
4533 while (b_pos > 0 && b_pos + word_shift > a_pos) {
4534 c->frac[--c_pos] = b->frac[--b_pos];
4535 }
4536 }
4537 if (b_pos == 0 && word_shift > a_pos) {
4538 while (word_shift-- > a_pos) {
4539 c->frac[--c_pos] = 0;
4540 }
4541 }
4542
4543 /* Just assign the last few digits of a to c because b has no */
4544 /* corresponding digits to be added. */
4545 b_pos_with_word_shift = b_pos + word_shift;
4546 while (a_pos > b_pos_with_word_shift) {
4547 c->frac[--c_pos] = a->frac[--a_pos];
4548 }
4549 carry = 0; /* set first carry be zero */
4550
4551 /* Now perform addition until every digits of b will be */
4552 /* exhausted. */
4553 while (b_pos > 0) {
4554 c->frac[--c_pos] = a->frac[--a_pos] + b->frac[--b_pos] + carry;
4555 if (c->frac[c_pos] >= BASE) {
4556 c->frac[c_pos] -= BASE;
4557 carry = 1;
4558 }
4559 else {
4560 carry = 0;
4561 }
4562 }
4563
4564 /* Just assign the first few digits of a with considering */
4565 /* the carry obtained so far because b has been exhausted. */
4566 while (a_pos > 0) {
4567 c->frac[--c_pos] = a->frac[--a_pos] + carry;
4568 if (c->frac[c_pos] >= BASE) {
4569 c->frac[c_pos] -= BASE;
4570 carry = 1;
4571 }
4572 else {
4573 carry = 0;
4574 }
4575 }
4576 if (c_pos) c->frac[c_pos - 1] += carry;
4577 goto Exit;
4578
4579 Assign_a:
4580 VpAsgn(c, a, 1);
4581 mrv = 0;
4582
4583 Exit:
4584
4585 #ifdef BIGDECIMAL_DEBUG
4586 if (gfDebug) {
4587 VPrint(stdout, "VpAddAbs exit: c=% \n", c);
4588 }
4589 #endif /* BIGDECIMAL_DEBUG */
4590 return mrv;
4591 }
4592
4593 /*
4594 * c = abs(a) - abs(b)
4595 */
4596 static BDIGIT
VpSubAbs(Real * a,Real * b,Real * c)4597 VpSubAbs(Real *a, Real *b, Real *c)
4598 {
4599 size_t word_shift;
4600 size_t ap;
4601 size_t bp;
4602 size_t cp;
4603 size_t a_pos;
4604 size_t b_pos, b_pos_with_word_shift;
4605 size_t c_pos;
4606 BDIGIT av, bv, borrow, mrv;
4607
4608 #ifdef BIGDECIMAL_DEBUG
4609 if (gfDebug) {
4610 VPrint(stdout, "VpSubAbs called: a = %\n", a);
4611 VPrint(stdout, " b = %\n", b);
4612 }
4613 #endif /* BIGDECIMAL_DEBUG */
4614
4615 word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
4616 a_pos = ap;
4617 b_pos = bp;
4618 c_pos = cp;
4619 if (word_shift == (size_t)-1L) return 0; /* Overflow */
4620 if (b_pos == (size_t)-1L) goto Assign_a;
4621
4622 if (av >= bv) {
4623 mrv = av - bv;
4624 borrow = 0;
4625 }
4626 else {
4627 mrv = 0;
4628 borrow = 1;
4629 }
4630
4631 /* Just assign the values which are the BASE subtracted by */
4632 /* each of the last few digits of the b because the a has no */
4633 /* corresponding digits to be subtracted. */
4634 if (b_pos + word_shift > a_pos) {
4635 while (b_pos > 0 && b_pos + word_shift > a_pos) {
4636 c->frac[--c_pos] = BASE - b->frac[--b_pos] - borrow;
4637 borrow = 1;
4638 }
4639 if (b_pos == 0) {
4640 while (word_shift > a_pos) {
4641 --word_shift;
4642 c->frac[--c_pos] = BASE - borrow;
4643 borrow = 1;
4644 }
4645 }
4646 }
4647 /* Just assign the last few digits of a to c because b has no */
4648 /* corresponding digits to subtract. */
4649
4650 b_pos_with_word_shift = b_pos + word_shift;
4651 while (a_pos > b_pos_with_word_shift) {
4652 c->frac[--c_pos] = a->frac[--a_pos];
4653 }
4654
4655 /* Now perform subtraction until every digits of b will be */
4656 /* exhausted. */
4657 while (b_pos > 0) {
4658 --c_pos;
4659 if (a->frac[--a_pos] < b->frac[--b_pos] + borrow) {
4660 c->frac[c_pos] = BASE + a->frac[a_pos] - b->frac[b_pos] - borrow;
4661 borrow = 1;
4662 }
4663 else {
4664 c->frac[c_pos] = a->frac[a_pos] - b->frac[b_pos] - borrow;
4665 borrow = 0;
4666 }
4667 }
4668
4669 /* Just assign the first few digits of a with considering */
4670 /* the borrow obtained so far because b has been exhausted. */
4671 while (a_pos > 0) {
4672 --c_pos;
4673 if (a->frac[--a_pos] < borrow) {
4674 c->frac[c_pos] = BASE + a->frac[a_pos] - borrow;
4675 borrow = 1;
4676 }
4677 else {
4678 c->frac[c_pos] = a->frac[a_pos] - borrow;
4679 borrow = 0;
4680 }
4681 }
4682 if (c_pos) c->frac[c_pos - 1] -= borrow;
4683 goto Exit;
4684
4685 Assign_a:
4686 VpAsgn(c, a, 1);
4687 mrv = 0;
4688
4689 Exit:
4690 #ifdef BIGDECIMAL_DEBUG
4691 if (gfDebug) {
4692 VPrint(stdout, "VpSubAbs exit: c=% \n", c);
4693 }
4694 #endif /* BIGDECIMAL_DEBUG */
4695 return mrv;
4696 }
4697
4698 /*
4699 * Note: If(av+bv)>= HALF_BASE,then 1 will be added to the least significant
4700 * digit of c(In case of addition).
4701 * ------------------------- figure of output -----------------------------------
4702 * a = xxxxxxxxxxx
4703 * b = xxxxxxxxxx
4704 * c =xxxxxxxxxxxxxxx
4705 * word_shift = | |
4706 * right_word = | | (Total digits in RHSV)
4707 * left_word = | | (Total digits in LHSV)
4708 * a_pos = |
4709 * b_pos = |
4710 * c_pos = |
4711 */
4712 static size_t
VpSetPTR(Real * a,Real * b,Real * c,size_t * a_pos,size_t * b_pos,size_t * c_pos,BDIGIT * av,BDIGIT * bv)4713 VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, BDIGIT *av, BDIGIT *bv)
4714 {
4715 size_t left_word, right_word, word_shift;
4716
4717 size_t const round_limit = (VpGetPrecLimit() + BASE_FIG - 1) / BASE_FIG;
4718
4719 assert(a->exponent >= b->exponent);
4720
4721 c->frac[0] = 0;
4722 *av = *bv = 0;
4723
4724 word_shift = (a->exponent - b->exponent);
4725 left_word = b->Prec + word_shift;
4726 right_word = Max(a->Prec, left_word);
4727 left_word = c->MaxPrec - 1; /* -1 ... prepare for round up */
4728
4729 /*
4730 * check if 'round' is needed.
4731 */
4732 if (right_word > left_word) { /* round ? */
4733 /*---------------------------------
4734 * Actual size of a = xxxxxxAxx
4735 * Actual size of b = xxxBxxxxx
4736 * Max. size of c = xxxxxx
4737 * Round off = |-----|
4738 * c_pos = |
4739 * right_word = |
4740 * a_pos = |
4741 */
4742 *c_pos = right_word = left_word + 1; /* Set resulting precision */
4743 /* be equal to that of c */
4744 if (a->Prec >= c->MaxPrec) {
4745 /*
4746 * a = xxxxxxAxxx
4747 * c = xxxxxx
4748 * a_pos = |
4749 */
4750 *a_pos = left_word;
4751 if (*a_pos <= round_limit) {
4752 *av = a->frac[*a_pos]; /* av is 'A' shown in above. */
4753 }
4754 }
4755 else {
4756 /*
4757 * a = xxxxxxx
4758 * c = xxxxxxxxxx
4759 * a_pos = |
4760 */
4761 *a_pos = a->Prec;
4762 }
4763 if (b->Prec + word_shift >= c->MaxPrec) {
4764 /*
4765 * a = xxxxxxxxx
4766 * b = xxxxxxxBxxx
4767 * c = xxxxxxxxxxx
4768 * b_pos = |
4769 */
4770 if (c->MaxPrec >= word_shift + 1) {
4771 *b_pos = c->MaxPrec - word_shift - 1;
4772 if (*b_pos + word_shift <= round_limit) {
4773 *bv = b->frac[*b_pos];
4774 }
4775 }
4776 else {
4777 *b_pos = -1L;
4778 }
4779 }
4780 else {
4781 /*
4782 * a = xxxxxxxxxxxxxxxx
4783 * b = xxxxxx
4784 * c = xxxxxxxxxxxxx
4785 * b_pos = |
4786 */
4787 *b_pos = b->Prec;
4788 }
4789 }
4790 else { /* The MaxPrec of c - 1 > The Prec of a + b */
4791 /*
4792 * a = xxxxxxx
4793 * b = xxxxxx
4794 * c = xxxxxxxxxxx
4795 * c_pos = |
4796 */
4797 *b_pos = b->Prec;
4798 *a_pos = a->Prec;
4799 *c_pos = right_word + 1;
4800 }
4801 c->Prec = *c_pos;
4802 c->exponent = a->exponent;
4803 if (!AddExponent(c, 1)) return (size_t)-1L;
4804 return word_shift;
4805 }
4806
4807 /*
4808 * Return number of significant digits
4809 * c = a * b , Where a = a0a1a2 ... an
4810 * b = b0b1b2 ... bm
4811 * c = c0c1c2 ... cl
4812 * a0 a1 ... an * bm
4813 * a0 a1 ... an * bm-1
4814 * . . .
4815 * . . .
4816 * a0 a1 .... an * b0
4817 * +_____________________________
4818 * c0 c1 c2 ...... cl
4819 * nc <---|
4820 * MaxAB |--------------------|
4821 */
4822 VP_EXPORT size_t
VpMult(Real * c,Real * a,Real * b)4823 VpMult(Real *c, Real *a, Real *b)
4824 {
4825 size_t MxIndA, MxIndB, MxIndAB, MxIndC;
4826 size_t ind_c, i, ii, nc;
4827 size_t ind_as, ind_ae, ind_bs;
4828 BDIGIT carry;
4829 BDIGIT_DBL s;
4830 Real *w;
4831
4832 #ifdef BIGDECIMAL_DEBUG
4833 if (gfDebug) {
4834 VPrint(stdout, "VpMult(Enter): a=% \n", a);
4835 VPrint(stdout, " b=% \n", b);
4836 }
4837 #endif /* BIGDECIMAL_DEBUG */
4838
4839 if (!VpIsDefOP(c, a, b, OP_SW_MULT)) return 0; /* No significant digit */
4840
4841 if (VpIsZero(a) || VpIsZero(b)) {
4842 /* at least a or b is zero */
4843 VpSetZero(c, VpGetSign(a) * VpGetSign(b));
4844 return 1; /* 0: 1 significant digit */
4845 }
4846
4847 if (VpIsOne(a)) {
4848 VpAsgn(c, b, VpGetSign(a));
4849 goto Exit;
4850 }
4851 if (VpIsOne(b)) {
4852 VpAsgn(c, a, VpGetSign(b));
4853 goto Exit;
4854 }
4855 if (b->Prec > a->Prec) {
4856 /* Adjust so that digits(a)>digits(b) */
4857 w = a;
4858 a = b;
4859 b = w;
4860 }
4861 w = NULL;
4862 MxIndA = a->Prec - 1;
4863 MxIndB = b->Prec - 1;
4864 MxIndC = c->MaxPrec - 1;
4865 MxIndAB = a->Prec + b->Prec - 1;
4866
4867 if (MxIndC < MxIndAB) { /* The Max. prec. of c < Prec(a)+Prec(b) */
4868 w = c;
4869 c = VpAlloc((size_t)((MxIndAB + 1) * BASE_FIG), "#0", 1, 1);
4870 MxIndC = MxIndAB;
4871 }
4872
4873 /* set LHSV c info */
4874
4875 c->exponent = a->exponent; /* set exponent */
4876 if (!AddExponent(c, b->exponent)) {
4877 if (w) VpFree(c);
4878 return 0;
4879 }
4880 VpSetSign(c, VpGetSign(a) * VpGetSign(b)); /* set sign */
4881 carry = 0;
4882 nc = ind_c = MxIndAB;
4883 memset(c->frac, 0, (nc + 1) * sizeof(BDIGIT)); /* Initialize c */
4884 c->Prec = nc + 1; /* set precision */
4885 for (nc = 0; nc < MxIndAB; ++nc, --ind_c) {
4886 if (nc < MxIndB) { /* The left triangle of the Fig. */
4887 ind_as = MxIndA - nc;
4888 ind_ae = MxIndA;
4889 ind_bs = MxIndB;
4890 }
4891 else if (nc <= MxIndA) { /* The middle rectangular of the Fig. */
4892 ind_as = MxIndA - nc;
4893 ind_ae = MxIndA - (nc - MxIndB);
4894 ind_bs = MxIndB;
4895 }
4896 else /* if (nc > MxIndA) */ { /* The right triangle of the Fig. */
4897 ind_as = 0;
4898 ind_ae = MxIndAB - nc - 1;
4899 ind_bs = MxIndB - (nc - MxIndA);
4900 }
4901
4902 for (i = ind_as; i <= ind_ae; ++i) {
4903 s = (BDIGIT_DBL)a->frac[i] * b->frac[ind_bs--];
4904 carry = (BDIGIT)(s / BASE);
4905 s -= (BDIGIT_DBL)carry * BASE;
4906 c->frac[ind_c] += (BDIGIT)s;
4907 if (c->frac[ind_c] >= BASE) {
4908 s = c->frac[ind_c] / BASE;
4909 carry += (BDIGIT)s;
4910 c->frac[ind_c] -= (BDIGIT)(s * BASE);
4911 }
4912 if (carry) {
4913 ii = ind_c;
4914 while (ii-- > 0) {
4915 c->frac[ii] += carry;
4916 if (c->frac[ii] >= BASE) {
4917 carry = c->frac[ii] / BASE;
4918 c->frac[ii] -= (carry * BASE);
4919 }
4920 else {
4921 break;
4922 }
4923 }
4924 }
4925 }
4926 }
4927 if (w != NULL) { /* free work variable */
4928 VpNmlz(c);
4929 VpAsgn(w, c, 1);
4930 VpFree(c);
4931 c = w;
4932 }
4933 else {
4934 VpLimitRound(c,0);
4935 }
4936
4937 Exit:
4938 #ifdef BIGDECIMAL_DEBUG
4939 if (gfDebug) {
4940 VPrint(stdout, "VpMult(c=a*b): c=% \n", c);
4941 VPrint(stdout, " a=% \n", a);
4942 VPrint(stdout, " b=% \n", b);
4943 }
4944 #endif /*BIGDECIMAL_DEBUG */
4945 return c->Prec*BASE_FIG;
4946 }
4947
4948 /*
4949 * c = a / b, remainder = r
4950 */
4951 VP_EXPORT size_t
VpDivd(Real * c,Real * r,Real * a,Real * b)4952 VpDivd(Real *c, Real *r, Real *a, Real *b)
4953 {
4954 size_t word_a, word_b, word_c, word_r;
4955 size_t i, n, ind_a, ind_b, ind_c, ind_r;
4956 size_t nLoop;
4957 BDIGIT_DBL q, b1, b1p1, b1b2, b1b2p1, r1r2;
4958 BDIGIT borrow, borrow1, borrow2;
4959 BDIGIT_DBL qb;
4960
4961 #ifdef BIGDECIMAL_DEBUG
4962 if (gfDebug) {
4963 VPrint(stdout, " VpDivd(c=a/b) a=% \n", a);
4964 VPrint(stdout, " b=% \n", b);
4965 }
4966 #endif /*BIGDECIMAL_DEBUG */
4967
4968 VpSetNaN(r);
4969 if (!VpIsDefOP(c, a, b, OP_SW_DIV)) goto Exit;
4970 if (VpIsZero(a) && VpIsZero(b)) {
4971 VpSetNaN(c);
4972 return VpException(VP_EXCEPTION_NaN, "Computation results to 'NaN'", 0);
4973 }
4974 if (VpIsZero(b)) {
4975 VpSetInf(c, VpGetSign(a) * VpGetSign(b));
4976 return VpException(VP_EXCEPTION_ZERODIVIDE, "Divide by zero", 0);
4977 }
4978 if (VpIsZero(a)) {
4979 /* numerator a is zero */
4980 VpSetZero(c, VpGetSign(a) * VpGetSign(b));
4981 VpSetZero(r, VpGetSign(a) * VpGetSign(b));
4982 goto Exit;
4983 }
4984 if (VpIsOne(b)) {
4985 /* divide by one */
4986 VpAsgn(c, a, VpGetSign(b));
4987 VpSetZero(r, VpGetSign(a));
4988 goto Exit;
4989 }
4990
4991 word_a = a->Prec;
4992 word_b = b->Prec;
4993 word_c = c->MaxPrec;
4994 word_r = r->MaxPrec;
4995
4996 ind_c = 0;
4997 ind_r = 1;
4998
4999 if (word_a >= word_r) goto space_error;
5000
5001 r->frac[0] = 0;
5002 while (ind_r <= word_a) {
5003 r->frac[ind_r] = a->frac[ind_r - 1];
5004 ++ind_r;
5005 }
5006
5007 while (ind_r < word_r) r->frac[ind_r++] = 0;
5008 while (ind_c < word_c) c->frac[ind_c++] = 0;
5009
5010 /* initial procedure */
5011 b1 = b1p1 = b->frac[0];
5012 if (b->Prec <= 1) {
5013 b1b2p1 = b1b2 = b1p1 * BASE;
5014 }
5015 else {
5016 b1p1 = b1 + 1;
5017 b1b2p1 = b1b2 = b1 * BASE + b->frac[1];
5018 if (b->Prec > 2) ++b1b2p1;
5019 }
5020
5021 /* */
5022 /* loop start */
5023 ind_c = word_r - 1;
5024 nLoop = Min(word_c,ind_c);
5025 ind_c = 1;
5026 while (ind_c < nLoop) {
5027 if (r->frac[ind_c] == 0) {
5028 ++ind_c;
5029 continue;
5030 }
5031 r1r2 = (BDIGIT_DBL)r->frac[ind_c] * BASE + r->frac[ind_c + 1];
5032 if (r1r2 == b1b2) {
5033 /* The first two word digits is the same */
5034 ind_b = 2;
5035 ind_a = ind_c + 2;
5036 while (ind_b < word_b) {
5037 if (r->frac[ind_a] < b->frac[ind_b]) goto div_b1p1;
5038 if (r->frac[ind_a] > b->frac[ind_b]) break;
5039 ++ind_a;
5040 ++ind_b;
5041 }
5042 /* The first few word digits of r and b is the same and */
5043 /* the first different word digit of w is greater than that */
5044 /* of b, so quotient is 1 and just subtract b from r. */
5045 borrow = 0; /* quotient=1, then just r-b */
5046 ind_b = b->Prec - 1;
5047 ind_r = ind_c + ind_b;
5048 if (ind_r >= word_r) goto space_error;
5049 n = ind_b;
5050 for (i = 0; i <= n; ++i) {
5051 if (r->frac[ind_r] < b->frac[ind_b] + borrow) {
5052 r->frac[ind_r] += (BASE - (b->frac[ind_b] + borrow));
5053 borrow = 1;
5054 }
5055 else {
5056 r->frac[ind_r] = r->frac[ind_r] - b->frac[ind_b] - borrow;
5057 borrow = 0;
5058 }
5059 --ind_r;
5060 --ind_b;
5061 }
5062 ++c->frac[ind_c];
5063 goto carry;
5064 }
5065 /* The first two word digits is not the same, */
5066 /* then compare magnitude, and divide actually. */
5067 if (r1r2 >= b1b2p1) {
5068 q = r1r2 / b1b2p1; /* q == (BDIGIT)q */
5069 c->frac[ind_c] += (BDIGIT)q;
5070 ind_r = b->Prec + ind_c - 1;
5071 goto sub_mult;
5072 }
5073
5074 div_b1p1:
5075 if (ind_c + 1 >= word_c) goto out_side;
5076 q = r1r2 / b1p1; /* q == (BDIGIT)q */
5077 c->frac[ind_c + 1] += (BDIGIT)q;
5078 ind_r = b->Prec + ind_c;
5079
5080 sub_mult:
5081 borrow1 = borrow2 = 0;
5082 ind_b = word_b - 1;
5083 if (ind_r >= word_r) goto space_error;
5084 n = ind_b;
5085 for (i = 0; i <= n; ++i) {
5086 /* now, perform r = r - q * b */
5087 qb = q * b->frac[ind_b];
5088 if (qb < BASE) borrow1 = 0;
5089 else {
5090 borrow1 = (BDIGIT)(qb / BASE);
5091 qb -= (BDIGIT_DBL)borrow1 * BASE; /* get qb < BASE */
5092 }
5093 if(r->frac[ind_r] < qb) {
5094 r->frac[ind_r] += (BDIGIT)(BASE - qb);
5095 borrow2 = borrow2 + borrow1 + 1;
5096 }
5097 else {
5098 r->frac[ind_r] -= (BDIGIT)qb;
5099 borrow2 += borrow1;
5100 }
5101 if (borrow2) {
5102 if(r->frac[ind_r - 1] < borrow2) {
5103 r->frac[ind_r - 1] += (BASE - borrow2);
5104 borrow2 = 1;
5105 }
5106 else {
5107 r->frac[ind_r - 1] -= borrow2;
5108 borrow2 = 0;
5109 }
5110 }
5111 --ind_r;
5112 --ind_b;
5113 }
5114
5115 r->frac[ind_r] -= borrow2;
5116 carry:
5117 ind_r = ind_c;
5118 while (c->frac[ind_r] >= BASE) {
5119 c->frac[ind_r] -= BASE;
5120 --ind_r;
5121 ++c->frac[ind_r];
5122 }
5123 }
5124 /* End of operation, now final arrangement */
5125 out_side:
5126 c->Prec = word_c;
5127 c->exponent = a->exponent;
5128 if (!AddExponent(c, 2)) return 0;
5129 if (!AddExponent(c, -(b->exponent))) return 0;
5130
5131 VpSetSign(c, VpGetSign(a) * VpGetSign(b));
5132 VpNmlz(c); /* normalize c */
5133 r->Prec = word_r;
5134 r->exponent = a->exponent;
5135 if (!AddExponent(r, 1)) return 0;
5136 VpSetSign(r, VpGetSign(a));
5137 VpNmlz(r); /* normalize r(remainder) */
5138 goto Exit;
5139
5140 space_error:
5141 #ifdef BIGDECIMAL_DEBUG
5142 if (gfDebug) {
5143 printf(" word_a=%"PRIuSIZE"\n", word_a);
5144 printf(" word_b=%"PRIuSIZE"\n", word_b);
5145 printf(" word_c=%"PRIuSIZE"\n", word_c);
5146 printf(" word_r=%"PRIuSIZE"\n", word_r);
5147 printf(" ind_r =%"PRIuSIZE"\n", ind_r);
5148 }
5149 #endif /* BIGDECIMAL_DEBUG */
5150 rb_bug("ERROR(VpDivd): space for remainder too small.");
5151
5152 Exit:
5153 #ifdef BIGDECIMAL_DEBUG
5154 if (gfDebug) {
5155 VPrint(stdout, " VpDivd(c=a/b), c=% \n", c);
5156 VPrint(stdout, " r=% \n", r);
5157 }
5158 #endif /* BIGDECIMAL_DEBUG */
5159 return c->Prec * BASE_FIG;
5160 }
5161
5162 /*
5163 * Input a = 00000xxxxxxxx En(5 preceding zeros)
5164 * Output a = xxxxxxxx En-5
5165 */
5166 static int
VpNmlz(Real * a)5167 VpNmlz(Real *a)
5168 {
5169 size_t ind_a, i;
5170
5171 if (!VpIsDef(a)) goto NoVal;
5172 if (VpIsZero(a)) goto NoVal;
5173
5174 ind_a = a->Prec;
5175 while (ind_a--) {
5176 if (a->frac[ind_a]) {
5177 a->Prec = ind_a + 1;
5178 i = 0;
5179 while (a->frac[i] == 0) ++i; /* skip the first few zeros */
5180 if (i) {
5181 a->Prec -= i;
5182 if (!AddExponent(a, -(SIGNED_VALUE)i)) return 0;
5183 memmove(&a->frac[0], &a->frac[i], a->Prec*sizeof(BDIGIT));
5184 }
5185 return 1;
5186 }
5187 }
5188 /* a is zero(no non-zero digit) */
5189 VpSetZero(a, VpGetSign(a));
5190 return 0;
5191
5192 NoVal:
5193 a->frac[0] = 0;
5194 a->Prec = 1;
5195 return 0;
5196 }
5197
5198 /*
5199 * VpComp = 0 ... if a=b,
5200 * Pos ... a>b,
5201 * Neg ... a<b.
5202 * 999 ... result undefined(NaN)
5203 */
5204 VP_EXPORT int
VpComp(Real * a,Real * b)5205 VpComp(Real *a, Real *b)
5206 {
5207 int val;
5208 size_t mx, ind;
5209 int e;
5210 val = 0;
5211 if (VpIsNaN(a) || VpIsNaN(b)) return 999;
5212 if (!VpIsDef(a)) {
5213 if (!VpIsDef(b)) e = a->sign - b->sign;
5214 else e = a->sign;
5215
5216 if (e > 0) return 1;
5217 else if (e < 0) return -1;
5218 else return 0;
5219 }
5220 if (!VpIsDef(b)) {
5221 e = -b->sign;
5222 if (e > 0) return 1;
5223 else return -1;
5224 }
5225 /* Zero check */
5226 if (VpIsZero(a)) {
5227 if (VpIsZero(b)) return 0; /* both zero */
5228 val = -VpGetSign(b);
5229 goto Exit;
5230 }
5231 if (VpIsZero(b)) {
5232 val = VpGetSign(a);
5233 goto Exit;
5234 }
5235
5236 /* compare sign */
5237 if (VpGetSign(a) > VpGetSign(b)) {
5238 val = 1; /* a>b */
5239 goto Exit;
5240 }
5241 if (VpGetSign(a) < VpGetSign(b)) {
5242 val = -1; /* a<b */
5243 goto Exit;
5244 }
5245
5246 /* a and b have same sign, && sign!=0,then compare exponent */
5247 if (a->exponent > b->exponent) {
5248 val = VpGetSign(a);
5249 goto Exit;
5250 }
5251 if (a->exponent < b->exponent) {
5252 val = -VpGetSign(b);
5253 goto Exit;
5254 }
5255
5256 /* a and b have same exponent, then compare their significand. */
5257 mx = (a->Prec < b->Prec) ? a->Prec : b->Prec;
5258 ind = 0;
5259 while (ind < mx) {
5260 if (a->frac[ind] > b->frac[ind]) {
5261 val = VpGetSign(a);
5262 goto Exit;
5263 }
5264 if (a->frac[ind] < b->frac[ind]) {
5265 val = -VpGetSign(b);
5266 goto Exit;
5267 }
5268 ++ind;
5269 }
5270 if (a->Prec > b->Prec) {
5271 val = VpGetSign(a);
5272 }
5273 else if (a->Prec < b->Prec) {
5274 val = -VpGetSign(b);
5275 }
5276
5277 Exit:
5278 if (val > 1) val = 1;
5279 else if (val < -1) val = -1;
5280
5281 #ifdef BIGDECIMAL_DEBUG
5282 if (gfDebug) {
5283 VPrint(stdout, " VpComp a=%\n", a);
5284 VPrint(stdout, " b=%\n", b);
5285 printf(" ans=%d\n", val);
5286 }
5287 #endif /* BIGDECIMAL_DEBUG */
5288 return (int)val;
5289 }
5290
5291 /*
5292 * cntl_chr ... ASCIIZ Character, print control characters
5293 * Available control codes:
5294 * % ... VP variable. To print '%', use '%%'.
5295 * \n ... new line
5296 * \b ... backspace
5297 * \t ... tab
5298 * Note: % must not appear more than once
5299 * a ... VP variable to be printed
5300 */
5301 #ifdef BIGDECIMAL_ENABLE_VPRINT
5302 static int
VPrint(FILE * fp,const char * cntl_chr,Real * a)5303 VPrint(FILE *fp, const char *cntl_chr, Real *a)
5304 {
5305 size_t i, j, nc, nd, ZeroSup, sep = 10;
5306 BDIGIT m, e, nn;
5307
5308 j = 0;
5309 nd = nc = 0; /* nd : number of digits in fraction part(every 10 digits, */
5310 /* nd<=10). */
5311 /* nc : number of characters printed */
5312 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
5313 while (*(cntl_chr + j)) {
5314 if (*(cntl_chr + j) == '%' && *(cntl_chr + j + 1) != '%') {
5315 nc = 0;
5316 if (VpIsNaN(a)) {
5317 fprintf(fp, SZ_NaN);
5318 nc += 8;
5319 }
5320 else if (VpIsPosInf(a)) {
5321 fprintf(fp, SZ_INF);
5322 nc += 8;
5323 }
5324 else if (VpIsNegInf(a)) {
5325 fprintf(fp, SZ_NINF);
5326 nc += 9;
5327 }
5328 else if (!VpIsZero(a)) {
5329 if (BIGDECIMAL_NEGATIVE_P(a)) {
5330 fprintf(fp, "-");
5331 ++nc;
5332 }
5333 nc += fprintf(fp, "0.");
5334 switch (*(cntl_chr + j + 1)) {
5335 default:
5336 break;
5337
5338 case '0': case 'z':
5339 ZeroSup = 0;
5340 ++j;
5341 sep = cntl_chr[j] == 'z' ? RMPD_COMPONENT_FIGURES : 10;
5342 break;
5343 }
5344 for (i = 0; i < a->Prec; ++i) {
5345 m = BASE1;
5346 e = a->frac[i];
5347 while (m) {
5348 nn = e / m;
5349 if (!ZeroSup || nn) {
5350 nc += fprintf(fp, "%lu", (unsigned long)nn); /* The leading zero(s) */
5351 /* as 0.00xx will not */
5352 /* be printed. */
5353 ++nd;
5354 ZeroSup = 0; /* Set to print succeeding zeros */
5355 }
5356 if (nd >= sep) { /* print ' ' after every 10 digits */
5357 nd = 0;
5358 nc += fprintf(fp, " ");
5359 }
5360 e = e - nn * m;
5361 m /= 10;
5362 }
5363 }
5364 nc += fprintf(fp, "E%"PRIdSIZE, VpExponent10(a));
5365 nc += fprintf(fp, " (%"PRIdVALUE", %lu, %lu)", a->exponent, a->Prec, a->MaxPrec);
5366 }
5367 else {
5368 nc += fprintf(fp, "0.0");
5369 }
5370 }
5371 else {
5372 ++nc;
5373 if (*(cntl_chr + j) == '\\') {
5374 switch (*(cntl_chr + j + 1)) {
5375 case 'n':
5376 fprintf(fp, "\n");
5377 ++j;
5378 break;
5379 case 't':
5380 fprintf(fp, "\t");
5381 ++j;
5382 break;
5383 case 'b':
5384 fprintf(fp, "\n");
5385 ++j;
5386 break;
5387 default:
5388 fprintf(fp, "%c", *(cntl_chr + j));
5389 break;
5390 }
5391 }
5392 else {
5393 fprintf(fp, "%c", *(cntl_chr + j));
5394 if (*(cntl_chr + j) == '%') ++j;
5395 }
5396 }
5397 j++;
5398 }
5399
5400 return (int)nc;
5401 }
5402 #endif
5403
5404 static void
VpFormatSt(char * psz,size_t fFmt)5405 VpFormatSt(char *psz, size_t fFmt)
5406 {
5407 size_t ie, i, nf = 0;
5408 char ch;
5409
5410 if (fFmt == 0) return;
5411
5412 ie = strlen(psz);
5413 for (i = 0; i < ie; ++i) {
5414 ch = psz[i];
5415 if (!ch) break;
5416 if (ISSPACE(ch) || ch=='-' || ch=='+') continue;
5417 if (ch == '.') { nf = 0; continue; }
5418 if (ch == 'E' || ch == 'e') break;
5419
5420 if (++nf > fFmt) {
5421 memmove(psz + i + 1, psz + i, ie - i + 1);
5422 ++ie;
5423 nf = 0;
5424 psz[i] = ' ';
5425 }
5426 }
5427 }
5428
5429 VP_EXPORT ssize_t
VpExponent10(Real * a)5430 VpExponent10(Real *a)
5431 {
5432 ssize_t ex;
5433 size_t n;
5434
5435 if (!VpHasVal(a)) return 0;
5436
5437 ex = a->exponent * (ssize_t)BASE_FIG;
5438 n = BASE1;
5439 while ((a->frac[0] / n) == 0) {
5440 --ex;
5441 n /= 10;
5442 }
5443 return ex;
5444 }
5445
5446 VP_EXPORT void
VpSzMantissa(Real * a,char * psz)5447 VpSzMantissa(Real *a,char *psz)
5448 {
5449 size_t i, n, ZeroSup;
5450 BDIGIT_DBL m, e, nn;
5451
5452 if (VpIsNaN(a)) {
5453 sprintf(psz, SZ_NaN);
5454 return;
5455 }
5456 if (VpIsPosInf(a)) {
5457 sprintf(psz, SZ_INF);
5458 return;
5459 }
5460 if (VpIsNegInf(a)) {
5461 sprintf(psz, SZ_NINF);
5462 return;
5463 }
5464
5465 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
5466 if (!VpIsZero(a)) {
5467 if (BIGDECIMAL_NEGATIVE_P(a)) *psz++ = '-';
5468 n = a->Prec;
5469 for (i = 0; i < n; ++i) {
5470 m = BASE1;
5471 e = a->frac[i];
5472 while (m) {
5473 nn = e / m;
5474 if (!ZeroSup || nn) {
5475 sprintf(psz, "%lu", (unsigned long)nn); /* The leading zero(s) */
5476 psz += strlen(psz);
5477 /* as 0.00xx will be ignored. */
5478 ZeroSup = 0; /* Set to print succeeding zeros */
5479 }
5480 e = e - nn * m;
5481 m /= 10;
5482 }
5483 }
5484 *psz = 0;
5485 while (psz[-1] == '0') *(--psz) = 0;
5486 }
5487 else {
5488 if (VpIsPosZero(a)) sprintf(psz, "0");
5489 else sprintf(psz, "-0");
5490 }
5491 }
5492
5493 VP_EXPORT int
VpToSpecialString(Real * a,char * psz,int fPlus)5494 VpToSpecialString(Real *a,char *psz,int fPlus)
5495 /* fPlus = 0: default, 1: set ' ' before digits, 2: set '+' before digits. */
5496 {
5497 if (VpIsNaN(a)) {
5498 sprintf(psz,SZ_NaN);
5499 return 1;
5500 }
5501
5502 if (VpIsPosInf(a)) {
5503 if (fPlus == 1) {
5504 *psz++ = ' ';
5505 }
5506 else if (fPlus == 2) {
5507 *psz++ = '+';
5508 }
5509 sprintf(psz, SZ_INF);
5510 return 1;
5511 }
5512 if (VpIsNegInf(a)) {
5513 sprintf(psz, SZ_NINF);
5514 return 1;
5515 }
5516 if (VpIsZero(a)) {
5517 if (VpIsPosZero(a)) {
5518 if (fPlus == 1) sprintf(psz, " 0.0");
5519 else if (fPlus == 2) sprintf(psz, "+0.0");
5520 else sprintf(psz, "0.0");
5521 }
5522 else sprintf(psz, "-0.0");
5523 return 1;
5524 }
5525 return 0;
5526 }
5527
5528 VP_EXPORT void
VpToString(Real * a,char * psz,size_t fFmt,int fPlus)5529 VpToString(Real *a, char *psz, size_t fFmt, int fPlus)
5530 /* fPlus = 0: default, 1: set ' ' before digits, 2: set '+' before digits. */
5531 {
5532 size_t i, n, ZeroSup;
5533 BDIGIT shift, m, e, nn;
5534 char *pszSav = psz;
5535 ssize_t ex;
5536
5537 if (VpToSpecialString(a, psz, fPlus)) return;
5538
5539 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
5540
5541 if (BIGDECIMAL_NEGATIVE_P(a)) *psz++ = '-';
5542 else if (fPlus == 1) *psz++ = ' ';
5543 else if (fPlus == 2) *psz++ = '+';
5544
5545 *psz++ = '0';
5546 *psz++ = '.';
5547 n = a->Prec;
5548 for (i = 0; i < n; ++i) {
5549 m = BASE1;
5550 e = a->frac[i];
5551 while (m) {
5552 nn = e / m;
5553 if (!ZeroSup || nn) {
5554 sprintf(psz, "%lu", (unsigned long)nn); /* The reading zero(s) */
5555 psz += strlen(psz);
5556 /* as 0.00xx will be ignored. */
5557 ZeroSup = 0; /* Set to print succeeding zeros */
5558 }
5559 e = e - nn * m;
5560 m /= 10;
5561 }
5562 }
5563 ex = a->exponent * (ssize_t)BASE_FIG;
5564 shift = BASE1;
5565 while (a->frac[0] / shift == 0) {
5566 --ex;
5567 shift /= 10;
5568 }
5569 while (psz[-1] == '0') {
5570 *(--psz) = 0;
5571 }
5572 sprintf(psz, "e%"PRIdSIZE, ex);
5573 if (fFmt) VpFormatSt(pszSav, fFmt);
5574 }
5575
5576 VP_EXPORT void
VpToFString(Real * a,char * psz,size_t fFmt,int fPlus)5577 VpToFString(Real *a, char *psz, size_t fFmt, int fPlus)
5578 /* fPlus = 0: default, 1: set ' ' before digits, 2: set '+' before digits. */
5579 {
5580 size_t i, n;
5581 BDIGIT m, e, nn;
5582 char *pszSav = psz;
5583 ssize_t ex;
5584
5585 if (VpToSpecialString(a, psz, fPlus)) return;
5586
5587 if (BIGDECIMAL_NEGATIVE_P(a)) *psz++ = '-';
5588 else if (fPlus == 1) *psz++ = ' ';
5589 else if (fPlus == 2) *psz++ = '+';
5590
5591 n = a->Prec;
5592 ex = a->exponent;
5593 if (ex <= 0) {
5594 *psz++ = '0';*psz++ = '.';
5595 while (ex < 0) {
5596 for (i=0; i < BASE_FIG; ++i) *psz++ = '0';
5597 ++ex;
5598 }
5599 ex = -1;
5600 }
5601
5602 for (i = 0; i < n; ++i) {
5603 --ex;
5604 if (i == 0 && ex >= 0) {
5605 sprintf(psz, "%lu", (unsigned long)a->frac[i]);
5606 psz += strlen(psz);
5607 }
5608 else {
5609 m = BASE1;
5610 e = a->frac[i];
5611 while (m) {
5612 nn = e / m;
5613 *psz++ = (char)(nn + '0');
5614 e = e - nn * m;
5615 m /= 10;
5616 }
5617 }
5618 if (ex == 0) *psz++ = '.';
5619 }
5620 while (--ex>=0) {
5621 m = BASE;
5622 while (m /= 10) *psz++ = '0';
5623 if (ex == 0) *psz++ = '.';
5624 }
5625 *psz = 0;
5626 while (psz[-1] == '0') *(--psz) = 0;
5627 if (psz[-1] == '.') sprintf(psz, "0");
5628 if (fFmt) VpFormatSt(pszSav, fFmt);
5629 }
5630
5631 /*
5632 * [Output]
5633 * a[] ... variable to be assigned the value.
5634 * [Input]
5635 * int_chr[] ... integer part(may include '+/-').
5636 * ni ... number of characters in int_chr[],not including '+/-'.
5637 * frac[] ... fraction part.
5638 * nf ... number of characters in frac[].
5639 * exp_chr[] ... exponent part(including '+/-').
5640 * ne ... number of characters in exp_chr[],not including '+/-'.
5641 */
5642 VP_EXPORT int
VpCtoV(Real * a,const char * int_chr,size_t ni,const char * frac,size_t nf,const char * exp_chr,size_t ne)5643 VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne)
5644 {
5645 size_t i, j, ind_a, ma, mi, me;
5646 SIGNED_VALUE e, es, eb, ef;
5647 int sign, signe, exponent_overflow;
5648
5649 /* get exponent part */
5650 e = 0;
5651 ma = a->MaxPrec;
5652 mi = ni;
5653 me = ne;
5654 signe = 1;
5655 exponent_overflow = 0;
5656 memset(a->frac, 0, ma * sizeof(BDIGIT));
5657 if (ne > 0) {
5658 i = 0;
5659 if (exp_chr[0] == '-') {
5660 signe = -1;
5661 ++i;
5662 ++me;
5663 }
5664 else if (exp_chr[0] == '+') {
5665 ++i;
5666 ++me;
5667 }
5668 while (i < me) {
5669 if (MUL_OVERFLOW_SIGNED_VALUE_P(e, (SIGNED_VALUE)BASE_FIG)) {
5670 es = e;
5671 goto exp_overflow;
5672 }
5673 es = e * (SIGNED_VALUE)BASE_FIG;
5674 if (MUL_OVERFLOW_SIGNED_VALUE_P(e, 10) ||
5675 SIGNED_VALUE_MAX - (exp_chr[i] - '0') < e * 10)
5676 goto exp_overflow;
5677 e = e * 10 + exp_chr[i] - '0';
5678 if (MUL_OVERFLOW_SIGNED_VALUE_P(e, (SIGNED_VALUE)BASE_FIG))
5679 goto exp_overflow;
5680 if (es > (SIGNED_VALUE)(e * BASE_FIG)) {
5681 exp_overflow:
5682 exponent_overflow = 1;
5683 e = es; /* keep sign */
5684 break;
5685 }
5686 ++i;
5687 }
5688 }
5689
5690 /* get integer part */
5691 i = 0;
5692 sign = 1;
5693 if (1 /*ni >= 0*/) {
5694 if (int_chr[0] == '-') {
5695 sign = -1;
5696 ++i;
5697 ++mi;
5698 }
5699 else if (int_chr[0] == '+') {
5700 ++i;
5701 ++mi;
5702 }
5703 }
5704
5705 e = signe * e; /* e: The value of exponent part. */
5706 e = e + ni; /* set actual exponent size. */
5707
5708 if (e > 0) signe = 1;
5709 else signe = -1;
5710
5711 /* Adjust the exponent so that it is the multiple of BASE_FIG. */
5712 j = 0;
5713 ef = 1;
5714 while (ef) {
5715 if (e >= 0) eb = e;
5716 else eb = -e;
5717 ef = eb / (SIGNED_VALUE)BASE_FIG;
5718 ef = eb - ef * (SIGNED_VALUE)BASE_FIG;
5719 if (ef) {
5720 ++j; /* Means to add one more preceding zero */
5721 ++e;
5722 }
5723 }
5724
5725 eb = e / (SIGNED_VALUE)BASE_FIG;
5726
5727 if (exponent_overflow) {
5728 int zero = 1;
5729 for ( ; i < mi && zero; i++) zero = int_chr[i] == '0';
5730 for (i = 0; i < nf && zero; i++) zero = frac[i] == '0';
5731 if (!zero && signe > 0) {
5732 VpSetInf(a, sign);
5733 VpException(VP_EXCEPTION_INFINITY, "exponent overflow",0);
5734 }
5735 else VpSetZero(a, sign);
5736 return 1;
5737 }
5738
5739 ind_a = 0;
5740 while (i < mi) {
5741 a->frac[ind_a] = 0;
5742 while (j < BASE_FIG && i < mi) {
5743 a->frac[ind_a] = a->frac[ind_a] * 10 + int_chr[i] - '0';
5744 ++j;
5745 ++i;
5746 }
5747 if (i < mi) {
5748 ++ind_a;
5749 if (ind_a >= ma) goto over_flow;
5750 j = 0;
5751 }
5752 }
5753
5754 /* get fraction part */
5755
5756 i = 0;
5757 while (i < nf) {
5758 while (j < BASE_FIG && i < nf) {
5759 a->frac[ind_a] = a->frac[ind_a] * 10 + frac[i] - '0';
5760 ++j;
5761 ++i;
5762 }
5763 if (i < nf) {
5764 ++ind_a;
5765 if (ind_a >= ma) goto over_flow;
5766 j = 0;
5767 }
5768 }
5769 goto Final;
5770
5771 over_flow:
5772 rb_warn("Conversion from String to BigDecimal overflow (last few digits discarded).");
5773
5774 Final:
5775 if (ind_a >= ma) ind_a = ma - 1;
5776 while (j < BASE_FIG) {
5777 a->frac[ind_a] = a->frac[ind_a] * 10;
5778 ++j;
5779 }
5780 a->Prec = ind_a + 1;
5781 a->exponent = eb;
5782 VpSetSign(a, sign);
5783 VpNmlz(a);
5784 return 1;
5785 }
5786
5787 /*
5788 * [Input]
5789 * *m ... Real
5790 * [Output]
5791 * *d ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig.
5792 * *e ... exponent of m.
5793 * DBLE_FIG ... Number of digits in a double variable.
5794 *
5795 * m -> d*10**e, 0<d<BASE
5796 * [Returns]
5797 * 0 ... Zero
5798 * 1 ... Normal
5799 * 2 ... Infinity
5800 * -1 ... NaN
5801 */
5802 VP_EXPORT int
VpVtoD(double * d,SIGNED_VALUE * e,Real * m)5803 VpVtoD(double *d, SIGNED_VALUE *e, Real *m)
5804 {
5805 size_t ind_m, mm, fig;
5806 double div;
5807 int f = 1;
5808
5809 if (VpIsNaN(m)) {
5810 *d = VpGetDoubleNaN();
5811 *e = 0;
5812 f = -1; /* NaN */
5813 goto Exit;
5814 }
5815 else if (VpIsPosZero(m)) {
5816 *d = 0.0;
5817 *e = 0;
5818 f = 0;
5819 goto Exit;
5820 }
5821 else if (VpIsNegZero(m)) {
5822 *d = VpGetDoubleNegZero();
5823 *e = 0;
5824 f = 0;
5825 goto Exit;
5826 }
5827 else if (VpIsPosInf(m)) {
5828 *d = VpGetDoublePosInf();
5829 *e = 0;
5830 f = 2;
5831 goto Exit;
5832 }
5833 else if (VpIsNegInf(m)) {
5834 *d = VpGetDoubleNegInf();
5835 *e = 0;
5836 f = 2;
5837 goto Exit;
5838 }
5839 /* Normal number */
5840 fig = (DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
5841 ind_m = 0;
5842 mm = Min(fig, m->Prec);
5843 *d = 0.0;
5844 div = 1.;
5845 while (ind_m < mm) {
5846 div /= (double)BASE;
5847 *d = *d + (double)m->frac[ind_m++] * div;
5848 }
5849 *e = m->exponent * (SIGNED_VALUE)BASE_FIG;
5850 *d *= VpGetSign(m);
5851
5852 Exit:
5853 #ifdef BIGDECIMAL_DEBUG
5854 if (gfDebug) {
5855 VPrint(stdout, " VpVtoD: m=%\n", m);
5856 printf(" d=%e * 10 **%ld\n", *d, *e);
5857 printf(" DBLE_FIG = %d\n", DBLE_FIG);
5858 }
5859 #endif /*BIGDECIMAL_DEBUG */
5860 return f;
5861 }
5862
5863 /*
5864 * m <- d
5865 */
5866 VP_EXPORT void
VpDtoV(Real * m,double d)5867 VpDtoV(Real *m, double d)
5868 {
5869 size_t ind_m, mm;
5870 SIGNED_VALUE ne;
5871 BDIGIT i;
5872 double val, val2;
5873
5874 if (isnan(d)) {
5875 VpSetNaN(m);
5876 goto Exit;
5877 }
5878 if (isinf(d)) {
5879 if (d > 0.0) VpSetPosInf(m);
5880 else VpSetNegInf(m);
5881 goto Exit;
5882 }
5883
5884 if (d == 0.0) {
5885 VpSetZero(m, 1);
5886 goto Exit;
5887 }
5888 val = (d > 0.) ? d : -d;
5889 ne = 0;
5890 if (val >= 1.0) {
5891 while (val >= 1.0) {
5892 val /= (double)BASE;
5893 ++ne;
5894 }
5895 }
5896 else {
5897 val2 = 1.0 / (double)BASE;
5898 while (val < val2) {
5899 val *= (double)BASE;
5900 --ne;
5901 }
5902 }
5903 /* Now val = 0.xxxxx*BASE**ne */
5904
5905 mm = m->MaxPrec;
5906 memset(m->frac, 0, mm * sizeof(BDIGIT));
5907 for (ind_m = 0; val > 0.0 && ind_m < mm; ind_m++) {
5908 val *= (double)BASE;
5909 i = (BDIGIT)val;
5910 val -= (double)i;
5911 m->frac[ind_m] = i;
5912 }
5913 if (ind_m >= mm) ind_m = mm - 1;
5914 VpSetSign(m, (d > 0.0) ? 1 : -1);
5915 m->Prec = ind_m + 1;
5916 m->exponent = ne;
5917
5918 VpInternalRound(m, 0, (m->Prec > 0) ? m->frac[m->Prec-1] : 0,
5919 (BDIGIT)(val*(double)BASE));
5920
5921 Exit:
5922 #ifdef BIGDECIMAL_DEBUG
5923 if (gfDebug) {
5924 printf("VpDtoV d=%30.30e\n", d);
5925 VPrint(stdout, " m=%\n", m);
5926 }
5927 #endif /* BIGDECIMAL_DEBUG */
5928 return;
5929 }
5930
5931 /*
5932 * m <- ival
5933 */
5934 #if 0 /* unused */
5935 VP_EXPORT void
5936 VpItoV(Real *m, SIGNED_VALUE ival)
5937 {
5938 size_t mm, ind_m;
5939 size_t val, v1, v2, v;
5940 int isign;
5941 SIGNED_VALUE ne;
5942
5943 if (ival == 0) {
5944 VpSetZero(m, 1);
5945 goto Exit;
5946 }
5947 isign = 1;
5948 val = ival;
5949 if (ival < 0) {
5950 isign = -1;
5951 val =(size_t)(-ival);
5952 }
5953 ne = 0;
5954 ind_m = 0;
5955 mm = m->MaxPrec;
5956 while (ind_m < mm) {
5957 m->frac[ind_m] = 0;
5958 ++ind_m;
5959 }
5960 ind_m = 0;
5961 while (val > 0) {
5962 if (val) {
5963 v1 = val;
5964 v2 = 1;
5965 while (v1 >= BASE) {
5966 v1 /= BASE;
5967 v2 *= BASE;
5968 }
5969 val = val - v2 * v1;
5970 v = v1;
5971 }
5972 else {
5973 v = 0;
5974 }
5975 m->frac[ind_m] = v;
5976 ++ind_m;
5977 ++ne;
5978 }
5979 m->Prec = ind_m - 1;
5980 m->exponent = ne;
5981 VpSetSign(m, isign);
5982 VpNmlz(m);
5983
5984 Exit:
5985 #ifdef BIGDECIMAL_DEBUG
5986 if (gfDebug) {
5987 printf(" VpItoV i=%d\n", ival);
5988 VPrint(stdout, " m=%\n", m);
5989 }
5990 #endif /* BIGDECIMAL_DEBUG */
5991 return;
5992 }
5993 #endif
5994
5995 /*
5996 * y = SQRT(x), y*y - x =>0
5997 */
5998 VP_EXPORT int
VpSqrt(Real * y,Real * x)5999 VpSqrt(Real *y, Real *x)
6000 {
6001 Real *f = NULL;
6002 Real *r = NULL;
6003 size_t y_prec;
6004 SIGNED_VALUE n, e;
6005 SIGNED_VALUE prec;
6006 ssize_t nr;
6007 double val;
6008
6009 /* Zero or +Infinity ? */
6010 if (VpIsZero(x) || VpIsPosInf(x)) {
6011 VpAsgn(y,x,1);
6012 goto Exit;
6013 }
6014
6015 /* Negative ? */
6016 if (BIGDECIMAL_NEGATIVE_P(x)) {
6017 VpSetNaN(y);
6018 return VpException(VP_EXCEPTION_OP, "sqrt of negative value", 0);
6019 }
6020
6021 /* NaN ? */
6022 if (VpIsNaN(x)) {
6023 VpSetNaN(y);
6024 return VpException(VP_EXCEPTION_OP, "sqrt of 'NaN'(Not a Number)", 0);
6025 }
6026
6027 /* One ? */
6028 if (VpIsOne(x)) {
6029 VpSetOne(y);
6030 goto Exit;
6031 }
6032
6033 n = (SIGNED_VALUE)y->MaxPrec;
6034 if (x->MaxPrec > (size_t)n) n = (ssize_t)x->MaxPrec;
6035
6036 /* allocate temporally variables */
6037 f = VpAlloc(y->MaxPrec * (BASE_FIG + 2), "#1", 1, 1);
6038 r = VpAlloc((n + n) * (BASE_FIG + 2), "#1", 1, 1);
6039
6040 nr = 0;
6041 y_prec = y->MaxPrec;
6042
6043 prec = x->exponent - (ssize_t)y_prec;
6044 if (x->exponent > 0)
6045 ++prec;
6046 else
6047 --prec;
6048
6049 VpVtoD(&val, &e, x); /* val <- x */
6050 e /= (SIGNED_VALUE)BASE_FIG;
6051 n = e / 2;
6052 if (e - n * 2 != 0) {
6053 val /= BASE;
6054 n = (e + 1) / 2;
6055 }
6056 VpDtoV(y, sqrt(val)); /* y <- sqrt(val) */
6057 y->exponent += n;
6058 n = (SIGNED_VALUE)((DBLE_FIG + BASE_FIG - 1) / BASE_FIG);
6059 y->MaxPrec = Min((size_t)n , y_prec);
6060 f->MaxPrec = y->MaxPrec + 1;
6061 n = (SIGNED_VALUE)(y_prec * BASE_FIG);
6062 if (n < (SIGNED_VALUE)maxnr) n = (SIGNED_VALUE)maxnr;
6063 do {
6064 y->MaxPrec *= 2;
6065 if (y->MaxPrec > y_prec) y->MaxPrec = y_prec;
6066 f->MaxPrec = y->MaxPrec;
6067 VpDivd(f, r, x, y); /* f = x/y */
6068 VpAddSub(r, f, y, -1); /* r = f - y */
6069 VpMult(f, VpPt5, r); /* f = 0.5*r */
6070 if (VpIsZero(f)) goto converge;
6071 VpAddSub(r, f, y, 1); /* r = y + f */
6072 VpAsgn(y, r, 1); /* y = r */
6073 } while (++nr < n);
6074
6075 #ifdef BIGDECIMAL_DEBUG
6076 if (gfDebug) {
6077 printf("ERROR(VpSqrt): did not converge within %ld iterations.\n", nr);
6078 }
6079 #endif /* BIGDECIMAL_DEBUG */
6080 y->MaxPrec = y_prec;
6081
6082 converge:
6083 VpChangeSign(y, 1);
6084 #ifdef BIGDECIMAL_DEBUG
6085 if (gfDebug) {
6086 VpMult(r, y, y);
6087 VpAddSub(f, x, r, -1);
6088 printf("VpSqrt: iterations = %"PRIdSIZE"\n", nr);
6089 VPrint(stdout, " y =% \n", y);
6090 VPrint(stdout, " x =% \n", x);
6091 VPrint(stdout, " x-y*y = % \n", f);
6092 }
6093 #endif /* BIGDECIMAL_DEBUG */
6094 y->MaxPrec = y_prec;
6095
6096 Exit:
6097 VpFree(f);
6098 VpFree(r);
6099 return 1;
6100 }
6101
6102 /*
6103 * Round relatively from the decimal point.
6104 * f: rounding mode
6105 * nf: digit location to round from the decimal point.
6106 */
6107 VP_EXPORT int
VpMidRound(Real * y,unsigned short f,ssize_t nf)6108 VpMidRound(Real *y, unsigned short f, ssize_t nf)
6109 {
6110 /* fracf: any positive digit under rounding position? */
6111 /* fracf_1further: any positive digits under one further than the rounding position? */
6112 /* exptoadd: number of digits needed to compensate negative nf */
6113 int fracf, fracf_1further;
6114 ssize_t n,i,ix,ioffset, exptoadd;
6115 BDIGIT v, shifter;
6116 BDIGIT div;
6117
6118 nf += y->exponent * (ssize_t)BASE_FIG;
6119 exptoadd=0;
6120 if (nf < 0) {
6121 /* rounding position too left(large). */
6122 if (f != VP_ROUND_CEIL && f != VP_ROUND_FLOOR) {
6123 VpSetZero(y, VpGetSign(y)); /* truncate everything */
6124 return 0;
6125 }
6126 exptoadd = -nf;
6127 nf = 0;
6128 }
6129
6130 ix = nf / (ssize_t)BASE_FIG;
6131 if ((size_t)ix >= y->Prec) return 0; /* rounding position too right(small). */
6132 v = y->frac[ix];
6133
6134 ioffset = nf - ix*(ssize_t)BASE_FIG;
6135 n = (ssize_t)BASE_FIG - ioffset - 1;
6136 for (shifter = 1, i = 0; i < n; ++i) shifter *= 10;
6137
6138 /* so the representation used (in y->frac) is an array of BDIGIT, where
6139 each BDIGIT contains a value between 0 and BASE-1, consisting of BASE_FIG
6140 decimal places.
6141
6142 (that numbers of decimal places are typed as ssize_t is somewhat confusing)
6143
6144 nf is now position (in decimal places) of the digit from the start of
6145 the array.
6146
6147 ix is the position (in BDIGITS) of the BDIGIT containing the decimal digit,
6148 from the start of the array.
6149
6150 v is the value of this BDIGIT
6151
6152 ioffset is the number of extra decimal places along of this decimal digit
6153 within v.
6154
6155 n is the number of decimal digits remaining within v after this decimal digit
6156 shifter is 10**n,
6157
6158 v % shifter are the remaining digits within v
6159 v % (shifter * 10) are the digit together with the remaining digits within v
6160 v / shifter are the digit's predecessors together with the digit
6161 div = v / shifter / 10 is just the digit's precessors
6162 (v / shifter) - div*10 is just the digit, which is what v ends up being reassigned to.
6163 */
6164
6165 fracf = (v % (shifter * 10) > 0);
6166 fracf_1further = ((v % shifter) > 0);
6167
6168 v /= shifter;
6169 div = v / 10;
6170 v = v - div*10;
6171 /* now v is just the digit required.
6172 now fracf is whether the digit or any of the remaining digits within v are non-zero
6173 now fracf_1further is whether any of the remaining digits within v are non-zero
6174 */
6175
6176 /* now check all the remaining BDIGITS for zero-ness a whole BDIGIT at a time.
6177 if we spot any non-zeroness, that means that we found a positive digit under
6178 rounding position, and we also found a positive digit under one further than
6179 the rounding position, so both searches (to see if any such non-zero digit exists)
6180 can stop */
6181
6182 for (i = ix + 1; (size_t)i < y->Prec; i++) {
6183 if (y->frac[i] % BASE) {
6184 fracf = fracf_1further = 1;
6185 break;
6186 }
6187 }
6188
6189 /* now fracf = does any positive digit exist under the rounding position?
6190 now fracf_1further = does any positive digit exist under one further than the
6191 rounding position?
6192 now v = the first digit under the rounding position */
6193
6194 /* drop digits after pointed digit */
6195 memset(y->frac + ix + 1, 0, (y->Prec - (ix + 1)) * sizeof(BDIGIT));
6196
6197 switch (f) {
6198 case VP_ROUND_DOWN: /* Truncate */
6199 break;
6200 case VP_ROUND_UP: /* Roundup */
6201 if (fracf) ++div;
6202 break;
6203 case VP_ROUND_HALF_UP:
6204 if (v>=5) ++div;
6205 break;
6206 case VP_ROUND_HALF_DOWN:
6207 if (v > 5 || (v == 5 && fracf_1further)) ++div;
6208 break;
6209 case VP_ROUND_CEIL:
6210 if (fracf && BIGDECIMAL_POSITIVE_P(y)) ++div;
6211 break;
6212 case VP_ROUND_FLOOR:
6213 if (fracf && BIGDECIMAL_NEGATIVE_P(y)) ++div;
6214 break;
6215 case VP_ROUND_HALF_EVEN: /* Banker's rounding */
6216 if (v > 5) ++div;
6217 else if (v == 5) {
6218 if (fracf_1further) {
6219 ++div;
6220 }
6221 else {
6222 if (ioffset == 0) {
6223 /* v is the first decimal digit of its BDIGIT;
6224 need to grab the previous BDIGIT if present
6225 to check for evenness of the previous decimal
6226 digit (which is same as that of the BDIGIT since
6227 base 10 has a factor of 2) */
6228 if (ix && (y->frac[ix-1] % 2)) ++div;
6229 }
6230 else {
6231 if (div % 2) ++div;
6232 }
6233 }
6234 }
6235 break;
6236 }
6237 for (i = 0; i <= n; ++i) div *= 10;
6238 if (div >= BASE) {
6239 if (ix) {
6240 y->frac[ix] = 0;
6241 VpRdup(y, ix);
6242 }
6243 else {
6244 short s = VpGetSign(y);
6245 SIGNED_VALUE e = y->exponent;
6246 VpSetOne(y);
6247 VpSetSign(y, s);
6248 y->exponent = e + 1;
6249 }
6250 }
6251 else {
6252 y->frac[ix] = div;
6253 VpNmlz(y);
6254 }
6255 if (exptoadd > 0) {
6256 y->exponent += (SIGNED_VALUE)(exptoadd / BASE_FIG);
6257 exptoadd %= (ssize_t)BASE_FIG;
6258 for (i = 0; i < exptoadd; i++) {
6259 y->frac[0] *= 10;
6260 if (y->frac[0] >= BASE) {
6261 y->frac[0] /= BASE;
6262 y->exponent++;
6263 }
6264 }
6265 }
6266 return 1;
6267 }
6268
6269 VP_EXPORT int
VpLeftRound(Real * y,unsigned short f,ssize_t nf)6270 VpLeftRound(Real *y, unsigned short f, ssize_t nf)
6271 /*
6272 * Round from the left hand side of the digits.
6273 */
6274 {
6275 BDIGIT v;
6276 if (!VpHasVal(y)) return 0; /* Unable to round */
6277 v = y->frac[0];
6278 nf -= VpExponent(y) * (ssize_t)BASE_FIG;
6279 while ((v /= 10) != 0) nf--;
6280 nf += (ssize_t)BASE_FIG-1;
6281 return VpMidRound(y, f, nf);
6282 }
6283
6284 VP_EXPORT int
VpActiveRound(Real * y,Real * x,unsigned short f,ssize_t nf)6285 VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t nf)
6286 {
6287 /* First,assign whole value in truncation mode */
6288 if (VpAsgn(y, x, 10) <= 1) return 0; /* Zero,NaN,or Infinity */
6289 return VpMidRound(y, f, nf);
6290 }
6291
6292 static int
VpLimitRound(Real * c,size_t ixDigit)6293 VpLimitRound(Real *c, size_t ixDigit)
6294 {
6295 size_t ix = VpGetPrecLimit();
6296 if (!VpNmlz(c)) return -1;
6297 if (!ix) return 0;
6298 if (!ixDigit) ixDigit = c->Prec-1;
6299 if ((ix + BASE_FIG - 1) / BASE_FIG > ixDigit + 1) return 0;
6300 return VpLeftRound(c, VpGetRoundMode(), (ssize_t)ix);
6301 }
6302
6303 /* If I understand correctly, this is only ever used to round off the final decimal
6304 digit of precision */
6305 static void
VpInternalRound(Real * c,size_t ixDigit,BDIGIT vPrev,BDIGIT v)6306 VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v)
6307 {
6308 int f = 0;
6309
6310 unsigned short const rounding_mode = VpGetRoundMode();
6311
6312 if (VpLimitRound(c, ixDigit)) return;
6313 if (!v) return;
6314
6315 v /= BASE1;
6316 switch (rounding_mode) {
6317 case VP_ROUND_DOWN:
6318 break;
6319 case VP_ROUND_UP:
6320 if (v) f = 1;
6321 break;
6322 case VP_ROUND_HALF_UP:
6323 if (v >= 5) f = 1;
6324 break;
6325 case VP_ROUND_HALF_DOWN:
6326 /* this is ok - because this is the last digit of precision,
6327 the case where v == 5 and some further digits are nonzero
6328 will never occur */
6329 if (v >= 6) f = 1;
6330 break;
6331 case VP_ROUND_CEIL:
6332 if (v && BIGDECIMAL_POSITIVE_P(c)) f = 1;
6333 break;
6334 case VP_ROUND_FLOOR:
6335 if (v && BIGDECIMAL_NEGATIVE_P(c)) f = 1;
6336 break;
6337 case VP_ROUND_HALF_EVEN: /* Banker's rounding */
6338 /* as per VP_ROUND_HALF_DOWN, because this is the last digit of precision,
6339 there is no case to worry about where v == 5 and some further digits are nonzero */
6340 if (v > 5) f = 1;
6341 else if (v == 5 && vPrev % 2) f = 1;
6342 break;
6343 }
6344 if (f) {
6345 VpRdup(c, ixDigit);
6346 VpNmlz(c);
6347 }
6348 }
6349
6350 /*
6351 * Rounds up m(plus one to final digit of m).
6352 */
6353 static int
VpRdup(Real * m,size_t ind_m)6354 VpRdup(Real *m, size_t ind_m)
6355 {
6356 BDIGIT carry;
6357
6358 if (!ind_m) ind_m = m->Prec;
6359
6360 carry = 1;
6361 while (carry > 0 && ind_m--) {
6362 m->frac[ind_m] += carry;
6363 if (m->frac[ind_m] >= BASE) m->frac[ind_m] -= BASE;
6364 else carry = 0;
6365 }
6366 if (carry > 0) { /* Overflow,count exponent and set fraction part be 1 */
6367 if (!AddExponent(m, 1)) return 0;
6368 m->Prec = m->frac[0] = 1;
6369 }
6370 else {
6371 VpNmlz(m);
6372 }
6373 return 1;
6374 }
6375
6376 /*
6377 * y = x - fix(x)
6378 */
6379 VP_EXPORT void
VpFrac(Real * y,Real * x)6380 VpFrac(Real *y, Real *x)
6381 {
6382 size_t my, ind_y, ind_x;
6383
6384 if (!VpHasVal(x)) {
6385 VpAsgn(y, x, 1);
6386 goto Exit;
6387 }
6388
6389 if (x->exponent > 0 && (size_t)x->exponent >= x->Prec) {
6390 VpSetZero(y, VpGetSign(x));
6391 goto Exit;
6392 }
6393 else if (x->exponent <= 0) {
6394 VpAsgn(y, x, 1);
6395 goto Exit;
6396 }
6397
6398 /* satisfy: x->exponent > 0 */
6399
6400 y->Prec = x->Prec - (size_t)x->exponent;
6401 y->Prec = Min(y->Prec, y->MaxPrec);
6402 y->exponent = 0;
6403 VpSetSign(y, VpGetSign(x));
6404 ind_y = 0;
6405 my = y->Prec;
6406 ind_x = x->exponent;
6407 while (ind_y < my) {
6408 y->frac[ind_y] = x->frac[ind_x];
6409 ++ind_y;
6410 ++ind_x;
6411 }
6412 VpNmlz(y);
6413
6414 Exit:
6415 #ifdef BIGDECIMAL_DEBUG
6416 if (gfDebug) {
6417 VPrint(stdout, "VpFrac y=%\n", y);
6418 VPrint(stdout, " x=%\n", x);
6419 }
6420 #endif /* BIGDECIMAL_DEBUG */
6421 return;
6422 }
6423
6424 /*
6425 * y = x ** n
6426 */
6427 VP_EXPORT int
VpPower(Real * y,Real * x,SIGNED_VALUE n)6428 VpPower(Real *y, Real *x, SIGNED_VALUE n)
6429 {
6430 size_t s, ss;
6431 ssize_t sign;
6432 Real *w1 = NULL;
6433 Real *w2 = NULL;
6434
6435 if (VpIsZero(x)) {
6436 if (n == 0) {
6437 VpSetOne(y);
6438 goto Exit;
6439 }
6440 sign = VpGetSign(x);
6441 if (n < 0) {
6442 n = -n;
6443 if (sign < 0) sign = (n % 2) ? -1 : 1;
6444 VpSetInf(y, sign);
6445 }
6446 else {
6447 if (sign < 0) sign = (n % 2) ? -1 : 1;
6448 VpSetZero(y,sign);
6449 }
6450 goto Exit;
6451 }
6452 if (VpIsNaN(x)) {
6453 VpSetNaN(y);
6454 goto Exit;
6455 }
6456 if (VpIsInf(x)) {
6457 if (n == 0) {
6458 VpSetOne(y);
6459 goto Exit;
6460 }
6461 if (n > 0) {
6462 VpSetInf(y, (n % 2 == 0 || VpIsPosInf(x)) ? 1 : -1);
6463 goto Exit;
6464 }
6465 VpSetZero(y, (n % 2 == 0 || VpIsPosInf(x)) ? 1 : -1);
6466 goto Exit;
6467 }
6468
6469 if (x->exponent == 1 && x->Prec == 1 && x->frac[0] == 1) {
6470 /* abs(x) = 1 */
6471 VpSetOne(y);
6472 if (BIGDECIMAL_POSITIVE_P(x)) goto Exit;
6473 if ((n % 2) == 0) goto Exit;
6474 VpSetSign(y, -1);
6475 goto Exit;
6476 }
6477
6478 if (n > 0) sign = 1;
6479 else if (n < 0) {
6480 sign = -1;
6481 n = -n;
6482 }
6483 else {
6484 VpSetOne(y);
6485 goto Exit;
6486 }
6487
6488 /* Allocate working variables */
6489
6490 w1 = VpAlloc((y->MaxPrec + 2) * BASE_FIG, "#0", 1, 1);
6491 w2 = VpAlloc((w1->MaxPrec * 2 + 1) * BASE_FIG, "#0", 1, 1);
6492 /* calculation start */
6493
6494 VpAsgn(y, x, 1);
6495 --n;
6496 while (n > 0) {
6497 VpAsgn(w1, x, 1);
6498 s = 1;
6499 while (ss = s, (s += s) <= (size_t)n) {
6500 VpMult(w2, w1, w1);
6501 VpAsgn(w1, w2, 1);
6502 }
6503 n -= (SIGNED_VALUE)ss;
6504 VpMult(w2, y, w1);
6505 VpAsgn(y, w2, 1);
6506 }
6507 if (sign < 0) {
6508 VpDivd(w1, w2, VpConstOne, y);
6509 VpAsgn(y, w1, 1);
6510 }
6511
6512 Exit:
6513 #ifdef BIGDECIMAL_DEBUG
6514 if (gfDebug) {
6515 VPrint(stdout, "VpPower y=%\n", y);
6516 VPrint(stdout, "VpPower x=%\n", x);
6517 printf(" n=%"PRIdVALUE"\n", n);
6518 }
6519 #endif /* BIGDECIMAL_DEBUG */
6520 VpFree(w2);
6521 VpFree(w1);
6522 return 1;
6523 }
6524
6525 #ifdef BIGDECIMAL_DEBUG
6526 int
VpVarCheck(Real * v)6527 VpVarCheck(Real * v)
6528 /*
6529 * Checks the validity of the Real variable v.
6530 * [Input]
6531 * v ... Real *, variable to be checked.
6532 * [Returns]
6533 * 0 ... correct v.
6534 * other ... error
6535 */
6536 {
6537 size_t i;
6538
6539 if (v->MaxPrec == 0) {
6540 printf("ERROR(VpVarCheck): Illegal Max. Precision(=%"PRIuSIZE")\n",
6541 v->MaxPrec);
6542 return 1;
6543 }
6544 if (v->Prec == 0 || v->Prec > v->MaxPrec) {
6545 printf("ERROR(VpVarCheck): Illegal Precision(=%"PRIuSIZE")\n", v->Prec);
6546 printf(" Max. Prec.=%"PRIuSIZE"\n", v->MaxPrec);
6547 return 2;
6548 }
6549 for (i = 0; i < v->Prec; ++i) {
6550 if (v->frac[i] >= BASE) {
6551 printf("ERROR(VpVarCheck): Illegal fraction\n");
6552 printf(" Frac[%"PRIuSIZE"]=%"PRIuBDIGIT"\n", i, v->frac[i]);
6553 printf(" Prec. =%"PRIuSIZE"\n", v->Prec);
6554 printf(" Exp. =%"PRIdVALUE"\n", v->exponent);
6555 printf(" BASE =%"PRIuBDIGIT"\n", BASE);
6556 return 3;
6557 }
6558 }
6559 return 0;
6560 }
6561 #endif /* BIGDECIMAL_DEBUG */
6562