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