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