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