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