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