Ruby  2.0.0p645(2015-04-13revision50299)
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
1192 BigDecimal_divide(Real **c, Real **res, Real **div, VALUE self, VALUE r)
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  ENTER(1);
2406  Real *pv = rb_check_typeddata(self, &BigDecimal_data_type);
2407  Real *x;
2408 
2409  GUARD_OBJ(x, BigDecimal_new(argc, argv));
2410  if (ToValue(x)) {
2411  pv = VpCopy(pv, x);
2412  }
2413  else {
2414  VpFree(pv);
2415  pv = x;
2416  }
2417  DATA_PTR(self) = pv;
2418  pv->obj = self;
2419  return self;
2420 }
2421 
2422 /* :nodoc:
2423  *
2424  * private method to dup and clone the provided BigDecimal +other+
2425  */
2426 static VALUE
2428 {
2429  Real *pv = rb_check_typeddata(self, &BigDecimal_data_type);
2430  Real *x = rb_check_typeddata(other, &BigDecimal_data_type);
2431 
2432  if (self != other) {
2433  DATA_PTR(self) = VpCopy(pv, x);
2434  }
2435  return self;
2436 }
2437 
2438 static Real *
2440 {
2441  size_t mf;
2442  VALUE nFig;
2443  VALUE iniValue;
2444 
2445  if (rb_scan_args(argc, argv, "11", &iniValue, &nFig) == 1) {
2446  mf = 0;
2447  }
2448  else {
2449  mf = GetPositiveInt(nFig);
2450  }
2451 
2452  switch (TYPE(iniValue)) {
2453  case T_DATA:
2454  if (is_kind_of_BigDecimal(iniValue)) {
2455  return DATA_PTR(iniValue);
2456  }
2457  break;
2458 
2459  case T_FIXNUM:
2460  /* fall through */
2461  case T_BIGNUM:
2462  return GetVpValue(iniValue, 1);
2463 
2464  case T_FLOAT:
2465  if (mf > DBL_DIG+1) {
2466  rb_raise(rb_eArgError, "precision too large.");
2467  }
2468  /* fall through */
2469  case T_RATIONAL:
2470  if (NIL_P(nFig)) {
2472  "can't omit precision for a %"PRIsVALUE".",
2473  RB_OBJ_CLASSNAME(iniValue));
2474  }
2475  return GetVpValueWithPrec(iniValue, mf, 1);
2476 
2477  case T_STRING:
2478  /* fall through */
2479  default:
2480  break;
2481  }
2482  StringValueCStr(iniValue);
2483  return VpAlloc(mf, RSTRING_PTR(iniValue));
2484 }
2485 
2486 static VALUE
2488 {
2489  ENTER(1);
2490  Real *pv;
2491 
2492  GUARD_OBJ(pv, BigDecimal_new(argc, argv));
2493  if (ToValue(pv)) pv = VpCopy(NULL, pv);
2494  pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv);
2495  return pv->obj;
2496 }
2497 
2498  /* call-seq:
2499  * BigDecimal.limit(digits)
2500  *
2501  * Limit the number of significant digits in newly created BigDecimal
2502  * numbers to the specified value. Rounding is performed as necessary,
2503  * as specified by BigDecimal.mode.
2504  *
2505  * A limit of 0, the default, means no upper limit.
2506  *
2507  * The limit specified by this method takes less priority over any limit
2508  * specified to instance methods such as ceil, floor, truncate, or round.
2509  */
2510 static VALUE
2512 {
2513  VALUE nFig;
2514  VALUE nCur = INT2NUM(VpGetPrecLimit());
2515 
2516  if(rb_scan_args(argc,argv,"01",&nFig)==1) {
2517  int nf;
2518  if(nFig==Qnil) return nCur;
2519  Check_Type(nFig, T_FIXNUM);
2520  nf = FIX2INT(nFig);
2521  if(nf<0) {
2522  rb_raise(rb_eArgError, "argument must be positive");
2523  }
2524  VpSetPrecLimit(nf);
2525  }
2526  return nCur;
2527 }
2528 
2529 /* Returns the sign of the value.
2530  *
2531  * Returns a positive value if > 0, a negative value if < 0, and a
2532  * zero if == 0.
2533  *
2534  * The specific value returned indicates the type and sign of the BigDecimal,
2535  * as follows:
2536  *
2537  * BigDecimal::SIGN_NaN:: value is Not a Number
2538  * BigDecimal::SIGN_POSITIVE_ZERO:: value is +0
2539  * BigDecimal::SIGN_NEGATIVE_ZERO:: value is -0
2540  * BigDecimal::SIGN_POSITIVE_INFINITE:: value is +infinity
2541  * BigDecimal::SIGN_NEGATIVE_INFINITE:: value is -infinity
2542  * BigDecimal::SIGN_POSITIVE_FINITE:: value is positive
2543  * BigDecimal::SIGN_NEGATIVE_FINITE:: value is negative
2544  */
2545 static VALUE
2547 { /* sign */
2548  int s = GetVpValue(self,1)->sign;
2549  return INT2FIX(s);
2550 }
2551 
2552 /*
2553  * call-seq: BigDecimal.save_exception_mode { ... }
2554  *
2555  * Excecute the provided block, but preserve the exception mode
2556  *
2557  * BigDecimal.save_exception_mode do
2558  * BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false)
2559  * BigDecimal.mode(BigDecimal::EXCEPTION_NaN, false)
2560  *
2561  * BigDecimal.new(BigDecimal('Infinity'))
2562  * BigDecimal.new(BigDecimal('-Infinity'))
2563  * BigDecimal(BigDecimal.new('NaN'))
2564  * end
2565  *
2566  * For use with the BigDecimal::EXCEPTION_*
2567  *
2568  * See BigDecimal.mode
2569  */
2570 static VALUE
2572 {
2573  unsigned short const exception_mode = VpGetException();
2574  int state;
2575  VALUE ret = rb_protect(rb_yield, Qnil, &state);
2576  VpSetException(exception_mode);
2577  if (state) rb_jump_tag(state);
2578  return ret;
2579 }
2580 
2581 /*
2582  * call-seq: BigDecimal.save_rounding_mode { ... }
2583  *
2584  * Excecute the provided block, but preserve the rounding mode
2585  *
2586  * BigDecimal.save_exception_mode do
2587  * BigDecimal.mode(BigDecimal::ROUND_MODE, :up)
2588  * puts BigDecimal.mode(BigDecimal::ROUND_MODE)
2589  * end
2590  *
2591  * For use with the BigDecimal::ROUND_*
2592  *
2593  * See BigDecimal.mode
2594  */
2595 static VALUE
2597 {
2598  unsigned short const round_mode = VpGetRoundMode();
2599  int state;
2600  VALUE ret = rb_protect(rb_yield, Qnil, &state);
2601  VpSetRoundMode(round_mode);
2602  if (state) rb_jump_tag(state);
2603  return ret;
2604 }
2605 
2606 /*
2607  * call-seq: BigDecimal.save_limit { ... }
2608  *
2609  * Excecute the provided block, but preserve the precision limit
2610  *
2611  * BigDecimal.limit(100)
2612  * puts BigDecimal.limit
2613  * BigDecimal.save_limit do
2614  * BigDecimal.limit(200)
2615  * puts BigDecimal.limit
2616  * end
2617  * puts BigDecimal.limit
2618  *
2619  */
2620 static VALUE
2622 {
2623  size_t const limit = VpGetPrecLimit();
2624  int state;
2625  VALUE ret = rb_protect(rb_yield, Qnil, &state);
2626  VpSetPrecLimit(limit);
2627  if (state) rb_jump_tag(state);
2628  return ret;
2629 }
2630 
2631 /* call-seq:
2632  * BigMath.exp(x, prec)
2633  *
2634  * Computes the value of e (the base of natural logarithms) raised to the
2635  * power of x, to the specified number of digits of precision.
2636  *
2637  * If x is infinity, returns Infinity.
2638  *
2639  * If x is NaN, returns NaN.
2640  */
2641 static VALUE
2642 BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
2643 {
2644  ssize_t prec, n, i;
2645  Real* vx = NULL;
2646  VALUE one, d, x1, y, z;
2647  int negative = 0;
2648  int infinite = 0;
2649  int nan = 0;
2650  double flo;
2651 
2652  prec = NUM2SSIZET(vprec);
2653  if (prec <= 0) {
2654  rb_raise(rb_eArgError, "Zero or negative precision for exp");
2655  }
2656 
2657  /* TODO: the following switch statement is almostly the same as one in the
2658  * BigDecimalCmp function. */
2659  switch (TYPE(x)) {
2660  case T_DATA:
2661  if (!is_kind_of_BigDecimal(x)) break;
2662  vx = DATA_PTR(x);
2663  negative = VpGetSign(vx) < 0;
2664  infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
2665  nan = VpIsNaN(vx);
2666  break;
2667 
2668  case T_FIXNUM:
2669  /* fall through */
2670  case T_BIGNUM:
2671  vx = GetVpValue(x, 0);
2672  break;
2673 
2674  case T_FLOAT:
2675  flo = RFLOAT_VALUE(x);
2676  negative = flo < 0;
2677  infinite = isinf(flo);
2678  nan = isnan(flo);
2679  if (!infinite && !nan) {
2680  vx = GetVpValueWithPrec(x, DBL_DIG+1, 0);
2681  }
2682  break;
2683 
2684  case T_RATIONAL:
2685  vx = GetVpValueWithPrec(x, prec, 0);
2686  break;
2687 
2688  default:
2689  break;
2690  }
2691  if (infinite) {
2692  if (negative) {
2693  return ToValue(GetVpValueWithPrec(INT2NUM(0), prec, 1));
2694  }
2695  else {
2696  Real* vy;
2697  vy = VpCreateRbObject(prec, "#0");
2698  RB_GC_GUARD(vy->obj);
2700  return ToValue(vy);
2701  }
2702  }
2703  else if (nan) {
2704  Real* vy;
2705  vy = VpCreateRbObject(prec, "#0");
2706  RB_GC_GUARD(vy->obj);
2707  VpSetNaN(vy);
2708  return ToValue(vy);
2709  }
2710  else if (vx == NULL) {
2712  }
2713  x = RB_GC_GUARD(vx->obj);
2714 
2715  n = prec + rmpd_double_figures();
2716  negative = VpGetSign(vx) < 0;
2717  if (negative) {
2718  VpSetSign(vx, 1);
2719  }
2720 
2721  RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
2722  RB_GC_GUARD(x1) = one;
2723  RB_GC_GUARD(y) = one;
2724  RB_GC_GUARD(d) = y;
2725  RB_GC_GUARD(z) = one;
2726  i = 0;
2727 
2728  while (!VpIsZero((Real*)DATA_PTR(d))) {
2729  VALUE argv[2];
2730  SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
2731  SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
2732  ssize_t m = n - vabs(ey - ed);
2733  if (m <= 0) {
2734  break;
2735  }
2736  else if ((size_t)m < rmpd_double_figures()) {
2737  m = rmpd_double_figures();
2738  }
2739 
2740  x1 = BigDecimal_mult2(x1, x, SSIZET2NUM(n));
2741  ++i;
2742  z = BigDecimal_mult(z, SSIZET2NUM(i));
2743  argv[0] = z;
2744  argv[1] = SSIZET2NUM(m);
2745  d = BigDecimal_div2(2, argv, x1);
2746  y = BigDecimal_add(y, d);
2747  }
2748 
2749  if (negative) {
2750  VALUE argv[2];
2751  argv[0] = y;
2752  argv[1] = vprec;
2753  return BigDecimal_div2(2, argv, one);
2754  }
2755  else {
2756  vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
2757  return BigDecimal_round(1, &vprec, y);
2758  }
2759 }
2760 
2761 /* call-seq:
2762  * BigMath.log(x, prec)
2763  *
2764  * Computes the natural logarithm of x to the specified number of digits of
2765  * precision.
2766  *
2767  * If x is zero or negative, raises Math::DomainError.
2768  *
2769  * If x is positive infinity, returns Infinity.
2770  *
2771  * If x is NaN, returns NaN.
2772  */
2773 static VALUE
2774 BigMath_s_log(VALUE klass, VALUE x, VALUE vprec)
2775 {
2776  ssize_t prec, n, i;
2777  SIGNED_VALUE expo;
2778  Real* vx = NULL;
2779  VALUE argv[2], vn, one, two, w, x2, y, d;
2780  int zero = 0;
2781  int negative = 0;
2782  int infinite = 0;
2783  int nan = 0;
2784  double flo;
2785  long fix;
2786 
2787  if (!is_integer(vprec)) {
2788  rb_raise(rb_eArgError, "precision must be an Integer");
2789  }
2790 
2791  prec = NUM2SSIZET(vprec);
2792  if (prec <= 0) {
2793  rb_raise(rb_eArgError, "Zero or negative precision for exp");
2794  }
2795 
2796  /* TODO: the following switch statement is almostly the same as one in the
2797  * BigDecimalCmp function. */
2798  switch (TYPE(x)) {
2799  case T_DATA:
2800  if (!is_kind_of_BigDecimal(x)) break;
2801  vx = DATA_PTR(x);
2802  zero = VpIsZero(vx);
2803  negative = VpGetSign(vx) < 0;
2804  infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
2805  nan = VpIsNaN(vx);
2806  break;
2807 
2808  case T_FIXNUM:
2809  fix = FIX2LONG(x);
2810  zero = fix == 0;
2811  negative = fix < 0;
2812  goto get_vp_value;
2813 
2814  case T_BIGNUM:
2815  zero = RBIGNUM_ZERO_P(x);
2816  negative = RBIGNUM_NEGATIVE_P(x);
2817 get_vp_value:
2818  if (zero || negative) break;
2819  vx = GetVpValue(x, 0);
2820  break;
2821 
2822  case T_FLOAT:
2823  flo = RFLOAT_VALUE(x);
2824  zero = flo == 0;
2825  negative = flo < 0;
2826  infinite = isinf(flo);
2827  nan = isnan(flo);
2828  if (!zero && !negative && !infinite && !nan) {
2829  vx = GetVpValueWithPrec(x, DBL_DIG+1, 1);
2830  }
2831  break;
2832 
2833  case T_RATIONAL:
2834  zero = RRATIONAL_ZERO_P(x);
2835  negative = RRATIONAL_NEGATIVE_P(x);
2836  if (zero || negative) break;
2837  vx = GetVpValueWithPrec(x, prec, 1);
2838  break;
2839 
2840  case T_COMPLEX:
2842  "Complex argument for BigMath.log");
2843 
2844  default:
2845  break;
2846  }
2847  if (infinite && !negative) {
2848  Real* vy;
2849  vy = VpCreateRbObject(prec, "#0");
2850  RB_GC_GUARD(vy->obj);
2852  return ToValue(vy);
2853  }
2854  else if (nan) {
2855  Real* vy;
2856  vy = VpCreateRbObject(prec, "#0");
2857  RB_GC_GUARD(vy->obj);
2858  VpSetNaN(vy);
2859  return ToValue(vy);
2860  }
2861  else if (zero || negative) {
2863  "Zero or negative argument for log");
2864  }
2865  else if (vx == NULL) {
2867  }
2868  x = ToValue(vx);
2869 
2870  RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
2871  RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2"));
2872 
2873  n = prec + rmpd_double_figures();
2874  RB_GC_GUARD(vn) = SSIZET2NUM(n);
2875  expo = VpExponent10(vx);
2876  if (expo < 0 || expo >= 3) {
2877  char buf[16];
2878  snprintf(buf, 16, "1E%"PRIdVALUE, -expo);
2879  x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
2880  }
2881  else {
2882  expo = 0;
2883  }
2884  w = BigDecimal_sub(x, one);
2885  argv[0] = BigDecimal_add(x, one);
2886  argv[1] = vn;
2887  x = BigDecimal_div2(2, argv, w);
2888  RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn);
2889  RB_GC_GUARD(y) = x;
2890  RB_GC_GUARD(d) = y;
2891  i = 1;
2892  while (!VpIsZero((Real*)DATA_PTR(d))) {
2893  SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
2894  SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
2895  ssize_t m = n - vabs(ey - ed);
2896  if (m <= 0) {
2897  break;
2898  }
2899  else if ((size_t)m < rmpd_double_figures()) {
2900  m = rmpd_double_figures();
2901  }
2902 
2903  x = BigDecimal_mult2(x2, x, vn);
2904  i += 2;
2905  argv[0] = SSIZET2NUM(i);
2906  argv[1] = SSIZET2NUM(m);
2907  d = BigDecimal_div2(2, argv, x);
2908  y = BigDecimal_add(y, d);
2909  }
2910 
2911  y = BigDecimal_mult(y, two);
2912  if (expo != 0) {
2913  VALUE log10, vexpo, dy;
2914  log10 = BigMath_s_log(klass, INT2FIX(10), vprec);
2915  vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1));
2916  dy = BigDecimal_mult(log10, vexpo);
2917  y = BigDecimal_add(y, dy);
2918  }
2919 
2920  return y;
2921 }
2922 
2923 /* Document-class: BigDecimal
2924  * BigDecimal provides arbitrary-precision floating point decimal arithmetic.
2925  *
2926  * Copyright (C) 2002 by Shigeo Kobayashi <shigeo@tinyforest.gr.jp>.
2927  *
2928  * You may distribute under the terms of either the GNU General Public
2929  * License or the Artistic License, as specified in the README file
2930  * of the BigDecimal distribution.
2931  *
2932  * Documented by mathew <meta@pobox.com>.
2933  *
2934  * = Introduction
2935  *
2936  * Ruby provides built-in support for arbitrary precision integer arithmetic.
2937  *
2938  * For example:
2939  *
2940  * 42**13 #=> 1265437718438866624512
2941  *
2942  * BigDecimal provides similar support for very large or very accurate floating
2943  * point numbers.
2944  *
2945  * Decimal arithmetic is also useful for general calculation, because it
2946  * provides the correct answers people expect--whereas normal binary floating
2947  * point arithmetic often introduces subtle errors because of the conversion
2948  * between base 10 and base 2.
2949  *
2950  * For example, try:
2951  *
2952  * sum = 0
2953  * for i in (1..10000)
2954  * sum = sum + 0.0001
2955  * end
2956  * print sum #=> 0.9999999999999062
2957  *
2958  * and contrast with the output from:
2959  *
2960  * require 'bigdecimal'
2961  *
2962  * sum = BigDecimal.new("0")
2963  * for i in (1..10000)
2964  * sum = sum + BigDecimal.new("0.0001")
2965  * end
2966  * print sum #=> 0.1E1
2967  *
2968  * Similarly:
2969  *
2970  * (BigDecimal.new("1.2") - BigDecimal("1.0")) == BigDecimal("0.2") #=> true
2971  *
2972  * (1.2 - 1.0) == 0.2 #=> false
2973  *
2974  * = Special features of accurate decimal arithmetic
2975  *
2976  * Because BigDecimal is more accurate than normal binary floating point
2977  * arithmetic, it requires some special values.
2978  *
2979  * == Infinity
2980  *
2981  * BigDecimal sometimes needs to return infinity, for example if you divide
2982  * a value by zero.
2983  *
2984  * BigDecimal.new("1.0") / BigDecimal.new("0.0") #=> infinity
2985  * BigDecimal.new("-1.0") / BigDecimal.new("0.0") #=> -infinity
2986  *
2987  * You can represent infinite numbers to BigDecimal using the strings
2988  * <code>'Infinity'</code>, <code>'+Infinity'</code> and
2989  * <code>'-Infinity'</code> (case-sensitive)
2990  *
2991  * == Not a Number
2992  *
2993  * When a computation results in an undefined value, the special value +NaN+
2994  * (for 'not a number') is returned.
2995  *
2996  * Example:
2997  *
2998  * BigDecimal.new("0.0") / BigDecimal.new("0.0") #=> NaN
2999  *
3000  * You can also create undefined values.
3001  *
3002  * NaN is never considered to be the same as any other value, even NaN itself:
3003  *
3004  * n = BigDecimal.new('NaN')
3005  * n == 0.0 #=> nil
3006  * n == n #=> nil
3007  *
3008  * == Positive and negative zero
3009  *
3010  * If a computation results in a value which is too small to be represented as
3011  * a BigDecimal within the currently specified limits of precision, zero must
3012  * be returned.
3013  *
3014  * If the value which is too small to be represented is negative, a BigDecimal
3015  * value of negative zero is returned.
3016  *
3017  * BigDecimal.new("1.0") / BigDecimal.new("-Infinity") #=> -0.0
3018  *
3019  * If the value is positive, a value of positive zero is returned.
3020  *
3021  * BigDecimal.new("1.0") / BigDecimal.new("Infinity") #=> 0.0
3022  *
3023  * (See BigDecimal.mode for how to specify limits of precision.)
3024  *
3025  * Note that +-0.0+ and +0.0+ are considered to be the same for the purposes of
3026  * comparison.
3027  *
3028  * Note also that in mathematics, there is no particular concept of negative
3029  * or positive zero; true mathematical zero has no sign.
3030  */
3031 void
3033 {
3034  VALUE arg;
3035 
3036  id_BigDecimal_exception_mode = rb_intern_const("BigDecimal.exception_mode");
3037  id_BigDecimal_rounding_mode = rb_intern_const("BigDecimal.rounding_mode");
3038  id_BigDecimal_precision_limit = rb_intern_const("BigDecimal.precision_limit");
3039 
3040  /* Initialize VP routines */
3041  VpInit(0UL);
3042 
3043  /* Class and method registration */
3044  rb_cBigDecimal = rb_define_class("BigDecimal",rb_cNumeric);
3046 
3047  /* Global function */
3049 
3050  /* Class methods */
3056 
3060 
3061  /* Constants definition */
3062 
3063  /*
3064  * Base value used in internal calculations. On a 32 bit system, BASE
3065  * is 10000, indicating that calculation is done in groups of 4 digits.
3066  * (If it were larger, BASE**2 wouldn't fit in 32 bits, so you couldn't
3067  * guarantee that two groups could always be multiplied together without
3068  * overflow.)
3069  */
3071 
3072  /* Exceptions */
3073 
3074  /*
3075  * 0xff: Determines whether overflow, underflow or zero divide result in
3076  * an exception being thrown. See BigDecimal.mode.
3077  */
3079 
3080  /*
3081  * 0x02: Determines what happens when the result of a computation is not a
3082  * number (NaN). See BigDecimal.mode.
3083  */
3085 
3086  /*
3087  * 0x01: Determines what happens when the result of a computation is
3088  * infinity. See BigDecimal.mode.
3089  */
3091 
3092  /*
3093  * 0x04: Determines what happens when the result of a computation is an
3094  * underflow (a result too small to be represented). See BigDecimal.mode.
3095  */
3097 
3098  /*
3099  * 0x01: Determines what happens when the result of a computation is an
3100  * overflow (a result too large to be represented). See BigDecimal.mode.
3101  */
3103 
3104  /*
3105  * 0x01: Determines what happens when a division by zero is performed.
3106  * See BigDecimal.mode.
3107  */
3109 
3110  /*
3111  * 0x100: Determines what happens when a result must be rounded in order to
3112  * fit in the appropriate number of significant digits. See
3113  * BigDecimal.mode.
3114  */
3116 
3117  /* 1: Indicates that values should be rounded away from zero. See
3118  * BigDecimal.mode.
3119  */
3121 
3122  /* 2: Indicates that values should be rounded towards zero. See
3123  * BigDecimal.mode.
3124  */
3126 
3127  /* 3: Indicates that digits >= 5 should be rounded up, others rounded down.
3128  * See BigDecimal.mode. */
3130 
3131  /* 4: Indicates that digits >= 6 should be rounded up, others rounded down.
3132  * See BigDecimal.mode.
3133  */
3135  /* 5: Round towards +infinity. See BigDecimal.mode. */
3137 
3138  /* 6: Round towards -infinity. See BigDecimal.mode. */
3140 
3141  /* 7: Round towards the even neighbor. See BigDecimal.mode. */
3143 
3144  /* 0: Indicates that a value is not a number. See BigDecimal.sign. */
3146 
3147  /* 1: Indicates that a value is +0. See BigDecimal.sign. */
3149 
3150  /* -1: Indicates that a value is -0. See BigDecimal.sign. */
3152 
3153  /* 2: Indicates that a value is positive and finite. See BigDecimal.sign. */
3155 
3156  /* -2: Indicates that a value is negative and finite. See BigDecimal.sign. */
3158 
3159  /* 3: Indicates that a value is positive and infinite. See BigDecimal.sign. */
3161 
3162  /* -3: Indicates that a value is negative and infinite. See BigDecimal.sign. */
3164 
3165  arg = rb_str_new2("+Infinity");
3166  /* Positive infinity value. */
3168  arg = rb_str_new2("NaN");
3169  /* 'Not a Number' value. */
3171 
3172 
3173  /* instance methods */
3177 
3199  /* rb_define_method(rb_cBigDecimal, "dup", BigDecimal_dup, 0); */
3229 
3230  /* mathematical functions */
3231  rb_mBigMath = rb_define_module("BigMath");
3234 
3235  id_up = rb_intern_const("up");
3236  id_down = rb_intern_const("down");
3237  id_truncate = rb_intern_const("truncate");
3238  id_half_up = rb_intern_const("half_up");
3239  id_default = rb_intern_const("default");
3240  id_half_down = rb_intern_const("half_down");
3241  id_half_even = rb_intern_const("half_even");
3242  id_banker = rb_intern_const("banker");
3243  id_ceiling = rb_intern_const("ceiling");
3244  id_ceil = rb_intern_const("ceil");
3245  id_floor = rb_intern_const("floor");
3246  id_to_r = rb_intern_const("to_r");
3247  id_eq = rb_intern_const("==");
3248 }
3249 
3250 /*
3251  *
3252  * ============================================================================
3253  *
3254  * vp_ routines begin from here.
3255  *
3256  * ============================================================================
3257  *
3258  */
3259 #ifdef BIGDECIMAL_DEBUG
3260 static int gfDebug = 1; /* Debug switch */
3261 #if 0
3262 static int gfCheckVal = 1; /* Value checking flag in VpNmlz() */
3263 #endif
3264 #endif /* BIGDECIMAL_DEBUG */
3265 
3266 static Real *VpConstOne; /* constant 1.0 */
3267 static Real *VpPt5; /* constant 0.5 */
3268 #define maxnr 100UL /* Maximum iterations for calcurating sqrt. */
3269  /* used in VpSqrt() */
3270 
3271 /* ETC */
3272 #define MemCmp(x,y,z) memcmp(x,y,z)
3273 #define StrCmp(x,y) strcmp(x,y)
3274 
3275 static int VpIsDefOP(Real *c,Real *a,Real *b,int sw);
3276 static int AddExponent(Real *a, SIGNED_VALUE n);
3277 static BDIGIT VpAddAbs(Real *a,Real *b,Real *c);
3278 static BDIGIT VpSubAbs(Real *a,Real *b,Real *c);
3279 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);
3280 static int VpNmlz(Real *a);
3281 static void VpFormatSt(char *psz, size_t fFmt);
3282 static int VpRdup(Real *m, size_t ind_m);
3283 
3284 #ifdef BIGDECIMAL_DEBUG
3285 static int gnAlloc=0; /* Memory allocation counter */
3286 #endif /* BIGDECIMAL_DEBUG */
3287 
3288 VP_EXPORT void *
3289 VpMemAlloc(size_t mb)
3290 {
3291  void *p = xmalloc(mb);
3292  if (!p) {
3293  VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
3294  }
3295  memset(p, 0, mb);
3296 #ifdef BIGDECIMAL_DEBUG
3297  gnAlloc++; /* Count allocation call */
3298 #endif /* BIGDECIMAL_DEBUG */
3299  return p;
3300 }
3301 
3302 VP_EXPORT void *
3303 VpMemRealloc(void *ptr, size_t mb)
3304 {
3305  void *p = xrealloc(ptr, mb);
3306  if (!p) {
3307  VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
3308  }
3309  return p;
3310 }
3311 
3312 VP_EXPORT void
3314 {
3315  if(pv != NULL) {
3316  xfree(pv);
3317 #ifdef BIGDECIMAL_DEBUG
3318  gnAlloc--; /* Decrement allocation count */
3319  if(gnAlloc==0) {
3320  printf(" *************** All memories allocated freed ****************");
3321  getchar();
3322  }
3323  if(gnAlloc<0) {
3324  printf(" ??????????? Too many memory free calls(%d) ?????????????\n",gnAlloc);
3325  getchar();
3326  }
3327 #endif /* BIGDECIMAL_DEBUG */
3328  }
3329 }
3330 
3331 /*
3332  * EXCEPTION Handling.
3333  */
3334 
3335 #define rmpd_set_thread_local_exception_mode(mode) \
3336  rb_thread_local_aset( \
3337  rb_thread_current(), \
3338  id_BigDecimal_exception_mode, \
3339  INT2FIX((int)(mode)) \
3340  )
3341 
3342 static unsigned short
3344 {
3345  VALUE const vmode = rb_thread_local_aref(
3348  );
3349 
3350  if (NIL_P(vmode)) {
3353  }
3354 
3355  return (unsigned short)FIX2UINT(vmode);
3356 }
3357 
3358 static void
3359 VpSetException(unsigned short f)
3360 {
3362 }
3363 
3364 /*
3365  * Precision limit.
3366  */
3367 
3368 #define rmpd_set_thread_local_precision_limit(limit) \
3369  rb_thread_local_aset( \
3370  rb_thread_current(), \
3371  id_BigDecimal_precision_limit, \
3372  SIZET2NUM(limit) \
3373  )
3374 #define RMPD_PRECISION_LIMIT_DEFAULT ((size_t)0)
3375 
3376 /* These 2 functions added at v1.1.7 */
3377 VP_EXPORT size_t
3379 {
3380  VALUE const vlimit = rb_thread_local_aref(
3383  );
3384 
3385  if (NIL_P(vlimit)) {
3388  }
3389 
3390  return NUM2SIZET(vlimit);
3391 }
3392 
3393 VP_EXPORT size_t
3395 {
3396  size_t const s = VpGetPrecLimit();
3398  return s;
3399 }
3400 
3401 /*
3402  * Rounding mode.
3403  */
3404 
3405 #define rmpd_set_thread_local_rounding_mode(mode) \
3406  rb_thread_local_aset( \
3407  rb_thread_current(), \
3408  id_BigDecimal_rounding_mode, \
3409  INT2FIX((int)(mode)) \
3410  )
3411 
3412 VP_EXPORT unsigned short
3414 {
3415  VALUE const vmode = rb_thread_local_aref(
3418  );
3419 
3420  if (NIL_P(vmode)) {
3423  }
3424 
3425  return (unsigned short)FIX2INT(vmode);
3426 }
3427 
3428 VP_EXPORT int
3429 VpIsRoundMode(unsigned short n)
3430 {
3431  switch (n) {
3432  case VP_ROUND_UP:
3433  case VP_ROUND_DOWN:
3434  case VP_ROUND_HALF_UP:
3435  case VP_ROUND_HALF_DOWN:
3436  case VP_ROUND_CEIL:
3437  case VP_ROUND_FLOOR:
3438  case VP_ROUND_HALF_EVEN:
3439  return 1;
3440 
3441  default:
3442  return 0;
3443  }
3444 }
3445 
3446 VP_EXPORT unsigned short
3447 VpSetRoundMode(unsigned short n)
3448 {
3449  if (VpIsRoundMode(n)) {
3451  return n;
3452  }
3453 
3454  return VpGetRoundMode();
3455 }
3456 
3457 /*
3458  * 0.0 & 1.0 generator
3459  * These gZero_..... and gOne_..... can be any name
3460  * referenced from nowhere except Zero() and One().
3461  * gZero_..... and gOne_..... must have global scope
3462  * (to let the compiler know they may be changed in outside
3463  * (... but not actually..)).
3464  */
3465 volatile const double gZero_ABCED9B1_CE73__00400511F31D = 0.0;
3466 volatile const double gOne_ABCED9B4_CE73__00400511F31D = 1.0;
3467 static double
3468 Zero(void)
3469 {
3471 }
3472 
3473 static double
3474 One(void)
3475 {
3477 }
3478 
3479 /*
3480  ----------------------------------------------------------------
3481  Value of sign in Real structure is reserved for future use.
3482  short sign;
3483  ==0 : NaN
3484  1 : Positive zero
3485  -1 : Negative zero
3486  2 : Positive number
3487  -2 : Negative number
3488  3 : Positive infinite number
3489  -3 : Negative infinite number
3490  ----------------------------------------------------------------
3491 */
3492 
3493 VP_EXPORT double
3494 VpGetDoubleNaN(void) /* Returns the value of NaN */
3495 {
3496  static double fNaN = 0.0;
3497  if(fNaN==0.0) fNaN = Zero()/Zero();
3498  return fNaN;
3499 }
3500 
3501 VP_EXPORT double
3502 VpGetDoublePosInf(void) /* Returns the value of +Infinity */
3503 {
3504  static double fInf = 0.0;
3505  if(fInf==0.0) fInf = One()/Zero();
3506  return fInf;
3507 }
3508 
3509 VP_EXPORT double
3510 VpGetDoubleNegInf(void) /* Returns the value of -Infinity */
3511 {
3512  static double fInf = 0.0;
3513  if(fInf==0.0) fInf = -(One()/Zero());
3514  return fInf;
3515 }
3516 
3517 VP_EXPORT double
3518 VpGetDoubleNegZero(void) /* Returns the value of -0 */
3519 {
3520  static double nzero = 1000.0;
3521  if(nzero!=0.0) nzero = (One()/VpGetDoubleNegInf());
3522  return nzero;
3523 }
3524 
3525 #if 0 /* unused */
3526 VP_EXPORT int
3527 VpIsNegDoubleZero(double v)
3528 {
3529  double z = VpGetDoubleNegZero();
3530  return MemCmp(&v,&z,sizeof(v))==0;
3531 }
3532 #endif
3533 
3534 VP_EXPORT int
3535 VpException(unsigned short f, const char *str,int always)
3536 {
3537  VALUE exc;
3538  int fatal=0;
3539  unsigned short const exception_mode = VpGetException();
3540 
3541  if(f==VP_EXCEPTION_OP || f==VP_EXCEPTION_MEMORY) always = 1;
3542 
3543  if (always || (exception_mode & f)) {
3544  switch(f)
3545  {
3546  /*
3547  case VP_EXCEPTION_OVERFLOW:
3548  */
3550  case VP_EXCEPTION_INFINITY:
3551  case VP_EXCEPTION_NaN:
3553  case VP_EXCEPTION_OP:
3554  exc = rb_eFloatDomainError;
3555  goto raise;
3556  case VP_EXCEPTION_MEMORY:
3557  fatal = 1;
3558  goto raise;
3559  default:
3560  fatal = 1;
3561  goto raise;
3562  }
3563  }
3564  return 0; /* 0 Means VpException() raised no exception */
3565 
3566 raise:
3567  if(fatal) rb_fatal("%s", str);
3568  else rb_raise(exc, "%s", str);
3569  return 0;
3570 }
3571 
3572 /* Throw exception or returns 0,when resulting c is Inf or NaN */
3573 /* sw=1:+ 2:- 3:* 4:/ */
3574 static int
3575 VpIsDefOP(Real *c,Real *a,Real *b,int sw)
3576 {
3577  if(VpIsNaN(a) || VpIsNaN(b)) {
3578  /* at least a or b is NaN */
3579  VpSetNaN(c);
3580  goto NaN;
3581  }
3582 
3583  if(VpIsInf(a)) {
3584  if(VpIsInf(b)) {
3585  switch(sw)
3586  {
3587  case 1: /* + */
3588  if(VpGetSign(a)==VpGetSign(b)) {
3589  VpSetInf(c,VpGetSign(a));
3590  goto Inf;
3591  } else {
3592  VpSetNaN(c);
3593  goto NaN;
3594  }
3595  case 2: /* - */
3596  if(VpGetSign(a)!=VpGetSign(b)) {
3597  VpSetInf(c,VpGetSign(a));
3598  goto Inf;
3599  } else {
3600  VpSetNaN(c);
3601  goto NaN;
3602  }
3603  break;
3604  case 3: /* * */
3605  VpSetInf(c,VpGetSign(a)*VpGetSign(b));
3606  goto Inf;
3607  break;
3608  case 4: /* / */
3609  VpSetNaN(c);
3610  goto NaN;
3611  }
3612  VpSetNaN(c);
3613  goto NaN;
3614  }
3615  /* Inf op Finite */
3616  switch(sw)
3617  {
3618  case 1: /* + */
3619  case 2: /* - */
3620  VpSetInf(c,VpGetSign(a));
3621  break;
3622  case 3: /* * */
3623  if(VpIsZero(b)) {
3624  VpSetNaN(c);
3625  goto NaN;
3626  }
3627  VpSetInf(c,VpGetSign(a)*VpGetSign(b));
3628  break;
3629  case 4: /* / */
3630  VpSetInf(c,VpGetSign(a)*VpGetSign(b));
3631  }
3632  goto Inf;
3633  }
3634 
3635  if(VpIsInf(b)) {
3636  switch(sw)
3637  {
3638  case 1: /* + */
3639  VpSetInf(c,VpGetSign(b));
3640  break;
3641  case 2: /* - */
3642  VpSetInf(c,-VpGetSign(b));
3643  break;
3644  case 3: /* * */
3645  if(VpIsZero(a)) {
3646  VpSetNaN(c);
3647  goto NaN;
3648  }
3649  VpSetInf(c,VpGetSign(a)*VpGetSign(b));
3650  break;
3651  case 4: /* / */
3652  VpSetZero(c,VpGetSign(a)*VpGetSign(b));
3653  }
3654  goto Inf;
3655  }
3656  return 1; /* Results OK */
3657 
3658 Inf:
3659  return VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0);
3660 NaN:
3661  return VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'",0);
3662 }
3663 
3664 /*
3665  ----------------------------------------------------------------
3666 */
3667 
3668 /*
3669  * returns number of chars needed to represent vp in specified format.
3670  */
3671 VP_EXPORT size_t
3672 VpNumOfChars(Real *vp,const char *pszFmt)
3673 {
3674  SIGNED_VALUE ex;
3675  size_t nc;
3676 
3677  if(vp == NULL) return BASE_FIG*2+6;
3678  if(!VpIsDef(vp)) return 32; /* not sure,may be OK */
3679 
3680  switch(*pszFmt)
3681  {
3682  case 'F':
3683  nc = BASE_FIG*(vp->Prec + 1)+2;
3684  ex = vp->exponent;
3685  if(ex < 0) {
3686  nc += BASE_FIG*(size_t)(-ex);
3687  }
3688  else {
3689  if((size_t)ex > vp->Prec) {
3690  nc += BASE_FIG*((size_t)ex - vp->Prec);
3691  }
3692  }
3693  break;
3694  case 'E':
3695  default:
3696  nc = BASE_FIG*(vp->Prec + 2)+6; /* 3: sign + exponent chars */
3697  }
3698  return nc;
3699 }
3700 
3701 /*
3702  * Initializer for Vp routines and constants used.
3703  * [Input]
3704  * BaseVal: Base value(assigned to BASE) for Vp calculation.
3705  * It must be the form BaseVal=10**n.(n=1,2,3,...)
3706  * If Base <= 0L,then the BASE will be calcurated so
3707  * that BASE is as large as possible satisfying the
3708  * relation MaxVal <= BASE*(BASE+1). Where the value
3709  * MaxVal is the largest value which can be represented
3710  * by one BDIGIT word in the computer used.
3711  *
3712  * [Returns]
3713  * 1+DBL_DIG ... OK
3714  */
3715 VP_EXPORT size_t
3716 VpInit(BDIGIT BaseVal)
3717 {
3718  /* Setup +/- Inf NaN -0 */
3719  VpGetDoubleNaN();
3723 
3724  /* Allocates Vp constants. */
3725  VpConstOne = VpAlloc(1UL, "1");
3726  VpPt5 = VpAlloc(1UL, ".5");
3727 
3728 #ifdef BIGDECIMAL_DEBUG
3729  gnAlloc = 0;
3730 #endif /* BIGDECIMAL_DEBUG */
3731 
3732 #ifdef BIGDECIMAL_DEBUG
3733  if(gfDebug) {
3734  printf("VpInit: BaseVal = %lu\n", BaseVal);
3735  printf(" BASE = %lu\n", BASE);
3736  printf(" HALF_BASE = %lu\n", HALF_BASE);
3737  printf(" BASE1 = %lu\n", BASE1);
3738  printf(" BASE_FIG = %u\n", BASE_FIG);
3739  printf(" DBLE_FIG = %d\n", DBLE_FIG);
3740  }
3741 #endif /* BIGDECIMAL_DEBUG */
3742 
3743  return rmpd_double_figures();
3744 }
3745 
3746 VP_EXPORT Real *
3747 VpOne(void)
3748 {
3749  return VpConstOne;
3750 }
3751 
3752 /* If exponent overflows,then raise exception or returns 0 */
3753 static int
3755 {
3756  SIGNED_VALUE e = a->exponent;
3757  SIGNED_VALUE m = e+n;
3758  SIGNED_VALUE eb, mb;
3759  if(e>0) {
3760  if(n>0) {
3761  mb = m*(SIGNED_VALUE)BASE_FIG;
3762  eb = e*(SIGNED_VALUE)BASE_FIG;
3763  if(mb<eb) goto overflow;
3764  }
3765  } else if(n<0) {
3766  mb = m*(SIGNED_VALUE)BASE_FIG;
3767  eb = e*(SIGNED_VALUE)BASE_FIG;
3768  if(mb>eb) goto underflow;
3769  }
3770  a->exponent = m;
3771  return 1;
3772 
3773 /* Overflow/Underflow ==> Raise exception or returns 0 */
3774 underflow:
3775  VpSetZero(a,VpGetSign(a));
3776  return VpException(VP_EXCEPTION_UNDERFLOW,"Exponent underflow",0);
3777 
3778 overflow:
3779  VpSetInf(a,VpGetSign(a));
3780  return VpException(VP_EXCEPTION_OVERFLOW,"Exponent overflow",0);
3781 }
3782 
3783 /*
3784  * Allocates variable.
3785  * [Input]
3786  * mx ... allocation unit, if zero then mx is determined by szVal.
3787  * The mx is the number of effective digits can to be stored.
3788  * szVal ... value assigned(char). If szVal==NULL,then zero is assumed.
3789  * If szVal[0]=='#' then Max. Prec. will not be considered(1.1.7),
3790  * full precision specified by szVal is allocated.
3791  *
3792  * [Returns]
3793  * Pointer to the newly allocated variable, or
3794  * NULL be returned if memory allocation is failed,or any error.
3795  */
3796 VP_EXPORT Real *
3797 VpAlloc(size_t mx, const char *szVal)
3798 {
3799  size_t i, ni, ipn, ipf, nf, ipe, ne, nalloc;
3800  char v,*psz;
3801  int sign=1;
3802  Real *vp = NULL;
3803  size_t mf = VpGetPrecLimit();
3804  VALUE buf;
3805 
3806  mx = (mx + BASE_FIG - 1) / BASE_FIG + 1; /* Determine allocation unit. */
3807  if (szVal) {
3808  while (ISSPACE(*szVal)) szVal++;
3809  if (*szVal != '#') {
3810  if (mf) {
3811  mf = (mf + BASE_FIG - 1) / BASE_FIG + 2; /* Needs 1 more for div */
3812  if (mx > mf) {
3813  mx = mf;
3814  }
3815  }
3816  }
3817  else {
3818  ++szVal;
3819  }
3820  }
3821  else {
3822  /* necessary to be able to store */
3823  /* at least mx digits. */
3824  /* szVal==NULL ==> allocate zero value. */
3825  vp = VpAllocReal(mx);
3826  /* xmalloc() alway returns(or throw interruption) */
3827  vp->MaxPrec = mx; /* set max precision */
3828  VpSetZero(vp,1); /* initialize vp to zero. */
3829  return vp;
3830  }
3831 
3832  /* Skip all '_' after digit: 2006-6-30 */
3833  ni = 0;
3834  buf = rb_str_tmp_new(strlen(szVal)+1);
3835  psz = RSTRING_PTR(buf);
3836  i = 0;
3837  ipn = 0;
3838  while ((psz[i]=szVal[ipn]) != 0) {
3839  if (ISDIGIT(psz[i])) ++ni;
3840  if (psz[i] == '_') {
3841  if (ni > 0) { ipn++; continue; }
3842  psz[i] = 0;
3843  break;
3844  }
3845  ++i;
3846  ++ipn;
3847  }
3848  /* Skip trailing spaces */
3849  while (--i > 0) {
3850  if (ISSPACE(psz[i])) psz[i] = 0;
3851  else break;
3852  }
3853  szVal = psz;
3854 
3855  /* Check on Inf & NaN */
3856  if (StrCmp(szVal, SZ_PINF) == 0 ||
3857  StrCmp(szVal, SZ_INF) == 0 ) {
3858  vp = VpAllocReal(1);
3859  vp->MaxPrec = 1; /* set max precision */
3860  VpSetPosInf(vp);
3861  return vp;
3862  }
3863  if (StrCmp(szVal, SZ_NINF) == 0) {
3864  vp = VpAllocReal(1);
3865  vp->MaxPrec = 1; /* set max precision */
3866  VpSetNegInf(vp);
3867  return vp;
3868  }
3869  if (StrCmp(szVal, SZ_NaN) == 0) {
3870  vp = VpAllocReal(1);
3871  vp->MaxPrec = 1; /* set max precision */
3872  VpSetNaN(vp);
3873  return vp;
3874  }
3875 
3876  /* check on number szVal[] */
3877  ipn = i = 0;
3878  if (szVal[i] == '-') { sign=-1; ++i; }
3879  else if (szVal[i] == '+') ++i;
3880  /* Skip digits */
3881  ni = 0; /* digits in mantissa */
3882  while ((v = szVal[i]) != 0) {
3883  if (!ISDIGIT(v)) break;
3884  ++i;
3885  ++ni;
3886  }
3887  nf = 0;
3888  ipf = 0;
3889  ipe = 0;
3890  ne = 0;
3891  if (v) {
3892  /* other than digit nor \0 */
3893  if (szVal[i] == '.') { /* xxx. */
3894  ++i;
3895  ipf = i;
3896  while ((v = szVal[i]) != 0) { /* get fraction part. */
3897  if (!ISDIGIT(v)) break;
3898  ++i;
3899  ++nf;
3900  }
3901  }
3902  ipe = 0; /* Exponent */
3903 
3904  switch (szVal[i]) {
3905  case '\0':
3906  break;
3907  case 'e': case 'E':
3908  case 'd': case 'D':
3909  ++i;
3910  ipe = i;
3911  v = szVal[i];
3912  if ((v == '-') || (v == '+')) ++i;
3913  while ((v=szVal[i]) != 0) {
3914  if (!ISDIGIT(v)) break;
3915  ++i;
3916  ++ne;
3917  }
3918  break;
3919  default:
3920  break;
3921  }
3922  }
3923  nalloc = (ni + nf + BASE_FIG - 1) / BASE_FIG + 1; /* set effective allocation */
3924  /* units for szVal[] */
3925  if (mx <= 0) mx = 1;
3926  nalloc = Max(nalloc, mx);
3927  mx = nalloc;
3928  vp = VpAllocReal(mx);
3929  /* xmalloc() alway returns(or throw interruption) */
3930  vp->MaxPrec = mx; /* set max precision */
3931  VpSetZero(vp, sign);
3932  VpCtoV(vp, &szVal[ipn], ni, &szVal[ipf], nf, &szVal[ipe], ne);
3933  rb_str_resize(buf, 0);
3934  return vp;
3935 }
3936 
3937 /*
3938  * Assignment(c=a).
3939  * [Input]
3940  * a ... RHSV
3941  * isw ... switch for assignment.
3942  * c = a when isw > 0
3943  * c = -a when isw < 0
3944  * if c->MaxPrec < a->Prec,then round operation
3945  * will be performed.
3946  * [Output]
3947  * c ... LHSV
3948  */
3949 VP_EXPORT size_t
3950 VpAsgn(Real *c, Real *a, int isw)
3951 {
3952  size_t n;
3953  if(VpIsNaN(a)) {
3954  VpSetNaN(c);
3955  return 0;
3956  }
3957  if(VpIsInf(a)) {
3958  VpSetInf(c,isw*VpGetSign(a));
3959  return 0;
3960  }
3961 
3962  /* check if the RHS is zero */
3963  if(!VpIsZero(a)) {
3964  c->exponent = a->exponent; /* store exponent */
3965  VpSetSign(c,(isw*VpGetSign(a))); /* set sign */
3966  n =(a->Prec < c->MaxPrec) ?(a->Prec) :(c->MaxPrec);
3967  c->Prec = n;
3968  memcpy(c->frac, a->frac, n * sizeof(BDIGIT));
3969  /* Needs round ? */
3970  if(isw!=10) {
3971  /* Not in ActiveRound */
3972  if(c->Prec < a->Prec) {
3973  VpInternalRound(c,n,(n>0)?a->frac[n-1]:0,a->frac[n]);
3974  } else {
3975  VpLimitRound(c,0);
3976  }
3977  }
3978  } else {
3979  /* The value of 'a' is zero. */
3980  VpSetZero(c,isw*VpGetSign(a));
3981  return 1;
3982  }
3983  return c->Prec*BASE_FIG;
3984 }
3985 
3986 /*
3987  * c = a + b when operation = 1 or 2
3988  * = a - b when operation = -1 or -2.
3989  * Returns number of significant digits of c
3990  */
3991 VP_EXPORT size_t
3992 VpAddSub(Real *c, Real *a, Real *b, int operation)
3993 {
3994  short sw, isw;
3995  Real *a_ptr, *b_ptr;
3996  size_t n, na, nb, i;
3997  BDIGIT mrv;
3998 
3999 #ifdef BIGDECIMAL_DEBUG
4000  if(gfDebug) {
4001  VPrint(stdout, "VpAddSub(enter) a=% \n", a);
4002  VPrint(stdout, " b=% \n", b);
4003  printf(" operation=%d\n", operation);
4004  }
4005 #endif /* BIGDECIMAL_DEBUG */
4006 
4007  if(!VpIsDefOP(c,a,b,(operation>0)?1:2)) return 0; /* No significant digits */
4008 
4009  /* check if a or b is zero */
4010  if(VpIsZero(a)) {
4011  /* a is zero,then assign b to c */
4012  if(!VpIsZero(b)) {
4013  VpAsgn(c, b, operation);
4014  } else {
4015  /* Both a and b are zero. */
4016  if(VpGetSign(a)<0 && operation*VpGetSign(b)<0) {
4017  /* -0 -0 */
4018  VpSetZero(c,-1);
4019  } else {
4020  VpSetZero(c,1);
4021  }
4022  return 1; /* 0: 1 significant digits */
4023  }
4024  return c->Prec*BASE_FIG;
4025  }
4026  if(VpIsZero(b)) {
4027  /* b is zero,then assign a to c. */
4028  VpAsgn(c, a, 1);
4029  return c->Prec*BASE_FIG;
4030  }
4031 
4032  if(operation < 0) sw = -1;
4033  else sw = 1;
4034 
4035  /* compare absolute value. As a result,|a_ptr|>=|b_ptr| */
4036  if(a->exponent > b->exponent) {
4037  a_ptr = a;
4038  b_ptr = b;
4039  } /* |a|>|b| */
4040  else if(a->exponent < b->exponent) {
4041  a_ptr = b;
4042  b_ptr = a;
4043  } /* |a|<|b| */
4044  else {
4045  /* Exponent part of a and b is the same,then compare fraction */
4046  /* part */
4047  na = a->Prec;
4048  nb = b->Prec;
4049  n = Min(na, nb);
4050  for(i=0;i < n; ++i) {
4051  if(a->frac[i] > b->frac[i]) {
4052  a_ptr = a;
4053  b_ptr = b;
4054  goto end_if;
4055  } else if(a->frac[i] < b->frac[i]) {
4056  a_ptr = b;
4057  b_ptr = a;
4058  goto end_if;
4059  }
4060  }
4061  if(na > nb) {
4062  a_ptr = a;
4063  b_ptr = b;
4064  goto end_if;
4065  } else if(na < nb) {
4066  a_ptr = b;
4067  b_ptr = a;
4068  goto end_if;
4069  }
4070  /* |a| == |b| */
4071  if(VpGetSign(a) + sw *VpGetSign(b) == 0) {
4072  VpSetZero(c,1); /* abs(a)=abs(b) and operation = '-' */
4073  return c->Prec*BASE_FIG;
4074  }
4075  a_ptr = a;
4076  b_ptr = b;
4077  }
4078 
4079 end_if:
4080  isw = VpGetSign(a) + sw *VpGetSign(b);
4081  /*
4082  * isw = 0 ...( 1)+(-1),( 1)-( 1),(-1)+(1),(-1)-(-1)
4083  * = 2 ...( 1)+( 1),( 1)-(-1)
4084  * =-2 ...(-1)+(-1),(-1)-( 1)
4085  * If isw==0, then c =(Sign a_ptr)(|a_ptr|-|b_ptr|)
4086  * else c =(Sign ofisw)(|a_ptr|+|b_ptr|)
4087  */
4088  if(isw) { /* addition */
4089  VpSetSign(c, 1);
4090  mrv = VpAddAbs(a_ptr, b_ptr, c);
4091  VpSetSign(c, isw / 2);
4092  } else { /* subtraction */
4093  VpSetSign(c, 1);
4094  mrv = VpSubAbs(a_ptr, b_ptr, c);
4095  if(a_ptr == a) {
4096  VpSetSign(c,VpGetSign(a));
4097  } else {
4098  VpSetSign(c,VpGetSign(a_ptr) * sw);
4099  }
4100  }
4101  VpInternalRound(c,0,(c->Prec>0)?c->frac[c->Prec-1]:0,mrv);
4102 
4103 #ifdef BIGDECIMAL_DEBUG
4104  if(gfDebug) {
4105  VPrint(stdout, "VpAddSub(result) c=% \n", c);
4106  VPrint(stdout, " a=% \n", a);
4107  VPrint(stdout, " b=% \n", b);
4108  printf(" operation=%d\n", operation);
4109  }
4110 #endif /* BIGDECIMAL_DEBUG */
4111  return c->Prec*BASE_FIG;
4112 }
4113 
4114 /*
4115  * Addition of two variable precisional variables
4116  * a and b assuming abs(a)>abs(b).
4117  * c = abs(a) + abs(b) ; where |a|>=|b|
4118  */
4119 static BDIGIT
4120 VpAddAbs(Real *a, Real *b, Real *c)
4121 {
4122  size_t word_shift;
4123  size_t ap;
4124  size_t bp;
4125  size_t cp;
4126  size_t a_pos;
4127  size_t b_pos, b_pos_with_word_shift;
4128  size_t c_pos;
4129  BDIGIT av, bv, carry, mrv;
4130 
4131 #ifdef BIGDECIMAL_DEBUG
4132  if(gfDebug) {
4133  VPrint(stdout, "VpAddAbs called: a = %\n", a);
4134  VPrint(stdout, " b = %\n", b);
4135  }
4136 #endif /* BIGDECIMAL_DEBUG */
4137 
4138  word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
4139  a_pos = ap;
4140  b_pos = bp;
4141  c_pos = cp;
4142  if(word_shift==(size_t)-1L) return 0; /* Overflow */
4143  if(b_pos == (size_t)-1L) goto Assign_a;
4144 
4145  mrv = av + bv; /* Most right val. Used for round. */
4146 
4147  /* Just assign the last few digits of b to c because a has no */
4148  /* corresponding digits to be added. */
4149  while(b_pos + word_shift > a_pos) {
4150  --c_pos;
4151  if(b_pos > 0) {
4152  c->frac[c_pos] = b->frac[--b_pos];
4153  } else {
4154  --word_shift;
4155  c->frac[c_pos] = 0;
4156  }
4157  }
4158 
4159  /* Just assign the last few digits of a to c because b has no */
4160  /* corresponding digits to be added. */
4161  b_pos_with_word_shift = b_pos + word_shift;
4162  while(a_pos > b_pos_with_word_shift) {
4163  c->frac[--c_pos] = a->frac[--a_pos];
4164  }
4165  carry = 0; /* set first carry be zero */
4166 
4167  /* Now perform addition until every digits of b will be */
4168  /* exhausted. */
4169  while(b_pos > 0) {
4170  c->frac[--c_pos] = a->frac[--a_pos] + b->frac[--b_pos] + carry;
4171  if(c->frac[c_pos] >= BASE) {
4172  c->frac[c_pos] -= BASE;
4173  carry = 1;
4174  } else {
4175  carry = 0;
4176  }
4177  }
4178 
4179  /* Just assign the first few digits of a with considering */
4180  /* the carry obtained so far because b has been exhausted. */
4181  while(a_pos > 0) {
4182  c->frac[--c_pos] = a->frac[--a_pos] + carry;
4183  if(c->frac[c_pos] >= BASE) {
4184  c->frac[c_pos] -= BASE;
4185  carry = 1;
4186  } else {
4187  carry = 0;
4188  }
4189  }
4190  if(c_pos) c->frac[c_pos - 1] += carry;
4191  goto Exit;
4192 
4193 Assign_a:
4194  VpAsgn(c, a, 1);
4195  mrv = 0;
4196 
4197 Exit:
4198 
4199 #ifdef BIGDECIMAL_DEBUG
4200  if(gfDebug) {
4201  VPrint(stdout, "VpAddAbs exit: c=% \n", c);
4202  }
4203 #endif /* BIGDECIMAL_DEBUG */
4204  return mrv;
4205 }
4206 
4207 /*
4208  * c = abs(a) - abs(b)
4209  */
4210 static BDIGIT
4211 VpSubAbs(Real *a, Real *b, Real *c)
4212 {
4213  size_t word_shift;
4214  size_t ap;
4215  size_t bp;
4216  size_t cp;
4217  size_t a_pos;
4218  size_t b_pos, b_pos_with_word_shift;
4219  size_t c_pos;
4220  BDIGIT av, bv, borrow, mrv;
4221 
4222 #ifdef BIGDECIMAL_DEBUG
4223  if(gfDebug) {
4224  VPrint(stdout, "VpSubAbs called: a = %\n", a);
4225  VPrint(stdout, " b = %\n", b);
4226  }
4227 #endif /* BIGDECIMAL_DEBUG */
4228 
4229  word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
4230  a_pos = ap;
4231  b_pos = bp;
4232  c_pos = cp;
4233  if(word_shift==(size_t)-1L) return 0; /* Overflow */
4234  if(b_pos == (size_t)-1L) goto Assign_a;
4235 
4236  if(av >= bv) {
4237  mrv = av - bv;
4238  borrow = 0;
4239  } else {
4240  mrv = 0;
4241  borrow = 1;
4242  }
4243 
4244  /* Just assign the values which are the BASE subtracted by */
4245  /* each of the last few digits of the b because the a has no */
4246  /* corresponding digits to be subtracted. */
4247  if(b_pos + word_shift > a_pos) {
4248  while(b_pos + word_shift > a_pos) {
4249  --c_pos;
4250  if(b_pos > 0) {
4251  c->frac[c_pos] = BASE - b->frac[--b_pos] - borrow;
4252  } else {
4253  --word_shift;
4254  c->frac[c_pos] = BASE - borrow;
4255  }
4256  borrow = 1;
4257  }
4258  }
4259  /* Just assign the last few digits of a to c because b has no */
4260  /* corresponding digits to subtract. */
4261 
4262  b_pos_with_word_shift = b_pos + word_shift;
4263  while(a_pos > b_pos_with_word_shift) {
4264  c->frac[--c_pos] = a->frac[--a_pos];
4265  }
4266 
4267  /* Now perform subtraction until every digits of b will be */
4268  /* exhausted. */
4269  while(b_pos > 0) {
4270  --c_pos;
4271  if(a->frac[--a_pos] < b->frac[--b_pos] + borrow) {
4272  c->frac[c_pos] = BASE + a->frac[a_pos] - b->frac[b_pos] - borrow;
4273  borrow = 1;
4274  } else {
4275  c->frac[c_pos] = a->frac[a_pos] - b->frac[b_pos] - borrow;
4276  borrow = 0;
4277  }
4278  }
4279 
4280  /* Just assign the first few digits of a with considering */
4281  /* the borrow obtained so far because b has been exhausted. */
4282  while(a_pos > 0) {
4283  --c_pos;
4284  if(a->frac[--a_pos] < borrow) {
4285  c->frac[c_pos] = BASE + a->frac[a_pos] - borrow;
4286  borrow = 1;
4287  } else {
4288  c->frac[c_pos] = a->frac[a_pos] - borrow;
4289  borrow = 0;
4290  }
4291  }
4292  if(c_pos) c->frac[c_pos - 1] -= borrow;
4293  goto Exit;
4294 
4295 Assign_a:
4296  VpAsgn(c, a, 1);
4297  mrv = 0;
4298 
4299 Exit:
4300 #ifdef BIGDECIMAL_DEBUG
4301  if(gfDebug) {
4302  VPrint(stdout, "VpSubAbs exit: c=% \n", c);
4303  }
4304 #endif /* BIGDECIMAL_DEBUG */
4305  return mrv;
4306 }
4307 
4308 /*
4309  * Note: If(av+bv)>= HALF_BASE,then 1 will be added to the least significant
4310  * digit of c(In case of addition).
4311  * ------------------------- figure of output -----------------------------------
4312  * a = xxxxxxxxxxx
4313  * b = xxxxxxxxxx
4314  * c =xxxxxxxxxxxxxxx
4315  * word_shift = | |
4316  * right_word = | | (Total digits in RHSV)
4317  * left_word = | | (Total digits in LHSV)
4318  * a_pos = |
4319  * b_pos = |
4320  * c_pos = |
4321  */
4322 static size_t
4323 VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, BDIGIT *av, BDIGIT *bv)
4324 {
4325  size_t left_word, right_word, word_shift;
4326  c->frac[0] = 0;
4327  *av = *bv = 0;
4328  word_shift =((a->exponent) -(b->exponent));
4329  left_word = b->Prec + word_shift;
4330  right_word = Max((a->Prec),left_word);
4331  left_word =(c->MaxPrec) - 1; /* -1 ... prepare for round up */
4332  /*
4333  * check if 'round' is needed.
4334  */
4335  if(right_word > left_word) { /* round ? */
4336  /*---------------------------------
4337  * Actual size of a = xxxxxxAxx
4338  * Actual size of b = xxxBxxxxx
4339  * Max. size of c = xxxxxx
4340  * Round off = |-----|
4341  * c_pos = |
4342  * right_word = |
4343  * a_pos = |
4344  */
4345  *c_pos = right_word = left_word + 1; /* Set resulting precision */
4346  /* be equal to that of c */
4347  if((a->Prec) >=(c->MaxPrec)) {
4348  /*
4349  * a = xxxxxxAxxx
4350  * c = xxxxxx
4351  * a_pos = |
4352  */
4353  *a_pos = left_word;
4354  *av = a->frac[*a_pos]; /* av is 'A' shown in above. */
4355  } else {
4356  /*
4357  * a = xxxxxxx
4358  * c = xxxxxxxxxx
4359  * a_pos = |
4360  */
4361  *a_pos = a->Prec;
4362  }
4363  if((b->Prec + word_shift) >= c->MaxPrec) {
4364  /*
4365  * a = xxxxxxxxx
4366  * b = xxxxxxxBxxx
4367  * c = xxxxxxxxxxx
4368  * b_pos = |
4369  */
4370  if(c->MaxPrec >=(word_shift + 1)) {
4371  *b_pos = c->MaxPrec - word_shift - 1;
4372  *bv = b->frac[*b_pos];
4373  } else {
4374  *b_pos = -1L;
4375  }
4376  } else {
4377  /*
4378  * a = xxxxxxxxxxxxxxxx
4379  * b = xxxxxx
4380  * c = xxxxxxxxxxxxx
4381  * b_pos = |
4382  */
4383  *b_pos = b->Prec;
4384  }
4385  } else { /* The MaxPrec of c - 1 > The Prec of a + b */
4386  /*
4387  * a = xxxxxxx
4388  * b = xxxxxx
4389  * c = xxxxxxxxxxx
4390  * c_pos = |
4391  */
4392  *b_pos = b->Prec;
4393  *a_pos = a->Prec;
4394  *c_pos = right_word + 1;
4395  }
4396  c->Prec = *c_pos;
4397  c->exponent = a->exponent;
4398  if(!AddExponent(c,1)) return (size_t)-1L;
4399  return word_shift;
4400 }
4401 
4402 /*
4403  * Return number og significant digits
4404  * c = a * b , Where a = a0a1a2 ... an
4405  * b = b0b1b2 ... bm
4406  * c = c0c1c2 ... cl
4407  * a0 a1 ... an * bm
4408  * a0 a1 ... an * bm-1
4409  * . . .
4410  * . . .
4411  * a0 a1 .... an * b0
4412  * +_____________________________
4413  * c0 c1 c2 ...... cl
4414  * nc <---|
4415  * MaxAB |--------------------|
4416  */
4417 VP_EXPORT size_t
4418 VpMult(Real *c, Real *a, Real *b)
4419 {
4420  size_t MxIndA, MxIndB, MxIndAB, MxIndC;
4421  size_t ind_c, i, ii, nc;
4422  size_t ind_as, ind_ae, ind_bs;
4423  BDIGIT carry;
4424  BDIGIT_DBL s;
4425  Real *w;
4426 
4427 #ifdef BIGDECIMAL_DEBUG
4428  if(gfDebug) {
4429  VPrint(stdout, "VpMult(Enter): a=% \n", a);
4430  VPrint(stdout, " b=% \n", b);
4431  }
4432 #endif /* BIGDECIMAL_DEBUG */
4433 
4434  if(!VpIsDefOP(c,a,b,3)) return 0; /* No significant digit */
4435 
4436  if(VpIsZero(a) || VpIsZero(b)) {
4437  /* at least a or b is zero */
4438  VpSetZero(c,VpGetSign(a)*VpGetSign(b));
4439  return 1; /* 0: 1 significant digit */
4440  }
4441 
4442  if(VpIsOne(a)) {
4443  VpAsgn(c, b, VpGetSign(a));
4444  goto Exit;
4445  }
4446  if(VpIsOne(b)) {
4447  VpAsgn(c, a, VpGetSign(b));
4448  goto Exit;
4449  }
4450  if((b->Prec) >(a->Prec)) {
4451  /* Adjust so that digits(a)>digits(b) */
4452  w = a;
4453  a = b;
4454  b = w;
4455  }
4456  w = NULL;
4457  MxIndA = a->Prec - 1;
4458  MxIndB = b->Prec - 1;
4459  MxIndC = c->MaxPrec - 1;
4460  MxIndAB = a->Prec + b->Prec - 1;
4461 
4462  if(MxIndC < MxIndAB) { /* The Max. prec. of c < Prec(a)+Prec(b) */
4463  w = c;
4464  c = VpAlloc((size_t)((MxIndAB + 1) * BASE_FIG), "#0");
4465  MxIndC = MxIndAB;
4466  }
4467 
4468  /* set LHSV c info */
4469 
4470  c->exponent = a->exponent; /* set exponent */
4471  if(!AddExponent(c,b->exponent)) {
4472  if(w) VpFree(c);
4473  return 0;
4474  }
4475  VpSetSign(c,VpGetSign(a)*VpGetSign(b)); /* set sign */
4476  carry = 0;
4477  nc = ind_c = MxIndAB;
4478  memset(c->frac, 0, (nc + 1) * sizeof(BDIGIT)); /* Initialize c */
4479  c->Prec = nc + 1; /* set precision */
4480  for(nc = 0; nc < MxIndAB; ++nc, --ind_c) {
4481  if(nc < MxIndB) { /* The left triangle of the Fig. */
4482  ind_as = MxIndA - nc;
4483  ind_ae = MxIndA;
4484  ind_bs = MxIndB;
4485  } else if(nc <= MxIndA) { /* The middle rectangular of the Fig. */
4486  ind_as = MxIndA - nc;
4487  ind_ae = MxIndA -(nc - MxIndB);
4488  ind_bs = MxIndB;
4489  } else if(nc > MxIndA) { /* The right triangle of the Fig. */
4490  ind_as = 0;
4491  ind_ae = MxIndAB - nc - 1;
4492  ind_bs = MxIndB -(nc - MxIndA);
4493  }
4494 
4495  for(i = ind_as; i <= ind_ae; ++i) {
4496  s = (BDIGIT_DBL)a->frac[i] * b->frac[ind_bs--];
4497  carry = (BDIGIT)(s / BASE);
4498  s -= (BDIGIT_DBL)carry * BASE;
4499  c->frac[ind_c] += (BDIGIT)s;
4500  if(c->frac[ind_c] >= BASE) {
4501  s = c->frac[ind_c] / BASE;
4502  carry += (BDIGIT)s;
4503  c->frac[ind_c] -= (BDIGIT)(s * BASE);
4504  }
4505  if(carry) {
4506  ii = ind_c;
4507  while(ii-- > 0) {
4508  c->frac[ii] += carry;
4509  if(c->frac[ii] >= BASE) {
4510  carry = c->frac[ii] / BASE;
4511  c->frac[ii] -= (carry * BASE);
4512  } else {
4513  break;
4514  }
4515  }
4516  }
4517  }
4518  }
4519  if(w != NULL) { /* free work variable */
4520  VpNmlz(c);
4521  VpAsgn(w, c, 1);
4522  VpFree(c);
4523  c = w;
4524  } else {
4525  VpLimitRound(c,0);
4526  }
4527 
4528 Exit:
4529 #ifdef BIGDECIMAL_DEBUG
4530  if(gfDebug) {
4531  VPrint(stdout, "VpMult(c=a*b): c=% \n", c);
4532  VPrint(stdout, " a=% \n", a);
4533  VPrint(stdout, " b=% \n", b);
4534  }
4535 #endif /*BIGDECIMAL_DEBUG */
4536  return c->Prec*BASE_FIG;
4537 }
4538 
4539 /*
4540  * c = a / b, remainder = r
4541  */
4542 VP_EXPORT size_t
4543 VpDivd(Real *c, Real *r, Real *a, Real *b)
4544 {
4545  size_t word_a, word_b, word_c, word_r;
4546  size_t i, n, ind_a, ind_b, ind_c, ind_r;
4547  size_t nLoop;
4548  BDIGIT_DBL q, b1, b1p1, b1b2, b1b2p1, r1r2;
4549  BDIGIT borrow, borrow1, borrow2;
4550  BDIGIT_DBL qb;
4551 
4552 #ifdef BIGDECIMAL_DEBUG
4553  if(gfDebug) {
4554  VPrint(stdout, " VpDivd(c=a/b) a=% \n", a);
4555  VPrint(stdout, " b=% \n", b);
4556  }
4557 #endif /*BIGDECIMAL_DEBUG */
4558 
4559  VpSetNaN(r);
4560  if(!VpIsDefOP(c,a,b,4)) goto Exit;
4561  if(VpIsZero(a)&&VpIsZero(b)) {
4562  VpSetNaN(c);
4563  return VpException(VP_EXCEPTION_NaN,"(VpDivd) 0/0 not defined(NaN)",0);
4564  }
4565  if(VpIsZero(b)) {
4566  VpSetInf(c,VpGetSign(a)*VpGetSign(b));
4567  return VpException(VP_EXCEPTION_ZERODIVIDE,"(VpDivd) Divide by zero",0);
4568  }
4569  if(VpIsZero(a)) {
4570  /* numerator a is zero */
4571  VpSetZero(c,VpGetSign(a)*VpGetSign(b));
4572  VpSetZero(r,VpGetSign(a)*VpGetSign(b));
4573  goto Exit;
4574  }
4575  if(VpIsOne(b)) {
4576  /* divide by one */
4577  VpAsgn(c, a, VpGetSign(b));
4578  VpSetZero(r,VpGetSign(a));
4579  goto Exit;
4580  }
4581 
4582  word_a = a->Prec;
4583  word_b = b->Prec;
4584  word_c = c->MaxPrec;
4585  word_r = r->MaxPrec;
4586 
4587  ind_c = 0;
4588  ind_r = 1;
4589 
4590  if(word_a >= word_r) goto space_error;
4591 
4592  r->frac[0] = 0;
4593  while(ind_r <= word_a) {
4594  r->frac[ind_r] = a->frac[ind_r - 1];
4595  ++ind_r;
4596  }
4597 
4598  while(ind_r < word_r) r->frac[ind_r++] = 0;
4599  while(ind_c < word_c) c->frac[ind_c++] = 0;
4600 
4601  /* initial procedure */
4602  b1 = b1p1 = b->frac[0];
4603  if(b->Prec <= 1) {
4604  b1b2p1 = b1b2 = b1p1 * BASE;
4605  } else {
4606  b1p1 = b1 + 1;
4607  b1b2p1 = b1b2 = b1 * BASE + b->frac[1];
4608  if(b->Prec > 2) ++b1b2p1;
4609  }
4610 
4611  /* */
4612  /* loop start */
4613  ind_c = word_r - 1;
4614  nLoop = Min(word_c,ind_c);
4615  ind_c = 1;
4616  while(ind_c < nLoop) {
4617  if(r->frac[ind_c] == 0) {
4618  ++ind_c;
4619  continue;
4620  }
4621  r1r2 = (BDIGIT_DBL)r->frac[ind_c] * BASE + r->frac[ind_c + 1];
4622  if(r1r2 == b1b2) {
4623  /* The first two word digits is the same */
4624  ind_b = 2;
4625  ind_a = ind_c + 2;
4626  while(ind_b < word_b) {
4627  if(r->frac[ind_a] < b->frac[ind_b]) goto div_b1p1;
4628  if(r->frac[ind_a] > b->frac[ind_b]) break;
4629  ++ind_a;
4630  ++ind_b;
4631  }
4632  /* The first few word digits of r and b is the same and */
4633  /* the first different word digit of w is greater than that */
4634  /* of b, so quotinet is 1 and just subtract b from r. */
4635  borrow = 0; /* quotient=1, then just r-b */
4636  ind_b = b->Prec - 1;
4637  ind_r = ind_c + ind_b;
4638  if(ind_r >= word_r) goto space_error;
4639  n = ind_b;
4640  for(i = 0; i <= n; ++i) {
4641  if(r->frac[ind_r] < b->frac[ind_b] + borrow) {
4642  r->frac[ind_r] += (BASE - (b->frac[ind_b] + borrow));
4643  borrow = 1;
4644  } else {
4645  r->frac[ind_r] = r->frac[ind_r] - b->frac[ind_b] - borrow;
4646  borrow = 0;
4647  }
4648  --ind_r;
4649  --ind_b;
4650  }
4651  ++c->frac[ind_c];
4652  goto carry;
4653  }
4654  /* The first two word digits is not the same, */
4655  /* then compare magnitude, and divide actually. */
4656  if(r1r2 >= b1b2p1) {
4657  q = r1r2 / b1b2p1; /* q == (BDIGIT)q */
4658  c->frac[ind_c] += (BDIGIT)q;
4659  ind_r = b->Prec + ind_c - 1;
4660  goto sub_mult;
4661  }
4662 
4663 div_b1p1:
4664  if(ind_c + 1 >= word_c) goto out_side;
4665  q = r1r2 / b1p1; /* q == (BDIGIT)q */
4666  c->frac[ind_c + 1] += (BDIGIT)q;
4667  ind_r = b->Prec + ind_c;
4668 
4669 sub_mult:
4670  borrow1 = borrow2 = 0;
4671  ind_b = word_b - 1;
4672  if(ind_r >= word_r) goto space_error;
4673  n = ind_b;
4674  for(i = 0; i <= n; ++i) {
4675  /* now, perform r = r - q * b */
4676  qb = q * b->frac[ind_b];
4677  if (qb < BASE) borrow1 = 0;
4678  else {
4679  borrow1 = (BDIGIT)(qb / BASE);
4680  qb -= (BDIGIT_DBL)borrow1 * BASE; /* get qb < BASE */
4681  }
4682  if(r->frac[ind_r] < qb) {
4683  r->frac[ind_r] += (BDIGIT)(BASE - qb);
4684  borrow2 = borrow2 + borrow1 + 1;
4685  } else {
4686  r->frac[ind_r] -= (BDIGIT)qb;
4687  borrow2 += borrow1;
4688  }
4689  if(borrow2) {
4690  if(r->frac[ind_r - 1] < borrow2) {
4691  r->frac[ind_r - 1] += (BASE - borrow2);
4692  borrow2 = 1;
4693  } else {
4694  r->frac[ind_r - 1] -= borrow2;
4695  borrow2 = 0;
4696  }
4697  }
4698  --ind_r;
4699  --ind_b;
4700  }
4701 
4702  r->frac[ind_r] -= borrow2;
4703 carry:
4704  ind_r = ind_c;
4705  while(c->frac[ind_r] >= BASE) {
4706  c->frac[ind_r] -= BASE;
4707  --ind_r;
4708  ++c->frac[ind_r];
4709  }
4710  }
4711  /* End of operation, now final arrangement */
4712 out_side:
4713  c->Prec = word_c;
4714  c->exponent = a->exponent;
4715  if(!AddExponent(c,2)) return 0;
4716  if(!AddExponent(c,-(b->exponent))) return 0;
4717 
4718  VpSetSign(c,VpGetSign(a)*VpGetSign(b));
4719  VpNmlz(c); /* normalize c */
4720  r->Prec = word_r;
4721  r->exponent = a->exponent;
4722  if(!AddExponent(r,1)) return 0;
4723  VpSetSign(r,VpGetSign(a));
4724  VpNmlz(r); /* normalize r(remainder) */
4725  goto Exit;
4726 
4727 space_error:
4728 #ifdef BIGDECIMAL_DEBUG
4729  if(gfDebug) {
4730  printf(" word_a=%lu\n", word_a);
4731  printf(" word_b=%lu\n", word_b);
4732  printf(" word_c=%lu\n", word_c);
4733  printf(" word_r=%lu\n", word_r);
4734  printf(" ind_r =%lu\n", ind_r);
4735  }
4736 #endif /* BIGDECIMAL_DEBUG */
4737  rb_bug("ERROR(VpDivd): space for remainder too small.");
4738 
4739 Exit:
4740 #ifdef BIGDECIMAL_DEBUG
4741  if(gfDebug) {
4742  VPrint(stdout, " VpDivd(c=a/b), c=% \n", c);
4743  VPrint(stdout, " r=% \n", r);
4744  }
4745 #endif /* BIGDECIMAL_DEBUG */
4746  return c->Prec*BASE_FIG;
4747 }
4748 
4749 /*
4750  * Input a = 00000xxxxxxxx En(5 preceeding zeros)
4751  * Output a = xxxxxxxx En-5
4752  */
4753 static int
4755 {
4756  size_t ind_a, i;
4757 
4758  if (!VpIsDef(a)) goto NoVal;
4759  if (VpIsZero(a)) goto NoVal;
4760 
4761  ind_a = a->Prec;
4762  while (ind_a--) {
4763  if (a->frac[ind_a]) {
4764  a->Prec = ind_a + 1;
4765  i = 0;
4766  while (a->frac[i] == 0) ++i; /* skip the first few zeros */
4767  if (i) {
4768  a->Prec -= i;
4769  if (!AddExponent(a, -(SIGNED_VALUE)i)) return 0;
4770  memmove(&a->frac[0], &a->frac[i], a->Prec*sizeof(BDIGIT));
4771  }
4772  return 1;
4773  }
4774  }
4775  /* a is zero(no non-zero digit) */
4776  VpSetZero(a, VpGetSign(a));
4777  return 0;
4778 
4779 NoVal:
4780  a->frac[0] = 0;
4781  a->Prec = 1;
4782  return 0;
4783 }
4784 
4785 /*
4786  * VpComp = 0 ... if a=b,
4787  * Pos ... a>b,
4788  * Neg ... a<b.
4789  * 999 ... result undefined(NaN)
4790  */
4791 VP_EXPORT int
4793 {
4794  int val;
4795  size_t mx, ind;
4796  int e;
4797  val = 0;
4798  if(VpIsNaN(a)||VpIsNaN(b)) return 999;
4799  if(!VpIsDef(a)) {
4800  if(!VpIsDef(b)) e = a->sign - b->sign;
4801  else e = a->sign;
4802  if(e>0) return 1;
4803  else if(e<0) return -1;
4804  else return 0;
4805  }
4806  if(!VpIsDef(b)) {
4807  e = -b->sign;
4808  if(e>0) return 1;
4809  else return -1;
4810  }
4811  /* Zero check */
4812  if(VpIsZero(a)) {
4813  if(VpIsZero(b)) return 0; /* both zero */
4814  val = -VpGetSign(b);
4815  goto Exit;
4816  }
4817  if(VpIsZero(b)) {
4818  val = VpGetSign(a);
4819  goto Exit;
4820  }
4821 
4822  /* compare sign */
4823  if(VpGetSign(a) > VpGetSign(b)) {
4824  val = 1; /* a>b */
4825  goto Exit;
4826  }
4827  if(VpGetSign(a) < VpGetSign(b)) {
4828  val = -1; /* a<b */
4829  goto Exit;
4830  }
4831 
4832  /* a and b have same sign, && signe!=0,then compare exponent */
4833  if((a->exponent) >(b->exponent)) {
4834  val = VpGetSign(a);
4835  goto Exit;
4836  }
4837  if((a->exponent) <(b->exponent)) {
4838  val = -VpGetSign(b);
4839  goto Exit;
4840  }
4841 
4842  /* a and b have same exponent, then compare significand. */
4843  mx =((a->Prec) <(b->Prec)) ?(a->Prec) :(b->Prec);
4844  ind = 0;
4845  while(ind < mx) {
4846  if((a->frac[ind]) >(b->frac[ind])) {
4847  val = VpGetSign(a);
4848  goto Exit;
4849  }
4850  if((a->frac[ind]) <(b->frac[ind])) {
4851  val = -VpGetSign(b);
4852  goto Exit;
4853  }
4854  ++ind;
4855  }
4856  if((a->Prec) >(b->Prec)) {
4857  val = VpGetSign(a);
4858  } else if((a->Prec) <(b->Prec)) {
4859  val = -VpGetSign(b);
4860  }
4861 
4862 Exit:
4863  if (val> 1) val = 1;
4864  else if(val<-1) val = -1;
4865 
4866 #ifdef BIGDECIMAL_DEBUG
4867  if(gfDebug) {
4868  VPrint(stdout, " VpComp a=%\n", a);
4869  VPrint(stdout, " b=%\n", b);
4870  printf(" ans=%d\n", val);
4871  }
4872 #endif /* BIGDECIMAL_DEBUG */
4873  return (int)val;
4874 }
4875 
4876 #ifdef BIGDECIMAL_ENABLE_VPRINT
4877 /*
4878  * cntl_chr ... ASCIIZ Character, print control characters
4879  * Available control codes:
4880  * % ... VP variable. To print '%', use '%%'.
4881  * \n ... new line
4882  * \b ... backspace
4883  * ... tab
4884  * Note: % must must not appear more than once
4885  * a ... VP variable to be printed
4886  */
4887 VP_EXPORT int
4888 VPrint(FILE *fp, const char *cntl_chr, Real *a)
4889 {
4890  size_t i, j, nc, nd, ZeroSup;
4891  BDIGIT m, e, nn;
4892 
4893  /* Check if NaN & Inf. */
4894  if(VpIsNaN(a)) {
4895  fprintf(fp,SZ_NaN);
4896  return 8;
4897  }
4898  if(VpIsPosInf(a)) {
4899  fprintf(fp,SZ_INF);
4900  return 8;
4901  }
4902  if(VpIsNegInf(a)) {
4903  fprintf(fp,SZ_NINF);
4904  return 9;
4905  }
4906  if(VpIsZero(a)) {
4907  fprintf(fp,"0.0");
4908  return 3;
4909  }
4910 
4911  j = 0;
4912  nd = nc = 0; /* nd : number of digits in fraction part(every 10 digits, */
4913  /* nd<=10). */
4914  /* nc : number of caracters printed */
4915  ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
4916  while(*(cntl_chr + j)) {
4917  if((*(cntl_chr + j) == '%') &&(*(cntl_chr + j + 1) != '%')) {
4918  nc = 0;
4919  if(!VpIsZero(a)) {
4920  if(VpGetSign(a) < 0) {
4921  fprintf(fp, "-");
4922  ++nc;
4923  }
4924  nc += fprintf(fp, "0.");
4925  for(i=0; i < a->Prec; ++i) {
4926  m = BASE1;
4927  e = a->frac[i];
4928  while(m) {
4929  nn = e / m;
4930  if((!ZeroSup) || nn) {
4931  nc += fprintf(fp, "%lu", (unsigned long)nn); /* The leading zero(s) */
4932  /* as 0.00xx will not */
4933  /* be printed. */
4934  ++nd;
4935  ZeroSup = 0; /* Set to print succeeding zeros */
4936  }
4937  if(nd >= 10) { /* print ' ' after every 10 digits */
4938  nd = 0;
4939  nc += fprintf(fp, " ");
4940  }
4941  e = e - nn * m;
4942  m /= 10;
4943  }
4944  }
4945  nc += fprintf(fp, "E%"PRIdSIZE, VpExponent10(a));
4946  } else {
4947  nc += fprintf(fp, "0.0");
4948  }
4949  } else {
4950  ++nc;
4951  if(*(cntl_chr + j) == '\\') {
4952  switch(*(cntl_chr + j + 1)) {
4953  case 'n':
4954  fprintf(fp, "\n");
4955  ++j;
4956  break;
4957  case 't':
4958  fprintf(fp, "\t");
4959  ++j;
4960  break;
4961  case 'b':
4962  fprintf(fp, "\n");
4963  ++j;
4964  break;
4965  default:
4966  fprintf(fp, "%c", *(cntl_chr + j));
4967  break;
4968  }
4969  } else {
4970  fprintf(fp, "%c", *(cntl_chr + j));
4971  if(*(cntl_chr + j) == '%') ++j;
4972  }
4973  }
4974  j++;
4975  }
4976  return (int)nc;
4977 }
4978 #endif /* BIGDECIMAL_ENABLE_VPRINT */
4979 
4980 static void
4981 VpFormatSt(char *psz, size_t fFmt)
4982 {
4983  size_t ie, i, nf = 0;
4984  char ch;
4985 
4986  if(fFmt<=0) return;
4987 
4988  ie = strlen(psz);
4989  for(i = 0; i < ie; ++i) {
4990  ch = psz[i];
4991  if(!ch) break;
4992  if(ISSPACE(ch) || ch=='-' || ch=='+') continue;
4993  if(ch == '.') { nf = 0;continue;}
4994  if(ch == 'E') break;
4995  nf++;
4996  if(nf > fFmt) {
4997  memmove(psz + i + 1, psz + i, ie - i + 1);
4998  ++ie;
4999  nf = 0;
5000  psz[i] = ' ';
5001  }
5002  }
5003 }
5004 
5005 VP_EXPORT ssize_t
5007 {
5008  ssize_t ex;
5009  size_t n;
5010 
5011  if (!VpHasVal(a)) return 0;
5012 
5013  ex = a->exponent * (ssize_t)BASE_FIG;
5014  n = BASE1;
5015  while ((a->frac[0] / n) == 0) {
5016  --ex;
5017  n /= 10;
5018  }
5019  return ex;
5020 }
5021 
5022 VP_EXPORT void
5023 VpSzMantissa(Real *a,char *psz)
5024 {
5025  size_t i, n, ZeroSup;
5026  BDIGIT_DBL m, e, nn;
5027 
5028  if(VpIsNaN(a)) {
5029  sprintf(psz,SZ_NaN);
5030  return;
5031  }
5032  if(VpIsPosInf(a)) {
5033  sprintf(psz,SZ_INF);
5034  return;
5035  }
5036  if(VpIsNegInf(a)) {
5037  sprintf(psz,SZ_NINF);
5038  return;
5039  }
5040 
5041  ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
5042  if(!VpIsZero(a)) {
5043  if(VpGetSign(a) < 0) *psz++ = '-';
5044  n = a->Prec;
5045  for (i=0; i < n; ++i) {
5046  m = BASE1;
5047  e = a->frac[i];
5048  while (m) {
5049  nn = e / m;
5050  if((!ZeroSup) || nn) {
5051  sprintf(psz, "%lu", (unsigned long)nn); /* The leading zero(s) */
5052  psz += strlen(psz);
5053  /* as 0.00xx will be ignored. */
5054  ZeroSup = 0; /* Set to print succeeding zeros */
5055  }
5056  e = e - nn * m;
5057  m /= 10;
5058  }
5059  }
5060  *psz = 0;
5061  while(psz[-1]=='0') *(--psz) = 0;
5062  } else {
5063  if(VpIsPosZero(a)) sprintf(psz, "0");
5064  else sprintf(psz, "-0");
5065  }
5066 }
5067 
5068 VP_EXPORT int
5069 VpToSpecialString(Real *a,char *psz,int fPlus)
5070 /* fPlus =0:default, =1: set ' ' before digits , =2: set '+' before digits. */
5071 {
5072  if(VpIsNaN(a)) {
5073  sprintf(psz,SZ_NaN);
5074  return 1;
5075  }
5076 
5077  if(VpIsPosInf(a)) {
5078  if(fPlus==1) {
5079  *psz++ = ' ';
5080  } else if(fPlus==2) {
5081  *psz++ = '+';
5082  }
5083  sprintf(psz,SZ_INF);
5084  return 1;
5085  }
5086  if(VpIsNegInf(a)) {
5087  sprintf(psz,SZ_NINF);
5088  return 1;
5089  }
5090  if(VpIsZero(a)) {
5091  if(VpIsPosZero(a)) {
5092  if(fPlus==1) sprintf(psz, " 0.0");
5093  else if(fPlus==2) sprintf(psz, "+0.0");
5094  else sprintf(psz, "0.0");
5095  } else sprintf(psz, "-0.0");
5096  return 1;
5097  }
5098  return 0;
5099 }
5100 
5101 VP_EXPORT void
5102 VpToString(Real *a, char *psz, size_t fFmt, int fPlus)
5103 /* fPlus =0:default, =1: set ' ' before digits , =2:set '+' before digits. */
5104 {
5105  size_t i, n, ZeroSup;
5106  BDIGIT shift, m, e, nn;
5107  char *pszSav = psz;
5108  ssize_t ex;
5109 
5110  if (VpToSpecialString(a, psz, fPlus)) return;
5111 
5112  ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */
5113 
5114  if (VpGetSign(a) < 0) *psz++ = '-';
5115  else if (fPlus == 1) *psz++ = ' ';
5116  else if (fPlus == 2) *psz++ = '+';
5117 
5118  *psz++ = '0';
5119  *psz++ = '.';
5120  n = a->Prec;
5121  for(i=0;i < n;++i) {
5122  m = BASE1;
5123  e = a->frac[i];
5124  while(m) {
5125  nn = e / m;
5126  if((!ZeroSup) || nn) {
5127  sprintf(psz, "%lu", (unsigned long)nn); /* The reading zero(s) */
5128  psz += strlen(psz);
5129  /* as 0.00xx will be ignored. */
5130  ZeroSup = 0; /* Set to print succeeding zeros */
5131  }
5132  e = e - nn * m;
5133  m /= 10;
5134  }
5135  }
5136  ex = a->exponent * (ssize_t)BASE_FIG;
5137  shift = BASE1;
5138  while(a->frac[0] / shift == 0) {
5139  --ex;
5140  shift /= 10;
5141  }
5142  while(psz[-1]=='0') *(--psz) = 0;
5143  sprintf(psz, "E%"PRIdSIZE, ex);
5144  if(fFmt) VpFormatSt(pszSav, fFmt);
5145 }
5146 
5147 VP_EXPORT void
5148 VpToFString(Real *a, char *psz, size_t fFmt, int fPlus)
5149 /* fPlus =0:default,=1: set ' ' before digits ,set '+' before digits. */
5150 {
5151  size_t i, n;
5152  BDIGIT m, e, nn;
5153  char *pszSav = psz;
5154  ssize_t ex;
5155 
5156  if(VpToSpecialString(a,psz,fPlus)) return;
5157 
5158  if(VpGetSign(a) < 0) *psz++ = '-';
5159  else if(fPlus==1) *psz++ = ' ';
5160  else if(fPlus==2) *psz++ = '+';
5161 
5162  n = a->Prec;
5163  ex = a->exponent;
5164  if(ex<=0) {
5165  *psz++ = '0';*psz++ = '.';
5166  while(ex<0) {
5167  for(i=0;i<BASE_FIG;++i) *psz++ = '0';
5168  ++ex;
5169  }
5170  ex = -1;
5171  }
5172 
5173  for(i=0;i < n;++i) {
5174  --ex;
5175  if(i==0 && ex >= 0) {
5176  sprintf(psz, "%lu", (unsigned long)a->frac[i]);
5177  psz += strlen(psz);
5178  } else {
5179  m = BASE1;
5180  e = a->frac[i];
5181  while(m) {
5182  nn = e / m;
5183  *psz++ = (char)(nn + '0');
5184  e = e - nn * m;
5185  m /= 10;
5186  }
5187  }
5188  if(ex == 0) *psz++ = '.';
5189  }
5190  while(--ex>=0) {
5191  m = BASE;
5192  while(m/=10) *psz++ = '0';
5193  if(ex == 0) *psz++ = '.';
5194  }
5195  *psz = 0;
5196  while(psz[-1]=='0') *(--psz) = 0;
5197  if(psz[-1]=='.') sprintf(psz, "0");
5198  if(fFmt) VpFormatSt(pszSav, fFmt);
5199 }
5200 
5201 /*
5202  * [Output]
5203  * a[] ... variable to be assigned the value.
5204  * [Input]
5205  * int_chr[] ... integer part(may include '+/-').
5206  * ni ... number of characters in int_chr[],not including '+/-'.
5207  * frac[] ... fraction part.
5208  * nf ... number of characters in frac[].
5209  * exp_chr[] ... exponent part(including '+/-').
5210  * ne ... number of characters in exp_chr[],not including '+/-'.
5211  */
5212 VP_EXPORT int
5213 VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne)
5214 {
5215  size_t i, j, ind_a, ma, mi, me;
5216  SIGNED_VALUE e, es, eb, ef;
5217  int sign, signe, exponent_overflow;
5218 
5219  /* get exponent part */
5220  e = 0;
5221  ma = a->MaxPrec;
5222  mi = ni;
5223  me = ne;
5224  signe = 1;
5225  exponent_overflow = 0;
5226  memset(a->frac, 0, ma * sizeof(BDIGIT));
5227  if (ne > 0) {
5228  i = 0;
5229  if (exp_chr[0] == '-') {
5230  signe = -1;
5231  ++i;
5232  ++me;
5233  }
5234  else if (exp_chr[0] == '+') {
5235  ++i;
5236  ++me;
5237  }
5238  while (i < me) {
5239  es = e * (SIGNED_VALUE)BASE_FIG;
5240  e = e * 10 + exp_chr[i] - '0';
5241  if (es > (SIGNED_VALUE)(e*BASE_FIG)) {
5242  exponent_overflow = 1;
5243  e = es; /* keep sign */
5244  break;
5245  }
5246  ++i;
5247  }
5248  }
5249 
5250  /* get integer part */
5251  i = 0;
5252  sign = 1;
5253  if(1 /*ni >= 0*/) {
5254  if(int_chr[0] == '-') {
5255  sign = -1;
5256  ++i;
5257  ++mi;
5258  } else if(int_chr[0] == '+') {
5259  ++i;
5260  ++mi;
5261  }
5262  }
5263 
5264  e = signe * e; /* e: The value of exponent part. */
5265  e = e + ni; /* set actual exponent size. */
5266 
5267  if (e > 0) signe = 1;
5268  else signe = -1;
5269 
5270  /* Adjust the exponent so that it is the multiple of BASE_FIG. */
5271  j = 0;
5272  ef = 1;
5273  while (ef) {
5274  if (e >= 0) eb = e;
5275  else eb = -e;
5276  ef = eb / (SIGNED_VALUE)BASE_FIG;
5277  ef = eb - ef * (SIGNED_VALUE)BASE_FIG;
5278  if (ef) {
5279  ++j; /* Means to add one more preceeding zero */
5280  ++e;
5281  }
5282  }
5283 
5284  eb = e / (SIGNED_VALUE)BASE_FIG;
5285 
5286  if (exponent_overflow) {
5287  int zero = 1;
5288  for ( ; i < mi && zero; i++) zero = int_chr[i] == '0';
5289  for (i = 0; i < nf && zero; i++) zero = frac[i] == '0';
5290  if (!zero && signe > 0) {
5291  VpSetInf(a, sign);
5292  VpException(VP_EXCEPTION_INFINITY, "exponent overflow",0);
5293  }
5294  else VpSetZero(a, sign);
5295  return 1;
5296  }
5297 
5298  ind_a = 0;
5299  while (i < mi) {
5300  a->frac[ind_a] = 0;
5301  while ((j < BASE_FIG) && (i < mi)) {
5302  a->frac[ind_a] = a->frac[ind_a] * 10 + int_chr[i] - '0';
5303  ++j;
5304  ++i;
5305  }
5306  if (i < mi) {
5307  ++ind_a;
5308  if (ind_a >= ma) goto over_flow;
5309  j = 0;
5310  }
5311  }
5312 
5313  /* get fraction part */
5314 
5315  i = 0;
5316  while(i < nf) {
5317  while((j < BASE_FIG) && (i < nf)) {
5318  a->frac[ind_a] = a->frac[ind_a] * 10 + frac[i] - '0';
5319  ++j;
5320  ++i;
5321  }
5322  if(i < nf) {
5323  ++ind_a;
5324  if(ind_a >= ma) goto over_flow;
5325  j = 0;
5326  }
5327  }
5328  goto Final;
5329 
5330 over_flow:
5331  rb_warn("Conversion from String to BigDecimal overflow (last few digits discarded).");
5332 
5333 Final:
5334  if (ind_a >= ma) ind_a = ma - 1;
5335  while (j < BASE_FIG) {
5336  a->frac[ind_a] = a->frac[ind_a] * 10;
5337  ++j;
5338  }
5339  a->Prec = ind_a + 1;
5340  a->exponent = eb;
5341  VpSetSign(a,sign);
5342  VpNmlz(a);
5343  return 1;
5344 }
5345 
5346 /*
5347  * [Input]
5348  * *m ... Real
5349  * [Output]
5350  * *d ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig.
5351  * *e ... exponent of m.
5352  * DBLE_FIG ... Number of digits in a double variable.
5353  *
5354  * m -> d*10**e, 0<d<BASE
5355  * [Returns]
5356  * 0 ... Zero
5357  * 1 ... Normal
5358  * 2 ... Infinity
5359  * -1 ... NaN
5360  */
5361 VP_EXPORT int
5362 VpVtoD(double *d, SIGNED_VALUE *e, Real *m)
5363 {
5364  size_t ind_m, mm, fig;
5365  double div;
5366  int f = 1;
5367 
5368  if(VpIsNaN(m)) {
5369  *d = VpGetDoubleNaN();
5370  *e = 0;
5371  f = -1; /* NaN */
5372  goto Exit;
5373  } else
5374  if(VpIsPosZero(m)) {
5375  *d = 0.0;
5376  *e = 0;
5377  f = 0;
5378  goto Exit;
5379  } else
5380  if(VpIsNegZero(m)) {
5381  *d = VpGetDoubleNegZero();
5382  *e = 0;
5383  f = 0;
5384  goto Exit;
5385  } else
5386  if(VpIsPosInf(m)) {
5387  *d = VpGetDoublePosInf();
5388  *e = 0;
5389  f = 2;
5390  goto Exit;
5391  } else
5392  if(VpIsNegInf(m)) {
5393  *d = VpGetDoubleNegInf();
5394  *e = 0;
5395  f = 2;
5396  goto Exit;
5397  }
5398  /* Normal number */
5399  fig =(DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
5400  ind_m = 0;
5401  mm = Min(fig,(m->Prec));
5402  *d = 0.0;
5403  div = 1.;
5404  while(ind_m < mm) {
5405  div /= (double)BASE;
5406  *d = *d + (double)m->frac[ind_m++] * div;
5407  }
5408  *e = m->exponent * (SIGNED_VALUE)BASE_FIG;
5409  *d *= VpGetSign(m);
5410 
5411 Exit:
5412 #ifdef BIGDECIMAL_DEBUG
5413  if(gfDebug) {
5414  VPrint(stdout, " VpVtoD: m=%\n", m);
5415  printf(" d=%e * 10 **%ld\n", *d, *e);
5416  printf(" DBLE_FIG = %d\n", DBLE_FIG);
5417  }
5418 #endif /*BIGDECIMAL_DEBUG */
5419  return f;
5420 }
5421 
5422 /*
5423  * m <- d
5424  */
5425 VP_EXPORT void
5426 VpDtoV(Real *m, double d)
5427 {
5428  size_t ind_m, mm;
5429  SIGNED_VALUE ne;
5430  BDIGIT i;
5431  double val, val2;
5432 
5433  if(isnan(d)) {
5434  VpSetNaN(m);
5435  goto Exit;
5436  }
5437  if(isinf(d)) {
5438  if(d>0.0) VpSetPosInf(m);
5439  else VpSetNegInf(m);
5440  goto Exit;
5441  }
5442 
5443  if(d == 0.0) {
5444  VpSetZero(m,1);
5445  goto Exit;
5446  }
5447  val =(d > 0.) ? d :(-d);
5448  ne = 0;
5449  if(val >= 1.0) {
5450  while(val >= 1.0) {
5451  val /= (double)BASE;
5452  ++ne;
5453  }
5454  } else {
5455  val2 = 1.0 /(double)BASE;
5456  while(val < val2) {
5457  val *= (double)BASE;
5458  --ne;
5459  }
5460  }
5461  /* Now val = 0.xxxxx*BASE**ne */
5462 
5463  mm = m->MaxPrec;
5464  memset(m->frac, 0, mm * sizeof(BDIGIT));
5465  for(ind_m = 0;val > 0.0 && ind_m < mm;ind_m++) {
5466  val *= (double)BASE;
5467  i = (BDIGIT)val;
5468  val -= (double)i;
5469  m->frac[ind_m] = i;
5470  }
5471  if(ind_m >= mm) ind_m = mm - 1;
5472  VpSetSign(m, (d > 0.0) ? 1 : -1);
5473  m->Prec = ind_m + 1;
5474  m->exponent = ne;
5475 
5476  VpInternalRound(m, 0, (m->Prec > 0) ? m->frac[m->Prec-1] : 0,
5477  (BDIGIT)(val*(double)BASE));
5478 
5479 Exit:
5480 #ifdef BIGDECIMAL_DEBUG
5481  if(gfDebug) {
5482  printf("VpDtoV d=%30.30e\n", d);
5483  VPrint(stdout, " m=%\n", m);
5484  }
5485 #endif /* BIGDECIMAL_DEBUG */
5486  return;
5487 }
5488 
5489 /*
5490  * m <- ival
5491  */
5492 #if 0 /* unused */
5493 VP_EXPORT void
5494 VpItoV(Real *m, SIGNED_VALUE ival)
5495 {
5496  size_t mm, ind_m;
5497  size_t val, v1, v2, v;
5498  int isign;
5499  SIGNED_VALUE ne;
5500 
5501  if(ival == 0) {
5502  VpSetZero(m,1);
5503  goto Exit;
5504  }
5505  isign = 1;
5506  val = ival;
5507  if(ival < 0) {
5508  isign = -1;
5509  val =(size_t)(-ival);
5510  }
5511  ne = 0;
5512  ind_m = 0;
5513  mm = m->MaxPrec;
5514  while(ind_m < mm) {
5515  m->frac[ind_m] = 0;
5516  ++ind_m;
5517  }
5518  ind_m = 0;
5519  while(val > 0) {
5520  if(val) {
5521  v1 = val;
5522  v2 = 1;
5523  while(v1 >= BASE) {
5524  v1 /= BASE;
5525  v2 *= BASE;
5526  }
5527  val = val - v2 * v1;
5528  v = v1;
5529  } else {
5530  v = 0;
5531  }
5532  m->frac[ind_m] = v;
5533  ++ind_m;
5534  ++ne;
5535  }
5536  m->Prec = ind_m - 1;
5537  m->exponent = ne;
5538  VpSetSign(m,isign);
5539  VpNmlz(m);
5540 
5541 Exit:
5542 #ifdef BIGDECIMAL_DEBUG
5543  if(gfDebug) {
5544  printf(" VpItoV i=%d\n", ival);
5545  VPrint(stdout, " m=%\n", m);
5546  }
5547 #endif /* BIGDECIMAL_DEBUG */
5548  return;
5549 }
5550 #endif
5551 
5552 /*
5553  * y = SQRT(x), y*y - x =>0
5554  */
5555 VP_EXPORT int
5557 {
5558  Real *f = NULL;
5559  Real *r = NULL;
5560  size_t y_prec;
5561  SIGNED_VALUE n, e;
5562  SIGNED_VALUE prec;
5563  ssize_t nr;
5564  double val;
5565 
5566  /* Zero, NaN or Infinity ? */
5567  if(!VpHasVal(x)) {
5568  if(VpIsZero(x)||VpGetSign(x)>0) {
5569  VpAsgn(y,x,1);
5570  goto Exit;
5571  }
5572  VpSetNaN(y);
5573  return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(NaN or negative value)",0);
5574  goto Exit;
5575  }
5576 
5577  /* Negative ? */
5578  if(VpGetSign(x) < 0) {
5579  VpSetNaN(y);
5580  return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(negative value)",0);
5581  }
5582 
5583  /* One ? */
5584  if(VpIsOne(x)) {
5585  VpSetOne(y);
5586  goto Exit;
5587  }
5588 
5589  n = (SIGNED_VALUE)y->MaxPrec;
5590  if (x->MaxPrec > (size_t)n) n = (ssize_t)x->MaxPrec;
5591  /* allocate temporally variables */
5592  f = VpAlloc(y->MaxPrec * (BASE_FIG + 2), "#1");
5593  r = VpAlloc((n + n) * (BASE_FIG + 2), "#1");
5594 
5595  nr = 0;
5596  y_prec = y->MaxPrec;
5597 
5598  prec = x->exponent - (ssize_t)y_prec;
5599  if (x->exponent > 0)
5600  ++prec;
5601  else
5602  --prec;
5603 
5604  VpVtoD(&val, &e, x); /* val <- x */
5605  e /= (SIGNED_VALUE)BASE_FIG;
5606  n = e / 2;
5607  if (e - n * 2 != 0) {
5608  val /= BASE;
5609  n = (e + 1) / 2;
5610  }
5611  VpDtoV(y, sqrt(val)); /* y <- sqrt(val) */
5612  y->exponent += n;
5613  n = (SIGNED_VALUE)((DBLE_FIG + BASE_FIG - 1) / BASE_FIG);
5614  y->MaxPrec = Min((size_t)n , y_prec);
5615  f->MaxPrec = y->MaxPrec + 1;
5616  n = (SIGNED_VALUE)(y_prec * BASE_FIG);
5617  if (n < (SIGNED_VALUE)maxnr) n = (SIGNED_VALUE)maxnr;
5618  do {
5619  y->MaxPrec *= 2;
5620  if (y->MaxPrec > y_prec) y->MaxPrec = y_prec;
5621  f->MaxPrec = y->MaxPrec;
5622  VpDivd(f, r, x, y); /* f = x/y */
5623  VpAddSub(r, f, y, -1); /* r = f - y */
5624  VpMult(f, VpPt5, r); /* f = 0.5*r */
5625  if(VpIsZero(f)) goto converge;
5626  VpAddSub(r, f, y, 1); /* r = y + f */
5627  VpAsgn(y, r, 1); /* y = r */
5628  if(f->exponent <= prec) goto converge;
5629  } while(++nr < n);
5630  /* */
5631 #ifdef BIGDECIMAL_DEBUG
5632  if(gfDebug) {
5633  printf("ERROR(VpSqrt): did not converge within %ld iterations.\n",
5634  nr);
5635  }
5636 #endif /* BIGDECIMAL_DEBUG */
5637  y->MaxPrec = y_prec;
5638 
5639 converge:
5640  VpChangeSign(y, 1);
5641 #ifdef BIGDECIMAL_DEBUG
5642  if(gfDebug) {
5643  VpMult(r, y, y);
5644  VpAddSub(f, x, r, -1);
5645  printf("VpSqrt: iterations = %"PRIdSIZE"\n", nr);
5646  VPrint(stdout, " y =% \n", y);
5647  VPrint(stdout, " x =% \n", x);
5648  VPrint(stdout, " x-y*y = % \n", f);
5649  }
5650 #endif /* BIGDECIMAL_DEBUG */
5651  y->MaxPrec = y_prec;
5652 
5653 Exit:
5654  VpFree(f);
5655  VpFree(r);
5656  return 1;
5657 }
5658 
5659 /*
5660  *
5661  * nf: digit position for operation.
5662  *
5663  */
5664 VP_EXPORT int
5665 VpMidRound(Real *y, unsigned short f, ssize_t nf)
5666 /*
5667  * Round reletively from the decimal point.
5668  * f: rounding mode
5669  * nf: digit location to round from the the decimal point.
5670  */
5671 {
5672  /* fracf: any positive digit under rounding position? */
5673  /* fracf_1further: any positive digits under one further than the rounding position? */
5674  /* exptoadd: number of digits needed to compensate negative nf */
5675  int fracf, fracf_1further;
5676  ssize_t n,i,ix,ioffset, exptoadd;
5677  BDIGIT v, shifter;
5678  BDIGIT div;
5679 
5680  nf += y->exponent * (ssize_t)BASE_FIG;
5681  exptoadd=0;
5682  if (nf < 0) {
5683  /* rounding position too left(large). */
5684  if((f!=VP_ROUND_CEIL) && (f!=VP_ROUND_FLOOR)) {
5685  VpSetZero(y,VpGetSign(y)); /* truncate everything */
5686  return 0;
5687  }
5688  exptoadd = -nf;
5689  nf = 0;
5690  }
5691 
5692  ix = nf / (ssize_t)BASE_FIG;
5693  if ((size_t)ix >= y->Prec) return 0; /* rounding position too right(small). */
5694  v = y->frac[ix];
5695 
5696  ioffset = nf - ix*(ssize_t)BASE_FIG;
5697  n = (ssize_t)BASE_FIG - ioffset - 1;
5698  for (shifter=1,i=0; i<n; ++i) shifter *= 10;
5699 
5700  /* so the representation used (in y->frac) is an array of BDIGIT, where
5701  each BDIGIT contains a value between 0 and BASE-1, consisting of BASE_FIG
5702  decimal places.
5703 
5704  (that numbers of decimal places are typed as ssize_t is somewhat confusing)
5705 
5706  nf is now position (in decimal places) of the digit from the start of
5707  the array.
5708  ix is the position (in BDIGITS) of the BDIGIT containing the decimal digit,
5709  from the start of the array.
5710  v is the value of this BDIGIT
5711  ioffset is the number of extra decimal places along of this decimal digit
5712  within v.
5713  n is the number of decimal digits remaining within v after this decimal digit
5714  shifter is 10**n,
5715  v % shifter are the remaining digits within v
5716  v % (shifter * 10) are the digit together with the remaining digits within v
5717  v / shifter are the digit's predecessors together with the digit
5718  div = v / shifter / 10 is just the digit's precessors
5719  (v / shifter) - div*10 is just the digit, which is what v ends up being reassigned to.
5720  */
5721 
5722  fracf = (v % (shifter * 10) > 0);
5723  fracf_1further = ((v % shifter) > 0);
5724 
5725  v /= shifter;
5726  div = v / 10;
5727  v = v - div*10;
5728  /* now v is just the digit required.
5729  now fracf is whether the digit or any of the remaining digits within v are non-zero
5730  now fracf_1further is whether any of the remaining digits within v are non-zero
5731  */
5732 
5733  /* now check all the remaining BDIGITS for zero-ness a whole BDIGIT at a time.
5734  if we spot any non-zeroness, that means that we foudn a positive digit under
5735  rounding position, and we also found a positive digit under one further than
5736  the rounding position, so both searches (to see if any such non-zero digit exists)
5737  can stop */
5738 
5739  for (i=ix+1; (size_t)i < y->Prec; i++) {
5740  if (y->frac[i] % BASE) {
5741  fracf = fracf_1further = 1;
5742  break;
5743  }
5744  }
5745 
5746  /* now fracf = does any positive digit exist under the rounding position?
5747  now fracf_1further = does any positive digit exist under one further than the
5748  rounding position?
5749  now v = the first digit under the rounding position */
5750 
5751  /* drop digits after pointed digit */
5752  memset(y->frac+ix+1, 0, (y->Prec - (ix+1)) * sizeof(BDIGIT));
5753 
5754  switch(f) {
5755  case VP_ROUND_DOWN: /* Truncate */
5756  break;
5757  case VP_ROUND_UP: /* Roundup */
5758  if (fracf) ++div;
5759  break;
5760  case VP_ROUND_HALF_UP:
5761  if (v>=5) ++div;
5762  break;
5763  case VP_ROUND_HALF_DOWN:
5764  if (v > 5 || (v == 5 && fracf_1further)) ++div;
5765  break;
5766  case VP_ROUND_CEIL:
5767  if (fracf && (VpGetSign(y)>0)) ++div;
5768  break;
5769  case VP_ROUND_FLOOR:
5770  if (fracf && (VpGetSign(y)<0)) ++div;
5771  break;
5772  case VP_ROUND_HALF_EVEN: /* Banker's rounding */
5773  if (v > 5) ++div;
5774  else if (v == 5) {
5775  if (fracf_1further) {
5776  ++div;
5777  }
5778  else {
5779  if (ioffset == 0) {
5780  /* v is the first decimal digit of its BDIGIT;
5781  need to grab the previous BDIGIT if present
5782  to check for evenness of the previous decimal
5783  digit (which is same as that of the BDIGIT since
5784  base 10 has a factor of 2) */
5785  if (ix && (y->frac[ix-1] % 2)) ++div;
5786  }
5787  else {
5788  if (div % 2) ++div;
5789  }
5790  }
5791  }
5792  break;
5793  }
5794  for (i=0; i<=n; ++i) div *= 10;
5795  if (div>=BASE) {
5796  if(ix) {
5797  y->frac[ix] = 0;
5798  VpRdup(y,ix);
5799  } else {
5800  short s = VpGetSign(y);
5801  SIGNED_VALUE e = y->exponent;
5802  VpSetOne(y);
5803  VpSetSign(y, s);
5804  y->exponent = e+1;
5805  }
5806  } else {
5807  y->frac[ix] = div;
5808  VpNmlz(y);
5809  }
5810  if (exptoadd > 0) {
5811  y->exponent += (SIGNED_VALUE)(exptoadd/BASE_FIG);
5812  exptoadd %= (ssize_t)BASE_FIG;
5813  for(i=0;i<exptoadd;i++) {
5814  y->frac[0] *= 10;
5815  if (y->frac[0] >= BASE) {
5816  y->frac[0] /= BASE;
5817  y->exponent++;
5818  }
5819  }
5820  }
5821  return 1;
5822 }
5823 
5824 VP_EXPORT int
5825 VpLeftRound(Real *y, unsigned short f, ssize_t nf)
5826 /*
5827  * Round from the left hand side of the digits.
5828  */
5829 {
5830  BDIGIT v;
5831  if (!VpHasVal(y)) return 0; /* Unable to round */
5832  v = y->frac[0];
5833  nf -= VpExponent(y)*(ssize_t)BASE_FIG;
5834  while ((v /= 10) != 0) nf--;
5835  nf += (ssize_t)BASE_FIG-1;
5836  return VpMidRound(y,f,nf);
5837 }
5838 
5839 VP_EXPORT int
5840 VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t nf)
5841 {
5842  /* First,assign whole value in truncation mode */
5843  if (VpAsgn(y, x, 10) <= 1) return 0; /* Zero,NaN,or Infinity */
5844  return VpMidRound(y,f,nf);
5845 }
5846 
5847 static int
5848 VpLimitRound(Real *c, size_t ixDigit)
5849 {
5850  size_t ix = VpGetPrecLimit();
5851  if(!VpNmlz(c)) return -1;
5852  if(!ix) return 0;
5853  if(!ixDigit) ixDigit = c->Prec-1;
5854  if((ix+BASE_FIG-1)/BASE_FIG > ixDigit+1) return 0;
5855  return VpLeftRound(c, VpGetRoundMode(), (ssize_t)ix);
5856 }
5857 
5858 /* If I understand correctly, this is only ever used to round off the final decimal
5859  digit of precision */
5860 static void
5861 VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v)
5862 {
5863  int f = 0;
5864 
5865  unsigned short const rounding_mode = VpGetRoundMode();
5866 
5867  if (VpLimitRound(c, ixDigit)) return;
5868  if (!v) return;
5869 
5870  v /= BASE1;
5871  switch (rounding_mode) {
5872  case VP_ROUND_DOWN:
5873  break;
5874  case VP_ROUND_UP:
5875  if (v) f = 1;
5876  break;
5877  case VP_ROUND_HALF_UP:
5878  if (v >= 5) f = 1;
5879  break;
5880  case VP_ROUND_HALF_DOWN:
5881  /* this is ok - because this is the last digit of precision,
5882  the case where v == 5 and some further digits are nonzero
5883  will never occur */
5884  if (v >= 6) f = 1;
5885  break;
5886  case VP_ROUND_CEIL:
5887  if (v && (VpGetSign(c) > 0)) f = 1;
5888  break;
5889  case VP_ROUND_FLOOR:
5890  if (v && (VpGetSign(c) < 0)) f = 1;
5891  break;
5892  case VP_ROUND_HALF_EVEN: /* Banker's rounding */
5893  /* as per VP_ROUND_HALF_DOWN, because this is the last digit of precision,
5894  there is no case to worry about where v == 5 and some further digits are nonzero */
5895  if (v > 5) f = 1;
5896  else if (v == 5 && vPrev % 2) f = 1;
5897  break;
5898  }
5899  if (f) {
5900  VpRdup(c, ixDigit);
5901  VpNmlz(c);
5902  }
5903 }
5904 
5905 /*
5906  * Rounds up m(plus one to final digit of m).
5907  */
5908 static int
5909 VpRdup(Real *m, size_t ind_m)
5910 {
5911  BDIGIT carry;
5912 
5913  if (!ind_m) ind_m = m->Prec;
5914 
5915  carry = 1;
5916  while (carry > 0 && (ind_m--)) {
5917  m->frac[ind_m] += carry;
5918  if (m->frac[ind_m] >= BASE) m->frac[ind_m] -= BASE;
5919  else carry = 0;
5920  }
5921  if(carry > 0) { /* Overflow,count exponent and set fraction part be 1 */
5922  if (!AddExponent(m, 1)) return 0;
5923  m->Prec = m->frac[0] = 1;
5924  } else {
5925  VpNmlz(m);
5926  }
5927  return 1;
5928 }
5929 
5930 /*
5931  * y = x - fix(x)
5932  */
5933 VP_EXPORT void
5935 {
5936  size_t my, ind_y, ind_x;
5937 
5938  if(!VpHasVal(x)) {
5939  VpAsgn(y,x,1);
5940  goto Exit;
5941  }
5942 
5943  if (x->exponent > 0 && (size_t)x->exponent >= x->Prec) {
5944  VpSetZero(y,VpGetSign(x));
5945  goto Exit;
5946  }
5947  else if(x->exponent <= 0) {
5948  VpAsgn(y, x, 1);
5949  goto Exit;
5950  }
5951 
5952  /* satisfy: x->exponent > 0 */
5953 
5954  y->Prec = x->Prec - (size_t)x->exponent;
5955  y->Prec = Min(y->Prec, y->MaxPrec);
5956  y->exponent = 0;
5957  VpSetSign(y,VpGetSign(x));
5958  ind_y = 0;
5959  my = y->Prec;
5960  ind_x = x->exponent;
5961  while(ind_y < my) {
5962  y->frac[ind_y] = x->frac[ind_x];
5963  ++ind_y;
5964  ++ind_x;
5965  }
5966  VpNmlz(y);
5967 
5968 Exit:
5969 #ifdef BIGDECIMAL_DEBUG
5970  if(gfDebug) {
5971  VPrint(stdout, "VpFrac y=%\n", y);
5972  VPrint(stdout, " x=%\n", x);
5973  }
5974 #endif /* BIGDECIMAL_DEBUG */
5975  return;
5976 }
5977 
5978 /*
5979  * y = x ** n
5980  */
5981 VP_EXPORT int
5983 {
5984  size_t s, ss;
5985  ssize_t sign;
5986  Real *w1 = NULL;
5987  Real *w2 = NULL;
5988 
5989  if(VpIsZero(x)) {
5990  if(n==0) {
5991  VpSetOne(y);
5992  goto Exit;
5993  }
5994  sign = VpGetSign(x);
5995  if(n<0) {
5996  n = -n;
5997  if(sign<0) sign = (n%2)?(-1):(1);
5998  VpSetInf (y,sign);
5999  } else {
6000  if(sign<0) sign = (n%2)?(-1):(1);
6001  VpSetZero(y,sign);
6002  }
6003  goto Exit;
6004  }
6005  if(VpIsNaN(x)) {
6006  VpSetNaN(y);
6007  goto Exit;
6008  }
6009  if(VpIsInf(x)) {
6010  if(n==0) {
6011  VpSetOne(y);
6012  goto Exit;
6013  }
6014  if(n>0) {
6015  VpSetInf(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1);
6016  goto Exit;
6017  }
6018  VpSetZero(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1);
6019  goto Exit;
6020  }
6021 
6022  if((x->exponent == 1) &&(x->Prec == 1) &&(x->frac[0] == 1)) {
6023  /* abs(x) = 1 */
6024  VpSetOne(y);
6025  if(VpGetSign(x) > 0) goto Exit;
6026  if((n % 2) == 0) goto Exit;
6027  VpSetSign(y, -1);
6028  goto Exit;
6029  }
6030 
6031  if(n > 0) sign = 1;
6032  else if(n < 0) {
6033  sign = -1;
6034  n = -n;
6035  } else {
6036  VpSetOne(y);
6037  goto Exit;
6038  }
6039 
6040  /* Allocate working variables */
6041 
6042  w1 = VpAlloc((y->MaxPrec + 2) * BASE_FIG, "#0");
6043  w2 = VpAlloc((w1->MaxPrec * 2 + 1) * BASE_FIG, "#0");
6044  /* calculation start */
6045 
6046  VpAsgn(y, x, 1);
6047  --n;
6048  while(n > 0) {
6049  VpAsgn(w1, x, 1);
6050  s = 1;
6051  while (ss = s, (s += s) <= (size_t)n) {
6052  VpMult(w2, w1, w1);
6053  VpAsgn(w1, w2, 1);
6054  }
6055  n -= (SIGNED_VALUE)ss;
6056  VpMult(w2, y, w1);
6057  VpAsgn(y, w2, 1);
6058  }
6059  if(sign < 0) {
6060  VpDivd(w1, w2, VpConstOne, y);
6061  VpAsgn(y, w1, 1);
6062  }
6063 
6064 Exit:
6065 #ifdef BIGDECIMAL_DEBUG
6066  if(gfDebug) {
6067  VPrint(stdout, "VpPower y=%\n", y);
6068  VPrint(stdout, "VpPower x=%\n", x);
6069  printf(" n=%d\n", n);
6070  }
6071 #endif /* BIGDECIMAL_DEBUG */
6072  VpFree(w2);
6073  VpFree(w1);
6074  return 1;
6075 }
6076 
6077 #ifdef BIGDECIMAL_DEBUG
6078 int
6079 VpVarCheck(Real * v)
6080 /*
6081  * Checks the validity of the Real variable v.
6082  * [Input]
6083  * v ... Real *, variable to be checked.
6084  * [Returns]
6085  * 0 ... correct v.
6086  * other ... error
6087  */
6088 {
6089  size_t i;
6090 
6091  if(v->MaxPrec <= 0) {
6092  printf("ERROR(VpVarCheck): Illegal Max. Precision(=%"PRIuSIZE")\n",
6093  v->MaxPrec);
6094  return 1;
6095  }
6096  if((v->Prec <= 0) ||((v->Prec) >(v->MaxPrec))) {
6097  printf("ERROR(VpVarCheck): Illegal Precision(=%"PRIuSIZE")\n", v->Prec);
6098  printf(" Max. Prec.=%"PRIuSIZE"\n", v->MaxPrec);
6099  return 2;
6100  }
6101  for(i = 0; i < v->Prec; ++i) {
6102  if((v->frac[i] >= BASE)) {
6103  printf("ERROR(VpVarCheck): Illegal fraction\n");
6104  printf(" Frac[%"PRIuSIZE"]=%lu\n", i, v->frac[i]);
6105  printf(" Prec. =%"PRIuSIZE"\n", v->Prec);
6106  printf(" Exp. =%"PRIdVALUE"\n", v->exponent);
6107  printf(" BASE =%lu\n", BASE);
6108  return 3;
6109  }
6110  }
6111  return 0;
6112 }
6113 #endif /* BIGDECIMAL_DEBUG */
VP_EXPORT size_t VpDivd(Real *c, Real *r, Real *a, Real *b)
Definition: bigdecimal.c:4543
VP_EXPORT void VpToString(Real *a, char *psz, size_t fFmt, int fPlus)
Definition: bigdecimal.c:5102
VP_EXPORT double VpGetDoublePosInf(void)
Definition: bigdecimal.c:3502
#define T_SYMBOL
Definition: ruby.h:502
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
static BDIGIT VpAddAbs(Real *a, Real *b, Real *c)
Definition: bigdecimal.c:4120
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:4323
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 VALUE BigDecimal_sign(VALUE self)
Definition: bigdecimal.c:2546
VP_EXPORT int VpException(unsigned short f, const char *str, int always)
Definition: bigdecimal.c:3535
void rb_bug(const char *fmt,...)
Definition: error.c:295
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:4981
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 *)
#define INT2NUM(x)
Definition: ruby.h:1178
int i
Definition: win32ole.c:784
#define VP_ROUND_HALF_UP
Definition: bigdecimal.h:110
#define T_FIXNUM
Definition: ruby.h:497
static VALUE BigDecimal_to_r(VALUE self)
Definition: bigdecimal.c:749
VALUE rb_Rational(VALUE, VALUE)
Definition: rational.c:1755
#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
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
if(dispIdMember==DISPID_VALUE)
Definition: win32ole.c:791
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:67
VALUE obj
Definition: bigdecimal.h:137
VP_EXPORT int VpIsRoundMode(unsigned short n)
Definition: bigdecimal.c:3429
static double Zero(void)
Definition: bigdecimal.c:3468
#define VP_ROUND_FLOOR
Definition: bigdecimal.h:113
#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:3343
#define VpIsDef(a)
Definition: bigdecimal.h:273
#define Qtrue
Definition: ruby.h:434
#define VP_SIGN_NEGATIVE_ZERO
Definition: bigdecimal.h:120
static VALUE rb_float_new(double d)
Definition: ruby.h:790
#define TypedData_Wrap_Struct(klass, data_type, sval)
Definition: ruby.h:1016
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
static VALUE ToValue(Real *p)
Definition: bigdecimal.c:176
static VALUE BigMath_s_log(VALUE, VALUE, VALUE)
Definition: bigdecimal.c:2774
static ID id_half_down
Definition: bigdecimal.c:53
VALUE rb_eTypeError
Definition: error.c:516
#define T_RATIONAL
Definition: ruby.h:503
#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:2511
#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:5909
RUBY_EXTERN VALUE rb_eMathDomainError
Definition: ruby.h:1498
#define SYM2ID(x)
Definition: ruby.h:364
st_index_t rb_memhash(const void *ptr, long len)
Definition: random.c:1422
static VALUE BigMath_s_exp(VALUE, VALUE, VALUE)
Definition: bigdecimal.c:2642
VP_EXPORT unsigned short VpGetRoundMode(void)
Definition: bigdecimal.c:3413
VP_EXPORT unsigned short VpSetRoundMode(unsigned short n)
Definition: bigdecimal.c:3447
VALUE rb_funcall(VALUE, ID, int,...)
Calls a method.
Definition: vm_eval.c:773
static VALUE BigDecimal_mod(VALUE self, VALUE r)
Definition: bigdecimal.c:1355
VP_EXPORT void VpDtoV(Real *m, double d)
Definition: bigdecimal.c:5426
VALUE rb_protect(VALUE(*proc)(VALUE), VALUE data, int *state)
Definition: eval.c:771
#define is_positive(x)
Definition: bigdecimal.c:2049
#define PRIxVALUE
Definition: ruby.h:145
static void cannot_be_coerced_into_BigDecimal(VALUE exc_class, VALUE v)
Definition: bigdecimal.c:191
#define Check_Type(v, t)
Definition: ruby.h:539
void rb_raise(VALUE exc, const char *fmt,...)
Definition: error.c:1788
VP_EXPORT void VpFrac(Real *y, Real *x)
Definition: bigdecimal.c:5934
#define VP_SIGN_POSITIVE_ZERO
Definition: bigdecimal.h:119
static VALUE BigDecimal_sqrt(VALUE self, VALUE nFig)
Definition: bigdecimal.c:1578
#define RB_GC_GUARD(v)
Definition: ruby.h:530
void rb_define_alloc_func(VALUE, rb_alloc_func_t)
static void VpSetException(unsigned short f)
Definition: bigdecimal.c:3359
#define DATA_PTR(dta)
Definition: ruby.h:985
#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
st_data_t st_index_t
Definition: st.h:63
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:5362
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 FIXNUM_P(f)
Definition: ruby.h:355
#define SSIZET2NUM(v)
Definition: ruby.h:271
static Real * VpConstOne
Definition: bigdecimal.c:3266
static VALUE BigDecimal_ceil(int argc, VALUE *argv, VALUE self)
Definition: bigdecimal.c:1794
VALUE rb_str_tmp_new(long)
static void BigDecimal_delete(void *pv)
Definition: bigdecimal.c:152
RUBY_EXTERN VALUE rb_eZeroDivError
Definition: ruby.h:1482
#define DoSomeOne(x, y, f)
Definition: bigdecimal.c:122
VP_EXPORT int VpComp(Real *a, Real *b)
Definition: bigdecimal.c:4792
#define BDIGIT
Definition: defines.h:93
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:5665
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:5213
VP_EXPORT int VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t nf)
Definition: bigdecimal.c:5840
RUBY_EXTERN void * memmove(void *, const void *, size_t)
Definition: memmove.c:7
#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
VALUE rb_thread_local_aref(VALUE, ID)
Definition: thread.c:2671
#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:64
static ID id_half_up
Definition: bigdecimal.c:51
Win32OLEIDispatch * p
Definition: win32ole.c:786
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
#define RB_TYPE_P(obj, type)
Definition: ruby.h:1537
static size_t BigDecimal_memsize(const void *ptr)
Definition: bigdecimal.c:158
#define RRATIONAL(obj)
Definition: ruby.h:1108
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
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:2596
VALUE rb_class_name(VALUE)
Definition: variable.c:383
static VALUE BigDecimal_nonzero(VALUE self)
Definition: bigdecimal.c:1046
static int VpLimitRound(Real *c, size_t ixDigit)
Definition: bigdecimal.c:5848
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 ne(x, y)
Definition: time.c:66
static int VpIsDefOP(Real *c, Real *a, Real *b, int sw)
Definition: bigdecimal.c:3575
#define SZ_NINF
Definition: bigdecimal.h:84
#define PRIdVALUE
Definition: ruby.h:142
#define BDIGIT_DBL_SIGNED
Definition: defines.h:96
#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:478
VALUE rb_str_cat2(VALUE, const char *)
Definition: string.c:1986
static VALUE BigDecimal_to_f(VALUE self)
Definition: bigdecimal.c:702
#define VP_EXPORT
Definition: bigdecimal.h:90
#define rmpd_set_thread_local_exception_mode(mode)
Definition: bigdecimal.c:3335
#define snprintf
Definition: subst.h:6
VALUE rb_thread_current(void)
Definition: thread.c:2358
#define HALF_BASE
Definition: bigdecimal.c:71
#define NIL_P(v)
Definition: ruby.h:446
VP_EXPORT double VpGetDoubleNaN(void)
Definition: bigdecimal.c:3494
VP_EXPORT size_t VpGetPrecLimit(void)
Definition: bigdecimal.c:3378
VALUE rb_define_class(const char *name, VALUE super)
Defines a top-level class.
Definition: class.c:488
#define VpIsOne(a)
Definition: bigdecimal.h:278
static VALUE BigDecimal_div(VALUE self, VALUE r)
Definition: bigdecimal.c:1243
VP_EXPORT void * VpMemAlloc(size_t mb)
Definition: bigdecimal.c:3289
void rb_define_const(VALUE, const char *, VALUE)
Definition: variable.c:2204
#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:4418
static double one(void)
Definition: isinf.c:52
#define ENTER(n)
Definition: bigdecimal.c:63
#define T_FLOAT
Definition: ruby.h:489
static int VpNmlz(Real *a)
Definition: bigdecimal.c:4754
#define TYPE(x)
Definition: ruby.h:513
int argc
Definition: ruby.c:130
volatile const double gOne_ABCED9B4_CE73__00400511F31D
Definition: bigdecimal.c:3466
static int is_kind_of_BigDecimal(VALUE const v)
Definition: bigdecimal.c:170
#define RBIGNUM_LEN(b)
Definition: ruby.h:1081
#define Qfalse
Definition: ruby.h:433
static VALUE BigDecimal_add2(VALUE self, VALUE b, VALUE n)
Definition: bigdecimal.c:1499
VP_EXPORT void VpToFString(Real *a, char *psz, size_t fFmt, int fPlus)
Definition: bigdecimal.c:5148
#define T_BIGNUM
Definition: ruby.h:495
static VALUE BigDecimal_initialize_copy(VALUE self, VALUE other)
Definition: bigdecimal.c:2427
#define MEMCPY(p1, p2, type, n)
Definition: ruby.h:1242
RUBY_EXTERN int isinf(double)
Definition: isinf.c:56
VP_EXPORT Real * VpOne(void)
Definition: bigdecimal.c:3747
VP_EXPORT void VpSzMantissa(Real *a, char *psz)
Definition: bigdecimal.c:5023
#define T_COMPLEX
Definition: ruby.h:504
#define VpDblFig()
Definition: bigdecimal.h:178
static VALUE int_chr(int argc, VALUE *argv, VALUE num)
Definition: numeric.c:2478
VALUE rb_big2str(VALUE x, int base)
Definition: bignum.c:1159
#define RBIGNUM_DIGITS(b)
Definition: ruby.h:1087
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:3313
void Init_bigdecimal(void)
Definition: bigdecimal.c:3032
VALUE rb_str_resize(VALUE, long)
Definition: string.c:1854
static VALUE BigDecimal_add(VALUE self, VALUE r)
Definition: bigdecimal.c:850
#define VpSetInf(a, s)
Definition: bigdecimal.h:276
#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
static void VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v)
Definition: bigdecimal.c:5861
volatile const double gZero_ABCED9B1_CE73__00400511F31D
Definition: bigdecimal.c:3465
#define RSTRING_LEN(str)
Definition: ruby.h:862
VALUE rb_yield(VALUE)
Definition: vm_eval.c:933
#define VpSetPosInf(a)
Definition: bigdecimal.h:274
int errno
#define T_DATA
Definition: ruby.h:500
static const rb_data_type_t BigDecimal_data_type
Definition: bigdecimal.c:164
#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:1842
#define rb_Rational1(x)
Definition: intern.h:158
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:3797
#define rmpd_set_thread_local_rounding_mode(mode)
Definition: bigdecimal.c:3405
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 PRIsVALUE
Definition: ruby.h:147
unsigned long ID
Definition: ruby.h:105
static VALUE BigDecimal_global_new(int argc, VALUE *argv, VALUE self)
Definition: bigdecimal.c:2487
#define Qnil
Definition: ruby.h:435
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
static int is_integer(VALUE x)
Definition: bigdecimal.c:2029
#define MemCmp(x, y, z)
Definition: bigdecimal.c:3272
unsigned long VALUE
Definition: ruby.h:104
#define VP_SIGN_NaN
Definition: bigdecimal.h:118
static ID id_BigDecimal_precision_limit
Definition: bigdecimal.c:46
#define StrCmp(x, y)
Definition: bigdecimal.c:3273
VP_EXPORT size_t VpInit(BDIGIT BaseVal)
Definition: bigdecimal.c:3716
static VALUE rmpd_power_by_big_decimal(Real const *x, Real const *exp, ssize_t const n)
Definition: bigdecimal.c:2117
#define FIX2INT(x)
Definition: ruby.h:624
#define BDIGIT_DBL
Definition: defines.h:95
VALUE rb_num_coerce_relop(VALUE, VALUE, ID)
Definition: numeric.c:299
#define maxnr
Definition: bigdecimal.c:3268
#define isnan(x)
Definition: win32.h:327
static void shift(struct cparse_params *v, long act, VALUE tok, VALUE val)
Definition: cparse.c:662
RUBY_EXTERN VALUE rb_cNumeric
Definition: ruby.h:1448
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
VALUE rb_str_dup(VALUE)
Definition: string.c:946
static VALUE BigDecimal_s_allocate(VALUE klass)
Definition: bigdecimal.c:2377
VP_EXPORT double VpGetDoubleNegInf(void)
Definition: bigdecimal.c:3510
void xfree(void *)
#define LONG2NUM(x)
Definition: ruby.h:1199
static VALUE BigDecimal_zero(VALUE self)
Definition: bigdecimal.c:1038
#define StringValueCStr(v)
Definition: ruby.h:548
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:3950
#define RSTRING_PTR(str)
Definition: ruby.h:866
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 RMPD_PRECISION_LIMIT_DEFAULT
Definition: bigdecimal.c:3374
static VALUE BigDecimal_gt(VALUE self, VALUE r)
Definition: bigdecimal.c:1111
#define rmpd_set_thread_local_precision_limit(limit)
Definition: bigdecimal.c:3368
static VALUE BigDecimal_mult(VALUE self, VALUE r)
Definition: bigdecimal.c:1165
static int is_negative(VALUE x)
Definition: bigdecimal.c:2035
#define RFLOAT_VALUE(v)
Definition: ruby.h:836
#define f
#define INT2FIX(i)
Definition: ruby.h:241
#define VpChangeSign(a, s)
Definition: bigdecimal.h:250
void * rb_check_typeddata(VALUE obj, const rb_data_type_t *data_type)
Definition: error.c:488
#define NUM2SSIZET(x)
Definition: ruby.h:671
VALUE rb_exc_new3(VALUE etype, VALUE str)
Definition: error.c:553
RUBY_EXTERN double round(double)
Definition: numeric.c:92
#define RB_OBJ_CLASSNAME(obj)
Definition: bigdecimal.c:111
#define xmalloc
Definition: defines.h:64
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:2439
#define NUM2SIZET(x)
Definition: ruby.h:670
VP_EXPORT int VpToSpecialString(Real *a, char *psz, int fPlus)
Definition: bigdecimal.c:5069
#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:5825
#define RARRAY_PTR(a)
Definition: ruby.h:904
VP_EXPORT size_t VpAddSub(Real *c, Real *a, Real *b, int operation)
Definition: bigdecimal.c:3992
#define VpIsZero(a)
Definition: bigdecimal.h:260
#define DBL_MIN_10_EXP
Definition: numeric.c:61
static double One(void)
Definition: bigdecimal.c:3474
static VALUE BigDecimal_power_op(VALUE self, VALUE exp)
Definition: bigdecimal.c:2371
static int rb_special_const_p(VALUE obj)
Definition: ruby.h:1560
#define VpSetOne(a)
Definition: bigdecimal.h:255
#define SZ_NaN
Definition: bigdecimal.h:81
#define RTEST(v)
Definition: ruby.h:445
static int is_even(VALUE x)
Definition: bigdecimal.c:2100
#define T_STRING
Definition: ruby.h:490
#define DBLE_FIG
Definition: bigdecimal.c:75
static VALUE BigDecimal_sub2(VALUE self, VALUE b, VALUE n)
Definition: bigdecimal.c:1516
#define PRIuSIZE
Definition: ruby.h:189
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:5982
static unsigned int hash(const char *str, unsigned int len)
Definition: lex.c:56
#define SafeStringValue(v)
Definition: ruby.h:552
#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
#define BASE1
Definition: bigdecimal.c:72
#define PRIdSIZE
Definition: ruby.h:186
static Real * VpPt5
Definition: bigdecimal.c:3267
#define assert(condition)
Definition: ossl.h:45
#define VpSetZero(a, s)
Definition: bigdecimal.h:263
#define VpSetNegInf(a)
Definition: bigdecimal.h:275
#define xrealloc
Definition: defines.h:67
#define VP_ROUND_HALF_EVEN
Definition: bigdecimal.h:114
VP_EXPORT size_t VpSetPrecLimit(size_t n)
Definition: bigdecimal.c:3394
RUBY_EXTERN VALUE rb_eFloatDomainError
Definition: ruby.h:1486
#define RRATIONAL_ZERO_P(x)
Definition: bigdecimal.c:97
#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:4211
#define SZ_PINF
Definition: bigdecimal.h:83
VALUE rb_inspect(VALUE)
Definition: object.c:411
#define FIX2UINT(x)
Definition: ruby.h:625
VP_EXPORT size_t VpNumOfChars(Real *vp, const char *pszFmt)
Definition: bigdecimal.c:3672
#define VpHasVal(a)
Definition: bigdecimal.h:277
#define rb_intern_const(str)
Definition: ruby.h:1332
#define RBIGNUM_NEGATIVE_P(b)
Definition: ruby.h:1076
#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:606
static unsigned short check_rounding_mode(VALUE const v)
Definition: bigdecimal.c:411
static VALUE BigDecimal_save_limit(VALUE self)
Definition: bigdecimal.c:2621
#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:2571
static int AddExponent(Real *a, SIGNED_VALUE n)
Definition: bigdecimal.c:3754
#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:3518
VALUE rb_mBigMath
Definition: bigdecimal.c:42
#define NULL
Definition: _sdbm.c:102
#define FIX2LONG(x)
Definition: ruby.h:353
#define Qundef
Definition: ruby.h:436
#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
VALUE rb_str_new2(const char *)
void rb_warn(const char *fmt,...)
Definition: error.c:221
VP_EXPORT void * VpMemRealloc(void *ptr, size_t mb)
Definition: bigdecimal.c:3303
#define bp()
Definition: vm_debug.h:27
VALUE rb_eArgError
Definition: error.c:517
static VALUE BigDecimal_IsInfinite(VALUE self)
Definition: bigdecimal.c:620
VP_EXPORT ssize_t VpExponent10(Real *a)
Definition: bigdecimal.c:5006
#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_num_coerce_cmp(VALUE, VALUE, ID)
Definition: numeric.c:291
static VALUE BigDecimal_to_s(int argc, VALUE *argv, VALUE self)
Definition: bigdecimal.c:1854
VALUE rb_str_new(const char *, long)
Definition: string.c:425
VALUE rb_obj_class(VALUE)
Definition: object.c:194
#define SIGNED_VALUE
Definition: ruby.h:106
#define GUARD_OBJ(p, y)
Definition: bigdecimal.c:66
VP_EXPORT int VpSqrt(Real *y, Real *x)
Definition: bigdecimal.c:5556