Ruby  2.0.0p247(2013-06-27revision41674)
ext/bigdecimal/bigdecimal.c
Go to the documentation of this file.
00001 /*
00002  *
00003  * Ruby BigDecimal(Variable decimal precision) extension library.
00004  *
00005  * Copyright(C) 2002 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp)
00006  *
00007  * You may distribute under the terms of either the GNU General Public
00008  * License or the Artistic License, as specified in the README file
00009  * of this BigDecimal distribution.
00010  *
00011  *  NOTE: Change log in this source removed to reduce source code size.
00012  *        See rev. 1.25 if needed.
00013  *
00014  */
00015 
00016 /* #define BIGDECIMAL_DEBUG 1 */
00017 #ifdef BIGDECIMAL_DEBUG
00018 # define BIGDECIMAL_ENABLE_VPRINT 1
00019 #endif
00020 #include "bigdecimal.h"
00021 
00022 #ifndef BIGDECIMAL_DEBUG
00023 # define NDEBUG
00024 #endif
00025 #include <assert.h>
00026 
00027 #include <ctype.h>
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <string.h>
00031 #include <errno.h>
00032 #include <math.h>
00033 #include "math.h"
00034 
00035 #ifdef HAVE_IEEEFP_H
00036 #include <ieeefp.h>
00037 #endif
00038 
00039 /* #define ENABLE_NUMERIC_STRING */
00040 
00041 VALUE rb_cBigDecimal;
00042 VALUE rb_mBigMath;
00043 
00044 static ID id_BigDecimal_exception_mode;
00045 static ID id_BigDecimal_rounding_mode;
00046 static ID id_BigDecimal_precision_limit;
00047 
00048 static ID id_up;
00049 static ID id_down;
00050 static ID id_truncate;
00051 static ID id_half_up;
00052 static ID id_default;
00053 static ID id_half_down;
00054 static ID id_half_even;
00055 static ID id_banker;
00056 static ID id_ceiling;
00057 static ID id_ceil;
00058 static ID id_floor;
00059 static ID id_to_r;
00060 static ID id_eq;
00061 
00062 /* MACRO's to guard objects from GC by keeping them in stack */
00063 #define ENTER(n) volatile VALUE RB_UNUSED_VAR(vStack[n]);int iStack=0
00064 #define PUSH(x)  vStack[iStack++] = (VALUE)(x);
00065 #define SAVE(p)  PUSH(p->obj);
00066 #define GUARD_OBJ(p,y) {p=y;SAVE(p);}
00067 
00068 #define BASE_FIG  RMPD_COMPONENT_FIGURES
00069 #define BASE      RMPD_BASE
00070 
00071 #define HALF_BASE (BASE/2)
00072 #define BASE1 (BASE/10)
00073 
00074 #ifndef DBLE_FIG
00075 #define DBLE_FIG (DBL_DIG+1)    /* figure of double */
00076 #endif
00077 
00078 #ifndef RBIGNUM_ZERO_P
00079 # define RBIGNUM_ZERO_P(x) (RBIGNUM_LEN(x) == 0 || \
00080                             (RBIGNUM_DIGITS(x)[0] == 0 && \
00081                              (RBIGNUM_LEN(x) == 1 || bigzero_p(x))))
00082 #endif
00083 
00084 static inline int
00085 bigzero_p(VALUE x)
00086 {
00087     long i;
00088     BDIGIT *ds = RBIGNUM_DIGITS(x);
00089 
00090     for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) {
00091         if (ds[i]) return 0;
00092     }
00093     return 1;
00094 }
00095 
00096 #ifndef RRATIONAL_ZERO_P
00097 # define RRATIONAL_ZERO_P(x) (FIXNUM_P(RRATIONAL(x)->num) && \
00098                               FIX2LONG(RRATIONAL(x)->num) == 0)
00099 #endif
00100 
00101 #ifndef RRATIONAL_NEGATIVE_P
00102 # define RRATIONAL_NEGATIVE_P(x) RTEST(rb_funcall((x), '<', 1, INT2FIX(0)))
00103 #endif
00104 
00105 /*
00106  * ================== Ruby Interface part ==========================
00107  */
00108 #define DoSomeOne(x,y,f) rb_num_coerce_bin(x,y,f)
00109 
00110 /*
00111  * Returns the BigDecimal version number.
00112  */
00113 static VALUE
00114 BigDecimal_version(VALUE self)
00115 {
00116     /*
00117      * 1.0.0: Ruby 1.8.0
00118      * 1.0.1: Ruby 1.8.1
00119      * 1.1.0: Ruby 1.9.3
00120     */
00121     return rb_str_new2("1.1.0");
00122 }
00123 
00124 /*
00125  *   VP routines used in BigDecimal part
00126  */
00127 static unsigned short VpGetException(void);
00128 static void  VpSetException(unsigned short f);
00129 static void  VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v);
00130 static int   VpLimitRound(Real *c, size_t ixDigit);
00131 static Real *VpCopy(Real *pv, Real const* const x);
00132 
00133 /*
00134  *  **** BigDecimal part ****
00135  */
00136 
00137 static void
00138 BigDecimal_delete(void *pv)
00139 {
00140     VpFree(pv);
00141 }
00142 
00143 static size_t
00144 BigDecimal_memsize(const void *ptr)
00145 {
00146     const Real *pv = ptr;
00147     return pv ? (sizeof(*pv) + pv->MaxPrec * sizeof(BDIGIT)) : 0;
00148 }
00149 
00150 static const rb_data_type_t BigDecimal_data_type = {
00151     "BigDecimal",
00152     {0, BigDecimal_delete, BigDecimal_memsize,},
00153 };
00154 
00155 static inline int
00156 is_kind_of_BigDecimal(VALUE const v)
00157 {
00158     return rb_typeddata_is_kind_of(v, &BigDecimal_data_type);
00159 }
00160 
00161 static VALUE
00162 ToValue(Real *p)
00163 {
00164     if(VpIsNaN(p)) {
00165         VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'(Not a Number)",0);
00166     } else if(VpIsPosInf(p)) {
00167         VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0);
00168     } else if(VpIsNegInf(p)) {
00169         VpException(VP_EXCEPTION_INFINITY,"Computation results to '-Infinity'",0);
00170     }
00171     return p->obj;
00172 }
00173 
00174 NORETURN(static void cannot_be_coerced_into_BigDecimal(VALUE, VALUE));
00175 
00176 static void
00177 cannot_be_coerced_into_BigDecimal(VALUE exc_class, VALUE v)
00178 {
00179     VALUE str;
00180 
00181     if (rb_special_const_p(v)) {
00182         str = rb_inspect(v);
00183     }
00184     else {
00185         str = rb_class_name(rb_obj_class(v));
00186     }
00187 
00188     str = rb_str_cat2(rb_str_dup(str), " can't be coerced into BigDecimal");
00189     rb_exc_raise(rb_exc_new3(exc_class, str));
00190 }
00191 
00192 static VALUE BigDecimal_div2(int, VALUE*, VALUE);
00193 
00194 static Real*
00195 GetVpValueWithPrec(VALUE v, long prec, int must)
00196 {
00197     Real *pv;
00198     VALUE num, bg, args[2];
00199     char szD[128];
00200     VALUE orig = Qundef;
00201 
00202 again:
00203     switch(TYPE(v))
00204     {
00205       case T_FLOAT:
00206         if (prec < 0) goto unable_to_coerce_without_prec;
00207         if (prec > DBL_DIG+1)goto SomeOneMayDoIt;
00208         v = rb_funcall(v, id_to_r, 0);
00209         goto again;
00210       case T_RATIONAL:
00211         if (prec < 0) goto unable_to_coerce_without_prec;
00212 
00213         if (orig == Qundef ? (orig = v, 1) : orig != v) {
00214             num = RRATIONAL(v)->num;
00215             pv = GetVpValueWithPrec(num, -1, must);
00216             if (pv == NULL) goto SomeOneMayDoIt;
00217 
00218             args[0] = RRATIONAL(v)->den;
00219             args[1] = LONG2NUM(prec);
00220             v = BigDecimal_div2(2, args, ToValue(pv));
00221             goto again;
00222         }
00223 
00224         v = orig;
00225         goto SomeOneMayDoIt;
00226 
00227       case T_DATA:
00228         if (is_kind_of_BigDecimal(v)) {
00229             pv = DATA_PTR(v);
00230             return pv;
00231         }
00232         else {
00233             goto SomeOneMayDoIt;
00234         }
00235         break;
00236 
00237       case T_FIXNUM:
00238         sprintf(szD, "%ld", FIX2LONG(v));
00239         return VpCreateRbObject(VpBaseFig() * 2 + 1, szD);
00240 
00241 #ifdef ENABLE_NUMERIC_STRING
00242       case T_STRING:
00243         SafeStringValue(v);
00244         return VpCreateRbObject(strlen(RSTRING_PTR(v)) + VpBaseFig() + 1,
00245                                 RSTRING_PTR(v));
00246 #endif /* ENABLE_NUMERIC_STRING */
00247 
00248       case T_BIGNUM:
00249         bg = rb_big2str(v, 10);
00250         return VpCreateRbObject(strlen(RSTRING_PTR(bg)) + VpBaseFig() + 1,
00251                                 RSTRING_PTR(bg));
00252       default:
00253         goto SomeOneMayDoIt;
00254     }
00255 
00256 SomeOneMayDoIt:
00257     if (must) {
00258         cannot_be_coerced_into_BigDecimal(rb_eTypeError, v);
00259     }
00260     return NULL; /* NULL means to coerce */
00261 
00262 unable_to_coerce_without_prec:
00263     if (must) {
00264         rb_raise(rb_eArgError,
00265                  "%s can't be coerced into BigDecimal without a precision",
00266                  rb_obj_classname(v));
00267     }
00268     return NULL;
00269 }
00270 
00271 static Real*
00272 GetVpValue(VALUE v, int must)
00273 {
00274     return GetVpValueWithPrec(v, -1, must);
00275 }
00276 
00277 /* call-seq:
00278  * BigDecimal.double_fig
00279  *
00280  * The BigDecimal.double_fig class method returns the number of digits a
00281  * Float number is allowed to have. The result depends upon the CPU and OS
00282  * in use.
00283  */
00284 static VALUE
00285 BigDecimal_double_fig(VALUE self)
00286 {
00287     return INT2FIX(VpDblFig());
00288 }
00289 
00290 /* call-seq:
00291  * precs
00292  *
00293  * Returns an Array of two Integer values.
00294  *
00295  * The first value is the current number of significant digits in the
00296  * BigDecimal. The second value is the maximum number of significant digits
00297  * for the BigDecimal.
00298  */
00299 static VALUE
00300 BigDecimal_prec(VALUE self)
00301 {
00302     ENTER(1);
00303     Real *p;
00304     VALUE obj;
00305 
00306     GUARD_OBJ(p,GetVpValue(self,1));
00307     obj = rb_assoc_new(INT2NUM(p->Prec*VpBaseFig()),
00308                        INT2NUM(p->MaxPrec*VpBaseFig()));
00309     return obj;
00310 }
00311 
00312 /*
00313  * call-seq: hash
00314  *
00315  * Creates a hash for this BigDecimal.
00316  *
00317  * Two BigDecimals with equal sign,
00318  * fractional part and exponent have the same hash.
00319  */
00320 static VALUE
00321 BigDecimal_hash(VALUE self)
00322 {
00323     ENTER(1);
00324     Real *p;
00325     st_index_t hash;
00326 
00327     GUARD_OBJ(p,GetVpValue(self,1));
00328     hash = (st_index_t)p->sign;
00329     /* hash!=2: the case for 0(1),NaN(0) or +-Infinity(3) is sign itself */
00330     if(hash == 2 || hash == (st_index_t)-2) {
00331         hash ^= rb_memhash(p->frac, sizeof(BDIGIT)*p->Prec);
00332         hash += p->exponent;
00333     }
00334     return INT2FIX(hash);
00335 }
00336 
00337 /*
00338  * call-seq: _dump
00339  *
00340  * Method used to provide marshalling support.
00341  *
00342  *      inf = BigDecimal.new('Infinity')
00343  *      => #<BigDecimal:1e16fa8,'Infinity',9(9)>
00344  *      BigDecimal._load(inf._dump)
00345  *      => #<BigDecimal:1df8dc8,'Infinity',9(9)>
00346  *
00347  * See the Marshal module.
00348  */
00349 static VALUE
00350 BigDecimal_dump(int argc, VALUE *argv, VALUE self)
00351 {
00352     ENTER(5);
00353     Real *vp;
00354     char *psz;
00355     VALUE dummy;
00356     volatile VALUE dump;
00357 
00358     rb_scan_args(argc, argv, "01", &dummy);
00359     GUARD_OBJ(vp,GetVpValue(self,1));
00360     dump = rb_str_new(0,VpNumOfChars(vp,"E")+50);
00361     psz = RSTRING_PTR(dump);
00362     sprintf(psz, "%"PRIuSIZE":", VpMaxPrec(vp)*VpBaseFig());
00363     VpToString(vp, psz+strlen(psz), 0, 0);
00364     rb_str_resize(dump, strlen(psz));
00365     return dump;
00366 }
00367 
00368 /*
00369  * Internal method used to provide marshalling support. See the Marshal module.
00370  */
00371 static VALUE
00372 BigDecimal_load(VALUE self, VALUE str)
00373 {
00374     ENTER(2);
00375     Real *pv;
00376     unsigned char *pch;
00377     unsigned char ch;
00378     unsigned long m=0;
00379 
00380     SafeStringValue(str);
00381     pch = (unsigned char *)RSTRING_PTR(str);
00382     /* First get max prec */
00383     while((*pch)!=(unsigned char)'\0' && (ch=*pch++)!=(unsigned char)':') {
00384         if(!ISDIGIT(ch)) {
00385             rb_raise(rb_eTypeError, "load failed: invalid character in the marshaled string");
00386         }
00387         m = m*10 + (unsigned long)(ch-'0');
00388     }
00389     if(m>VpBaseFig()) m -= VpBaseFig();
00390     GUARD_OBJ(pv,VpNewRbClass(m,(char *)pch,self));
00391     m /= VpBaseFig();
00392     if(m && pv->MaxPrec>m) pv->MaxPrec = m+1;
00393     return ToValue(pv);
00394 }
00395 
00396 static unsigned short
00397 check_rounding_mode(VALUE const v)
00398 {
00399     unsigned short sw;
00400     ID id;
00401     switch (TYPE(v)) {
00402       case T_SYMBOL:
00403         id = SYM2ID(v);
00404         if (id == id_up)
00405             return VP_ROUND_UP;
00406         if (id == id_down || id == id_truncate)
00407             return VP_ROUND_DOWN;
00408         if (id == id_half_up || id == id_default)
00409             return VP_ROUND_HALF_UP;
00410         if (id == id_half_down)
00411             return VP_ROUND_HALF_DOWN;
00412         if (id == id_half_even || id == id_banker)
00413             return VP_ROUND_HALF_EVEN;
00414         if (id == id_ceiling || id == id_ceil)
00415             return VP_ROUND_CEIL;
00416         if (id == id_floor)
00417             return VP_ROUND_FLOOR;
00418         rb_raise(rb_eArgError, "invalid rounding mode");
00419 
00420       default:
00421         break;
00422     }
00423 
00424     Check_Type(v, T_FIXNUM);
00425     sw = (unsigned short)FIX2UINT(v);
00426     if (!VpIsRoundMode(sw)) {
00427         rb_raise(rb_eArgError, "invalid rounding mode");
00428     }
00429     return sw;
00430 }
00431 
00432 /* call-seq:
00433  * BigDecimal.mode(mode, value)
00434  *
00435  * Controls handling of arithmetic exceptions and rounding. If no value
00436  * is supplied, the current value is returned.
00437  *
00438  * Six values of the mode parameter control the handling of arithmetic
00439  * exceptions:
00440  *
00441  * BigDecimal::EXCEPTION_NaN
00442  * BigDecimal::EXCEPTION_INFINITY
00443  * BigDecimal::EXCEPTION_UNDERFLOW
00444  * BigDecimal::EXCEPTION_OVERFLOW
00445  * BigDecimal::EXCEPTION_ZERODIVIDE
00446  * BigDecimal::EXCEPTION_ALL
00447  *
00448  * For each mode parameter above, if the value set is false, computation
00449  * continues after an arithmetic exception of the appropriate type.
00450  * When computation continues, results are as follows:
00451  *
00452  * EXCEPTION_NaN:: NaN
00453  * EXCEPTION_INFINITY:: +infinity or -infinity
00454  * EXCEPTION_UNDERFLOW:: 0
00455  * EXCEPTION_OVERFLOW:: +infinity or -infinity
00456  * EXCEPTION_ZERODIVIDE:: +infinity or -infinity
00457  *
00458  * One value of the mode parameter controls the rounding of numeric values:
00459  * BigDecimal::ROUND_MODE. The values it can take are:
00460  *
00461  * ROUND_UP, :up:: round away from zero
00462  * ROUND_DOWN, :down, :truncate:: round towards zero (truncate)
00463  * ROUND_HALF_UP, :half_up, :default:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round away from zero. (default)
00464  * ROUND_HALF_DOWN, :half_down:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards zero.
00465  * 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)
00466  * ROUND_CEILING, :ceiling, :ceil:: round towards positive infinity (ceil)
00467  * ROUND_FLOOR, :floor:: round towards negative infinity (floor)
00468  *
00469  */
00470 static VALUE
00471 BigDecimal_mode(int argc, VALUE *argv, VALUE self)
00472 {
00473     VALUE which;
00474     VALUE val;
00475     unsigned long f,fo;
00476 
00477     if(rb_scan_args(argc,argv,"11",&which,&val)==1) val = Qnil;
00478 
00479     Check_Type(which, T_FIXNUM);
00480     f = (unsigned long)FIX2INT(which);
00481 
00482     if(f&VP_EXCEPTION_ALL) {
00483         /* Exception mode setting */
00484         fo = VpGetException();
00485         if(val==Qnil) return INT2FIX(fo);
00486         if(val!=Qfalse && val!=Qtrue) {
00487             rb_raise(rb_eArgError, "second argument must be true or false");
00488             return Qnil; /* Not reached */
00489         }
00490         if(f&VP_EXCEPTION_INFINITY) {
00491             VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_INFINITY):
00492                            (fo&(~VP_EXCEPTION_INFINITY))));
00493         }
00494         fo = VpGetException();
00495         if(f&VP_EXCEPTION_NaN) {
00496             VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_NaN):
00497                            (fo&(~VP_EXCEPTION_NaN))));
00498         }
00499         fo = VpGetException();
00500         if(f&VP_EXCEPTION_UNDERFLOW) {
00501             VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_UNDERFLOW):
00502                            (fo&(~VP_EXCEPTION_UNDERFLOW))));
00503         }
00504         fo = VpGetException();
00505         if(f&VP_EXCEPTION_ZERODIVIDE) {
00506             VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_ZERODIVIDE):
00507                            (fo&(~VP_EXCEPTION_ZERODIVIDE))));
00508         }
00509         fo = VpGetException();
00510         return INT2FIX(fo);
00511     }
00512     if (VP_ROUND_MODE == f) {
00513         /* Rounding mode setting */
00514         unsigned short sw;
00515         fo = VpGetRoundMode();
00516         if (NIL_P(val)) return INT2FIX(fo);
00517         sw = check_rounding_mode(val);
00518         fo = VpSetRoundMode(sw);
00519         return INT2FIX(fo);
00520     }
00521     rb_raise(rb_eTypeError, "first argument for BigDecimal#mode invalid");
00522     return Qnil;
00523 }
00524 
00525 static size_t
00526 GetAddSubPrec(Real *a, Real *b)
00527 {
00528     size_t mxs;
00529     size_t mx = a->Prec;
00530     SIGNED_VALUE d;
00531 
00532     if(!VpIsDef(a) || !VpIsDef(b)) return (size_t)-1L;
00533     if(mx < b->Prec) mx = b->Prec;
00534     if(a->exponent!=b->exponent) {
00535         mxs = mx;
00536         d = a->exponent - b->exponent;
00537         if (d < 0) d = -d;
00538         mx = mx + (size_t)d;
00539         if (mx<mxs) {
00540             return VpException(VP_EXCEPTION_INFINITY,"Exponent overflow",0);
00541         }
00542     }
00543     return mx;
00544 }
00545 
00546 static SIGNED_VALUE
00547 GetPositiveInt(VALUE v)
00548 {
00549     SIGNED_VALUE n;
00550     Check_Type(v, T_FIXNUM);
00551     n = FIX2INT(v);
00552     if (n < 0) {
00553         rb_raise(rb_eArgError, "argument must be positive");
00554     }
00555     return n;
00556 }
00557 
00558 VP_EXPORT Real *
00559 VpNewRbClass(size_t mx, const char *str, VALUE klass)
00560 {
00561     Real *pv = VpAlloc(mx,str);
00562     pv->obj = TypedData_Wrap_Struct(klass, &BigDecimal_data_type, pv);
00563     return pv;
00564 }
00565 
00566 VP_EXPORT Real *
00567 VpCreateRbObject(size_t mx, const char *str)
00568 {
00569     Real *pv = VpAlloc(mx,str);
00570     pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv);
00571     return pv;
00572 }
00573 
00574 #define VpAllocReal(prec) (Real *)VpMemAlloc(offsetof(Real, frac) + (prec) * sizeof(BDIGIT))
00575 #define VpReallocReal(ptr, prec) (Real *)VpMemRealloc((ptr), offsetof(Real, frac) + (prec) * sizeof(BDIGIT))
00576 
00577 static Real *
00578 VpCopy(Real *pv, Real const* const x)
00579 {
00580     assert(x != NULL);
00581 
00582     pv = VpReallocReal(pv, x->MaxPrec);
00583     pv->MaxPrec = x->MaxPrec;
00584     pv->Prec = x->Prec;
00585     pv->exponent = x->exponent;
00586     pv->sign = x->sign;
00587     pv->flag = x->flag;
00588     MEMCPY(pv->frac, x->frac, BDIGIT, pv->MaxPrec);
00589 
00590     return pv;
00591 }
00592 
00593 /* Returns True if the value is Not a Number */
00594 static VALUE
00595 BigDecimal_IsNaN(VALUE self)
00596 {
00597     Real *p = GetVpValue(self,1);
00598     if(VpIsNaN(p))  return Qtrue;
00599     return Qfalse;
00600 }
00601 
00602 /* Returns nil, -1, or +1 depending on whether the value is finite,
00603  * -infinity, or +infinity.
00604  */
00605 static VALUE
00606 BigDecimal_IsInfinite(VALUE self)
00607 {
00608     Real *p = GetVpValue(self,1);
00609     if(VpIsPosInf(p)) return INT2FIX(1);
00610     if(VpIsNegInf(p)) return INT2FIX(-1);
00611     return Qnil;
00612 }
00613 
00614 /* Returns True if the value is finite (not NaN or infinite) */
00615 static VALUE
00616 BigDecimal_IsFinite(VALUE self)
00617 {
00618     Real *p = GetVpValue(self,1);
00619     if(VpIsNaN(p)) return Qfalse;
00620     if(VpIsInf(p)) return Qfalse;
00621     return Qtrue;
00622 }
00623 
00624 static void
00625 BigDecimal_check_num(Real *p)
00626 {
00627     if(VpIsNaN(p)) {
00628        VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'(Not a Number)",1);
00629     } else if(VpIsPosInf(p)) {
00630        VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",1);
00631     } else if(VpIsNegInf(p)) {
00632        VpException(VP_EXCEPTION_INFINITY,"Computation results to '-Infinity'",1);
00633     }
00634 }
00635 
00636 static VALUE BigDecimal_split(VALUE self);
00637 
00638 /* Returns the value as an integer (Fixnum or Bignum).
00639  *
00640  * If the BigNumber is infinity or NaN, raises FloatDomainError.
00641  */
00642 static VALUE
00643 BigDecimal_to_i(VALUE self)
00644 {
00645     ENTER(5);
00646     ssize_t e, nf;
00647     Real *p;
00648 
00649     GUARD_OBJ(p,GetVpValue(self,1));
00650     BigDecimal_check_num(p);
00651 
00652     e = VpExponent10(p);
00653     if(e<=0) return INT2FIX(0);
00654     nf = VpBaseFig();
00655     if(e<=nf) {
00656         return LONG2NUM((long)(VpGetSign(p)*(BDIGIT_DBL_SIGNED)p->frac[0]));
00657     }
00658     else {
00659         VALUE a = BigDecimal_split(self);
00660         VALUE digits = RARRAY_PTR(a)[1];
00661         VALUE numerator = rb_funcall(digits, rb_intern("to_i"), 0);
00662         VALUE ret;
00663         ssize_t dpower = e - (ssize_t)RSTRING_LEN(digits);
00664 
00665         if (VpGetSign(p) < 0) {
00666             numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1));
00667         }
00668         if (dpower < 0) {
00669             ret = rb_funcall(numerator, rb_intern("div"), 1,
00670                               rb_funcall(INT2FIX(10), rb_intern("**"), 1,
00671                                          INT2FIX(-dpower)));
00672         }
00673         else
00674             ret = rb_funcall(numerator, '*', 1,
00675                              rb_funcall(INT2FIX(10), rb_intern("**"), 1,
00676                                         INT2FIX(dpower)));
00677         if (RB_TYPE_P(ret, T_FLOAT))
00678             rb_raise(rb_eFloatDomainError, "Infinity");
00679         return ret;
00680     }
00681 }
00682 
00683 /* Returns a new Float object having approximately the same value as the
00684  * BigDecimal number. Normal accuracy limits and built-in errors of binary
00685  * Float arithmetic apply.
00686  */
00687 static VALUE
00688 BigDecimal_to_f(VALUE self)
00689 {
00690     ENTER(1);
00691     Real *p;
00692     double d;
00693     SIGNED_VALUE e;
00694     char *buf;
00695     volatile VALUE str;
00696 
00697     GUARD_OBJ(p, GetVpValue(self, 1));
00698     if (VpVtoD(&d, &e, p) != 1)
00699         return rb_float_new(d);
00700     if (e > (SIGNED_VALUE)(DBL_MAX_10_EXP+BASE_FIG))
00701         goto overflow;
00702     if (e < (SIGNED_VALUE)(DBL_MIN_10_EXP-BASE_FIG))
00703         goto underflow;
00704 
00705     str = rb_str_new(0, VpNumOfChars(p,"E"));
00706     buf = RSTRING_PTR(str);
00707     VpToString(p, buf, 0, 0);
00708     errno = 0;
00709     d = strtod(buf, 0);
00710     if (errno == ERANGE) {
00711         if (d == 0.0) goto underflow;
00712         if (fabs(d) >= HUGE_VAL) goto overflow;
00713     }
00714     return rb_float_new(d);
00715 
00716 overflow:
00717     VpException(VP_EXCEPTION_OVERFLOW, "BigDecimal to Float conversion", 0);
00718     if (p->sign >= 0)
00719         return rb_float_new(VpGetDoublePosInf());
00720     else
00721         return rb_float_new(VpGetDoubleNegInf());
00722 
00723 underflow:
00724     VpException(VP_EXCEPTION_UNDERFLOW, "BigDecimal to Float conversion", 0);
00725     if (p->sign >= 0)
00726         return rb_float_new(0.0);
00727     else
00728         return rb_float_new(-0.0);
00729 }
00730 
00731 
00732 /* Converts a BigDecimal to a Rational.
00733  */
00734 static VALUE
00735 BigDecimal_to_r(VALUE self)
00736 {
00737     Real *p;
00738     ssize_t sign, power, denomi_power;
00739     VALUE a, digits, numerator;
00740 
00741     p = GetVpValue(self,1);
00742     BigDecimal_check_num(p);
00743 
00744     sign = VpGetSign(p);
00745     power = VpExponent10(p);
00746     a = BigDecimal_split(self);
00747     digits = RARRAY_PTR(a)[1];
00748     denomi_power = power - RSTRING_LEN(digits);
00749     numerator = rb_funcall(digits, rb_intern("to_i"), 0);
00750 
00751     if (sign < 0) {
00752         numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1));
00753     }
00754     if (denomi_power < 0) {
00755         return rb_Rational(numerator,
00756                            rb_funcall(INT2FIX(10), rb_intern("**"), 1,
00757                                       INT2FIX(-denomi_power)));
00758     }
00759     else {
00760         return rb_Rational1(rb_funcall(numerator, '*', 1,
00761                                        rb_funcall(INT2FIX(10), rb_intern("**"), 1,
00762                                                   INT2FIX(denomi_power))));
00763     }
00764 }
00765 
00766 /* The coerce method provides support for Ruby type coercion. It is not
00767  * enabled by default.
00768  *
00769  * This means that binary operations like + * / or - can often be performed
00770  * on a BigDecimal and an object of another type, if the other object can
00771  * be coerced into a BigDecimal value.
00772  *
00773  * e.g.
00774  * a = BigDecimal.new("1.0")
00775  * b = a / 2.0  -> 0.5
00776  *
00777  * Note that coercing a String to a BigDecimal is not supported by default;
00778  * it requires a special compile-time option when building Ruby.
00779  */
00780 static VALUE
00781 BigDecimal_coerce(VALUE self, VALUE other)
00782 {
00783     ENTER(2);
00784     VALUE obj;
00785     Real *b;
00786 
00787     if (RB_TYPE_P(other, T_FLOAT)) {
00788         obj = rb_assoc_new(other, BigDecimal_to_f(self));
00789     }
00790     else {
00791         if (RB_TYPE_P(other, T_RATIONAL)) {
00792             Real* pv = DATA_PTR(self);
00793             GUARD_OBJ(b, GetVpValueWithPrec(other, pv->Prec*VpBaseFig(), 1));
00794         }
00795         else {
00796             GUARD_OBJ(b, GetVpValue(other, 1));
00797         }
00798         obj = rb_assoc_new(b->obj, self);
00799     }
00800 
00801     return obj;
00802 }
00803 
00804 /*
00805  * call-seq: +@
00806  *
00807  * Return self.
00808  *
00809  * e.g.
00810  *   b = +a  # b == a
00811  */
00812 static VALUE
00813 BigDecimal_uplus(VALUE self)
00814 {
00815     return self;
00816 }
00817 
00818  /*
00819   * Document-method: BigDecimal#add
00820   * Document-method: BigDecimal#+
00821   *
00822   * call-seq:
00823   * add(value, digits)
00824   *
00825   * Add the specified value.
00826   *
00827   * e.g.
00828   *   c = a.add(b,n)
00829   *   c = a + b
00830   *
00831   * digits:: If specified and less than the number of significant digits of the
00832   * result, the result is rounded to that number of digits, according to
00833   * BigDecimal.mode.
00834   */
00835 static VALUE
00836 BigDecimal_add(VALUE self, VALUE r)
00837 {
00838     ENTER(5);
00839     Real *c, *a, *b;
00840     size_t mx;
00841 
00842     GUARD_OBJ(a, GetVpValue(self, 1));
00843     if (RB_TYPE_P(r, T_FLOAT)) {
00844         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
00845     }
00846     else if (RB_TYPE_P(r, T_RATIONAL)) {
00847         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
00848     }
00849     else {
00850         b = GetVpValue(r,0);
00851     }
00852 
00853     if (!b) return DoSomeOne(self,r,'+');
00854     SAVE(b);
00855 
00856     if (VpIsNaN(b)) return b->obj;
00857     if (VpIsNaN(a)) return a->obj;
00858 
00859     mx = GetAddSubPrec(a, b);
00860     if (mx == (size_t)-1L) {
00861         GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
00862         VpAddSub(c, a, b, 1);
00863     }
00864     else {
00865         GUARD_OBJ(c, VpCreateRbObject(mx * (VpBaseFig() + 1), "0"));
00866         if(!mx) {
00867             VpSetInf(c, VpGetSign(a));
00868         }
00869         else {
00870             VpAddSub(c, a, b, 1);
00871         }
00872     }
00873     return ToValue(c);
00874 }
00875 
00876  /* call-seq:
00877   * sub(value, digits)
00878   *
00879   * Subtract the specified value.
00880   *
00881   * e.g.
00882   *   c = a.sub(b,n)
00883   *   c = a - b
00884   *
00885   * digits:: If specified and less than the number of significant digits of the
00886   * result, the result is rounded to that number of digits, according to
00887   * BigDecimal.mode.
00888   */
00889 static VALUE
00890 BigDecimal_sub(VALUE self, VALUE r)
00891 {
00892     ENTER(5);
00893     Real *c, *a, *b;
00894     size_t mx;
00895 
00896     GUARD_OBJ(a,GetVpValue(self,1));
00897     if (RB_TYPE_P(r, T_FLOAT)) {
00898         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
00899     }
00900     else if (RB_TYPE_P(r, T_RATIONAL)) {
00901         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
00902     }
00903     else {
00904         b = GetVpValue(r,0);
00905     }
00906 
00907     if(!b) return DoSomeOne(self,r,'-');
00908     SAVE(b);
00909 
00910     if(VpIsNaN(b)) return b->obj;
00911     if(VpIsNaN(a)) return a->obj;
00912 
00913     mx = GetAddSubPrec(a,b);
00914     if (mx == (size_t)-1L) {
00915         GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
00916         VpAddSub(c, a, b, -1);
00917     } else {
00918         GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
00919         if(!mx) {
00920             VpSetInf(c,VpGetSign(a));
00921         } else {
00922             VpAddSub(c, a, b, -1);
00923         }
00924     }
00925     return ToValue(c);
00926 }
00927 
00928 static VALUE
00929 BigDecimalCmp(VALUE self, VALUE r,char op)
00930 {
00931     ENTER(5);
00932     SIGNED_VALUE e;
00933     Real *a, *b=0;
00934     GUARD_OBJ(a,GetVpValue(self,1));
00935     switch (TYPE(r)) {
00936     case T_DATA:
00937         if (!is_kind_of_BigDecimal(r)) break;
00938         /* fall through */
00939     case T_FIXNUM:
00940         /* fall through */
00941     case T_BIGNUM:
00942         GUARD_OBJ(b, GetVpValue(r,0));
00943         break;
00944 
00945     case T_FLOAT:
00946         GUARD_OBJ(b, GetVpValueWithPrec(r, DBL_DIG+1, 0));
00947         break;
00948 
00949     case T_RATIONAL:
00950         GUARD_OBJ(b, GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 0));
00951         break;
00952 
00953     default:
00954         break;
00955     }
00956     if (b == NULL) {
00957         ID f = 0;
00958 
00959         switch (op) {
00960         case '*':
00961             return rb_num_coerce_cmp(self, r, rb_intern("<=>"));
00962 
00963         case '=':
00964             return RTEST(rb_num_coerce_cmp(self, r, rb_intern("=="))) ? Qtrue : Qfalse;
00965 
00966         case 'G':
00967             f = rb_intern(">=");
00968             break;
00969 
00970         case 'L':
00971             f = rb_intern("<=");
00972             break;
00973 
00974         case '>':
00975             /* fall through */
00976         case '<':
00977             f = (ID)op;
00978             break;
00979 
00980         default:
00981             break;
00982         }
00983         return rb_num_coerce_relop(self, r, f);
00984     }
00985     SAVE(b);
00986     e = VpComp(a, b);
00987     if (e == 999)
00988         return (op == '*') ? Qnil : Qfalse;
00989     switch (op) {
00990     case '*':
00991         return   INT2FIX(e); /* any op */
00992 
00993     case '=':
00994         if(e==0) return Qtrue;
00995         return Qfalse;
00996 
00997     case 'G':
00998         if(e>=0) return Qtrue;
00999         return Qfalse;
01000 
01001     case '>':
01002         if(e> 0) return Qtrue;
01003         return Qfalse;
01004 
01005     case 'L':
01006         if(e<=0) return Qtrue;
01007         return Qfalse;
01008 
01009     case '<':
01010         if(e< 0) return Qtrue;
01011         return Qfalse;
01012 
01013     default:
01014         break;
01015     }
01016 
01017     rb_bug("Undefined operation in BigDecimalCmp()");
01018 
01019     UNREACHABLE;
01020 }
01021 
01022 /* Returns True if the value is zero. */
01023 static VALUE
01024 BigDecimal_zero(VALUE self)
01025 {
01026     Real *a = GetVpValue(self,1);
01027     return VpIsZero(a) ? Qtrue : Qfalse;
01028 }
01029 
01030 /* Returns self if the value is non-zero, nil otherwise. */
01031 static VALUE
01032 BigDecimal_nonzero(VALUE self)
01033 {
01034     Real *a = GetVpValue(self,1);
01035     return VpIsZero(a) ? Qnil : self;
01036 }
01037 
01038 /* The comparison operator.
01039  * a <=> b is 0 if a == b, 1 if a > b, -1 if a < b.
01040  */
01041 static VALUE
01042 BigDecimal_comp(VALUE self, VALUE r)
01043 {
01044     return BigDecimalCmp(self, r, '*');
01045 }
01046 
01047 /*
01048  * Tests for value equality; returns true if the values are equal.
01049  *
01050  * The == and === operators and the eql? method have the same implementation
01051  * for BigDecimal.
01052  *
01053  * Values may be coerced to perform the comparison:
01054  *
01055  * BigDecimal.new('1.0') == 1.0  -> true
01056  */
01057 static VALUE
01058 BigDecimal_eq(VALUE self, VALUE r)
01059 {
01060     return BigDecimalCmp(self, r, '=');
01061 }
01062 
01063 /* call-seq:
01064  * a < b
01065  *
01066  * Returns true if a is less than b.
01067  *
01068  * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
01069  */
01070 static VALUE
01071 BigDecimal_lt(VALUE self, VALUE r)
01072 {
01073     return BigDecimalCmp(self, r, '<');
01074 }
01075 
01076 /* call-seq:
01077  * a <= b
01078  *
01079  * Returns true if a is less than or equal to b.
01080  *
01081  * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
01082  */
01083 static VALUE
01084 BigDecimal_le(VALUE self, VALUE r)
01085 {
01086     return BigDecimalCmp(self, r, 'L');
01087 }
01088 
01089 /* call-seq:
01090  * a > b
01091  *
01092  * Returns true if a is greater than b.
01093  *
01094  * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce).
01095  */
01096 static VALUE
01097 BigDecimal_gt(VALUE self, VALUE r)
01098 {
01099     return BigDecimalCmp(self, r, '>');
01100 }
01101 
01102 /* call-seq:
01103  * a >= b
01104  *
01105  * Returns true if a is greater than or equal to b.
01106  *
01107  * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce)
01108  */
01109 static VALUE
01110 BigDecimal_ge(VALUE self, VALUE r)
01111 {
01112     return BigDecimalCmp(self, r, 'G');
01113 }
01114 
01115 /*
01116  * call-seq: -@
01117  *
01118  * Return the negation of self.
01119  *
01120  * e.g.
01121  *   b = -a
01122  *   b == a * -1
01123  */
01124 static VALUE
01125 BigDecimal_neg(VALUE self)
01126 {
01127     ENTER(5);
01128     Real *c, *a;
01129     GUARD_OBJ(a,GetVpValue(self,1));
01130     GUARD_OBJ(c,VpCreateRbObject(a->Prec *(VpBaseFig() + 1), "0"));
01131     VpAsgn(c, a, -1);
01132     return ToValue(c);
01133 }
01134 
01135  /*
01136   * Document-method: BigDecimal#mult
01137   *
01138   * call-seq: mult(value, digits)
01139   *
01140   * Multiply by the specified value.
01141   *
01142   * e.g.
01143   *   c = a.mult(b,n)
01144   *   c = a * b
01145   *
01146   * digits:: If specified and less than the number of significant digits of the
01147   * result, the result is rounded to that number of digits, according to
01148   * BigDecimal.mode.
01149   */
01150 static VALUE
01151 BigDecimal_mult(VALUE self, VALUE r)
01152 {
01153     ENTER(5);
01154     Real *c, *a, *b;
01155     size_t mx;
01156 
01157     GUARD_OBJ(a,GetVpValue(self,1));
01158     if (RB_TYPE_P(r, T_FLOAT)) {
01159         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
01160     }
01161     else if (RB_TYPE_P(r, T_RATIONAL)) {
01162         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
01163     }
01164     else {
01165         b = GetVpValue(r,0);
01166     }
01167 
01168     if(!b) return DoSomeOne(self,r,'*');
01169     SAVE(b);
01170 
01171     mx = a->Prec + b->Prec;
01172     GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
01173     VpMult(c, a, b);
01174     return ToValue(c);
01175 }
01176 
01177 static VALUE
01178 BigDecimal_divide(Real **c, Real **res, Real **div, VALUE self, VALUE r)
01179 /* For c = self.div(r): with round operation */
01180 {
01181     ENTER(5);
01182     Real *a, *b;
01183     size_t mx;
01184 
01185     GUARD_OBJ(a,GetVpValue(self,1));
01186     if (RB_TYPE_P(r, T_FLOAT)) {
01187         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
01188     }
01189     else if (RB_TYPE_P(r, T_RATIONAL)) {
01190         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
01191     }
01192     else {
01193         b = GetVpValue(r,0);
01194     }
01195 
01196     if(!b) return DoSomeOne(self,r,'/');
01197     SAVE(b);
01198 
01199     *div = b;
01200     mx = a->Prec + vabs(a->exponent);
01201     if(mx<b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent);
01202     mx =(mx + 1) * VpBaseFig();
01203     GUARD_OBJ((*c),VpCreateRbObject(mx, "#0"));
01204     GUARD_OBJ((*res),VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
01205     VpDivd(*c, *res, a, b);
01206     return (VALUE)0;
01207 }
01208 
01209  /* call-seq:
01210   * div(value, digits)
01211   * quo(value)
01212   *
01213   * Divide by the specified value.
01214   *
01215   * e.g.
01216   *   c = a.div(b,n)
01217   *
01218   * digits:: If specified and less than the number of significant digits of the
01219   * result, the result is rounded to that number of digits, according to
01220   * BigDecimal.mode.
01221   *
01222   * If digits is 0, the result is the same as the / operator. If not, the
01223   * result is an integer BigDecimal, by analogy with Float#div.
01224   *
01225   * The alias quo is provided since <code>div(value, 0)</code> is the same as
01226   * computing the quotient; see BigDecimal#divmod.
01227   */
01228 static VALUE
01229 BigDecimal_div(VALUE self, VALUE r)
01230 /* For c = self/r: with round operation */
01231 {
01232     ENTER(5);
01233     Real *c=NULL, *res=NULL, *div = NULL;
01234     r = BigDecimal_divide(&c, &res, &div, self, r);
01235     if(r!=(VALUE)0) return r; /* coerced by other */
01236     SAVE(c);SAVE(res);SAVE(div);
01237     /* a/b = c + r/b */
01238     /* c xxxxx
01239        r 00000yyyyy  ==> (y/b)*BASE >= HALF_BASE
01240      */
01241     /* Round */
01242     if(VpHasVal(div)) { /* frac[0] must be zero for NaN,INF,Zero */
01243         VpInternalRound(c, 0, c->frac[c->Prec-1], (BDIGIT)(VpBaseVal()*(BDIGIT_DBL)res->frac[0]/div->frac[0]));
01244     }
01245     return ToValue(c);
01246 }
01247 
01248 /*
01249  * %: mod = a%b = a - (a.to_f/b).floor * b
01250  * div = (a.to_f/b).floor
01251  */
01252 static VALUE
01253 BigDecimal_DoDivmod(VALUE self, VALUE r, Real **div, Real **mod)
01254 {
01255     ENTER(8);
01256     Real *c=NULL, *d=NULL, *res=NULL;
01257     Real *a, *b;
01258     size_t mx;
01259 
01260     GUARD_OBJ(a,GetVpValue(self,1));
01261     if (RB_TYPE_P(r, T_FLOAT)) {
01262         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
01263     }
01264     else if (RB_TYPE_P(r, T_RATIONAL)) {
01265         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
01266     }
01267     else {
01268         b = GetVpValue(r,0);
01269     }
01270 
01271     if(!b) return Qfalse;
01272     SAVE(b);
01273 
01274     if(VpIsNaN(a) || VpIsNaN(b)) goto NaN;
01275     if(VpIsInf(a) && VpIsInf(b)) goto NaN;
01276     if(VpIsZero(b)) {
01277         rb_raise(rb_eZeroDivError, "divided by 0");
01278     }
01279     if(VpIsInf(a)) {
01280        GUARD_OBJ(d,VpCreateRbObject(1, "0"));
01281        VpSetInf(d, (SIGNED_VALUE)(VpGetSign(a) == VpGetSign(b) ? 1 : -1));
01282        GUARD_OBJ(c,VpCreateRbObject(1, "NaN"));
01283        *div = d;
01284        *mod = c;
01285        return Qtrue;
01286     }
01287     if(VpIsInf(b)) {
01288        GUARD_OBJ(d,VpCreateRbObject(1, "0"));
01289        *div = d;
01290        *mod = a;
01291        return Qtrue;
01292     }
01293     if(VpIsZero(a)) {
01294        GUARD_OBJ(c,VpCreateRbObject(1, "0"));
01295        GUARD_OBJ(d,VpCreateRbObject(1, "0"));
01296        *div = d;
01297        *mod = c;
01298        return Qtrue;
01299     }
01300 
01301     mx = a->Prec + vabs(a->exponent);
01302     if(mx<b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent);
01303     mx =(mx + 1) * VpBaseFig();
01304     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01305     GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
01306     VpDivd(c, res, a, b);
01307     mx = c->Prec *(VpBaseFig() + 1);
01308     GUARD_OBJ(d,VpCreateRbObject(mx, "0"));
01309     VpActiveRound(d,c,VP_ROUND_DOWN,0);
01310     VpMult(res,d,b);
01311     VpAddSub(c,a,res,-1);
01312     if(!VpIsZero(c) && (VpGetSign(a)*VpGetSign(b)<0)) {
01313         VpAddSub(res,d,VpOne(),-1);
01314         GUARD_OBJ(d,VpCreateRbObject(GetAddSubPrec(c, b)*(VpBaseFig() + 1), "0"));
01315         VpAddSub(d  ,c,b,       1);
01316         *div = res;
01317         *mod = d;
01318     } else {
01319         *div = d;
01320         *mod = c;
01321     }
01322     return Qtrue;
01323 
01324 NaN:
01325     GUARD_OBJ(c,VpCreateRbObject(1, "NaN"));
01326     GUARD_OBJ(d,VpCreateRbObject(1, "NaN"));
01327     *div = d;
01328     *mod = c;
01329     return Qtrue;
01330 }
01331 
01332 /* call-seq:
01333  * a % b
01334  * a.modulo(b)
01335  *
01336  * Returns the modulus from dividing by b.
01337  *
01338  * See BigDecimal#divmod.
01339  */
01340 static VALUE
01341 BigDecimal_mod(VALUE self, VALUE r) /* %: a%b = a - (a.to_f/b).floor * b */
01342 {
01343     ENTER(3);
01344     Real *div=NULL, *mod=NULL;
01345 
01346     if(BigDecimal_DoDivmod(self,r,&div,&mod)) {
01347         SAVE(div); SAVE(mod);
01348         return ToValue(mod);
01349     }
01350     return DoSomeOne(self,r,'%');
01351 }
01352 
01353 static VALUE
01354 BigDecimal_divremain(VALUE self, VALUE r, Real **dv, Real **rv)
01355 {
01356     ENTER(10);
01357     size_t mx;
01358     Real *a=NULL, *b=NULL, *c=NULL, *res=NULL, *d=NULL, *rr=NULL, *ff=NULL;
01359     Real *f=NULL;
01360 
01361     GUARD_OBJ(a,GetVpValue(self,1));
01362     if (RB_TYPE_P(r, T_FLOAT)) {
01363         b = GetVpValueWithPrec(r, DBL_DIG+1, 1);
01364     }
01365     else if (RB_TYPE_P(r, T_RATIONAL)) {
01366         b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1);
01367     }
01368     else {
01369         b = GetVpValue(r,0);
01370     }
01371 
01372     if(!b) return DoSomeOne(self,r,rb_intern("remainder"));
01373     SAVE(b);
01374 
01375     mx  =(a->MaxPrec + b->MaxPrec) *VpBaseFig();
01376     GUARD_OBJ(c  ,VpCreateRbObject(mx, "0"));
01377     GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
01378     GUARD_OBJ(rr ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
01379     GUARD_OBJ(ff ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
01380 
01381     VpDivd(c, res, a, b);
01382 
01383     mx = c->Prec *(VpBaseFig() + 1);
01384 
01385     GUARD_OBJ(d,VpCreateRbObject(mx, "0"));
01386     GUARD_OBJ(f,VpCreateRbObject(mx, "0"));
01387 
01388     VpActiveRound(d,c,VP_ROUND_DOWN,0); /* 0: round off */
01389 
01390     VpFrac(f, c);
01391     VpMult(rr,f,b);
01392     VpAddSub(ff,res,rr,1);
01393 
01394     *dv = d;
01395     *rv = ff;
01396     return (VALUE)0;
01397 }
01398 
01399 /* Returns the remainder from dividing by the value.
01400  *
01401  * x.remainder(y) means x-y*(x/y).truncate
01402  */
01403 static VALUE
01404 BigDecimal_remainder(VALUE self, VALUE r) /* remainder */
01405 {
01406     VALUE  f;
01407     Real  *d,*rv=0;
01408     f = BigDecimal_divremain(self,r,&d,&rv);
01409     if(f!=(VALUE)0) return f;
01410     return ToValue(rv);
01411 }
01412 
01413 /* Divides by the specified value, and returns the quotient and modulus
01414  * as BigDecimal numbers. The quotient is rounded towards negative infinity.
01415  *
01416  * For example:
01417  *
01418  * require 'bigdecimal'
01419  *
01420  * a = BigDecimal.new("42")
01421  * b = BigDecimal.new("9")
01422  *
01423  * q,m = a.divmod(b)
01424  *
01425  * c = q * b + m
01426  *
01427  * a == c  -> true
01428  *
01429  * The quotient q is (a/b).floor, and the modulus is the amount that must be
01430  * added to q * b to get a.
01431  */
01432 static VALUE
01433 BigDecimal_divmod(VALUE self, VALUE r)
01434 {
01435     ENTER(5);
01436     Real *div=NULL, *mod=NULL;
01437 
01438     if(BigDecimal_DoDivmod(self,r,&div,&mod)) {
01439         SAVE(div); SAVE(mod);
01440         return rb_assoc_new(ToValue(div), ToValue(mod));
01441     }
01442     return DoSomeOne(self,r,rb_intern("divmod"));
01443 }
01444 
01445 /*
01446  * See BigDecimal#quo
01447  */
01448 static VALUE
01449 BigDecimal_div2(int argc, VALUE *argv, VALUE self)
01450 {
01451     ENTER(5);
01452     VALUE b,n;
01453     int na = rb_scan_args(argc,argv,"11",&b,&n);
01454     if(na==1) { /* div in Float sense */
01455        Real *div=NULL;
01456        Real *mod;
01457        if(BigDecimal_DoDivmod(self,b,&div,&mod)) {
01458           return BigDecimal_to_i(ToValue(div));
01459        }
01460        return DoSomeOne(self,b,rb_intern("div"));
01461     } else {    /* div in BigDecimal sense */
01462        SIGNED_VALUE ix = GetPositiveInt(n);
01463        if (ix == 0) return BigDecimal_div(self, b);
01464        else {
01465           Real *res=NULL;
01466           Real *av=NULL, *bv=NULL, *cv=NULL;
01467           size_t mx = (ix+VpBaseFig()*2);
01468           size_t pl = VpSetPrecLimit(0);
01469 
01470           GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
01471           GUARD_OBJ(av,GetVpValue(self,1));
01472           GUARD_OBJ(bv,GetVpValue(b,1));
01473           mx = av->Prec + bv->Prec + 2;
01474           if(mx <= cv->MaxPrec) mx = cv->MaxPrec+1;
01475           GUARD_OBJ(res,VpCreateRbObject((mx * 2  + 2)*VpBaseFig(), "#0"));
01476           VpDivd(cv,res,av,bv);
01477           VpSetPrecLimit(pl);
01478           VpLeftRound(cv,VpGetRoundMode(),ix);
01479           return ToValue(cv);
01480        }
01481     }
01482 }
01483 
01484 static VALUE
01485 BigDecimal_add2(VALUE self, VALUE b, VALUE n)
01486 {
01487     ENTER(2);
01488     Real   *cv;
01489     SIGNED_VALUE mx = GetPositiveInt(n);
01490     if (mx == 0) return BigDecimal_add(self, b);
01491     else {
01492        size_t pl = VpSetPrecLimit(0);
01493        VALUE   c = BigDecimal_add(self,b);
01494        VpSetPrecLimit(pl);
01495        GUARD_OBJ(cv,GetVpValue(c,1));
01496        VpLeftRound(cv,VpGetRoundMode(),mx);
01497        return ToValue(cv);
01498     }
01499 }
01500 
01501 static VALUE
01502 BigDecimal_sub2(VALUE self, VALUE b, VALUE n)
01503 {
01504     ENTER(2);
01505     Real *cv;
01506     SIGNED_VALUE mx = GetPositiveInt(n);
01507     if (mx == 0) return BigDecimal_sub(self, b);
01508     else {
01509        size_t pl = VpSetPrecLimit(0);
01510        VALUE   c = BigDecimal_sub(self,b);
01511        VpSetPrecLimit(pl);
01512        GUARD_OBJ(cv,GetVpValue(c,1));
01513        VpLeftRound(cv,VpGetRoundMode(),mx);
01514        return ToValue(cv);
01515     }
01516 }
01517 
01518 static VALUE
01519 BigDecimal_mult2(VALUE self, VALUE b, VALUE n)
01520 {
01521     ENTER(2);
01522     Real *cv;
01523     SIGNED_VALUE mx = GetPositiveInt(n);
01524     if (mx == 0) return BigDecimal_mult(self, b);
01525     else {
01526        size_t pl = VpSetPrecLimit(0);
01527        VALUE   c = BigDecimal_mult(self,b);
01528        VpSetPrecLimit(pl);
01529        GUARD_OBJ(cv,GetVpValue(c,1));
01530        VpLeftRound(cv,VpGetRoundMode(),mx);
01531        return ToValue(cv);
01532     }
01533 }
01534 
01535 /* Returns the absolute value.
01536  *
01537  * BigDecimal('5').abs -> 5
01538  *
01539  * BigDecimal('-3').abs -> 3
01540  */
01541 static VALUE
01542 BigDecimal_abs(VALUE self)
01543 {
01544     ENTER(5);
01545     Real *c, *a;
01546     size_t mx;
01547 
01548     GUARD_OBJ(a,GetVpValue(self,1));
01549     mx = a->Prec *(VpBaseFig() + 1);
01550     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01551     VpAsgn(c, a, 1);
01552     VpChangeSign(c, 1);
01553     return ToValue(c);
01554 }
01555 
01556 /* call-seq:
01557  * sqrt(n)
01558  *
01559  * Returns the square root of the value.
01560  *
01561  * Result has at least n significant digits.
01562  */
01563 static VALUE
01564 BigDecimal_sqrt(VALUE self, VALUE nFig)
01565 {
01566     ENTER(5);
01567     Real *c, *a;
01568     size_t mx, n;
01569 
01570     GUARD_OBJ(a,GetVpValue(self,1));
01571     mx = a->Prec *(VpBaseFig() + 1);
01572 
01573     n = GetPositiveInt(nFig) + VpDblFig() + 1;
01574     if(mx <= n) mx = n;
01575     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01576     VpSqrt(c, a);
01577     return ToValue(c);
01578 }
01579 
01580 /* Return the integer part of the number.
01581  */
01582 static VALUE
01583 BigDecimal_fix(VALUE self)
01584 {
01585     ENTER(5);
01586     Real *c, *a;
01587     size_t mx;
01588 
01589     GUARD_OBJ(a,GetVpValue(self,1));
01590     mx = a->Prec *(VpBaseFig() + 1);
01591     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01592     VpActiveRound(c,a,VP_ROUND_DOWN,0); /* 0: round off */
01593     return ToValue(c);
01594 }
01595 
01596 /* call-seq:
01597  * round(n, mode)
01598  *
01599  * Round to the nearest 1 (by default), returning the result as a BigDecimal.
01600  *
01601  *      BigDecimal('3.14159').round #=> 3
01602  *      BigDecimal('8.7').round #=> 9
01603  *
01604  * If n is specified and positive, the fractional part of the result has no
01605  * more than that many digits.
01606  *
01607  * If n is specified and negative, at least that many digits to the left of the
01608  * decimal point will be 0 in the result.
01609  *
01610  *      BigDecimal('3.14159').round(3) #=> 3.142
01611  *      BigDecimal('13345.234').round(-2) #=> 13300.0
01612  *
01613  * The value of the optional mode argument can be used to determine how
01614  * rounding is performed; see BigDecimal.mode.
01615  */
01616 static VALUE
01617 BigDecimal_round(int argc, VALUE *argv, VALUE self)
01618 {
01619     ENTER(5);
01620     Real   *c, *a;
01621     int    iLoc = 0;
01622     VALUE  vLoc;
01623     VALUE  vRound;
01624     size_t mx, pl;
01625 
01626     unsigned short sw = VpGetRoundMode();
01627 
01628     switch (rb_scan_args(argc, argv, "02", &vLoc, &vRound)) {
01629     case 0:
01630         iLoc = 0;
01631         break;
01632     case 1:
01633         Check_Type(vLoc, T_FIXNUM);
01634         iLoc = FIX2INT(vLoc);
01635         break;
01636     case 2:
01637         Check_Type(vLoc, T_FIXNUM);
01638         iLoc = FIX2INT(vLoc);
01639         sw = check_rounding_mode(vRound);
01640         break;
01641     }
01642 
01643     pl = VpSetPrecLimit(0);
01644     GUARD_OBJ(a,GetVpValue(self,1));
01645     mx = a->Prec *(VpBaseFig() + 1);
01646     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01647     VpSetPrecLimit(pl);
01648     VpActiveRound(c,a,sw,iLoc);
01649     if (argc == 0) {
01650         return BigDecimal_to_i(ToValue(c));
01651     }
01652     return ToValue(c);
01653 }
01654 
01655 /* call-seq:
01656  * truncate(n)
01657  *
01658  * Truncate to the nearest 1, returning the result as a BigDecimal.
01659  *
01660  *      BigDecimal('3.14159').truncate #=> 3
01661  *      BigDecimal('8.7').truncate #=> 8
01662  *
01663  * If n is specified and positive, the fractional part of the result has no
01664  * more than that many digits.
01665  *
01666  * If n is specified and negative, at least that many digits to the left of the
01667  * decimal point will be 0 in the result.
01668  *
01669  *      BigDecimal('3.14159').truncate(3) #=> 3.141
01670  *      BigDecimal('13345.234').truncate(-2) #=> 13300.0
01671  */
01672 static VALUE
01673 BigDecimal_truncate(int argc, VALUE *argv, VALUE self)
01674 {
01675     ENTER(5);
01676     Real *c, *a;
01677     int iLoc;
01678     VALUE vLoc;
01679     size_t mx, pl = VpSetPrecLimit(0);
01680 
01681     if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
01682         iLoc = 0;
01683     } else {
01684         Check_Type(vLoc, T_FIXNUM);
01685         iLoc = FIX2INT(vLoc);
01686     }
01687 
01688     GUARD_OBJ(a,GetVpValue(self,1));
01689     mx = a->Prec *(VpBaseFig() + 1);
01690     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01691     VpSetPrecLimit(pl);
01692     VpActiveRound(c,a,VP_ROUND_DOWN,iLoc); /* 0: truncate */
01693     if (argc == 0) {
01694         return BigDecimal_to_i(ToValue(c));
01695     }
01696     return ToValue(c);
01697 }
01698 
01699 /* Return the fractional part of the number.
01700  */
01701 static VALUE
01702 BigDecimal_frac(VALUE self)
01703 {
01704     ENTER(5);
01705     Real *c, *a;
01706     size_t mx;
01707 
01708     GUARD_OBJ(a,GetVpValue(self,1));
01709     mx = a->Prec *(VpBaseFig() + 1);
01710     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01711     VpFrac(c, a);
01712     return ToValue(c);
01713 }
01714 
01715 /* call-seq:
01716  * floor(n)
01717  *
01718  * Return the largest integer less than or equal to the value, as a BigDecimal.
01719  *
01720  *      BigDecimal('3.14159').floor #=> 3
01721  *      BigDecimal('-9.1').floor #=> -10
01722  *
01723  * If n is specified and positive, the fractional part of the result has no
01724  * more than that many digits.
01725  *
01726  * If n is specified and negative, at least that
01727  * many digits to the left of the decimal point will be 0 in the result.
01728  *
01729  *      BigDecimal('3.14159').floor(3) #=> 3.141
01730  *      BigDecimal('13345.234').floor(-2) #=> 13300.0
01731  */
01732 static VALUE
01733 BigDecimal_floor(int argc, VALUE *argv, VALUE self)
01734 {
01735     ENTER(5);
01736     Real *c, *a;
01737     int iLoc;
01738     VALUE vLoc;
01739     size_t mx, pl = VpSetPrecLimit(0);
01740 
01741     if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
01742         iLoc = 0;
01743     } else {
01744         Check_Type(vLoc, T_FIXNUM);
01745         iLoc = FIX2INT(vLoc);
01746     }
01747 
01748     GUARD_OBJ(a,GetVpValue(self,1));
01749     mx = a->Prec *(VpBaseFig() + 1);
01750     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01751     VpSetPrecLimit(pl);
01752     VpActiveRound(c,a,VP_ROUND_FLOOR,iLoc);
01753 #ifdef BIGDECIMAL_DEBUG
01754     VPrint(stderr, "floor: c=%\n", c);
01755 #endif
01756     if (argc == 0) {
01757         return BigDecimal_to_i(ToValue(c));
01758     }
01759     return ToValue(c);
01760 }
01761 
01762 /* call-seq:
01763  * ceil(n)
01764  *
01765  * Return the smallest integer greater than or equal to the value, as a BigDecimal.
01766  *
01767  *      BigDecimal('3.14159').ceil #=> 4
01768  *      BigDecimal('-9.1').ceil #=> -9
01769  *
01770  * If n is specified and positive, the fractional part of the result has no
01771  * more than that many digits.
01772  *
01773  * If n is specified and negative, at least that
01774  * many digits to the left of the decimal point will be 0 in the result.
01775  *
01776  *      BigDecimal('3.14159').ceil(3) #=> 3.142
01777  *      BigDecimal('13345.234').ceil(-2) #=> 13400.0
01778  */
01779 static VALUE
01780 BigDecimal_ceil(int argc, VALUE *argv, VALUE self)
01781 {
01782     ENTER(5);
01783     Real *c, *a;
01784     int iLoc;
01785     VALUE vLoc;
01786     size_t mx, pl = VpSetPrecLimit(0);
01787 
01788     if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
01789         iLoc = 0;
01790     } else {
01791         Check_Type(vLoc, T_FIXNUM);
01792         iLoc = FIX2INT(vLoc);
01793     }
01794 
01795     GUARD_OBJ(a,GetVpValue(self,1));
01796     mx = a->Prec *(VpBaseFig() + 1);
01797     GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
01798     VpSetPrecLimit(pl);
01799     VpActiveRound(c,a,VP_ROUND_CEIL,iLoc);
01800     if (argc == 0) {
01801         return BigDecimal_to_i(ToValue(c));
01802     }
01803     return ToValue(c);
01804 }
01805 
01806 /* call-seq:
01807  * to_s(s)
01808  *
01809  * Converts the value to a string.
01810  *
01811  * The default format looks like  0.xxxxEnn.
01812  *
01813  * The optional parameter s consists of either an integer; or an optional '+'
01814  * or ' ', followed by an optional number, followed by an optional 'E' or 'F'.
01815  *
01816  * If there is a '+' at the start of s, positive values are returned with
01817  * a leading '+'.
01818  *
01819  * A space at the start of s returns positive values with a leading space.
01820  *
01821  * If s contains a number, a space is inserted after each group of that many
01822  * fractional digits.
01823  *
01824  * If s ends with an 'E', engineering notation (0.xxxxEnn) is used.
01825  *
01826  * If s ends with an 'F', conventional floating point notation is used.
01827  *
01828  * Examples:
01829  *
01830  *      BigDecimal.new('-123.45678901234567890').to_s('5F')
01831  *          #=> '-123.45678 90123 45678 9'
01832  *
01833  *      BigDecimal.new('123.45678901234567890').to_s('+8F')
01834  *          #=> '+123.45678901 23456789'
01835  *
01836  *      BigDecimal.new('123.45678901234567890').to_s(' F')
01837  *          #=> ' 123.4567890123456789'
01838  */
01839 static VALUE
01840 BigDecimal_to_s(int argc, VALUE *argv, VALUE self)
01841 {
01842     ENTER(5);
01843     int   fmt = 0;   /* 0:E format */
01844     int   fPlus = 0; /* =0:default,=1: set ' ' before digits ,set '+' before digits. */
01845     Real  *vp;
01846     volatile VALUE str;
01847     char  *psz;
01848     char   ch;
01849     size_t nc, mc = 0;
01850     VALUE  f;
01851 
01852     GUARD_OBJ(vp,GetVpValue(self,1));
01853 
01854     if (rb_scan_args(argc,argv,"01",&f)==1) {
01855         if (RB_TYPE_P(f, T_STRING)) {
01856             SafeStringValue(f);
01857             psz = RSTRING_PTR(f);
01858             if (*psz == ' ') {
01859                 fPlus = 1;
01860                 psz++;
01861             }
01862             else if (*psz == '+') {
01863                 fPlus = 2;
01864                 psz++;
01865             }
01866             while ((ch = *psz++) != 0) {
01867                 if (ISSPACE(ch)) {
01868                     continue;
01869                 }
01870                 if (!ISDIGIT(ch)) {
01871                     if (ch == 'F' || ch == 'f') {
01872                         fmt = 1; /* F format */
01873                     }
01874                     break;
01875                 }
01876                 mc = mc * 10 + ch - '0';
01877             }
01878         }
01879         else {
01880             mc = (size_t)GetPositiveInt(f);
01881         }
01882     }
01883     if (fmt) {
01884         nc = VpNumOfChars(vp, "F");
01885     }
01886     else {
01887         nc = VpNumOfChars(vp, "E");
01888     }
01889     if (mc > 0) {
01890         nc += (nc + mc - 1) / mc + 1;
01891     }
01892 
01893     str = rb_str_new(0, nc);
01894     psz = RSTRING_PTR(str);
01895 
01896     if (fmt) {
01897         VpToFString(vp, psz, mc, fPlus);
01898     }
01899     else {
01900         VpToString (vp, psz, mc, fPlus);
01901     }
01902     rb_str_resize(str, strlen(psz));
01903     return str;
01904 }
01905 
01906 /* Splits a BigDecimal number into four parts, returned as an array of values.
01907  *
01908  * The first value represents the sign of the BigDecimal, and is -1 or 1, or 0
01909  * if the BigDecimal is Not a Number.
01910  *
01911  * The second value is a string representing the significant digits of the
01912  * BigDecimal, with no leading zeros.
01913  *
01914  * The third value is the base used for arithmetic (currently always 10) as an
01915  * Integer.
01916  *
01917  * The fourth value is an Integer exponent.
01918  *
01919  * If the BigDecimal can be represented as 0.xxxxxx*10**n, then xxxxxx is the
01920  * string of significant digits with no leading zeros, and n is the exponent.
01921  *
01922  * From these values, you can translate a BigDecimal to a float as follows:
01923  *
01924  *   sign, significant_digits, base, exponent = a.split
01925  *   f = sign * "0.#{significant_digits}".to_f * (base ** exponent)
01926  *
01927  * (Note that the to_f method is provided as a more convenient way to translate
01928  * a BigDecimal to a Float.)
01929  */
01930 static VALUE
01931 BigDecimal_split(VALUE self)
01932 {
01933     ENTER(5);
01934     Real *vp;
01935     VALUE obj,str;
01936     ssize_t e, s;
01937     char *psz1;
01938 
01939     GUARD_OBJ(vp,GetVpValue(self,1));
01940     str = rb_str_new(0, VpNumOfChars(vp,"E"));
01941     psz1 = RSTRING_PTR(str);
01942     VpSzMantissa(vp,psz1);
01943     s = 1;
01944     if(psz1[0]=='-') {
01945         size_t len = strlen(psz1+1);
01946 
01947         memmove(psz1, psz1+1, len);
01948         psz1[len] = '\0';
01949         s = -1;
01950     }
01951     if(psz1[0]=='N') s=0; /* NaN */
01952     e = VpExponent10(vp);
01953     obj  = rb_ary_new2(4);
01954     rb_ary_push(obj, INT2FIX(s));
01955     rb_ary_push(obj, str);
01956     rb_str_resize(str, strlen(psz1));
01957     rb_ary_push(obj, INT2FIX(10));
01958     rb_ary_push(obj, INT2NUM(e));
01959     return obj;
01960 }
01961 
01962 /* Returns the exponent of the BigDecimal number, as an Integer.
01963  *
01964  * If the number can be represented as 0.xxxxxx*10**n where xxxxxx is a string
01965  * of digits with no leading zeros, then n is the exponent.
01966  */
01967 static VALUE
01968 BigDecimal_exponent(VALUE self)
01969 {
01970     ssize_t e = VpExponent10(GetVpValue(self, 1));
01971     return INT2NUM(e);
01972 }
01973 
01974 /* Returns debugging information about the value as a string of comma-separated
01975  * values in angle brackets with a leading #:
01976  *
01977  * BigDecimal.new("1234.5678").inspect ->
01978  * "#<BigDecimal:b7ea1130,'0.12345678E4',8(12)>"
01979  *
01980  * The first part is the address, the second is the value as a string, and
01981  * the final part ss(mm) is the current number of significant digits and the
01982  * maximum number of significant digits, respectively.
01983  */
01984 static VALUE
01985 BigDecimal_inspect(VALUE self)
01986 {
01987     ENTER(5);
01988     Real *vp;
01989     volatile VALUE obj;
01990     size_t nc;
01991     char *psz, *tmp;
01992 
01993     GUARD_OBJ(vp,GetVpValue(self,1));
01994     nc = VpNumOfChars(vp,"E");
01995     nc +=(nc + 9) / 10;
01996 
01997     obj = rb_str_new(0, nc+256);
01998     psz = RSTRING_PTR(obj);
01999     sprintf(psz,"#<BigDecimal:%"PRIxVALUE",'",self);
02000     tmp = psz + strlen(psz);
02001     VpToString(vp, tmp, 10, 0);
02002     tmp += strlen(tmp);
02003     sprintf(tmp, "',%"PRIuSIZE"(%"PRIuSIZE")>", VpPrec(vp)*VpBaseFig(), VpMaxPrec(vp)*VpBaseFig());
02004     rb_str_resize(obj, strlen(psz));
02005     return obj;
02006 }
02007 
02008 static VALUE BigMath_s_exp(VALUE, VALUE, VALUE);
02009 static VALUE BigMath_s_log(VALUE, VALUE, VALUE);
02010 
02011 #define BigMath_exp(x, n) BigMath_s_exp(rb_mBigMath, (x), (n))
02012 #define BigMath_log(x, n) BigMath_s_log(rb_mBigMath, (x), (n))
02013 
02014 inline static int
02015 is_integer(VALUE x)
02016 {
02017     return (RB_TYPE_P(x, T_FIXNUM) || RB_TYPE_P(x, T_BIGNUM));
02018 }
02019 
02020 inline static int
02021 is_negative(VALUE x)
02022 {
02023     if (FIXNUM_P(x)) {
02024         return FIX2LONG(x) < 0;
02025     }
02026     else if (RB_TYPE_P(x, T_BIGNUM)) {
02027         return RBIGNUM_NEGATIVE_P(x);
02028     }
02029     else if (RB_TYPE_P(x, T_FLOAT)) {
02030         return RFLOAT_VALUE(x) < 0.0;
02031     }
02032     return RTEST(rb_funcall(x, '<', 1, INT2FIX(0)));
02033 }
02034 
02035 #define is_positive(x) (!is_negative(x))
02036 
02037 inline static int
02038 is_zero(VALUE x)
02039 {
02040     VALUE num;
02041 
02042     switch (TYPE(x)) {
02043       case T_FIXNUM:
02044         return FIX2LONG(x) == 0;
02045 
02046       case T_BIGNUM:
02047         return Qfalse;
02048 
02049       case T_RATIONAL:
02050         num = RRATIONAL(x)->num;
02051         return FIXNUM_P(num) && FIX2LONG(num) == 0;
02052 
02053       default:
02054         break;
02055     }
02056 
02057     return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(0)));
02058 }
02059 
02060 inline static int
02061 is_one(VALUE x)
02062 {
02063     VALUE num, den;
02064 
02065     switch (TYPE(x)) {
02066       case T_FIXNUM:
02067         return FIX2LONG(x) == 1;
02068 
02069       case T_BIGNUM:
02070         return Qfalse;
02071 
02072       case T_RATIONAL:
02073         num = RRATIONAL(x)->num;
02074         den = RRATIONAL(x)->den;
02075         return FIXNUM_P(den) && FIX2LONG(den) == 1 &&
02076                FIXNUM_P(num) && FIX2LONG(num) == 1;
02077 
02078       default:
02079         break;
02080     }
02081 
02082     return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(1)));
02083 }
02084 
02085 inline static int
02086 is_even(VALUE x)
02087 {
02088     switch (TYPE(x)) {
02089       case T_FIXNUM:
02090         return (FIX2LONG(x) % 2) == 0;
02091 
02092       case T_BIGNUM:
02093         return (RBIGNUM_DIGITS(x)[0] % 2) == 0;
02094 
02095       default:
02096         break;
02097     }
02098 
02099     return 0;
02100 }
02101 
02102 static VALUE
02103 rmpd_power_by_big_decimal(Real const* x, Real const* exp, ssize_t const n)
02104 {
02105     VALUE log_x, multiplied, y;
02106     volatile VALUE obj = exp->obj;
02107 
02108     if (VpIsZero(exp)) {
02109         return ToValue(VpCreateRbObject(n, "1"));
02110     }
02111 
02112     log_x = BigMath_log(x->obj, SSIZET2NUM(n+1));
02113     multiplied = BigDecimal_mult2(exp->obj, log_x, SSIZET2NUM(n+1));
02114     y = BigMath_exp(multiplied, SSIZET2NUM(n));
02115     RB_GC_GUARD(obj);
02116 
02117     return y;
02118 }
02119 
02120 /* call-seq:
02121  * power(n)
02122  * power(n, prec)
02123  *
02124  * Returns the value raised to the power of n.
02125  *
02126  * Note that n must be an Integer.
02127  *
02128  * Also available as the operator **
02129  */
02130 static VALUE
02131 BigDecimal_power(int argc, VALUE*argv, VALUE self)
02132 {
02133     ENTER(5);
02134     VALUE vexp, prec;
02135     Real* exp = NULL;
02136     Real *x, *y;
02137     ssize_t mp, ma, n;
02138     SIGNED_VALUE int_exp;
02139     double d;
02140 
02141     rb_scan_args(argc, argv, "11", &vexp, &prec);
02142 
02143     GUARD_OBJ(x, GetVpValue(self, 1));
02144     n = NIL_P(prec) ? (ssize_t)(x->Prec*VpBaseFig()) : NUM2SSIZET(prec);
02145 
02146     if (VpIsNaN(x)) {
02147         y = VpCreateRbObject(n, "0#");
02148         RB_GC_GUARD(y->obj);
02149         VpSetNaN(y);
02150         return ToValue(y);
02151     }
02152 
02153 retry:
02154     switch (TYPE(vexp)) {
02155       case T_FIXNUM:
02156         break;
02157 
02158       case T_BIGNUM:
02159         break;
02160 
02161       case T_FLOAT:
02162         d = RFLOAT_VALUE(vexp);
02163         if (d == round(d)) {
02164             vexp = LL2NUM((LONG_LONG)round(d));
02165             goto retry;
02166         }
02167         exp = GetVpValueWithPrec(vexp, DBL_DIG+1, 1);
02168         break;
02169 
02170       case T_RATIONAL:
02171         if (is_zero(RRATIONAL(vexp)->num)) {
02172             if (is_positive(vexp)) {
02173                 vexp = INT2FIX(0);
02174                 goto retry;
02175             }
02176         }
02177         else if (is_one(RRATIONAL(vexp)->den)) {
02178             vexp = RRATIONAL(vexp)->num;
02179             goto retry;
02180         }
02181         exp = GetVpValueWithPrec(vexp, n, 1);
02182         break;
02183 
02184       case T_DATA:
02185         if (is_kind_of_BigDecimal(vexp)) {
02186             VALUE zero = INT2FIX(0);
02187             VALUE rounded = BigDecimal_round(1, &zero, vexp);
02188             if (RTEST(BigDecimal_eq(vexp, rounded))) {
02189                 vexp = BigDecimal_to_i(vexp);
02190                 goto retry;
02191             }
02192             exp = DATA_PTR(vexp);
02193             break;
02194         }
02195         /* fall through */
02196       default:
02197         rb_raise(rb_eTypeError,
02198                  "wrong argument type %s (expected scalar Numeric)",
02199                  rb_obj_classname(vexp));
02200     }
02201 
02202     if (VpIsZero(x)) {
02203         if (is_negative(vexp)) {
02204             y = VpCreateRbObject(n, "#0");
02205             RB_GC_GUARD(y->obj);
02206             if (VpGetSign(x) < 0) {
02207                 if (is_integer(vexp)) {
02208                     if (is_even(vexp)) {
02209                         /* (-0) ** (-even_integer)  -> Infinity */
02210                         VpSetPosInf(y);
02211                     }
02212                     else {
02213                         /* (-0) ** (-odd_integer)  -> -Infinity */
02214                         VpSetNegInf(y);
02215                     }
02216                 }
02217                 else {
02218                     /* (-0) ** (-non_integer)  -> Infinity */
02219                     VpSetPosInf(y);
02220                 }
02221             }
02222             else {
02223                 /* (+0) ** (-num)  -> Infinity */
02224                 VpSetPosInf(y);
02225             }
02226             return ToValue(y);
02227         }
02228         else if (is_zero(vexp)) {
02229             return ToValue(VpCreateRbObject(n, "1"));
02230         }
02231         else {
02232             return ToValue(VpCreateRbObject(n, "0"));
02233         }
02234     }
02235 
02236     if (is_zero(vexp)) {
02237         return ToValue(VpCreateRbObject(n, "1"));
02238     }
02239     else if (is_one(vexp)) {
02240         return self;
02241     }
02242 
02243     if (VpIsInf(x)) {
02244         if (is_negative(vexp)) {
02245             if (VpGetSign(x) < 0) {
02246                 if (is_integer(vexp)) {
02247                     if (is_even(vexp)) {
02248                         /* (-Infinity) ** (-even_integer) -> +0 */
02249                         return ToValue(VpCreateRbObject(n, "0"));
02250                     }
02251                     else {
02252                         /* (-Infinity) ** (-odd_integer) -> -0 */
02253                         return ToValue(VpCreateRbObject(n, "-0"));
02254                     }
02255                 }
02256                 else {
02257                     /* (-Infinity) ** (-non_integer) -> -0 */
02258                     return ToValue(VpCreateRbObject(n, "-0"));
02259                 }
02260             }
02261             else {
02262                 return ToValue(VpCreateRbObject(n, "0"));
02263             }
02264         }
02265         else {
02266             y = VpCreateRbObject(n, "0#");
02267             if (VpGetSign(x) < 0) {
02268                 if (is_integer(vexp)) {
02269                     if (is_even(vexp)) {
02270                         VpSetPosInf(y);
02271                     }
02272                     else {
02273                         VpSetNegInf(y);
02274                     }
02275                 }
02276                 else {
02277                     /* TODO: support complex */
02278                     rb_raise(rb_eMathDomainError,
02279                              "a non-integral exponent for a negative base");
02280                 }
02281             }
02282             else {
02283                 VpSetPosInf(y);
02284             }
02285             return ToValue(y);
02286         }
02287     }
02288 
02289     if (exp != NULL) {
02290         return rmpd_power_by_big_decimal(x, exp, n);
02291     }
02292     else if (RB_TYPE_P(vexp, T_BIGNUM)) {
02293         VALUE abs_value = BigDecimal_abs(self);
02294         if (is_one(abs_value)) {
02295             return ToValue(VpCreateRbObject(n, "1"));
02296         }
02297         else if (RTEST(rb_funcall(abs_value, '<', 1, INT2FIX(1)))) {
02298             if (is_negative(vexp)) {
02299                 y = VpCreateRbObject(n, "0#");
02300                 if (is_even(vexp)) {
02301                     VpSetInf(y, VpGetSign(x));
02302                 }
02303                 else {
02304                     VpSetInf(y, -VpGetSign(x));
02305                 }
02306                 return ToValue(y);
02307             }
02308             else if (VpGetSign(x) < 0 && is_even(vexp)) {
02309                 return ToValue(VpCreateRbObject(n, "-0"));
02310             }
02311             else {
02312                 return ToValue(VpCreateRbObject(n, "0"));
02313             }
02314         }
02315         else {
02316             if (is_positive(vexp)) {
02317                 y = VpCreateRbObject(n, "0#");
02318                 if (is_even(vexp)) {
02319                     VpSetInf(y, VpGetSign(x));
02320                 }
02321                 else {
02322                     VpSetInf(y, -VpGetSign(x));
02323                 }
02324                 return ToValue(y);
02325             }
02326             else if (VpGetSign(x) < 0 && is_even(vexp)) {
02327                 return ToValue(VpCreateRbObject(n, "-0"));
02328             }
02329             else {
02330                 return ToValue(VpCreateRbObject(n, "0"));
02331             }
02332         }
02333     }
02334 
02335     int_exp = FIX2INT(vexp);
02336     ma = int_exp;
02337     if (ma < 0)  ma = -ma;
02338     if (ma == 0) ma = 1;
02339 
02340     if (VpIsDef(x)) {
02341         mp = x->Prec * (VpBaseFig() + 1);
02342         GUARD_OBJ(y, VpCreateRbObject(mp * (ma + 1), "0"));
02343     }
02344     else {
02345         GUARD_OBJ(y, VpCreateRbObject(1, "0"));
02346     }
02347     VpPower(y, x, int_exp);
02348     return ToValue(y);
02349 }
02350 
02351 /* call-seq:
02352  *   big_decimal ** exp  -> big_decimal
02353  *
02354  * It is a synonym of BigDecimal#power(exp).
02355  */
02356 static VALUE
02357 BigDecimal_power_op(VALUE self, VALUE exp)
02358 {
02359     return BigDecimal_power(1, &exp, self);
02360 }
02361 
02362 static VALUE
02363 BigDecimal_s_allocate(VALUE klass)
02364 {
02365     return VpNewRbClass(0, NULL, klass)->obj;
02366 }
02367 
02368 static Real *BigDecimal_new(int argc, VALUE *argv);
02369 
02370 /* call-seq:
02371  *   new(initial, digits)
02372  *
02373  * Create a new BigDecimal object.
02374  *
02375  * initial:: The initial value, as an Integer, a Float, a Rational,
02376  *           a BigDecimal, or a String.
02377  *
02378  *           If it is a String, spaces are ignored and unrecognized characters
02379  *           terminate the value.
02380  *
02381  * digits:: The number of significant digits, as a Fixnum. If omitted or 0,
02382  *          the number of significant digits is determined from the initial
02383  *          value.
02384  *
02385  * The actual number of significant digits used in computation is usually
02386  * larger than the specified number.
02387  */
02388 static VALUE
02389 BigDecimal_initialize(int argc, VALUE *argv, VALUE self)
02390 {
02391     Real *pv = rb_check_typeddata(self, &BigDecimal_data_type);
02392     Real *x = BigDecimal_new(argc, argv);
02393 
02394     if (ToValue(x)) {
02395         pv = VpCopy(pv, x);
02396     }
02397     else {
02398         VpFree(pv);
02399         pv = x;
02400     }
02401     DATA_PTR(self) = pv;
02402     pv->obj = self;
02403     return self;
02404 }
02405 
02406 /* :nodoc:
02407  *
02408  * private method to dup and clone the provided BigDecimal +other+
02409  */
02410 static VALUE
02411 BigDecimal_initialize_copy(VALUE self, VALUE other)
02412 {
02413     Real *pv = rb_check_typeddata(self, &BigDecimal_data_type);
02414     Real *x = rb_check_typeddata(other, &BigDecimal_data_type);
02415 
02416     DATA_PTR(self) = VpCopy(pv, x);
02417     return self;
02418 }
02419 
02420 static Real *
02421 BigDecimal_new(int argc, VALUE *argv)
02422 {
02423     size_t mf;
02424     VALUE  nFig;
02425     VALUE  iniValue;
02426 
02427     if (rb_scan_args(argc, argv, "11", &iniValue, &nFig) == 1) {
02428         mf = 0;
02429     }
02430     else {
02431         mf = GetPositiveInt(nFig);
02432     }
02433 
02434     switch (TYPE(iniValue)) {
02435       case T_DATA:
02436         if (is_kind_of_BigDecimal(iniValue)) {
02437             return DATA_PTR(iniValue);
02438         }
02439         break;
02440 
02441       case T_FIXNUM:
02442         /* fall through */
02443       case T_BIGNUM:
02444         return GetVpValue(iniValue, 1);
02445 
02446       case T_FLOAT:
02447         if (mf > DBL_DIG+1) {
02448             rb_raise(rb_eArgError, "precision too large.");
02449         }
02450         /* fall through */
02451       case T_RATIONAL:
02452         if (NIL_P(nFig)) {
02453             rb_raise(rb_eArgError,
02454                      "can't omit precision for a %s.",
02455                      rb_class2name(CLASS_OF(iniValue)));
02456         }
02457         return GetVpValueWithPrec(iniValue, mf, 1);
02458 
02459       case T_STRING:
02460         /* fall through */
02461       default:
02462         break;
02463     }
02464     StringValueCStr(iniValue);
02465     return VpAlloc(mf, RSTRING_PTR(iniValue));
02466 }
02467 
02468 static VALUE
02469 BigDecimal_global_new(int argc, VALUE *argv, VALUE self)
02470 {
02471     Real *pv = BigDecimal_new(argc, argv);
02472     if (ToValue(pv)) pv = VpCopy(NULL, pv);
02473     pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv);
02474     return pv->obj;
02475 }
02476 
02477  /* call-seq:
02478   * BigDecimal.limit(digits)
02479   *
02480   * Limit the number of significant digits in newly created BigDecimal
02481   * numbers to the specified value. Rounding is performed as necessary,
02482   * as specified by BigDecimal.mode.
02483   *
02484   * A limit of 0, the default, means no upper limit.
02485   *
02486   * The limit specified by this method takes less priority over any limit
02487   * specified to instance methods such as ceil, floor, truncate, or round.
02488   */
02489 static VALUE
02490 BigDecimal_limit(int argc, VALUE *argv, VALUE self)
02491 {
02492     VALUE  nFig;
02493     VALUE  nCur = INT2NUM(VpGetPrecLimit());
02494 
02495     if(rb_scan_args(argc,argv,"01",&nFig)==1) {
02496         int nf;
02497         if(nFig==Qnil) return nCur;
02498         Check_Type(nFig, T_FIXNUM);
02499         nf = FIX2INT(nFig);
02500         if(nf<0) {
02501             rb_raise(rb_eArgError, "argument must be positive");
02502         }
02503         VpSetPrecLimit(nf);
02504     }
02505     return nCur;
02506 }
02507 
02508 /* Returns the sign of the value.
02509  *
02510  * Returns a positive value if > 0, a negative value if < 0, and a
02511  * zero if == 0.
02512  *
02513  * The specific value returned indicates the type and sign of the BigDecimal,
02514  * as follows:
02515  *
02516  * BigDecimal::SIGN_NaN:: value is Not a Number
02517  * BigDecimal::SIGN_POSITIVE_ZERO:: value is +0
02518  * BigDecimal::SIGN_NEGATIVE_ZERO:: value is -0
02519  * BigDecimal::SIGN_POSITIVE_INFINITE:: value is +infinity
02520  * BigDecimal::SIGN_NEGATIVE_INFINITE:: value is -infinity
02521  * BigDecimal::SIGN_POSITIVE_FINITE:: value is positive
02522  * BigDecimal::SIGN_NEGATIVE_FINITE:: value is negative
02523  */
02524 static VALUE
02525 BigDecimal_sign(VALUE self)
02526 { /* sign */
02527     int s = GetVpValue(self,1)->sign;
02528     return INT2FIX(s);
02529 }
02530 
02531 /*
02532  * call-seq: BigDecimal.save_exception_mode { ... }
02533  *
02534  * Excecute the provided block, but preserve the exception mode
02535  *
02536  *     BigDecimal.save_exception_mode do
02537  *       BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false)
02538  *       BigDecimal.mode(BigDecimal::EXCEPTION_NaN, false)
02539  *
02540  *       BigDecimal.new(BigDecimal('Infinity'))
02541  *       BigDecimal.new(BigDecimal('-Infinity'))
02542  *       BigDecimal(BigDecimal.new('NaN'))
02543  *     end
02544  *
02545  * For use with the BigDecimal::EXCEPTION_*
02546  *
02547  * See BigDecimal.mode
02548  */
02549 static VALUE
02550 BigDecimal_save_exception_mode(VALUE self)
02551 {
02552     unsigned short const exception_mode = VpGetException();
02553     int state;
02554     VALUE ret = rb_protect(rb_yield, Qnil, &state);
02555     VpSetException(exception_mode);
02556     if (state) rb_jump_tag(state);
02557     return ret;
02558 }
02559 
02560 /*
02561  * call-seq: BigDecimal.save_rounding_mode { ... }
02562  *
02563  * Excecute the provided block, but preserve the rounding mode
02564  *
02565  *     BigDecimal.save_exception_mode do
02566  *       BigDecimal.mode(BigDecimal::ROUND_MODE, :up)
02567  *       puts BigDecimal.mode(BigDecimal::ROUND_MODE)
02568  *     end
02569  *
02570  * For use with the BigDecimal::ROUND_*
02571  *
02572  * See BigDecimal.mode
02573  */
02574 static VALUE
02575 BigDecimal_save_rounding_mode(VALUE self)
02576 {
02577     unsigned short const round_mode = VpGetRoundMode();
02578     int state;
02579     VALUE ret = rb_protect(rb_yield, Qnil, &state);
02580     VpSetRoundMode(round_mode);
02581     if (state) rb_jump_tag(state);
02582     return ret;
02583 }
02584 
02585 /*
02586  * call-seq: BigDecimal.save_limit { ... }
02587  *
02588  * Excecute the provided block, but preserve the precision limit
02589  *
02590  *      BigDecimal.limit(100)
02591  *      puts BigDecimal.limit
02592  *      BigDecimal.save_limit do
02593  *          BigDecimal.limit(200)
02594  *          puts BigDecimal.limit
02595  *      end
02596  *      puts BigDecimal.limit
02597  *
02598  */
02599 static VALUE
02600 BigDecimal_save_limit(VALUE self)
02601 {
02602     size_t const limit = VpGetPrecLimit();
02603     int state;
02604     VALUE ret = rb_protect(rb_yield, Qnil, &state);
02605     VpSetPrecLimit(limit);
02606     if (state) rb_jump_tag(state);
02607     return ret;
02608 }
02609 
02610 /* call-seq:
02611  * BigMath.exp(x, prec)
02612  *
02613  * Computes the value of e (the base of natural logarithms) raised to the
02614  * power of x, to the specified number of digits of precision.
02615  *
02616  * If x is infinity, returns Infinity.
02617  *
02618  * If x is NaN, returns NaN.
02619  */
02620 static VALUE
02621 BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
02622 {
02623     ssize_t prec, n, i;
02624     Real* vx = NULL;
02625     VALUE one, d, x1, y, z;
02626     int negative = 0;
02627     int infinite = 0;
02628     int nan = 0;
02629     double flo;
02630 
02631     prec = NUM2SSIZET(vprec);
02632     if (prec <= 0) {
02633         rb_raise(rb_eArgError, "Zero or negative precision for exp");
02634     }
02635 
02636     /* TODO: the following switch statement is almostly the same as one in the
02637      *       BigDecimalCmp function. */
02638     switch (TYPE(x)) {
02639       case T_DATA:
02640           if (!is_kind_of_BigDecimal(x)) break;
02641           vx = DATA_PTR(x);
02642           negative = VpGetSign(vx) < 0;
02643           infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
02644           nan = VpIsNaN(vx);
02645           break;
02646 
02647       case T_FIXNUM:
02648           /* fall through */
02649       case T_BIGNUM:
02650           vx = GetVpValue(x, 0);
02651           break;
02652 
02653       case T_FLOAT:
02654         flo = RFLOAT_VALUE(x);
02655         negative = flo < 0;
02656         infinite = isinf(flo);
02657         nan = isnan(flo);
02658         if (!infinite && !nan) {
02659             vx = GetVpValueWithPrec(x, DBL_DIG+1, 0);
02660         }
02661         break;
02662 
02663       case T_RATIONAL:
02664         vx = GetVpValueWithPrec(x, prec, 0);
02665         break;
02666 
02667       default:
02668         break;
02669     }
02670     if (infinite) {
02671         if (negative) {
02672             return ToValue(GetVpValueWithPrec(INT2NUM(0), prec, 1));
02673         }
02674         else {
02675             Real* vy;
02676             vy = VpCreateRbObject(prec, "#0");
02677             RB_GC_GUARD(vy->obj);
02678             VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
02679             return ToValue(vy);
02680         }
02681     }
02682     else if (nan) {
02683         Real* vy;
02684         vy = VpCreateRbObject(prec, "#0");
02685         RB_GC_GUARD(vy->obj);
02686         VpSetNaN(vy);
02687         return ToValue(vy);
02688     }
02689     else if (vx == NULL) {
02690         cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
02691     }
02692     RB_GC_GUARD(vx->obj);
02693 
02694     n = prec + rmpd_double_figures();
02695     negative = VpGetSign(vx) < 0;
02696     if (negative) {
02697         VpSetSign(vx, 1);
02698     }
02699 
02700     RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
02701     RB_GC_GUARD(x1) = one;
02702     RB_GC_GUARD(y)  = one;
02703     RB_GC_GUARD(d)  = y;
02704     RB_GC_GUARD(z)  = one;
02705     i  = 0;
02706 
02707     while (!VpIsZero((Real*)DATA_PTR(d))) {
02708         VALUE argv[2];
02709         SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
02710         SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
02711         ssize_t m = n - vabs(ey - ed);
02712         if (m <= 0) {
02713             break;
02714         }
02715         else if ((size_t)m < rmpd_double_figures()) {
02716             m = rmpd_double_figures();
02717         }
02718 
02719         x1 = BigDecimal_mult2(x1, x, SSIZET2NUM(n));
02720         ++i;
02721         z = BigDecimal_mult(z, SSIZET2NUM(i));
02722         argv[0] = z;
02723         argv[1] = SSIZET2NUM(m);
02724         d = BigDecimal_div2(2, argv, x1);
02725         y = BigDecimal_add(y, d);
02726     }
02727 
02728     if (negative) {
02729         VALUE argv[2];
02730         argv[0] = y;
02731         argv[1] = vprec;
02732         return BigDecimal_div2(2, argv, one);
02733     }
02734     else {
02735         vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
02736         return BigDecimal_round(1, &vprec, y);
02737     }
02738 }
02739 
02740 /* call-seq:
02741  * BigMath.log(x, prec)
02742  *
02743  * Computes the natural logarithm of x to the specified number of digits of
02744  * precision.
02745  *
02746  * If x is zero or negative, raises Math::DomainError.
02747  *
02748  * If x is positive infinity, returns Infinity.
02749  *
02750  * If x is NaN, returns NaN.
02751  */
02752 static VALUE
02753 BigMath_s_log(VALUE klass, VALUE x, VALUE vprec)
02754 {
02755     ssize_t prec, n, i;
02756     SIGNED_VALUE expo;
02757     Real* vx = NULL;
02758     VALUE argv[2], vn, one, two, w, x2, y, d;
02759     int zero = 0;
02760     int negative = 0;
02761     int infinite = 0;
02762     int nan = 0;
02763     double flo;
02764     long fix;
02765 
02766     if (!is_integer(vprec)) {
02767         rb_raise(rb_eArgError, "precision must be an Integer");
02768     }
02769 
02770     prec = NUM2SSIZET(vprec);
02771     if (prec <= 0) {
02772         rb_raise(rb_eArgError, "Zero or negative precision for exp");
02773     }
02774 
02775     /* TODO: the following switch statement is almostly the same as one in the
02776      *       BigDecimalCmp function. */
02777     switch (TYPE(x)) {
02778       case T_DATA:
02779           if (!is_kind_of_BigDecimal(x)) break;
02780           vx = DATA_PTR(x);
02781           zero = VpIsZero(vx);
02782           negative = VpGetSign(vx) < 0;
02783           infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
02784           nan = VpIsNaN(vx);
02785           break;
02786 
02787       case T_FIXNUM:
02788         fix = FIX2LONG(x);
02789         zero = fix == 0;
02790         negative = fix < 0;
02791         goto get_vp_value;
02792 
02793       case T_BIGNUM:
02794         zero = RBIGNUM_ZERO_P(x);
02795         negative = RBIGNUM_NEGATIVE_P(x);
02796 get_vp_value:
02797         if (zero || negative) break;
02798         vx = GetVpValue(x, 0);
02799         break;
02800 
02801       case T_FLOAT:
02802         flo = RFLOAT_VALUE(x);
02803         zero = flo == 0;
02804         negative = flo < 0;
02805         infinite = isinf(flo);
02806         nan = isnan(flo);
02807         if (!zero && !negative && !infinite && !nan) {
02808             vx = GetVpValueWithPrec(x, DBL_DIG+1, 1);
02809         }
02810         break;
02811 
02812       case T_RATIONAL:
02813         zero = RRATIONAL_ZERO_P(x);
02814         negative = RRATIONAL_NEGATIVE_P(x);
02815         if (zero || negative) break;
02816         vx = GetVpValueWithPrec(x, prec, 1);
02817         break;
02818 
02819       case T_COMPLEX:
02820         rb_raise(rb_eMathDomainError,
02821                  "Complex argument for BigMath.log");
02822 
02823       default:
02824         break;
02825     }
02826     if (infinite && !negative) {
02827         Real* vy;
02828         vy = VpCreateRbObject(prec, "#0");
02829         RB_GC_GUARD(vy->obj);
02830         VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
02831         return ToValue(vy);
02832     }
02833     else if (nan) {
02834         Real* vy;
02835         vy = VpCreateRbObject(prec, "#0");
02836         RB_GC_GUARD(vy->obj);
02837         VpSetNaN(vy);
02838         return ToValue(vy);
02839     }
02840     else if (zero || negative) {
02841         rb_raise(rb_eMathDomainError,
02842                  "Zero or negative argument for log");
02843     }
02844     else if (vx == NULL) {
02845         cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
02846     }
02847     x = ToValue(vx);
02848 
02849     RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
02850     RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2"));
02851 
02852     n = prec + rmpd_double_figures();
02853     RB_GC_GUARD(vn) = SSIZET2NUM(n);
02854     expo = VpExponent10(vx);
02855     if (expo < 0 || expo >= 3) {
02856         char buf[16];
02857         snprintf(buf, 16, "1E%"PRIdVALUE, -expo);
02858         x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
02859     }
02860     else {
02861         expo = 0;
02862     }
02863     w = BigDecimal_sub(x, one);
02864     argv[0] = BigDecimal_add(x, one);
02865     argv[1] = vn;
02866     x = BigDecimal_div2(2, argv, w);
02867     RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn);
02868     RB_GC_GUARD(y)  = x;
02869     RB_GC_GUARD(d)  = y;
02870     i = 1;
02871     while (!VpIsZero((Real*)DATA_PTR(d))) {
02872         SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
02873         SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
02874         ssize_t m = n - vabs(ey - ed);
02875         if (m <= 0) {
02876             break;
02877         }
02878         else if ((size_t)m < rmpd_double_figures()) {
02879             m = rmpd_double_figures();
02880         }
02881 
02882         x = BigDecimal_mult2(x2, x, vn);
02883         i += 2;
02884         argv[0] = SSIZET2NUM(i);
02885         argv[1] = SSIZET2NUM(m);
02886         d = BigDecimal_div2(2, argv, x);
02887         y = BigDecimal_add(y, d);
02888     }
02889 
02890     y = BigDecimal_mult(y, two);
02891     if (expo != 0) {
02892         VALUE log10, vexpo, dy;
02893         log10 = BigMath_s_log(klass, INT2FIX(10), vprec);
02894         vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1));
02895         dy = BigDecimal_mult(log10, vexpo);
02896         y = BigDecimal_add(y, dy);
02897     }
02898 
02899     return y;
02900 }
02901 
02902 /* Document-class: BigDecimal
02903  * BigDecimal provides arbitrary-precision floating point decimal arithmetic.
02904  *
02905  * Copyright (C) 2002 by Shigeo Kobayashi <shigeo@tinyforest.gr.jp>.
02906  *
02907  * You may distribute under the terms of either the GNU General Public
02908  * License or the Artistic License, as specified in the README file
02909  * of the BigDecimal distribution.
02910  *
02911  * Documented by mathew <meta@pobox.com>.
02912  *
02913  * = Introduction
02914  *
02915  * Ruby provides built-in support for arbitrary precision integer arithmetic.
02916  *
02917  * For example:
02918  *
02919  *      42**13  #=>   1265437718438866624512
02920  *
02921  * BigDecimal provides similar support for very large or very accurate floating
02922  * point numbers.
02923  *
02924  * Decimal arithmetic is also useful for general calculation, because it
02925  * provides the correct answers people expect--whereas normal binary floating
02926  * point arithmetic often introduces subtle errors because of the conversion
02927  * between base 10 and base 2.
02928  *
02929  * For example, try:
02930  *
02931  *   sum = 0
02932  *   for i in (1..10000)
02933  *     sum = sum + 0.0001
02934  *   end
02935  *   print sum #=> 0.9999999999999062
02936  *
02937  * and contrast with the output from:
02938  *
02939  *   require 'bigdecimal'
02940  *
02941  *   sum = BigDecimal.new("0")
02942  *   for i in (1..10000)
02943  *     sum = sum + BigDecimal.new("0.0001")
02944  *   end
02945  *   print sum #=> 0.1E1
02946  *
02947  * Similarly:
02948  *
02949  *      (BigDecimal.new("1.2") - BigDecimal("1.0")) == BigDecimal("0.2") #=> true
02950  *
02951  *      (1.2 - 1.0) == 0.2 #=> false
02952  *
02953  * = Special features of accurate decimal arithmetic
02954  *
02955  * Because BigDecimal is more accurate than normal binary floating point
02956  * arithmetic, it requires some special values.
02957  *
02958  * == Infinity
02959  *
02960  * BigDecimal sometimes needs to return infinity, for example if you divide
02961  * a value by zero.
02962  *
02963  *      BigDecimal.new("1.0") / BigDecimal.new("0.0")  #=> infinity
02964  *      BigDecimal.new("-1.0") / BigDecimal.new("0.0")  #=> -infinity
02965  *
02966  * You can represent infinite numbers to BigDecimal using the strings
02967  * <code>'Infinity'</code>, <code>'+Infinity'</code> and
02968  * <code>'-Infinity'</code> (case-sensitive)
02969  *
02970  * == Not a Number
02971  *
02972  * When a computation results in an undefined value, the special value +NaN+
02973  * (for 'not a number') is returned.
02974  *
02975  * Example:
02976  *
02977  *      BigDecimal.new("0.0") / BigDecimal.new("0.0") #=> NaN
02978  *
02979  * You can also create undefined values.
02980  *
02981  * NaN is never considered to be the same as any other value, even NaN itself:
02982  *
02983  *      n = BigDecimal.new('NaN')
02984  *      n == 0.0 #=> nil
02985  *      n == n #=> nil
02986  *
02987  * == Positive and negative zero
02988  *
02989  * If a computation results in a value which is too small to be represented as
02990  * a BigDecimal within the currently specified limits of precision, zero must
02991  * be returned.
02992  *
02993  * If the value which is too small to be represented is negative, a BigDecimal
02994  * value of negative zero is returned.
02995  *
02996  *      BigDecimal.new("1.0") / BigDecimal.new("-Infinity") #=> -0.0
02997  *
02998  * If the value is positive, a value of positive zero is returned.
02999  *
03000  *      BigDecimal.new("1.0") / BigDecimal.new("Infinity") #=> 0.0
03001  *
03002  * (See BigDecimal.mode for how to specify limits of precision.)
03003  *
03004  * Note that +-0.0+ and +0.0+ are considered to be the same for the purposes of
03005  * comparison.
03006  *
03007  * Note also that in mathematics, there is no particular concept of negative
03008  * or positive zero; true mathematical zero has no sign.
03009  */
03010 void
03011 Init_bigdecimal(void)
03012 {
03013     VALUE arg;
03014 
03015     id_BigDecimal_exception_mode = rb_intern_const("BigDecimal.exception_mode");
03016     id_BigDecimal_rounding_mode = rb_intern_const("BigDecimal.rounding_mode");
03017     id_BigDecimal_precision_limit = rb_intern_const("BigDecimal.precision_limit");
03018 
03019     /* Initialize VP routines */
03020     VpInit(0UL);
03021 
03022     /* Class and method registration */
03023     rb_cBigDecimal = rb_define_class("BigDecimal",rb_cNumeric);
03024     rb_define_alloc_func(rb_cBigDecimal, BigDecimal_s_allocate);
03025 
03026     /* Global function */
03027     rb_define_global_function("BigDecimal", BigDecimal_global_new, -1);
03028 
03029     /* Class methods */
03030     rb_define_singleton_method(rb_cBigDecimal, "mode", BigDecimal_mode, -1);
03031     rb_define_singleton_method(rb_cBigDecimal, "limit", BigDecimal_limit, -1);
03032     rb_define_singleton_method(rb_cBigDecimal, "double_fig", BigDecimal_double_fig, 0);
03033     rb_define_singleton_method(rb_cBigDecimal, "_load", BigDecimal_load, 1);
03034     rb_define_singleton_method(rb_cBigDecimal, "ver", BigDecimal_version, 0);
03035 
03036     rb_define_singleton_method(rb_cBigDecimal, "save_exception_mode", BigDecimal_save_exception_mode, 0);
03037     rb_define_singleton_method(rb_cBigDecimal, "save_rounding_mode", BigDecimal_save_rounding_mode, 0);
03038     rb_define_singleton_method(rb_cBigDecimal, "save_limit", BigDecimal_save_limit, 0);
03039 
03040     /* Constants definition */
03041 
03042     /*
03043      * Base value used in internal calculations.  On a 32 bit system, BASE
03044      * is 10000, indicating that calculation is done in groups of 4 digits.
03045      * (If it were larger, BASE**2 wouldn't fit in 32 bits, so you couldn't
03046      * guarantee that two groups could always be multiplied together without
03047      * overflow.)
03048      */
03049     rb_define_const(rb_cBigDecimal, "BASE", INT2FIX((SIGNED_VALUE)VpBaseVal()));
03050 
03051     /* Exceptions */
03052 
03053     /*
03054      * 0xff: Determines whether overflow, underflow or zero divide result in
03055      * an exception being thrown. See BigDecimal.mode.
03056      */
03057     rb_define_const(rb_cBigDecimal, "EXCEPTION_ALL",INT2FIX(VP_EXCEPTION_ALL));
03058 
03059     /*
03060      * 0x02: Determines what happens when the result of a computation is not a
03061      * number (NaN). See BigDecimal.mode.
03062      */
03063     rb_define_const(rb_cBigDecimal, "EXCEPTION_NaN",INT2FIX(VP_EXCEPTION_NaN));
03064 
03065     /*
03066      * 0x01: Determines what happens when the result of a computation is
03067      * infinity.  See BigDecimal.mode.
03068      */
03069     rb_define_const(rb_cBigDecimal, "EXCEPTION_INFINITY",INT2FIX(VP_EXCEPTION_INFINITY));
03070 
03071     /*
03072      * 0x04: Determines what happens when the result of a computation is an
03073      * underflow (a result too small to be represented). See BigDecimal.mode.
03074      */
03075     rb_define_const(rb_cBigDecimal, "EXCEPTION_UNDERFLOW",INT2FIX(VP_EXCEPTION_UNDERFLOW));
03076 
03077     /*
03078      * 0x01: Determines what happens when the result of a computation is an
03079      * overflow (a result too large to be represented). See BigDecimal.mode.
03080      */
03081     rb_define_const(rb_cBigDecimal, "EXCEPTION_OVERFLOW",INT2FIX(VP_EXCEPTION_OVERFLOW));
03082 
03083     /*
03084      * 0x01: Determines what happens when a division by zero is performed.
03085      * See BigDecimal.mode.
03086      */
03087     rb_define_const(rb_cBigDecimal, "EXCEPTION_ZERODIVIDE",INT2FIX(VP_EXCEPTION_ZERODIVIDE));
03088 
03089     /*
03090      * 0x100: Determines what happens when a result must be rounded in order to
03091      * fit in the appropriate number of significant digits. See
03092      * BigDecimal.mode.
03093      */
03094     rb_define_const(rb_cBigDecimal, "ROUND_MODE",INT2FIX(VP_ROUND_MODE));
03095 
03096     /* 1: Indicates that values should be rounded away from zero. See
03097      * BigDecimal.mode.
03098      */
03099     rb_define_const(rb_cBigDecimal, "ROUND_UP",INT2FIX(VP_ROUND_UP));
03100 
03101     /* 2: Indicates that values should be rounded towards zero. See
03102      * BigDecimal.mode.
03103      */
03104     rb_define_const(rb_cBigDecimal, "ROUND_DOWN",INT2FIX(VP_ROUND_DOWN));
03105 
03106     /* 3: Indicates that digits >= 5 should be rounded up, others rounded down.
03107      * See BigDecimal.mode. */
03108     rb_define_const(rb_cBigDecimal, "ROUND_HALF_UP",INT2FIX(VP_ROUND_HALF_UP));
03109 
03110     /* 4: Indicates that digits >= 6 should be rounded up, others rounded down.
03111      * See BigDecimal.mode.
03112      */
03113     rb_define_const(rb_cBigDecimal, "ROUND_HALF_DOWN",INT2FIX(VP_ROUND_HALF_DOWN));
03114     /* 5: Round towards +infinity. See BigDecimal.mode. */
03115     rb_define_const(rb_cBigDecimal, "ROUND_CEILING",INT2FIX(VP_ROUND_CEIL));
03116 
03117     /* 6: Round towards -infinity. See BigDecimal.mode. */
03118     rb_define_const(rb_cBigDecimal, "ROUND_FLOOR",INT2FIX(VP_ROUND_FLOOR));
03119 
03120     /* 7: Round towards the even neighbor. See BigDecimal.mode. */
03121     rb_define_const(rb_cBigDecimal, "ROUND_HALF_EVEN",INT2FIX(VP_ROUND_HALF_EVEN));
03122 
03123     /* 0: Indicates that a value is not a number. See BigDecimal.sign. */
03124     rb_define_const(rb_cBigDecimal, "SIGN_NaN",INT2FIX(VP_SIGN_NaN));
03125 
03126     /* 1: Indicates that a value is +0. See BigDecimal.sign. */
03127     rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_ZERO",INT2FIX(VP_SIGN_POSITIVE_ZERO));
03128 
03129     /* -1: Indicates that a value is -0. See BigDecimal.sign. */
03130     rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_ZERO",INT2FIX(VP_SIGN_NEGATIVE_ZERO));
03131 
03132     /* 2: Indicates that a value is positive and finite. See BigDecimal.sign. */
03133     rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_FINITE",INT2FIX(VP_SIGN_POSITIVE_FINITE));
03134 
03135     /* -2: Indicates that a value is negative and finite. See BigDecimal.sign. */
03136     rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_FINITE",INT2FIX(VP_SIGN_NEGATIVE_FINITE));
03137 
03138     /* 3: Indicates that a value is positive and infinite. See BigDecimal.sign. */
03139     rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_INFINITE",INT2FIX(VP_SIGN_POSITIVE_INFINITE));
03140 
03141     /* -3: Indicates that a value is negative and infinite. See BigDecimal.sign. */
03142     rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_INFINITE",INT2FIX(VP_SIGN_NEGATIVE_INFINITE));
03143 
03144     arg = rb_str_new2("+Infinity");
03145     /* Positive infinity value. */
03146     rb_define_const(rb_cBigDecimal, "INFINITY", BigDecimal_global_new(1, &arg, rb_cBigDecimal));
03147     arg = rb_str_new2("NaN");
03148     /* 'Not a Number' value. */
03149     rb_define_const(rb_cBigDecimal, "NAN", BigDecimal_global_new(1, &arg, rb_cBigDecimal));
03150 
03151 
03152     /* instance methods */
03153     rb_define_method(rb_cBigDecimal, "initialize", BigDecimal_initialize, -1);
03154     rb_define_method(rb_cBigDecimal, "initialize_copy", BigDecimal_initialize_copy, 1);
03155     rb_define_method(rb_cBigDecimal, "precs", BigDecimal_prec, 0);
03156 
03157     rb_define_method(rb_cBigDecimal, "add", BigDecimal_add2, 2);
03158     rb_define_method(rb_cBigDecimal, "sub", BigDecimal_sub2, 2);
03159     rb_define_method(rb_cBigDecimal, "mult", BigDecimal_mult2, 2);
03160     rb_define_method(rb_cBigDecimal, "div", BigDecimal_div2, -1);
03161     rb_define_method(rb_cBigDecimal, "hash", BigDecimal_hash, 0);
03162     rb_define_method(rb_cBigDecimal, "to_s", BigDecimal_to_s, -1);
03163     rb_define_method(rb_cBigDecimal, "to_i", BigDecimal_to_i, 0);
03164     rb_define_method(rb_cBigDecimal, "to_int", BigDecimal_to_i, 0);
03165     rb_define_method(rb_cBigDecimal, "to_r", BigDecimal_to_r, 0);
03166     rb_define_method(rb_cBigDecimal, "split", BigDecimal_split, 0);
03167     rb_define_method(rb_cBigDecimal, "+", BigDecimal_add, 1);
03168     rb_define_method(rb_cBigDecimal, "-", BigDecimal_sub, 1);
03169     rb_define_method(rb_cBigDecimal, "+@", BigDecimal_uplus, 0);
03170     rb_define_method(rb_cBigDecimal, "-@", BigDecimal_neg, 0);
03171     rb_define_method(rb_cBigDecimal, "*", BigDecimal_mult, 1);
03172     rb_define_method(rb_cBigDecimal, "/", BigDecimal_div, 1);
03173     rb_define_method(rb_cBigDecimal, "quo", BigDecimal_div, 1);
03174     rb_define_method(rb_cBigDecimal, "%", BigDecimal_mod, 1);
03175     rb_define_method(rb_cBigDecimal, "modulo", BigDecimal_mod, 1);
03176     rb_define_method(rb_cBigDecimal, "remainder", BigDecimal_remainder, 1);
03177     rb_define_method(rb_cBigDecimal, "divmod", BigDecimal_divmod, 1);
03178     /* rb_define_method(rb_cBigDecimal, "dup", BigDecimal_dup, 0); */
03179     rb_define_method(rb_cBigDecimal, "to_f", BigDecimal_to_f, 0);
03180     rb_define_method(rb_cBigDecimal, "abs", BigDecimal_abs, 0);
03181     rb_define_method(rb_cBigDecimal, "sqrt", BigDecimal_sqrt, 1);
03182     rb_define_method(rb_cBigDecimal, "fix", BigDecimal_fix, 0);
03183     rb_define_method(rb_cBigDecimal, "round", BigDecimal_round, -1);
03184     rb_define_method(rb_cBigDecimal, "frac", BigDecimal_frac, 0);
03185     rb_define_method(rb_cBigDecimal, "floor", BigDecimal_floor, -1);
03186     rb_define_method(rb_cBigDecimal, "ceil", BigDecimal_ceil, -1);
03187     rb_define_method(rb_cBigDecimal, "power", BigDecimal_power, -1);
03188     rb_define_method(rb_cBigDecimal, "**", BigDecimal_power_op, 1);
03189     rb_define_method(rb_cBigDecimal, "<=>", BigDecimal_comp, 1);
03190     rb_define_method(rb_cBigDecimal, "==", BigDecimal_eq, 1);
03191     rb_define_method(rb_cBigDecimal, "===", BigDecimal_eq, 1);
03192     rb_define_method(rb_cBigDecimal, "eql?", BigDecimal_eq, 1);
03193     rb_define_method(rb_cBigDecimal, "<", BigDecimal_lt, 1);
03194     rb_define_method(rb_cBigDecimal, "<=", BigDecimal_le, 1);
03195     rb_define_method(rb_cBigDecimal, ">", BigDecimal_gt, 1);
03196     rb_define_method(rb_cBigDecimal, ">=", BigDecimal_ge, 1);
03197     rb_define_method(rb_cBigDecimal, "zero?", BigDecimal_zero, 0);
03198     rb_define_method(rb_cBigDecimal, "nonzero?", BigDecimal_nonzero, 0);
03199     rb_define_method(rb_cBigDecimal, "coerce", BigDecimal_coerce, 1);
03200     rb_define_method(rb_cBigDecimal, "inspect", BigDecimal_inspect, 0);
03201     rb_define_method(rb_cBigDecimal, "exponent", BigDecimal_exponent, 0);
03202     rb_define_method(rb_cBigDecimal, "sign", BigDecimal_sign, 0);
03203     rb_define_method(rb_cBigDecimal, "nan?",      BigDecimal_IsNaN, 0);
03204     rb_define_method(rb_cBigDecimal, "infinite?", BigDecimal_IsInfinite, 0);
03205     rb_define_method(rb_cBigDecimal, "finite?",   BigDecimal_IsFinite, 0);
03206     rb_define_method(rb_cBigDecimal, "truncate",  BigDecimal_truncate, -1);
03207     rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1);
03208 
03209     /* mathematical functions */
03210     rb_mBigMath = rb_define_module("BigMath");
03211     rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2);
03212     rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2);
03213 
03214     id_up = rb_intern_const("up");
03215     id_down = rb_intern_const("down");
03216     id_truncate = rb_intern_const("truncate");
03217     id_half_up = rb_intern_const("half_up");
03218     id_default = rb_intern_const("default");
03219     id_half_down = rb_intern_const("half_down");
03220     id_half_even = rb_intern_const("half_even");
03221     id_banker = rb_intern_const("banker");
03222     id_ceiling = rb_intern_const("ceiling");
03223     id_ceil = rb_intern_const("ceil");
03224     id_floor = rb_intern_const("floor");
03225     id_to_r = rb_intern_const("to_r");
03226     id_eq = rb_intern_const("==");
03227 }
03228 
03229 /*
03230  *
03231  *  ============================================================================
03232  *
03233  *  vp_ routines begin from here.
03234  *
03235  *  ============================================================================
03236  *
03237  */
03238 #ifdef BIGDECIMAL_DEBUG
03239 static int gfDebug = 1;         /* Debug switch */
03240 #if 0
03241 static int gfCheckVal = 1;      /* Value checking flag in VpNmlz()  */
03242 #endif
03243 #endif /* BIGDECIMAL_DEBUG */
03244 
03245 static Real *VpConstOne;    /* constant 1.0 */
03246 static Real *VpPt5;        /* constant 0.5 */
03247 #define maxnr 100UL    /* Maximum iterations for calcurating sqrt. */
03248                 /* used in VpSqrt() */
03249 
03250 /* ETC */
03251 #define MemCmp(x,y,z) memcmp(x,y,z)
03252 #define StrCmp(x,y)   strcmp(x,y)
03253 
03254 static int VpIsDefOP(Real *c,Real *a,Real *b,int sw);
03255 static int AddExponent(Real *a, SIGNED_VALUE n);
03256 static BDIGIT VpAddAbs(Real *a,Real *b,Real *c);
03257 static BDIGIT VpSubAbs(Real *a,Real *b,Real *c);
03258 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);
03259 static int VpNmlz(Real *a);
03260 static void VpFormatSt(char *psz, size_t fFmt);
03261 static int VpRdup(Real *m, size_t ind_m);
03262 
03263 #ifdef BIGDECIMAL_DEBUG
03264 static int gnAlloc=0; /* Memory allocation counter */
03265 #endif /* BIGDECIMAL_DEBUG */
03266 
03267 VP_EXPORT void *
03268 VpMemAlloc(size_t mb)
03269 {
03270     void *p = xmalloc(mb);
03271     if (!p) {
03272         VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
03273     }
03274     memset(p, 0, mb);
03275 #ifdef BIGDECIMAL_DEBUG
03276     gnAlloc++; /* Count allocation call */
03277 #endif /* BIGDECIMAL_DEBUG */
03278     return p;
03279 }
03280 
03281 VP_EXPORT void *
03282 VpMemRealloc(void *ptr, size_t mb)
03283 {
03284     void *p = xrealloc(ptr, mb);
03285     if (!p) {
03286         VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1);
03287     }
03288     return p;
03289 }
03290 
03291 VP_EXPORT void
03292 VpFree(Real *pv)
03293 {
03294     if(pv != NULL) {
03295         xfree(pv);
03296 #ifdef BIGDECIMAL_DEBUG
03297         gnAlloc--; /* Decrement allocation count */
03298         if(gnAlloc==0) {
03299             printf(" *************** All memories allocated freed ****************");
03300             getchar();
03301         }
03302         if(gnAlloc<0) {
03303             printf(" ??????????? Too many memory free calls(%d) ?????????????\n",gnAlloc);
03304             getchar();
03305         }
03306 #endif /* BIGDECIMAL_DEBUG */
03307     }
03308 }
03309 
03310 /*
03311  * EXCEPTION Handling.
03312  */
03313 
03314 #define rmpd_set_thread_local_exception_mode(mode) \
03315     rb_thread_local_aset( \
03316         rb_thread_current(), \
03317         id_BigDecimal_exception_mode, \
03318         INT2FIX((int)(mode)) \
03319     )
03320 
03321 static unsigned short
03322 VpGetException (void)
03323 {
03324     VALUE const vmode = rb_thread_local_aref(
03325         rb_thread_current(),
03326         id_BigDecimal_exception_mode
03327     );
03328 
03329     if (NIL_P(vmode)) {
03330         rmpd_set_thread_local_exception_mode(RMPD_EXCEPTION_MODE_DEFAULT);
03331         return RMPD_EXCEPTION_MODE_DEFAULT;
03332     }
03333 
03334     return (unsigned short)FIX2UINT(vmode);
03335 }
03336 
03337 static void
03338 VpSetException(unsigned short f)
03339 {
03340     rmpd_set_thread_local_exception_mode(f);
03341 }
03342 
03343 /*
03344  * Precision limit.
03345  */
03346 
03347 #define rmpd_set_thread_local_precision_limit(limit) \
03348     rb_thread_local_aset( \
03349         rb_thread_current(), \
03350         id_BigDecimal_precision_limit, \
03351         SIZET2NUM(limit) \
03352     )
03353 #define RMPD_PRECISION_LIMIT_DEFAULT ((size_t)0)
03354 
03355 /* These 2 functions added at v1.1.7 */
03356 VP_EXPORT size_t
03357 VpGetPrecLimit(void)
03358 {
03359     VALUE const vlimit = rb_thread_local_aref(
03360         rb_thread_current(),
03361         id_BigDecimal_precision_limit
03362     );
03363 
03364     if (NIL_P(vlimit)) {
03365         rmpd_set_thread_local_precision_limit(RMPD_PRECISION_LIMIT_DEFAULT);
03366         return RMPD_PRECISION_LIMIT_DEFAULT;
03367     }
03368 
03369     return NUM2SIZET(vlimit);
03370 }
03371 
03372 VP_EXPORT size_t
03373 VpSetPrecLimit(size_t n)
03374 {
03375     size_t const s = VpGetPrecLimit();
03376     rmpd_set_thread_local_precision_limit(n);
03377     return s;
03378 }
03379 
03380 /*
03381  * Rounding mode.
03382  */
03383 
03384 #define rmpd_set_thread_local_rounding_mode(mode) \
03385     rb_thread_local_aset( \
03386         rb_thread_current(), \
03387         id_BigDecimal_rounding_mode, \
03388         INT2FIX((int)(mode)) \
03389     )
03390 
03391 VP_EXPORT unsigned short
03392 VpGetRoundMode(void)
03393 {
03394     VALUE const vmode = rb_thread_local_aref(
03395         rb_thread_current(),
03396         id_BigDecimal_rounding_mode
03397     );
03398 
03399     if (NIL_P(vmode)) {
03400         rmpd_set_thread_local_rounding_mode(RMPD_ROUNDING_MODE_DEFAULT);
03401         return RMPD_ROUNDING_MODE_DEFAULT;
03402     }
03403 
03404     return (unsigned short)FIX2INT(vmode);
03405 }
03406 
03407 VP_EXPORT int
03408 VpIsRoundMode(unsigned short n)
03409 {
03410     switch (n) {
03411       case VP_ROUND_UP:
03412       case VP_ROUND_DOWN:
03413       case VP_ROUND_HALF_UP:
03414       case VP_ROUND_HALF_DOWN:
03415       case VP_ROUND_CEIL:
03416       case VP_ROUND_FLOOR:
03417       case VP_ROUND_HALF_EVEN:
03418         return 1;
03419 
03420       default:
03421         return 0;
03422     }
03423 }
03424 
03425 VP_EXPORT unsigned short
03426 VpSetRoundMode(unsigned short n)
03427 {
03428     if (VpIsRoundMode(n)) {
03429         rmpd_set_thread_local_rounding_mode(n);
03430         return n;
03431     }
03432 
03433     return VpGetRoundMode();
03434 }
03435 
03436 /*
03437  *  0.0 & 1.0 generator
03438  *    These gZero_..... and gOne_..... can be any name
03439  *    referenced from nowhere except Zero() and One().
03440  *    gZero_..... and gOne_..... must have global scope
03441  *    (to let the compiler know they may be changed in outside
03442  *    (... but not actually..)).
03443  */
03444 volatile const double gZero_ABCED9B1_CE73__00400511F31D = 0.0;
03445 volatile const double gOne_ABCED9B4_CE73__00400511F31D  = 1.0;
03446 static double
03447 Zero(void)
03448 {
03449     return gZero_ABCED9B1_CE73__00400511F31D;
03450 }
03451 
03452 static double
03453 One(void)
03454 {
03455     return gOne_ABCED9B4_CE73__00400511F31D;
03456 }
03457 
03458 /*
03459   ----------------------------------------------------------------
03460   Value of sign in Real structure is reserved for future use.
03461   short sign;
03462                     ==0 : NaN
03463                       1 : Positive zero
03464                      -1 : Negative zero
03465                       2 : Positive number
03466                      -2 : Negative number
03467                       3 : Positive infinite number
03468                      -3 : Negative infinite number
03469   ----------------------------------------------------------------
03470 */
03471 
03472 VP_EXPORT double
03473 VpGetDoubleNaN(void) /* Returns the value of NaN */
03474 {
03475     static double fNaN = 0.0;
03476     if(fNaN==0.0) fNaN = Zero()/Zero();
03477     return fNaN;
03478 }
03479 
03480 VP_EXPORT double
03481 VpGetDoublePosInf(void) /* Returns the value of +Infinity */
03482 {
03483     static double fInf = 0.0;
03484     if(fInf==0.0) fInf = One()/Zero();
03485     return fInf;
03486 }
03487 
03488 VP_EXPORT double
03489 VpGetDoubleNegInf(void) /* Returns the value of -Infinity */
03490 {
03491     static double fInf = 0.0;
03492     if(fInf==0.0) fInf = -(One()/Zero());
03493     return fInf;
03494 }
03495 
03496 VP_EXPORT double
03497 VpGetDoubleNegZero(void) /* Returns the value of -0 */
03498 {
03499     static double nzero = 1000.0;
03500     if(nzero!=0.0) nzero = (One()/VpGetDoubleNegInf());
03501     return nzero;
03502 }
03503 
03504 #if 0  /* unused */
03505 VP_EXPORT int
03506 VpIsNegDoubleZero(double v)
03507 {
03508     double z = VpGetDoubleNegZero();
03509     return MemCmp(&v,&z,sizeof(v))==0;
03510 }
03511 #endif
03512 
03513 VP_EXPORT int
03514 VpException(unsigned short f, const char *str,int always)
03515 {
03516     VALUE exc;
03517     int   fatal=0;
03518     unsigned short const exception_mode = VpGetException();
03519 
03520     if(f==VP_EXCEPTION_OP || f==VP_EXCEPTION_MEMORY) always = 1;
03521 
03522     if (always || (exception_mode & f)) {
03523         switch(f)
03524         {
03525         /*
03526         case VP_EXCEPTION_OVERFLOW:
03527         */
03528         case VP_EXCEPTION_ZERODIVIDE:
03529         case VP_EXCEPTION_INFINITY:
03530         case VP_EXCEPTION_NaN:
03531         case VP_EXCEPTION_UNDERFLOW:
03532         case VP_EXCEPTION_OP:
03533              exc = rb_eFloatDomainError;
03534              goto raise;
03535         case VP_EXCEPTION_MEMORY:
03536              fatal = 1;
03537              goto raise;
03538         default:
03539              fatal = 1;
03540              goto raise;
03541         }
03542     }
03543     return 0; /* 0 Means VpException() raised no exception */
03544 
03545 raise:
03546     if(fatal) rb_fatal("%s", str);
03547     else   rb_raise(exc, "%s", str);
03548     return 0;
03549 }
03550 
03551 /* Throw exception or returns 0,when resulting c is Inf or NaN */
03552 /*  sw=1:+ 2:- 3:* 4:/ */
03553 static int
03554 VpIsDefOP(Real *c,Real *a,Real *b,int sw)
03555 {
03556     if(VpIsNaN(a) || VpIsNaN(b)) {
03557         /* at least a or b is NaN */
03558         VpSetNaN(c);
03559         goto NaN;
03560     }
03561 
03562     if(VpIsInf(a)) {
03563         if(VpIsInf(b)) {
03564             switch(sw)
03565             {
03566             case 1: /* + */
03567                 if(VpGetSign(a)==VpGetSign(b)) {
03568                     VpSetInf(c,VpGetSign(a));
03569                     goto Inf;
03570                 } else {
03571                     VpSetNaN(c);
03572                     goto NaN;
03573                 }
03574             case 2: /* - */
03575                 if(VpGetSign(a)!=VpGetSign(b)) {
03576                     VpSetInf(c,VpGetSign(a));
03577                     goto Inf;
03578                 } else {
03579                     VpSetNaN(c);
03580                     goto NaN;
03581                 }
03582                 break;
03583             case 3: /* * */
03584                 VpSetInf(c,VpGetSign(a)*VpGetSign(b));
03585                 goto Inf;
03586                 break;
03587             case 4: /* / */
03588                 VpSetNaN(c);
03589                 goto NaN;
03590             }
03591             VpSetNaN(c);
03592             goto NaN;
03593         }
03594         /* Inf op Finite */
03595         switch(sw)
03596         {
03597         case 1: /* + */
03598         case 2: /* - */
03599                 VpSetInf(c,VpGetSign(a));
03600                 break;
03601         case 3: /* * */
03602                 if(VpIsZero(b)) {
03603                     VpSetNaN(c);
03604                     goto NaN;
03605                 }
03606                 VpSetInf(c,VpGetSign(a)*VpGetSign(b));
03607                 break;
03608         case 4: /* / */
03609                 VpSetInf(c,VpGetSign(a)*VpGetSign(b));
03610         }
03611         goto Inf;
03612     }
03613 
03614     if(VpIsInf(b)) {
03615         switch(sw)
03616         {
03617         case 1: /* + */
03618                 VpSetInf(c,VpGetSign(b));
03619                 break;
03620         case 2: /* - */
03621                 VpSetInf(c,-VpGetSign(b));
03622                 break;
03623         case 3: /* * */
03624                 if(VpIsZero(a)) {
03625                     VpSetNaN(c);
03626                     goto NaN;
03627                 }
03628                 VpSetInf(c,VpGetSign(a)*VpGetSign(b));
03629                 break;
03630         case 4: /* / */
03631                 VpSetZero(c,VpGetSign(a)*VpGetSign(b));
03632         }
03633         goto Inf;
03634     }
03635     return 1; /* Results OK */
03636 
03637 Inf:
03638     return VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0);
03639 NaN:
03640     return VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'",0);
03641 }
03642 
03643 /*
03644   ----------------------------------------------------------------
03645 */
03646 
03647 /*
03648  *    returns number of chars needed to represent vp in specified format.
03649  */
03650 VP_EXPORT size_t
03651 VpNumOfChars(Real *vp,const char *pszFmt)
03652 {
03653     SIGNED_VALUE  ex;
03654     size_t nc;
03655 
03656     if(vp == NULL)   return BASE_FIG*2+6;
03657     if(!VpIsDef(vp)) return 32; /* not sure,may be OK */
03658 
03659     switch(*pszFmt)
03660     {
03661     case 'F':
03662          nc = BASE_FIG*(vp->Prec + 1)+2;
03663          ex = vp->exponent;
03664          if(ex < 0) {
03665              nc += BASE_FIG*(size_t)(-ex);
03666          }
03667          else {
03668              if((size_t)ex > vp->Prec) {
03669                  nc += BASE_FIG*((size_t)ex - vp->Prec);
03670              }
03671          }
03672          break;
03673     case 'E':
03674     default:
03675          nc = BASE_FIG*(vp->Prec + 2)+6; /* 3: sign + exponent chars */
03676     }
03677     return nc;
03678 }
03679 
03680 /*
03681  * Initializer for Vp routines and constants used.
03682  * [Input]
03683  *   BaseVal: Base value(assigned to BASE) for Vp calculation.
03684  *   It must be the form BaseVal=10**n.(n=1,2,3,...)
03685  *   If Base <= 0L,then the BASE will be calcurated so
03686  *   that BASE is as large as possible satisfying the
03687  *   relation MaxVal <= BASE*(BASE+1). Where the value
03688  *   MaxVal is the largest value which can be represented
03689  *   by one BDIGIT word in the computer used.
03690  *
03691  * [Returns]
03692  * 1+DBL_DIG   ... OK
03693  */
03694 VP_EXPORT size_t
03695 VpInit(BDIGIT BaseVal)
03696 {
03697     /* Setup +/- Inf  NaN -0 */
03698     VpGetDoubleNaN();
03699     VpGetDoublePosInf();
03700     VpGetDoubleNegInf();
03701     VpGetDoubleNegZero();
03702 
03703     /* Allocates Vp constants. */
03704     VpConstOne = VpAlloc(1UL, "1");
03705     VpPt5 = VpAlloc(1UL, ".5");
03706 
03707 #ifdef BIGDECIMAL_DEBUG
03708     gnAlloc = 0;
03709 #endif /* BIGDECIMAL_DEBUG */
03710 
03711 #ifdef BIGDECIMAL_DEBUG
03712     if(gfDebug) {
03713         printf("VpInit: BaseVal   = %lu\n", BaseVal);
03714         printf("  BASE   = %lu\n", BASE);
03715         printf("  HALF_BASE = %lu\n", HALF_BASE);
03716         printf("  BASE1  = %lu\n", BASE1);
03717         printf("  BASE_FIG  = %u\n", BASE_FIG);
03718         printf("  DBLE_FIG  = %d\n", DBLE_FIG);
03719     }
03720 #endif /* BIGDECIMAL_DEBUG */
03721 
03722     return rmpd_double_figures();
03723 }
03724 
03725 VP_EXPORT Real *
03726 VpOne(void)
03727 {
03728     return VpConstOne;
03729 }
03730 
03731 /* If exponent overflows,then raise exception or returns 0 */
03732 static int
03733 AddExponent(Real *a, SIGNED_VALUE n)
03734 {
03735     SIGNED_VALUE e = a->exponent;
03736     SIGNED_VALUE m = e+n;
03737     SIGNED_VALUE eb, mb;
03738     if(e>0) {
03739         if(n>0) {
03740             mb = m*(SIGNED_VALUE)BASE_FIG;
03741             eb = e*(SIGNED_VALUE)BASE_FIG;
03742             if(mb<eb) goto overflow;
03743         }
03744     } else if(n<0) {
03745         mb = m*(SIGNED_VALUE)BASE_FIG;
03746         eb = e*(SIGNED_VALUE)BASE_FIG;
03747         if(mb>eb) goto underflow;
03748     }
03749     a->exponent = m;
03750     return 1;
03751 
03752 /* Overflow/Underflow ==> Raise exception or returns 0 */
03753 underflow:
03754     VpSetZero(a,VpGetSign(a));
03755     return VpException(VP_EXCEPTION_UNDERFLOW,"Exponent underflow",0);
03756 
03757 overflow:
03758     VpSetInf(a,VpGetSign(a));
03759     return VpException(VP_EXCEPTION_OVERFLOW,"Exponent overflow",0);
03760 }
03761 
03762 /*
03763  * Allocates variable.
03764  * [Input]
03765  *   mx ... allocation unit, if zero then mx is determined by szVal.
03766  *    The mx is the number of effective digits can to be stored.
03767  *   szVal ... value assigned(char). If szVal==NULL,then zero is assumed.
03768  *            If szVal[0]=='#' then Max. Prec. will not be considered(1.1.7),
03769  *            full precision specified by szVal is allocated.
03770  *
03771  * [Returns]
03772  *   Pointer to the newly allocated variable, or
03773  *   NULL be returned if memory allocation is failed,or any error.
03774  */
03775 VP_EXPORT Real *
03776 VpAlloc(size_t mx, const char *szVal)
03777 {
03778     size_t i, ni, ipn, ipf, nf, ipe, ne, nalloc;
03779     char v,*psz;
03780     int  sign=1;
03781     Real *vp = NULL;
03782     size_t mf = VpGetPrecLimit();
03783     VALUE buf;
03784 
03785     mx = (mx + BASE_FIG - 1) / BASE_FIG + 1;    /* Determine allocation unit. */
03786     if (szVal) {
03787         while (ISSPACE(*szVal)) szVal++;
03788         if (*szVal != '#') {
03789              if (mf) {
03790                 mf = (mf + BASE_FIG - 1) / BASE_FIG + 2; /* Needs 1 more for div */
03791                 if (mx > mf) {
03792                     mx = mf;
03793                 }
03794             }
03795         }
03796         else {
03797             ++szVal;
03798         }
03799     }
03800     else {
03801        /* necessary to be able to store */
03802        /* at least mx digits. */
03803        /* szVal==NULL ==> allocate zero value. */
03804        vp = VpAllocReal(mx);
03805        /* xmalloc() alway returns(or throw interruption) */
03806        vp->MaxPrec = mx;    /* set max precision */
03807        VpSetZero(vp,1);    /* initialize vp to zero. */
03808        return vp;
03809     }
03810 
03811     /* Skip all '_' after digit: 2006-6-30 */
03812     ni = 0;
03813     buf = rb_str_tmp_new(strlen(szVal)+1);
03814     psz = RSTRING_PTR(buf);
03815     i   = 0;
03816     ipn = 0;
03817     while ((psz[i]=szVal[ipn]) != 0) {
03818         if (ISDIGIT(psz[i])) ++ni;
03819         if (psz[i] == '_') {
03820             if (ni > 0) { ipn++; continue; }
03821             psz[i] = 0;
03822             break;
03823         }
03824         ++i;
03825         ++ipn;
03826     }
03827     /* Skip trailing spaces */
03828     while (--i > 0) {
03829         if (ISSPACE(psz[i])) psz[i] = 0;
03830         else break;
03831     }
03832     szVal = psz;
03833 
03834     /* Check on Inf & NaN */
03835     if (StrCmp(szVal, SZ_PINF) == 0 ||
03836         StrCmp(szVal, SZ_INF)  == 0 ) {
03837         vp = VpAllocReal(1);
03838         vp->MaxPrec = 1;    /* set max precision */
03839         VpSetPosInf(vp);
03840         return vp;
03841     }
03842     if (StrCmp(szVal, SZ_NINF) == 0) {
03843         vp = VpAllocReal(1);
03844         vp->MaxPrec = 1;    /* set max precision */
03845         VpSetNegInf(vp);
03846         return vp;
03847     }
03848     if (StrCmp(szVal, SZ_NaN) == 0) {
03849         vp = VpAllocReal(1);
03850         vp->MaxPrec = 1;    /* set max precision */
03851         VpSetNaN(vp);
03852         return vp;
03853     }
03854 
03855     /* check on number szVal[] */
03856     ipn = i = 0;
03857     if      (szVal[i] == '-') { sign=-1; ++i; }
03858     else if (szVal[i] == '+')            ++i;
03859     /* Skip digits */
03860     ni = 0;            /* digits in mantissa */
03861     while ((v = szVal[i]) != 0) {
03862         if (!ISDIGIT(v)) break;
03863         ++i;
03864         ++ni;
03865     }
03866     nf  = 0;
03867     ipf = 0;
03868     ipe = 0;
03869     ne  = 0;
03870     if (v) {
03871         /* other than digit nor \0 */
03872         if (szVal[i] == '.') {    /* xxx. */
03873             ++i;
03874             ipf = i;
03875             while ((v = szVal[i]) != 0) {    /* get fraction part. */
03876                 if (!ISDIGIT(v)) break;
03877                 ++i;
03878                 ++nf;
03879             }
03880         }
03881         ipe = 0;        /* Exponent */
03882 
03883         switch (szVal[i]) {
03884         case '\0':
03885             break;
03886         case 'e': case 'E':
03887         case 'd': case 'D':
03888             ++i;
03889             ipe = i;
03890             v = szVal[i];
03891             if ((v == '-') || (v == '+')) ++i;
03892             while ((v=szVal[i]) != 0) {
03893                 if (!ISDIGIT(v)) break;
03894                 ++i;
03895                 ++ne;
03896             }
03897             break;
03898         default:
03899             break;
03900         }
03901     }
03902     nalloc = (ni + nf + BASE_FIG - 1) / BASE_FIG + 1;    /* set effective allocation  */
03903     /* units for szVal[]  */
03904     if (mx <= 0) mx = 1;
03905     nalloc = Max(nalloc, mx);
03906     mx = nalloc;
03907     vp = VpAllocReal(mx);
03908     /* xmalloc() alway returns(or throw interruption) */
03909     vp->MaxPrec = mx;        /* set max precision */
03910     VpSetZero(vp, sign);
03911     VpCtoV(vp, &szVal[ipn], ni, &szVal[ipf], nf, &szVal[ipe], ne);
03912     rb_str_resize(buf, 0);
03913     return vp;
03914 }
03915 
03916 /*
03917  * Assignment(c=a).
03918  * [Input]
03919  *   a   ... RHSV
03920  *   isw ... switch for assignment.
03921  *    c = a  when isw > 0
03922  *    c = -a when isw < 0
03923  *    if c->MaxPrec < a->Prec,then round operation
03924  *    will be performed.
03925  * [Output]
03926  *  c  ... LHSV
03927  */
03928 VP_EXPORT size_t
03929 VpAsgn(Real *c, Real *a, int isw)
03930 {
03931     size_t n;
03932     if(VpIsNaN(a)) {
03933         VpSetNaN(c);
03934         return 0;
03935     }
03936     if(VpIsInf(a)) {
03937         VpSetInf(c,isw*VpGetSign(a));
03938         return 0;
03939     }
03940 
03941     /* check if the RHS is zero */
03942     if(!VpIsZero(a)) {
03943         c->exponent = a->exponent;    /* store  exponent */
03944         VpSetSign(c,(isw*VpGetSign(a)));    /* set sign */
03945         n =(a->Prec < c->MaxPrec) ?(a->Prec) :(c->MaxPrec);
03946         c->Prec = n;
03947         memcpy(c->frac, a->frac, n * sizeof(BDIGIT));
03948         /* Needs round ? */
03949         if(isw!=10) {
03950             /* Not in ActiveRound */
03951             if(c->Prec < a->Prec) {
03952                 VpInternalRound(c,n,(n>0)?a->frac[n-1]:0,a->frac[n]);
03953             } else {
03954                 VpLimitRound(c,0);
03955             }
03956         }
03957     } else {
03958         /* The value of 'a' is zero.  */
03959         VpSetZero(c,isw*VpGetSign(a));
03960         return 1;
03961     }
03962     return c->Prec*BASE_FIG;
03963 }
03964 
03965 /*
03966  *   c = a + b  when operation =  1 or 2
03967  *  = a - b  when operation = -1 or -2.
03968  *   Returns number of significant digits of c
03969  */
03970 VP_EXPORT size_t
03971 VpAddSub(Real *c, Real *a, Real *b, int operation)
03972 {
03973     short sw, isw;
03974     Real *a_ptr, *b_ptr;
03975     size_t n, na, nb, i;
03976     BDIGIT mrv;
03977 
03978 #ifdef BIGDECIMAL_DEBUG
03979     if(gfDebug) {
03980         VPrint(stdout, "VpAddSub(enter) a=% \n", a);
03981         VPrint(stdout, "     b=% \n", b);
03982         printf(" operation=%d\n", operation);
03983     }
03984 #endif /* BIGDECIMAL_DEBUG */
03985 
03986     if(!VpIsDefOP(c,a,b,(operation>0)?1:2)) return 0; /* No significant digits */
03987 
03988     /* check if a or b is zero  */
03989     if(VpIsZero(a)) {
03990         /* a is zero,then assign b to c */
03991         if(!VpIsZero(b)) {
03992             VpAsgn(c, b, operation);
03993         } else {
03994             /* Both a and b are zero. */
03995             if(VpGetSign(a)<0 && operation*VpGetSign(b)<0) {
03996                 /* -0 -0 */
03997                 VpSetZero(c,-1);
03998             } else {
03999                 VpSetZero(c,1);
04000             }
04001             return 1; /* 0: 1 significant digits */
04002         }
04003         return c->Prec*BASE_FIG;
04004     }
04005     if(VpIsZero(b)) {
04006         /* b is zero,then assign a to c. */
04007         VpAsgn(c, a, 1);
04008         return c->Prec*BASE_FIG;
04009     }
04010 
04011     if(operation < 0) sw = -1;
04012     else              sw =  1;
04013 
04014     /* compare absolute value. As a result,|a_ptr|>=|b_ptr| */
04015     if(a->exponent > b->exponent) {
04016         a_ptr = a;
04017         b_ptr = b;
04018     }         /* |a|>|b| */
04019     else if(a->exponent < b->exponent) {
04020         a_ptr = b;
04021         b_ptr = a;
04022     }                /* |a|<|b| */
04023     else {
04024         /* Exponent part of a and b is the same,then compare fraction */
04025         /* part */
04026         na = a->Prec;
04027         nb = b->Prec;
04028         n = Min(na, nb);
04029         for(i=0;i < n; ++i) {
04030             if(a->frac[i] > b->frac[i]) {
04031                 a_ptr = a;
04032                 b_ptr = b;
04033                 goto end_if;
04034             } else if(a->frac[i] < b->frac[i]) {
04035                 a_ptr = b;
04036                 b_ptr = a;
04037                 goto end_if;
04038             }
04039         }
04040         if(na > nb) {
04041          a_ptr = a;
04042             b_ptr = b;
04043             goto end_if;
04044         } else if(na < nb) {
04045             a_ptr = b;
04046             b_ptr = a;
04047             goto end_if;
04048         }
04049         /* |a| == |b| */
04050         if(VpGetSign(a) + sw *VpGetSign(b) == 0) {
04051             VpSetZero(c,1);        /* abs(a)=abs(b) and operation = '-'  */
04052             return c->Prec*BASE_FIG;
04053         }
04054         a_ptr = a;
04055         b_ptr = b;
04056     }
04057 
04058 end_if:
04059     isw = VpGetSign(a) + sw *VpGetSign(b);
04060     /*
04061      *  isw = 0 ...( 1)+(-1),( 1)-( 1),(-1)+(1),(-1)-(-1)
04062      *      = 2 ...( 1)+( 1),( 1)-(-1)
04063      *      =-2 ...(-1)+(-1),(-1)-( 1)
04064      *   If isw==0, then c =(Sign a_ptr)(|a_ptr|-|b_ptr|)
04065      *              else c =(Sign ofisw)(|a_ptr|+|b_ptr|)
04066     */
04067     if(isw) {            /* addition */
04068         VpSetSign(c, 1);
04069         mrv = VpAddAbs(a_ptr, b_ptr, c);
04070         VpSetSign(c, isw / 2);
04071     } else {            /* subtraction */
04072         VpSetSign(c, 1);
04073         mrv = VpSubAbs(a_ptr, b_ptr, c);
04074         if(a_ptr == a) {
04075             VpSetSign(c,VpGetSign(a));
04076         } else    {
04077             VpSetSign(c,VpGetSign(a_ptr) * sw);
04078         }
04079     }
04080     VpInternalRound(c,0,(c->Prec>0)?c->frac[c->Prec-1]:0,mrv);
04081 
04082 #ifdef BIGDECIMAL_DEBUG
04083     if(gfDebug) {
04084         VPrint(stdout, "VpAddSub(result) c=% \n", c);
04085         VPrint(stdout, "     a=% \n", a);
04086         VPrint(stdout, "     b=% \n", b);
04087         printf(" operation=%d\n", operation);
04088     }
04089 #endif /* BIGDECIMAL_DEBUG */
04090     return c->Prec*BASE_FIG;
04091 }
04092 
04093 /*
04094  * Addition of two variable precisional variables
04095  * a and b assuming abs(a)>abs(b).
04096  *   c = abs(a) + abs(b) ; where |a|>=|b|
04097  */
04098 static BDIGIT
04099 VpAddAbs(Real *a, Real *b, Real *c)
04100 {
04101     size_t word_shift;
04102     size_t ap;
04103     size_t bp;
04104     size_t cp;
04105     size_t a_pos;
04106     size_t b_pos, b_pos_with_word_shift;
04107     size_t c_pos;
04108     BDIGIT av, bv, carry, mrv;
04109 
04110 #ifdef BIGDECIMAL_DEBUG
04111     if(gfDebug) {
04112         VPrint(stdout, "VpAddAbs called: a = %\n", a);
04113         VPrint(stdout, "     b = %\n", b);
04114     }
04115 #endif /* BIGDECIMAL_DEBUG */
04116 
04117     word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
04118     a_pos = ap;
04119     b_pos = bp;
04120     c_pos = cp;
04121     if(word_shift==(size_t)-1L) return 0; /* Overflow */
04122     if(b_pos == (size_t)-1L) goto Assign_a;
04123 
04124     mrv = av + bv; /* Most right val. Used for round. */
04125 
04126     /* Just assign the last few digits of b to c because a has no  */
04127     /* corresponding digits to be added. */
04128     while(b_pos + word_shift > a_pos) {
04129         --c_pos;
04130         if(b_pos > 0) {
04131             c->frac[c_pos] = b->frac[--b_pos];
04132         } else {
04133             --word_shift;
04134             c->frac[c_pos] = 0;
04135         }
04136     }
04137 
04138     /* Just assign the last few digits of a to c because b has no */
04139     /* corresponding digits to be added. */
04140     b_pos_with_word_shift = b_pos + word_shift;
04141     while(a_pos > b_pos_with_word_shift) {
04142         c->frac[--c_pos] = a->frac[--a_pos];
04143     }
04144     carry = 0;    /* set first carry be zero */
04145 
04146     /* Now perform addition until every digits of b will be */
04147     /* exhausted. */
04148     while(b_pos > 0) {
04149         c->frac[--c_pos] = a->frac[--a_pos] + b->frac[--b_pos] + carry;
04150         if(c->frac[c_pos] >= BASE) {
04151             c->frac[c_pos] -= BASE;
04152             carry = 1;
04153         } else {
04154             carry = 0;
04155         }
04156     }
04157 
04158     /* Just assign the first few digits of a with considering */
04159     /* the carry obtained so far because b has been exhausted. */
04160     while(a_pos > 0) {
04161         c->frac[--c_pos] = a->frac[--a_pos] + carry;
04162         if(c->frac[c_pos] >= BASE) {
04163             c->frac[c_pos] -= BASE;
04164             carry = 1;
04165         } else {
04166             carry = 0;
04167         }
04168     }
04169     if(c_pos) c->frac[c_pos - 1] += carry;
04170     goto Exit;
04171 
04172 Assign_a:
04173     VpAsgn(c, a, 1);
04174     mrv = 0;
04175 
04176 Exit:
04177 
04178 #ifdef BIGDECIMAL_DEBUG
04179     if(gfDebug) {
04180         VPrint(stdout, "VpAddAbs exit: c=% \n", c);
04181     }
04182 #endif /* BIGDECIMAL_DEBUG */
04183     return mrv;
04184 }
04185 
04186 /*
04187  * c = abs(a) - abs(b)
04188  */
04189 static BDIGIT
04190 VpSubAbs(Real *a, Real *b, Real *c)
04191 {
04192     size_t word_shift;
04193     size_t ap;
04194     size_t bp;
04195     size_t cp;
04196     size_t a_pos;
04197     size_t b_pos, b_pos_with_word_shift;
04198     size_t c_pos;
04199     BDIGIT av, bv, borrow, mrv;
04200 
04201 #ifdef BIGDECIMAL_DEBUG
04202     if(gfDebug) {
04203         VPrint(stdout, "VpSubAbs called: a = %\n", a);
04204         VPrint(stdout, "     b = %\n", b);
04205     }
04206 #endif /* BIGDECIMAL_DEBUG */
04207 
04208     word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
04209     a_pos = ap;
04210     b_pos = bp;
04211     c_pos = cp;
04212     if(word_shift==(size_t)-1L) return 0; /* Overflow */
04213     if(b_pos == (size_t)-1L) goto Assign_a;
04214 
04215     if(av >= bv) {
04216         mrv = av - bv;
04217         borrow = 0;
04218     } else {
04219         mrv    = 0;
04220         borrow = 1;
04221     }
04222 
04223     /* Just assign the values which are the BASE subtracted by   */
04224     /* each of the last few digits of the b because the a has no */
04225     /* corresponding digits to be subtracted. */
04226     if(b_pos + word_shift > a_pos) {
04227         while(b_pos + word_shift > a_pos) {
04228             --c_pos;
04229             if(b_pos > 0) {
04230                 c->frac[c_pos] = BASE - b->frac[--b_pos] - borrow;
04231             } else {
04232                 --word_shift;
04233                 c->frac[c_pos] = BASE - borrow;
04234             }
04235             borrow = 1;
04236         }
04237     }
04238     /* Just assign the last few digits of a to c because b has no */
04239     /* corresponding digits to subtract. */
04240 
04241     b_pos_with_word_shift = b_pos + word_shift;
04242     while(a_pos > b_pos_with_word_shift) {
04243         c->frac[--c_pos] = a->frac[--a_pos];
04244     }
04245 
04246     /* Now perform subtraction until every digits of b will be */
04247     /* exhausted. */
04248     while(b_pos > 0) {
04249         --c_pos;
04250         if(a->frac[--a_pos] < b->frac[--b_pos] + borrow) {
04251             c->frac[c_pos] = BASE + a->frac[a_pos] - b->frac[b_pos] - borrow;
04252             borrow = 1;
04253         } else {
04254             c->frac[c_pos] = a->frac[a_pos] - b->frac[b_pos] - borrow;
04255             borrow = 0;
04256         }
04257     }
04258 
04259     /* Just assign the first few digits of a with considering */
04260     /* the borrow obtained so far because b has been exhausted. */
04261     while(a_pos > 0) {
04262         --c_pos;
04263         if(a->frac[--a_pos] < borrow) {
04264             c->frac[c_pos] = BASE + a->frac[a_pos] - borrow;
04265             borrow = 1;
04266         } else {
04267             c->frac[c_pos] = a->frac[a_pos] - borrow;
04268             borrow = 0;
04269         }
04270     }
04271     if(c_pos) c->frac[c_pos - 1] -= borrow;
04272     goto Exit;
04273 
04274 Assign_a:
04275     VpAsgn(c, a, 1);
04276     mrv = 0;
04277 
04278 Exit:
04279 #ifdef BIGDECIMAL_DEBUG
04280     if(gfDebug) {
04281         VPrint(stdout, "VpSubAbs exit: c=% \n", c);
04282     }
04283 #endif /* BIGDECIMAL_DEBUG */
04284     return mrv;
04285 }
04286 
04287 /*
04288  * Note: If(av+bv)>= HALF_BASE,then 1 will be added to the least significant
04289  *    digit of c(In case of addition).
04290  * ------------------------- figure of output -----------------------------------
04291  *      a =  xxxxxxxxxxx
04292  *      b =    xxxxxxxxxx
04293  *      c =xxxxxxxxxxxxxxx
04294  *      word_shift =  |   |
04295  *      right_word =  |    | (Total digits in RHSV)
04296  *      left_word  = |   |   (Total digits in LHSV)
04297  *      a_pos      =    |
04298  *      b_pos      =     |
04299  *      c_pos      =      |
04300  */
04301 static size_t
04302 VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, BDIGIT *av, BDIGIT *bv)
04303 {
04304     size_t left_word, right_word, word_shift;
04305     c->frac[0] = 0;
04306     *av = *bv = 0;
04307     word_shift =((a->exponent) -(b->exponent));
04308     left_word = b->Prec + word_shift;
04309     right_word = Max((a->Prec),left_word);
04310     left_word =(c->MaxPrec) - 1;    /* -1 ... prepare for round up */
04311     /*
04312      * check if 'round' is needed.
04313      */
04314     if(right_word > left_word) {    /* round ? */
04315         /*---------------------------------
04316          *  Actual size of a = xxxxxxAxx
04317          *  Actual size of b = xxxBxxxxx
04318          *  Max. size of   c = xxxxxx
04319          *  Round off        =   |-----|
04320          *  c_pos            =   |
04321          *  right_word       =   |
04322          *  a_pos            =    |
04323          */
04324         *c_pos = right_word = left_word + 1;    /* Set resulting precision */
04325         /* be equal to that of c */
04326         if((a->Prec) >=(c->MaxPrec)) {
04327             /*
04328              *   a =  xxxxxxAxxx
04329              *   c =  xxxxxx
04330              *   a_pos =    |
04331              */
04332             *a_pos = left_word;
04333             *av = a->frac[*a_pos];    /* av is 'A' shown in above. */
04334         } else {
04335             /*
04336              *   a = xxxxxxx
04337              *   c = xxxxxxxxxx
04338              *  a_pos =     |
04339              */
04340             *a_pos = a->Prec;
04341         }
04342         if((b->Prec + word_shift) >= c->MaxPrec) {
04343             /*
04344              *   a = xxxxxxxxx
04345              *   b =  xxxxxxxBxxx
04346              *   c = xxxxxxxxxxx
04347              *  b_pos =   |
04348              */
04349             if(c->MaxPrec >=(word_shift + 1)) {
04350                 *b_pos = c->MaxPrec - word_shift - 1;
04351                 *bv = b->frac[*b_pos];
04352             } else {
04353                 *b_pos = -1L;
04354             }
04355         } else {
04356             /*
04357              *   a = xxxxxxxxxxxxxxxx
04358              *   b =  xxxxxx
04359              *   c = xxxxxxxxxxxxx
04360              *  b_pos =     |
04361              */
04362             *b_pos = b->Prec;
04363         }
04364     } else {            /* The MaxPrec of c - 1 > The Prec of a + b  */
04365         /*
04366          *    a =   xxxxxxx
04367          *    b =   xxxxxx
04368          *    c = xxxxxxxxxxx
04369          *   c_pos =   |
04370          */
04371         *b_pos = b->Prec;
04372         *a_pos = a->Prec;
04373         *c_pos = right_word + 1;
04374     }
04375     c->Prec = *c_pos;
04376     c->exponent = a->exponent;
04377     if(!AddExponent(c,1)) return (size_t)-1L;
04378     return word_shift;
04379 }
04380 
04381 /*
04382  * Return number og significant digits
04383  *       c = a * b , Where a = a0a1a2 ... an
04384  *             b = b0b1b2 ... bm
04385  *             c = c0c1c2 ... cl
04386  *          a0 a1 ... an   * bm
04387  *       a0 a1 ... an   * bm-1
04388  *         .   .    .
04389  *       .   .   .
04390  *        a0 a1 .... an    * b0
04391  *      +_____________________________
04392  *     c0 c1 c2  ......  cl
04393  *     nc      <---|
04394  *     MaxAB |--------------------|
04395  */
04396 VP_EXPORT size_t
04397 VpMult(Real *c, Real *a, Real *b)
04398 {
04399     size_t MxIndA, MxIndB, MxIndAB, MxIndC;
04400     size_t ind_c, i, ii, nc;
04401     size_t ind_as, ind_ae, ind_bs;
04402     BDIGIT carry;
04403     BDIGIT_DBL s;
04404     Real *w;
04405 
04406 #ifdef BIGDECIMAL_DEBUG
04407     if(gfDebug) {
04408         VPrint(stdout, "VpMult(Enter): a=% \n", a);
04409         VPrint(stdout, "      b=% \n", b);
04410     }
04411 #endif /* BIGDECIMAL_DEBUG */
04412 
04413     if(!VpIsDefOP(c,a,b,3)) return 0; /* No significant digit */
04414 
04415     if(VpIsZero(a) || VpIsZero(b)) {
04416         /* at least a or b is zero */
04417         VpSetZero(c,VpGetSign(a)*VpGetSign(b));
04418         return 1; /* 0: 1 significant digit */
04419     }
04420 
04421     if(VpIsOne(a)) {
04422         VpAsgn(c, b, VpGetSign(a));
04423         goto Exit;
04424     }
04425     if(VpIsOne(b)) {
04426         VpAsgn(c, a, VpGetSign(b));
04427         goto Exit;
04428     }
04429     if((b->Prec) >(a->Prec)) {
04430         /* Adjust so that digits(a)>digits(b) */
04431         w = a;
04432         a = b;
04433         b = w;
04434     }
04435     w = NULL;
04436     MxIndA = a->Prec - 1;
04437     MxIndB = b->Prec - 1;
04438     MxIndC = c->MaxPrec - 1;
04439     MxIndAB = a->Prec + b->Prec - 1;
04440 
04441     if(MxIndC < MxIndAB) {    /* The Max. prec. of c < Prec(a)+Prec(b) */
04442         w = c;
04443         c = VpAlloc((size_t)((MxIndAB + 1) * BASE_FIG), "#0");
04444         MxIndC = MxIndAB;
04445     }
04446 
04447     /* set LHSV c info */
04448 
04449     c->exponent = a->exponent;    /* set exponent */
04450     if(!AddExponent(c,b->exponent)) {
04451         if(w) VpFree(c);
04452         return 0;
04453     }
04454     VpSetSign(c,VpGetSign(a)*VpGetSign(b));    /* set sign  */
04455     carry = 0;
04456     nc = ind_c = MxIndAB;
04457     memset(c->frac, 0, (nc + 1) * sizeof(BDIGIT));        /* Initialize c  */
04458     c->Prec = nc + 1;        /* set precision */
04459     for(nc = 0; nc < MxIndAB; ++nc, --ind_c) {
04460         if(nc < MxIndB) {    /* The left triangle of the Fig. */
04461             ind_as = MxIndA - nc;
04462             ind_ae = MxIndA;
04463             ind_bs = MxIndB;
04464         } else if(nc <= MxIndA) {    /* The middle rectangular of the Fig. */
04465             ind_as = MxIndA - nc;
04466             ind_ae = MxIndA -(nc - MxIndB);
04467             ind_bs = MxIndB;
04468         } else if(nc > MxIndA) {    /*  The right triangle of the Fig. */
04469             ind_as = 0;
04470             ind_ae = MxIndAB - nc - 1;
04471             ind_bs = MxIndB -(nc - MxIndA);
04472         }
04473 
04474         for(i = ind_as; i <= ind_ae; ++i) {
04475             s = (BDIGIT_DBL)a->frac[i] * b->frac[ind_bs--];
04476             carry = (BDIGIT)(s / BASE);
04477             s -= (BDIGIT_DBL)carry * BASE;
04478             c->frac[ind_c] += (BDIGIT)s;
04479             if(c->frac[ind_c] >= BASE) {
04480                 s = c->frac[ind_c] / BASE;
04481                 carry += (BDIGIT)s;
04482                 c->frac[ind_c] -= (BDIGIT)(s * BASE);
04483             }
04484             if(carry) {
04485                 ii = ind_c;
04486                 while(ii-- > 0) {
04487                     c->frac[ii] += carry;
04488                     if(c->frac[ii] >= BASE) {
04489                         carry = c->frac[ii] / BASE;
04490                         c->frac[ii] -= (carry * BASE);
04491                     } else {
04492                         break;
04493                     }
04494                 }
04495             }
04496         }
04497     }
04498     if(w != NULL) {        /* free work variable */
04499         VpNmlz(c);
04500         VpAsgn(w, c, 1);
04501         VpFree(c);
04502         c = w;
04503     } else {
04504         VpLimitRound(c,0);
04505     }
04506 
04507 Exit:
04508 #ifdef BIGDECIMAL_DEBUG
04509     if(gfDebug) {
04510         VPrint(stdout, "VpMult(c=a*b): c=% \n", c);
04511         VPrint(stdout, "      a=% \n", a);
04512         VPrint(stdout, "      b=% \n", b);
04513     }
04514 #endif /*BIGDECIMAL_DEBUG */
04515     return c->Prec*BASE_FIG;
04516 }
04517 
04518 /*
04519  *   c = a / b,  remainder = r
04520  */
04521 VP_EXPORT size_t
04522 VpDivd(Real *c, Real *r, Real *a, Real *b)
04523 {
04524     size_t word_a, word_b, word_c, word_r;
04525     size_t i, n, ind_a, ind_b, ind_c, ind_r;
04526     size_t nLoop;
04527     BDIGIT_DBL q, b1, b1p1, b1b2, b1b2p1, r1r2;
04528     BDIGIT borrow, borrow1, borrow2;
04529     BDIGIT_DBL qb;
04530 
04531 #ifdef BIGDECIMAL_DEBUG
04532     if(gfDebug) {
04533         VPrint(stdout, " VpDivd(c=a/b)  a=% \n", a);
04534         VPrint(stdout, "    b=% \n", b);
04535     }
04536 #endif /*BIGDECIMAL_DEBUG */
04537 
04538     VpSetNaN(r);
04539     if(!VpIsDefOP(c,a,b,4)) goto Exit;
04540     if(VpIsZero(a)&&VpIsZero(b)) {
04541         VpSetNaN(c);
04542         return VpException(VP_EXCEPTION_NaN,"(VpDivd) 0/0 not defined(NaN)",0);
04543     }
04544     if(VpIsZero(b)) {
04545         VpSetInf(c,VpGetSign(a)*VpGetSign(b));
04546         return VpException(VP_EXCEPTION_ZERODIVIDE,"(VpDivd) Divide by zero",0);
04547     }
04548     if(VpIsZero(a)) {
04549         /* numerator a is zero  */
04550         VpSetZero(c,VpGetSign(a)*VpGetSign(b));
04551         VpSetZero(r,VpGetSign(a)*VpGetSign(b));
04552         goto Exit;
04553     }
04554     if(VpIsOne(b)) {
04555         /* divide by one  */
04556         VpAsgn(c, a, VpGetSign(b));
04557         VpSetZero(r,VpGetSign(a));
04558         goto Exit;
04559     }
04560 
04561     word_a = a->Prec;
04562     word_b = b->Prec;
04563     word_c = c->MaxPrec;
04564     word_r = r->MaxPrec;
04565 
04566     ind_c = 0;
04567     ind_r = 1;
04568 
04569     if(word_a >= word_r) goto space_error;
04570 
04571     r->frac[0] = 0;
04572     while(ind_r <= word_a) {
04573         r->frac[ind_r] = a->frac[ind_r - 1];
04574         ++ind_r;
04575     }
04576 
04577     while(ind_r < word_r) r->frac[ind_r++] = 0;
04578     while(ind_c < word_c) c->frac[ind_c++] = 0;
04579 
04580     /* initial procedure */
04581     b1 = b1p1 = b->frac[0];
04582     if(b->Prec <= 1) {
04583         b1b2p1 = b1b2 = b1p1 * BASE;
04584     } else {
04585         b1p1 = b1 + 1;
04586         b1b2p1 = b1b2 = b1 * BASE + b->frac[1];
04587         if(b->Prec > 2) ++b1b2p1;
04588     }
04589 
04590     /* */
04591     /* loop start */
04592     ind_c = word_r - 1;
04593     nLoop = Min(word_c,ind_c);
04594     ind_c = 1;
04595     while(ind_c < nLoop) {
04596         if(r->frac[ind_c] == 0) {
04597             ++ind_c;
04598             continue;
04599         }
04600         r1r2 = (BDIGIT_DBL)r->frac[ind_c] * BASE + r->frac[ind_c + 1];
04601         if(r1r2 == b1b2) {
04602             /* The first two word digits is the same */
04603             ind_b = 2;
04604             ind_a = ind_c + 2;
04605             while(ind_b < word_b) {
04606                 if(r->frac[ind_a] < b->frac[ind_b]) goto div_b1p1;
04607                 if(r->frac[ind_a] > b->frac[ind_b]) break;
04608                 ++ind_a;
04609                 ++ind_b;
04610             }
04611             /* The first few word digits of r and b is the same and */
04612             /* the first different word digit of w is greater than that */
04613             /* of b, so quotinet is 1 and just subtract b from r. */
04614             borrow = 0;        /* quotient=1, then just r-b */
04615             ind_b = b->Prec - 1;
04616             ind_r = ind_c + ind_b;
04617             if(ind_r >= word_r) goto space_error;
04618             n = ind_b;
04619             for(i = 0; i <= n; ++i) {
04620                 if(r->frac[ind_r] < b->frac[ind_b] + borrow) {
04621                     r->frac[ind_r] += (BASE - (b->frac[ind_b] + borrow));
04622                     borrow = 1;
04623                 } else {
04624                     r->frac[ind_r] = r->frac[ind_r] - b->frac[ind_b] - borrow;
04625                     borrow = 0;
04626                 }
04627                 --ind_r;
04628                 --ind_b;
04629             }
04630             ++c->frac[ind_c];
04631             goto carry;
04632         }
04633         /* The first two word digits is not the same, */
04634         /* then compare magnitude, and divide actually. */
04635         if(r1r2 >= b1b2p1) {
04636             q = r1r2 / b1b2p1;  /* q == (BDIGIT)q  */
04637             c->frac[ind_c] += (BDIGIT)q;
04638             ind_r = b->Prec + ind_c - 1;
04639             goto sub_mult;
04640         }
04641 
04642 div_b1p1:
04643         if(ind_c + 1 >= word_c) goto out_side;
04644         q = r1r2 / b1p1;  /* q == (BDIGIT)q */
04645         c->frac[ind_c + 1] += (BDIGIT)q;
04646         ind_r = b->Prec + ind_c;
04647 
04648 sub_mult:
04649         borrow1 = borrow2 = 0;
04650         ind_b = word_b - 1;
04651         if(ind_r >= word_r) goto space_error;
04652         n = ind_b;
04653         for(i = 0; i <= n; ++i) {
04654             /* now, perform r = r - q * b */
04655             qb = q * b->frac[ind_b];
04656             if (qb < BASE) borrow1 = 0;
04657             else {
04658                 borrow1 = (BDIGIT)(qb / BASE);
04659                 qb -= (BDIGIT_DBL)borrow1 * BASE;       /* get qb < BASE */
04660             }
04661             if(r->frac[ind_r] < qb) {
04662                 r->frac[ind_r] += (BDIGIT)(BASE - qb);
04663                 borrow2 = borrow2 + borrow1 + 1;
04664             } else {
04665                 r->frac[ind_r] -= (BDIGIT)qb;
04666                 borrow2 += borrow1;
04667             }
04668             if(borrow2) {
04669                 if(r->frac[ind_r - 1] < borrow2) {
04670                     r->frac[ind_r - 1] += (BASE - borrow2);
04671                     borrow2 = 1;
04672                 } else {
04673                     r->frac[ind_r - 1] -= borrow2;
04674                     borrow2 = 0;
04675                 }
04676             }
04677             --ind_r;
04678             --ind_b;
04679         }
04680 
04681         r->frac[ind_r] -= borrow2;
04682 carry:
04683         ind_r = ind_c;
04684         while(c->frac[ind_r] >= BASE) {
04685             c->frac[ind_r] -= BASE;
04686             --ind_r;
04687             ++c->frac[ind_r];
04688         }
04689     }
04690     /* End of operation, now final arrangement */
04691 out_side:
04692     c->Prec = word_c;
04693     c->exponent = a->exponent;
04694     if(!AddExponent(c,2))   return 0;
04695     if(!AddExponent(c,-(b->exponent))) return 0;
04696 
04697     VpSetSign(c,VpGetSign(a)*VpGetSign(b));
04698     VpNmlz(c);            /* normalize c */
04699     r->Prec = word_r;
04700     r->exponent = a->exponent;
04701     if(!AddExponent(r,1)) return 0;
04702     VpSetSign(r,VpGetSign(a));
04703     VpNmlz(r);            /* normalize r(remainder) */
04704     goto Exit;
04705 
04706 space_error:
04707 #ifdef BIGDECIMAL_DEBUG
04708     if(gfDebug) {
04709         printf("   word_a=%lu\n", word_a);
04710         printf("   word_b=%lu\n", word_b);
04711         printf("   word_c=%lu\n", word_c);
04712         printf("   word_r=%lu\n", word_r);
04713         printf("   ind_r =%lu\n", ind_r);
04714     }
04715 #endif /* BIGDECIMAL_DEBUG */
04716     rb_bug("ERROR(VpDivd): space for remainder too small.");
04717 
04718 Exit:
04719 #ifdef BIGDECIMAL_DEBUG
04720     if(gfDebug) {
04721         VPrint(stdout, " VpDivd(c=a/b), c=% \n", c);
04722         VPrint(stdout, "    r=% \n", r);
04723     }
04724 #endif /* BIGDECIMAL_DEBUG */
04725     return c->Prec*BASE_FIG;
04726 }
04727 
04728 /*
04729  *  Input  a = 00000xxxxxxxx En(5 preceeding zeros)
04730  *  Output a = xxxxxxxx En-5
04731  */
04732 static int
04733 VpNmlz(Real *a)
04734 {
04735     size_t ind_a, i;
04736 
04737     if (!VpIsDef(a)) goto NoVal;
04738     if (VpIsZero(a)) goto NoVal;
04739 
04740     ind_a = a->Prec;
04741     while (ind_a--) {
04742         if (a->frac[ind_a]) {
04743             a->Prec = ind_a + 1;
04744             i = 0;
04745             while (a->frac[i] == 0) ++i;        /* skip the first few zeros */
04746             if (i) {
04747                 a->Prec -= i;
04748                 if (!AddExponent(a, -(SIGNED_VALUE)i)) return 0;
04749                 memmove(&a->frac[0], &a->frac[i], a->Prec*sizeof(BDIGIT));
04750             }
04751             return 1;
04752         }
04753     }
04754     /* a is zero(no non-zero digit) */
04755     VpSetZero(a, VpGetSign(a));
04756     return 0;
04757 
04758 NoVal:
04759     a->frac[0] = 0;
04760     a->Prec = 1;
04761     return 0;
04762 }
04763 
04764 /*
04765  *  VpComp = 0  ... if a=b,
04766  *   Pos  ... a>b,
04767  *   Neg  ... a<b.
04768  *   999  ... result undefined(NaN)
04769  */
04770 VP_EXPORT int
04771 VpComp(Real *a, Real *b)
04772 {
04773     int val;
04774     size_t mx, ind;
04775     int e;
04776     val = 0;
04777     if(VpIsNaN(a)||VpIsNaN(b)) return 999;
04778     if(!VpIsDef(a)) {
04779         if(!VpIsDef(b)) e = a->sign - b->sign;
04780         else             e = a->sign;
04781         if(e>0)   return  1;
04782         else if(e<0) return -1;
04783         else   return  0;
04784     }
04785     if(!VpIsDef(b)) {
04786         e = -b->sign;
04787         if(e>0) return  1;
04788         else return -1;
04789     }
04790     /* Zero check */
04791     if(VpIsZero(a)) {
04792         if(VpIsZero(b))      return 0; /* both zero */
04793         val = -VpGetSign(b);
04794         goto Exit;
04795     }
04796     if(VpIsZero(b)) {
04797         val = VpGetSign(a);
04798         goto Exit;
04799     }
04800 
04801     /* compare sign */
04802     if(VpGetSign(a) > VpGetSign(b)) {
04803         val = 1;        /* a>b */
04804         goto Exit;
04805     }
04806     if(VpGetSign(a) < VpGetSign(b)) {
04807         val = -1;        /* a<b */
04808         goto Exit;
04809     }
04810 
04811     /* a and b have same sign, && signe!=0,then compare exponent */
04812     if((a->exponent) >(b->exponent)) {
04813         val = VpGetSign(a);
04814         goto Exit;
04815     }
04816     if((a->exponent) <(b->exponent)) {
04817         val = -VpGetSign(b);
04818         goto Exit;
04819     }
04820 
04821     /* a and b have same exponent, then compare significand. */
04822     mx =((a->Prec) <(b->Prec)) ?(a->Prec) :(b->Prec);
04823     ind = 0;
04824     while(ind < mx) {
04825         if((a->frac[ind]) >(b->frac[ind])) {
04826             val = VpGetSign(a);
04827          goto Exit;
04828         }
04829         if((a->frac[ind]) <(b->frac[ind])) {
04830             val = -VpGetSign(b);
04831             goto Exit;
04832         }
04833         ++ind;
04834     }
04835     if((a->Prec) >(b->Prec)) {
04836         val = VpGetSign(a);
04837     } else if((a->Prec) <(b->Prec)) {
04838         val = -VpGetSign(b);
04839     }
04840 
04841 Exit:
04842     if  (val> 1) val =  1;
04843     else if(val<-1) val = -1;
04844 
04845 #ifdef BIGDECIMAL_DEBUG
04846     if(gfDebug) {
04847         VPrint(stdout, " VpComp a=%\n", a);
04848         VPrint(stdout, "  b=%\n", b);
04849         printf("  ans=%d\n", val);
04850     }
04851 #endif /* BIGDECIMAL_DEBUG */
04852     return (int)val;
04853 }
04854 
04855 #ifdef BIGDECIMAL_ENABLE_VPRINT
04856 /*
04857  *    cntl_chr ... ASCIIZ Character, print control characters
04858  *     Available control codes:
04859  *      %  ... VP variable. To print '%', use '%%'.
04860  *      \n ... new line
04861  *      \b ... backspace
04862  *           ... tab
04863  *     Note: % must must not appear more than once
04864  *    a  ... VP variable to be printed
04865  */
04866 VP_EXPORT int
04867 VPrint(FILE *fp, const char *cntl_chr, Real *a)
04868 {
04869     size_t i, j, nc, nd, ZeroSup;
04870     BDIGIT m, e, nn;
04871 
04872     /* Check if NaN & Inf. */
04873     if(VpIsNaN(a)) {
04874         fprintf(fp,SZ_NaN);
04875         return 8;
04876     }
04877     if(VpIsPosInf(a)) {
04878         fprintf(fp,SZ_INF);
04879         return 8;
04880     }
04881     if(VpIsNegInf(a)) {
04882         fprintf(fp,SZ_NINF);
04883         return 9;
04884     }
04885     if(VpIsZero(a)) {
04886         fprintf(fp,"0.0");
04887         return 3;
04888     }
04889 
04890     j = 0;
04891     nd = nc = 0;        /*  nd : number of digits in fraction part(every 10 digits, */
04892     /*    nd<=10). */
04893     /*  nc : number of caracters printed  */
04894     ZeroSup = 1;        /* Flag not to print the leading zeros as 0.00xxxxEnn */
04895     while(*(cntl_chr + j)) {
04896         if((*(cntl_chr + j) == '%') &&(*(cntl_chr + j + 1) != '%')) {
04897          nc = 0;
04898          if(!VpIsZero(a)) {
04899                 if(VpGetSign(a) < 0) {
04900                     fprintf(fp, "-");
04901                     ++nc;
04902                 }
04903                 nc += fprintf(fp, "0.");
04904                 for(i=0; i < a->Prec; ++i) {
04905                     m = BASE1;
04906                     e = a->frac[i];
04907                     while(m) {
04908                         nn = e / m;
04909                         if((!ZeroSup) || nn) {
04910                             nc += fprintf(fp, "%lu", (unsigned long)nn);    /* The leading zero(s) */
04911                             /* as 0.00xx will not */
04912                             /* be printed. */
04913                             ++nd;
04914                             ZeroSup = 0;    /* Set to print succeeding zeros */
04915                         }
04916                         if(nd >= 10) {    /* print ' ' after every 10 digits */
04917                             nd = 0;
04918                             nc += fprintf(fp, " ");
04919                         }
04920                         e = e - nn * m;
04921                         m /= 10;
04922                     }
04923                 }
04924                 nc += fprintf(fp, "E%"PRIdSIZE, VpExponent10(a));
04925             } else {
04926                 nc += fprintf(fp, "0.0");
04927             }
04928         } else {
04929             ++nc;
04930             if(*(cntl_chr + j) == '\\') {
04931                 switch(*(cntl_chr + j + 1)) {
04932                 case 'n':
04933                     fprintf(fp, "\n");
04934                     ++j;
04935                     break;
04936                 case 't':
04937                     fprintf(fp, "\t");
04938                     ++j;
04939                  break;
04940                 case 'b':
04941                     fprintf(fp, "\n");
04942                     ++j;
04943                     break;
04944                 default:
04945                     fprintf(fp, "%c", *(cntl_chr + j));
04946                     break;
04947                 }
04948             } else {
04949                 fprintf(fp, "%c", *(cntl_chr + j));
04950                 if(*(cntl_chr + j) == '%') ++j;
04951             }
04952         }
04953         j++;
04954     }
04955     return (int)nc;
04956 }
04957 #endif /* BIGDECIMAL_ENABLE_VPRINT */
04958 
04959 static void
04960 VpFormatSt(char *psz, size_t fFmt)
04961 {
04962     size_t ie, i, nf = 0;
04963     char ch;
04964 
04965     if(fFmt<=0) return;
04966 
04967     ie = strlen(psz);
04968     for(i = 0; i < ie; ++i) {
04969         ch = psz[i];
04970         if(!ch) break;
04971         if(ISSPACE(ch) || ch=='-' || ch=='+') continue;
04972         if(ch == '.')                { nf = 0;continue;}
04973         if(ch == 'E') break;
04974         nf++;
04975         if(nf > fFmt) {
04976             memmove(psz + i + 1, psz + i, ie - i + 1);
04977             ++ie;
04978             nf = 0;
04979             psz[i] = ' ';
04980         }
04981     }
04982 }
04983 
04984 VP_EXPORT ssize_t
04985 VpExponent10(Real *a)
04986 {
04987     ssize_t ex;
04988     size_t n;
04989 
04990     if (!VpHasVal(a)) return 0;
04991 
04992     ex = a->exponent * (ssize_t)BASE_FIG;
04993     n = BASE1;
04994     while ((a->frac[0] / n) == 0) {
04995          --ex;
04996          n /= 10;
04997     }
04998     return ex;
04999 }
05000 
05001 VP_EXPORT void
05002 VpSzMantissa(Real *a,char *psz)
05003 {
05004     size_t i, n, ZeroSup;
05005     BDIGIT_DBL m, e, nn;
05006 
05007     if(VpIsNaN(a)) {
05008         sprintf(psz,SZ_NaN);
05009         return;
05010     }
05011     if(VpIsPosInf(a)) {
05012         sprintf(psz,SZ_INF);
05013         return;
05014     }
05015     if(VpIsNegInf(a)) {
05016         sprintf(psz,SZ_NINF);
05017         return;
05018     }
05019 
05020     ZeroSup = 1;        /* Flag not to print the leading zeros as 0.00xxxxEnn */
05021     if(!VpIsZero(a)) {
05022         if(VpGetSign(a) < 0) *psz++ = '-';
05023         n = a->Prec;
05024         for (i=0; i < n; ++i) {
05025             m = BASE1;
05026             e = a->frac[i];
05027             while (m) {
05028                 nn = e / m;
05029                 if((!ZeroSup) || nn) {
05030                     sprintf(psz, "%lu", (unsigned long)nn);    /* The leading zero(s) */
05031                     psz += strlen(psz);
05032                     /* as 0.00xx will be ignored. */
05033                     ZeroSup = 0;    /* Set to print succeeding zeros */
05034                 }
05035                 e = e - nn * m;
05036                 m /= 10;
05037             }
05038         }
05039         *psz = 0;
05040         while(psz[-1]=='0') *(--psz) = 0;
05041     } else {
05042         if(VpIsPosZero(a)) sprintf(psz, "0");
05043         else      sprintf(psz, "-0");
05044     }
05045 }
05046 
05047 VP_EXPORT int
05048 VpToSpecialString(Real *a,char *psz,int fPlus)
05049 /* fPlus =0:default, =1: set ' ' before digits , =2: set '+' before digits. */
05050 {
05051     if(VpIsNaN(a)) {
05052         sprintf(psz,SZ_NaN);
05053         return 1;
05054     }
05055 
05056     if(VpIsPosInf(a)) {
05057         if(fPlus==1) {
05058            *psz++ = ' ';
05059         } else if(fPlus==2) {
05060            *psz++ = '+';
05061         }
05062         sprintf(psz,SZ_INF);
05063         return 1;
05064     }
05065     if(VpIsNegInf(a)) {
05066         sprintf(psz,SZ_NINF);
05067         return 1;
05068     }
05069     if(VpIsZero(a)) {
05070         if(VpIsPosZero(a)) {
05071             if(fPlus==1)      sprintf(psz, " 0.0");
05072             else if(fPlus==2) sprintf(psz, "+0.0");
05073             else              sprintf(psz, "0.0");
05074         } else    sprintf(psz, "-0.0");
05075         return 1;
05076     }
05077     return 0;
05078 }
05079 
05080 VP_EXPORT void
05081 VpToString(Real *a, char *psz, size_t fFmt, int fPlus)
05082 /* fPlus =0:default, =1: set ' ' before digits , =2:set '+' before digits. */
05083 {
05084     size_t i, n, ZeroSup;
05085     BDIGIT shift, m, e, nn;
05086     char *pszSav = psz;
05087     ssize_t ex;
05088 
05089     if (VpToSpecialString(a, psz, fPlus)) return;
05090 
05091     ZeroSup = 1;    /* Flag not to print the leading zeros as 0.00xxxxEnn */
05092 
05093     if (VpGetSign(a) < 0) *psz++ = '-';
05094     else if (fPlus == 1)  *psz++ = ' ';
05095     else if (fPlus == 2)  *psz++ = '+';
05096 
05097     *psz++ = '0';
05098     *psz++ = '.';
05099     n = a->Prec;
05100     for(i=0;i < n;++i) {
05101         m = BASE1;
05102         e = a->frac[i];
05103         while(m) {
05104             nn = e / m;
05105             if((!ZeroSup) || nn) {
05106                 sprintf(psz, "%lu", (unsigned long)nn);    /* The reading zero(s) */
05107                 psz += strlen(psz);
05108                 /* as 0.00xx will be ignored. */
05109                 ZeroSup = 0;    /* Set to print succeeding zeros */
05110             }
05111             e = e - nn * m;
05112             m /= 10;
05113         }
05114     }
05115     ex = a->exponent * (ssize_t)BASE_FIG;
05116     shift = BASE1;
05117     while(a->frac[0] / shift == 0) {
05118         --ex;
05119         shift /= 10;
05120     }
05121     while(psz[-1]=='0') *(--psz) = 0;
05122     sprintf(psz, "E%"PRIdSIZE, ex);
05123     if(fFmt) VpFormatSt(pszSav, fFmt);
05124 }
05125 
05126 VP_EXPORT void
05127 VpToFString(Real *a, char *psz, size_t fFmt, int fPlus)
05128 /* fPlus =0:default,=1: set ' ' before digits ,set '+' before digits. */
05129 {
05130     size_t i, n;
05131     BDIGIT m, e, nn;
05132     char *pszSav = psz;
05133     ssize_t ex;
05134 
05135     if(VpToSpecialString(a,psz,fPlus)) return;
05136 
05137     if(VpGetSign(a) < 0) *psz++ = '-';
05138     else if(fPlus==1)    *psz++ = ' ';
05139     else if(fPlus==2)    *psz++ = '+';
05140 
05141     n  = a->Prec;
05142     ex = a->exponent;
05143     if(ex<=0) {
05144        *psz++ = '0';*psz++ = '.';
05145        while(ex<0) {
05146           for(i=0;i<BASE_FIG;++i) *psz++ = '0';
05147           ++ex;
05148        }
05149        ex = -1;
05150     }
05151 
05152     for(i=0;i < n;++i) {
05153        --ex;
05154        if(i==0 && ex >= 0) {
05155            sprintf(psz, "%lu", (unsigned long)a->frac[i]);
05156            psz += strlen(psz);
05157        } else {
05158            m = BASE1;
05159            e = a->frac[i];
05160            while(m) {
05161                nn = e / m;
05162                *psz++ = (char)(nn + '0');
05163                e = e - nn * m;
05164                m /= 10;
05165            }
05166        }
05167        if(ex == 0) *psz++ = '.';
05168     }
05169     while(--ex>=0) {
05170        m = BASE;
05171        while(m/=10) *psz++ = '0';
05172        if(ex == 0) *psz++ = '.';
05173     }
05174     *psz = 0;
05175     while(psz[-1]=='0') *(--psz) = 0;
05176     if(psz[-1]=='.') sprintf(psz, "0");
05177     if(fFmt) VpFormatSt(pszSav, fFmt);
05178 }
05179 
05180 /*
05181  *  [Output]
05182  *   a[]  ... variable to be assigned the value.
05183  *  [Input]
05184  *   int_chr[]  ... integer part(may include '+/-').
05185  *   ni   ... number of characters in int_chr[],not including '+/-'.
05186  *   frac[]  ... fraction part.
05187  *   nf   ... number of characters in frac[].
05188  *   exp_chr[]  ... exponent part(including '+/-').
05189  *   ne   ... number of characters in exp_chr[],not including '+/-'.
05190  */
05191 VP_EXPORT int
05192 VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne)
05193 {
05194     size_t i, j, ind_a, ma, mi, me;
05195     SIGNED_VALUE e, es, eb, ef;
05196     int  sign, signe, exponent_overflow;
05197 
05198     /* get exponent part */
05199     e = 0;
05200     ma = a->MaxPrec;
05201     mi = ni;
05202     me = ne;
05203     signe = 1;
05204     exponent_overflow = 0;
05205     memset(a->frac, 0, ma * sizeof(BDIGIT));
05206     if (ne > 0) {
05207         i = 0;
05208         if (exp_chr[0] == '-') {
05209             signe = -1;
05210             ++i;
05211             ++me;
05212         }
05213         else if (exp_chr[0] == '+') {
05214             ++i;
05215             ++me;
05216         }
05217         while (i < me) {
05218             es = e * (SIGNED_VALUE)BASE_FIG;
05219             e = e * 10 + exp_chr[i] - '0';
05220             if (es > (SIGNED_VALUE)(e*BASE_FIG)) {
05221                 exponent_overflow = 1;
05222                 e = es; /* keep sign */
05223                 break;
05224             }
05225             ++i;
05226         }
05227     }
05228 
05229     /* get integer part */
05230     i = 0;
05231     sign = 1;
05232     if(1 /*ni >= 0*/) {
05233         if(int_chr[0] == '-') {
05234             sign = -1;
05235             ++i;
05236             ++mi;
05237         } else if(int_chr[0] == '+') {
05238             ++i;
05239             ++mi;
05240         }
05241     }
05242 
05243     e = signe * e;        /* e: The value of exponent part. */
05244     e = e + ni;        /* set actual exponent size. */
05245 
05246     if (e > 0) signe = 1;
05247     else       signe = -1;
05248 
05249     /* Adjust the exponent so that it is the multiple of BASE_FIG. */
05250     j = 0;
05251     ef = 1;
05252     while (ef) {
05253         if (e >= 0) eb =  e;
05254         else        eb = -e;
05255         ef = eb / (SIGNED_VALUE)BASE_FIG;
05256         ef = eb - ef * (SIGNED_VALUE)BASE_FIG;
05257         if (ef) {
05258             ++j;        /* Means to add one more preceeding zero */
05259             ++e;
05260         }
05261     }
05262 
05263     eb = e / (SIGNED_VALUE)BASE_FIG;
05264 
05265     if (exponent_overflow) {
05266         int zero = 1;
05267         for (     ; i < mi && zero; i++) zero = int_chr[i] == '0';
05268         for (i = 0; i < nf && zero; i++) zero = frac[i] == '0';
05269         if (!zero && signe > 0) {
05270             VpSetInf(a, sign);
05271             VpException(VP_EXCEPTION_INFINITY, "exponent overflow",0);
05272         }
05273         else VpSetZero(a, sign);
05274         return 1;
05275     }
05276 
05277     ind_a = 0;
05278     while (i < mi) {
05279         a->frac[ind_a] = 0;
05280         while ((j < BASE_FIG) && (i < mi)) {
05281             a->frac[ind_a] = a->frac[ind_a] * 10 + int_chr[i] - '0';
05282             ++j;
05283             ++i;
05284         }
05285         if (i < mi) {
05286             ++ind_a;
05287             if (ind_a >= ma) goto over_flow;
05288             j = 0;
05289         }
05290     }
05291 
05292     /* get fraction part */
05293 
05294     i = 0;
05295     while(i < nf) {
05296         while((j < BASE_FIG) && (i < nf)) {
05297             a->frac[ind_a] = a->frac[ind_a] * 10 + frac[i] - '0';
05298             ++j;
05299             ++i;
05300         }
05301         if(i < nf) {
05302             ++ind_a;
05303             if(ind_a >= ma) goto over_flow;
05304             j = 0;
05305         }
05306     }
05307     goto Final;
05308 
05309 over_flow:
05310     rb_warn("Conversion from String to BigDecimal overflow (last few digits discarded).");
05311 
05312 Final:
05313     if (ind_a >= ma) ind_a = ma - 1;
05314     while (j < BASE_FIG) {
05315         a->frac[ind_a] = a->frac[ind_a] * 10;
05316         ++j;
05317     }
05318     a->Prec = ind_a + 1;
05319     a->exponent = eb;
05320     VpSetSign(a,sign);
05321     VpNmlz(a);
05322     return 1;
05323 }
05324 
05325 /*
05326  * [Input]
05327  *   *m  ... Real
05328  * [Output]
05329  *   *d  ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig.
05330  *   *e  ... exponent of m.
05331  * DBLE_FIG ... Number of digits in a double variable.
05332  *
05333  *  m -> d*10**e, 0<d<BASE
05334  * [Returns]
05335  *   0 ... Zero
05336  *   1 ... Normal
05337  *   2 ... Infinity
05338  *  -1 ... NaN
05339  */
05340 VP_EXPORT int
05341 VpVtoD(double *d, SIGNED_VALUE *e, Real *m)
05342 {
05343     size_t ind_m, mm, fig;
05344     double div;
05345     int    f = 1;
05346 
05347     if(VpIsNaN(m)) {
05348         *d = VpGetDoubleNaN();
05349         *e = 0;
05350         f = -1; /* NaN */
05351         goto Exit;
05352     } else
05353     if(VpIsPosZero(m)) {
05354         *d = 0.0;
05355         *e = 0;
05356         f  = 0;
05357         goto Exit;
05358     } else
05359     if(VpIsNegZero(m)) {
05360         *d = VpGetDoubleNegZero();
05361         *e = 0;
05362         f  = 0;
05363         goto Exit;
05364     } else
05365     if(VpIsPosInf(m)) {
05366         *d = VpGetDoublePosInf();
05367         *e = 0;
05368         f  = 2;
05369         goto Exit;
05370     } else
05371     if(VpIsNegInf(m)) {
05372         *d = VpGetDoubleNegInf();
05373         *e = 0;
05374         f  = 2;
05375         goto Exit;
05376     }
05377     /* Normal number */
05378     fig =(DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
05379     ind_m = 0;
05380     mm = Min(fig,(m->Prec));
05381     *d = 0.0;
05382     div = 1.;
05383     while(ind_m < mm) {
05384         div /= (double)BASE;
05385         *d = *d + (double)m->frac[ind_m++] * div;
05386     }
05387     *e = m->exponent * (SIGNED_VALUE)BASE_FIG;
05388     *d *= VpGetSign(m);
05389 
05390 Exit:
05391 #ifdef BIGDECIMAL_DEBUG
05392     if(gfDebug) {
05393         VPrint(stdout, " VpVtoD: m=%\n", m);
05394         printf("   d=%e * 10 **%ld\n", *d, *e);
05395         printf("   DBLE_FIG = %d\n", DBLE_FIG);
05396     }
05397 #endif /*BIGDECIMAL_DEBUG */
05398     return f;
05399 }
05400 
05401 /*
05402  * m <- d
05403  */
05404 VP_EXPORT void
05405 VpDtoV(Real *m, double d)
05406 {
05407     size_t ind_m, mm;
05408     SIGNED_VALUE ne;
05409     BDIGIT i;
05410     double  val, val2;
05411 
05412     if(isnan(d)) {
05413         VpSetNaN(m);
05414         goto Exit;
05415     }
05416     if(isinf(d)) {
05417         if(d>0.0) VpSetPosInf(m);
05418         else   VpSetNegInf(m);
05419         goto Exit;
05420     }
05421 
05422     if(d == 0.0) {
05423         VpSetZero(m,1);
05424         goto Exit;
05425     }
05426     val =(d > 0.) ? d :(-d);
05427     ne = 0;
05428     if(val >= 1.0) {
05429         while(val >= 1.0) {
05430             val /= (double)BASE;
05431             ++ne;
05432         }
05433     } else {
05434         val2 = 1.0 /(double)BASE;
05435         while(val < val2) {
05436             val *= (double)BASE;
05437             --ne;
05438         }
05439     }
05440     /* Now val = 0.xxxxx*BASE**ne */
05441 
05442     mm = m->MaxPrec;
05443     memset(m->frac, 0, mm * sizeof(BDIGIT));
05444     for(ind_m = 0;val > 0.0 && ind_m < mm;ind_m++) {
05445         val *= (double)BASE;
05446         i = (BDIGIT)val;
05447         val -= (double)i;
05448         m->frac[ind_m] = i;
05449     }
05450     if(ind_m >= mm) ind_m = mm - 1;
05451     VpSetSign(m, (d > 0.0) ? 1 : -1);
05452     m->Prec = ind_m + 1;
05453     m->exponent = ne;
05454 
05455     VpInternalRound(m, 0, (m->Prec > 0) ? m->frac[m->Prec-1] : 0,
05456                       (BDIGIT)(val*(double)BASE));
05457 
05458 Exit:
05459 #ifdef BIGDECIMAL_DEBUG
05460     if(gfDebug) {
05461         printf("VpDtoV d=%30.30e\n", d);
05462         VPrint(stdout, "  m=%\n", m);
05463     }
05464 #endif /* BIGDECIMAL_DEBUG */
05465     return;
05466 }
05467 
05468 /*
05469  *  m <- ival
05470  */
05471 #if 0  /* unused */
05472 VP_EXPORT void
05473 VpItoV(Real *m, SIGNED_VALUE ival)
05474 {
05475     size_t mm, ind_m;
05476     size_t val, v1, v2, v;
05477     int isign;
05478     SIGNED_VALUE ne;
05479 
05480     if(ival == 0) {
05481         VpSetZero(m,1);
05482         goto Exit;
05483     }
05484     isign = 1;
05485     val = ival;
05486     if(ival < 0) {
05487         isign = -1;
05488         val =(size_t)(-ival);
05489     }
05490     ne = 0;
05491     ind_m = 0;
05492     mm = m->MaxPrec;
05493     while(ind_m < mm) {
05494         m->frac[ind_m] = 0;
05495         ++ind_m;
05496     }
05497     ind_m = 0;
05498     while(val > 0) {
05499         if(val) {
05500          v1 = val;
05501          v2 = 1;
05502             while(v1 >= BASE) {
05503                 v1 /= BASE;
05504                 v2 *= BASE;
05505             }
05506             val = val - v2 * v1;
05507             v = v1;
05508         } else {
05509             v = 0;
05510         }
05511         m->frac[ind_m] = v;
05512         ++ind_m;
05513         ++ne;
05514     }
05515     m->Prec = ind_m - 1;
05516     m->exponent = ne;
05517     VpSetSign(m,isign);
05518     VpNmlz(m);
05519 
05520 Exit:
05521 #ifdef BIGDECIMAL_DEBUG
05522     if(gfDebug) {
05523         printf(" VpItoV i=%d\n", ival);
05524         VPrint(stdout, "  m=%\n", m);
05525     }
05526 #endif /* BIGDECIMAL_DEBUG */
05527     return;
05528 }
05529 #endif
05530 
05531 /*
05532  * y = SQRT(x),  y*y - x =>0
05533  */
05534 VP_EXPORT int
05535 VpSqrt(Real *y, Real *x)
05536 {
05537     Real *f = NULL;
05538     Real *r = NULL;
05539     size_t y_prec;
05540     SIGNED_VALUE n, e;
05541     SIGNED_VALUE prec;
05542     ssize_t nr;
05543     double val;
05544 
05545     /* Zero, NaN or Infinity ? */
05546     if(!VpHasVal(x)) {
05547         if(VpIsZero(x)||VpGetSign(x)>0) {
05548             VpAsgn(y,x,1);
05549             goto Exit;
05550         }
05551         VpSetNaN(y);
05552         return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(NaN or negative value)",0);
05553         goto Exit;
05554     }
05555 
05556      /* Negative ? */
05557     if(VpGetSign(x) < 0) {
05558         VpSetNaN(y);
05559         return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(negative value)",0);
05560     }
05561 
05562     /* One ? */
05563     if(VpIsOne(x)) {
05564         VpSetOne(y);
05565         goto Exit;
05566     }
05567 
05568     n = (SIGNED_VALUE)y->MaxPrec;
05569     if (x->MaxPrec > (size_t)n) n = (ssize_t)x->MaxPrec;
05570     /* allocate temporally variables  */
05571     f = VpAlloc(y->MaxPrec * (BASE_FIG + 2), "#1");
05572     r = VpAlloc((n + n) * (BASE_FIG + 2), "#1");
05573 
05574     nr = 0;
05575     y_prec = y->MaxPrec;
05576 
05577     prec = x->exponent - (ssize_t)y_prec;
05578     if (x->exponent > 0)
05579         ++prec;
05580     else
05581         --prec;
05582 
05583     VpVtoD(&val, &e, x);    /* val <- x  */
05584     e /= (SIGNED_VALUE)BASE_FIG;
05585     n = e / 2;
05586     if (e - n * 2 != 0) {
05587         val /= BASE;
05588         n = (e + 1) / 2;
05589     }
05590     VpDtoV(y, sqrt(val));    /* y <- sqrt(val) */
05591     y->exponent += n;
05592     n = (SIGNED_VALUE)((DBLE_FIG + BASE_FIG - 1) / BASE_FIG);
05593     y->MaxPrec = Min((size_t)n , y_prec);
05594     f->MaxPrec = y->MaxPrec + 1;
05595     n = (SIGNED_VALUE)(y_prec * BASE_FIG);
05596     if (n < (SIGNED_VALUE)maxnr) n = (SIGNED_VALUE)maxnr;
05597     do {
05598         y->MaxPrec *= 2;
05599         if (y->MaxPrec > y_prec) y->MaxPrec = y_prec;
05600         f->MaxPrec = y->MaxPrec;
05601         VpDivd(f, r, x, y);     /* f = x/y    */
05602         VpAddSub(r, f, y, -1);  /* r = f - y  */
05603         VpMult(f, VpPt5, r);    /* f = 0.5*r  */
05604         if(VpIsZero(f))         goto converge;
05605         VpAddSub(r, f, y, 1);   /* r = y + f  */
05606         VpAsgn(y, r, 1);        /* y = r      */
05607         if(f->exponent <= prec) goto converge;
05608     } while(++nr < n);
05609     /* */
05610 #ifdef BIGDECIMAL_DEBUG
05611     if(gfDebug) {
05612         printf("ERROR(VpSqrt): did not converge within %ld iterations.\n",
05613             nr);
05614     }
05615 #endif /* BIGDECIMAL_DEBUG */
05616     y->MaxPrec = y_prec;
05617 
05618 converge:
05619     VpChangeSign(y, 1);
05620 #ifdef BIGDECIMAL_DEBUG
05621     if(gfDebug) {
05622         VpMult(r, y, y);
05623         VpAddSub(f, x, r, -1);
05624         printf("VpSqrt: iterations = %"PRIdSIZE"\n", nr);
05625         VPrint(stdout, "  y =% \n", y);
05626         VPrint(stdout, "  x =% \n", x);
05627         VPrint(stdout, "  x-y*y = % \n", f);
05628     }
05629 #endif /* BIGDECIMAL_DEBUG */
05630     y->MaxPrec = y_prec;
05631 
05632 Exit:
05633     VpFree(f);
05634     VpFree(r);
05635     return 1;
05636 }
05637 
05638 /*
05639  *
05640  * nf: digit position for operation.
05641  *
05642  */
05643 VP_EXPORT int
05644 VpMidRound(Real *y, unsigned short f, ssize_t nf)
05645 /*
05646  * Round reletively from the decimal point.
05647  *    f: rounding mode
05648  *   nf: digit location to round from the the decimal point.
05649  */
05650 {
05651     /* fracf: any positive digit under rounding position? */
05652     /* fracf_1further: any positive digits under one further than the rounding position? */
05653     /* exptoadd: number of digits needed to compensate negative nf */
05654     int fracf, fracf_1further;
05655     ssize_t n,i,ix,ioffset, exptoadd;
05656     BDIGIT v, shifter;
05657     BDIGIT div;
05658 
05659     nf += y->exponent * (ssize_t)BASE_FIG;
05660     exptoadd=0;
05661     if (nf < 0) {
05662         /* rounding position too left(large). */
05663         if((f!=VP_ROUND_CEIL) && (f!=VP_ROUND_FLOOR)) {
05664             VpSetZero(y,VpGetSign(y)); /* truncate everything */
05665             return 0;
05666         }
05667         exptoadd = -nf;
05668         nf = 0;
05669     }
05670 
05671     ix = nf / (ssize_t)BASE_FIG;
05672     if ((size_t)ix >= y->Prec) return 0;  /* rounding position too right(small). */
05673     v = y->frac[ix];
05674 
05675     ioffset = nf - ix*(ssize_t)BASE_FIG;
05676     n = (ssize_t)BASE_FIG - ioffset - 1;
05677     for (shifter=1,i=0; i<n; ++i) shifter *= 10;
05678 
05679     /* so the representation used (in y->frac) is an array of BDIGIT, where
05680        each BDIGIT contains a value between 0 and BASE-1, consisting of BASE_FIG
05681        decimal places.
05682 
05683        (that numbers of decimal places are typed as ssize_t is somewhat confusing)
05684 
05685        nf is now position (in decimal places) of the digit from the start of
05686           the array.
05687        ix is the position (in BDIGITS) of the BDIGIT containing the decimal digit,
05688           from the start of the array.
05689        v is the value of this BDIGIT
05690        ioffset is the number of extra decimal places along of this decimal digit
05691           within v.
05692        n is the number of decimal digits remaining within v after this decimal digit
05693        shifter is 10**n,
05694        v % shifter are the remaining digits within v
05695        v % (shifter * 10) are the digit together with the remaining digits within v
05696        v / shifter are the digit's predecessors together with the digit
05697        div = v / shifter / 10 is just the digit's precessors
05698        (v / shifter) - div*10 is just the digit, which is what v ends up being reassigned to.
05699     */
05700 
05701     fracf = (v % (shifter * 10) > 0);
05702     fracf_1further = ((v % shifter) > 0);
05703 
05704     v /= shifter;
05705     div = v / 10;
05706     v = v - div*10;
05707     /* now v is just the digit required.
05708        now fracf is whether the digit or any of the remaining digits within v are non-zero
05709        now fracf_1further is whether any of the remaining digits within v are non-zero
05710     */
05711 
05712     /* now check all the remaining BDIGITS for zero-ness a whole BDIGIT at a time.
05713        if we spot any non-zeroness, that means that we foudn a positive digit under
05714        rounding position, and we also found a positive digit under one further than
05715        the rounding position, so both searches (to see if any such non-zero digit exists)
05716        can stop */
05717 
05718     for (i=ix+1; (size_t)i < y->Prec; i++) {
05719         if (y->frac[i] % BASE) {
05720             fracf = fracf_1further = 1;
05721             break;
05722         }
05723     }
05724 
05725     /* now fracf = does any positive digit exist under the rounding position?
05726        now fracf_1further = does any positive digit exist under one further than the
05727        rounding position?
05728        now v = the first digit under the rounding position */
05729 
05730     /* drop digits after pointed digit */
05731     memset(y->frac+ix+1, 0, (y->Prec - (ix+1)) * sizeof(BDIGIT));
05732 
05733     switch(f) {
05734     case VP_ROUND_DOWN: /* Truncate */
05735          break;
05736     case VP_ROUND_UP:   /* Roundup */
05737         if (fracf) ++div;
05738         break;
05739     case VP_ROUND_HALF_UP:
05740         if (v>=5) ++div;
05741         break;
05742     case VP_ROUND_HALF_DOWN:
05743         if (v > 5 || (v == 5 && fracf_1further)) ++div;
05744         break;
05745     case VP_ROUND_CEIL:
05746         if (fracf && (VpGetSign(y)>0)) ++div;
05747         break;
05748     case VP_ROUND_FLOOR:
05749         if (fracf && (VpGetSign(y)<0)) ++div;
05750         break;
05751     case VP_ROUND_HALF_EVEN: /* Banker's rounding */
05752         if (v > 5) ++div;
05753         else if (v == 5) {
05754             if (fracf_1further) {
05755               ++div;
05756             }
05757             else {
05758                 if (ioffset == 0) {
05759                     /* v is the first decimal digit of its BDIGIT;
05760                        need to grab the previous BDIGIT if present
05761                        to check for evenness of the previous decimal
05762                        digit (which is same as that of the BDIGIT since
05763                        base 10 has a factor of 2) */
05764                     if (ix && (y->frac[ix-1] % 2)) ++div;
05765                 }
05766                 else {
05767                     if (div % 2) ++div;
05768                 }
05769             }
05770         }
05771         break;
05772     }
05773     for (i=0; i<=n; ++i) div *= 10;
05774     if (div>=BASE) {
05775         if(ix) {
05776             y->frac[ix] = 0;
05777             VpRdup(y,ix);
05778         } else {
05779             short s = VpGetSign(y);
05780             SIGNED_VALUE e = y->exponent;
05781             VpSetOne(y);
05782             VpSetSign(y, s);
05783             y->exponent = e+1;
05784         }
05785     } else {
05786         y->frac[ix] = div;
05787         VpNmlz(y);
05788     }
05789     if (exptoadd > 0) {
05790         y->exponent += (SIGNED_VALUE)(exptoadd/BASE_FIG);
05791         exptoadd %= (ssize_t)BASE_FIG;
05792         for(i=0;i<exptoadd;i++) {
05793             y->frac[0] *= 10;
05794             if (y->frac[0] >= BASE) {
05795                 y->frac[0] /= BASE;
05796                 y->exponent++;
05797             }
05798         }
05799     }
05800     return 1;
05801 }
05802 
05803 VP_EXPORT int
05804 VpLeftRound(Real *y, unsigned short f, ssize_t nf)
05805 /*
05806  * Round from the left hand side of the digits.
05807  */
05808 {
05809     BDIGIT v;
05810     if (!VpHasVal(y)) return 0; /* Unable to round */
05811     v = y->frac[0];
05812     nf -= VpExponent(y)*(ssize_t)BASE_FIG;
05813     while ((v /= 10) != 0) nf--;
05814     nf += (ssize_t)BASE_FIG-1;
05815     return VpMidRound(y,f,nf);
05816 }
05817 
05818 VP_EXPORT int
05819 VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t nf)
05820 {
05821     /* First,assign whole value in truncation mode */
05822     if (VpAsgn(y, x, 10) <= 1) return 0; /* Zero,NaN,or Infinity */
05823     return VpMidRound(y,f,nf);
05824 }
05825 
05826 static int
05827 VpLimitRound(Real *c, size_t ixDigit)
05828 {
05829     size_t ix = VpGetPrecLimit();
05830     if(!VpNmlz(c))    return -1;
05831     if(!ix)           return 0;
05832     if(!ixDigit) ixDigit = c->Prec-1;
05833     if((ix+BASE_FIG-1)/BASE_FIG > ixDigit+1) return 0;
05834     return VpLeftRound(c, VpGetRoundMode(), (ssize_t)ix);
05835 }
05836 
05837 /* If I understand correctly, this is only ever used to round off the final decimal
05838    digit of precision */
05839 static void
05840 VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v)
05841 {
05842     int f = 0;
05843 
05844     unsigned short const rounding_mode = VpGetRoundMode();
05845 
05846     if (VpLimitRound(c, ixDigit)) return;
05847     if (!v) return;
05848 
05849     v /= BASE1;
05850     switch (rounding_mode) {
05851     case VP_ROUND_DOWN:
05852         break;
05853     case VP_ROUND_UP:
05854         if (v) f = 1;
05855         break;
05856     case VP_ROUND_HALF_UP:
05857         if (v >= 5) f = 1;
05858         break;
05859     case VP_ROUND_HALF_DOWN:
05860         /* this is ok - because this is the last digit of precision,
05861            the case where v == 5 and some further digits are nonzero
05862            will never occur */
05863         if (v >= 6) f = 1;
05864         break;
05865     case VP_ROUND_CEIL:
05866         if (v && (VpGetSign(c) > 0)) f = 1;
05867         break;
05868     case VP_ROUND_FLOOR:
05869         if (v && (VpGetSign(c) < 0)) f = 1;
05870         break;
05871     case VP_ROUND_HALF_EVEN:  /* Banker's rounding */
05872         /* as per VP_ROUND_HALF_DOWN, because this is the last digit of precision,
05873            there is no case to worry about where v == 5 and some further digits are nonzero */
05874         if (v > 5) f = 1;
05875         else if (v == 5 && vPrev % 2) f = 1;
05876         break;
05877     }
05878     if (f) {
05879         VpRdup(c, ixDigit);
05880         VpNmlz(c);
05881     }
05882 }
05883 
05884 /*
05885  *  Rounds up m(plus one to final digit of m).
05886  */
05887 static int
05888 VpRdup(Real *m, size_t ind_m)
05889 {
05890     BDIGIT carry;
05891 
05892     if (!ind_m) ind_m = m->Prec;
05893 
05894     carry = 1;
05895     while (carry > 0 && (ind_m--)) {
05896         m->frac[ind_m] += carry;
05897         if (m->frac[ind_m] >= BASE) m->frac[ind_m] -= BASE;
05898         else                        carry = 0;
05899     }
05900     if(carry > 0) {        /* Overflow,count exponent and set fraction part be 1  */
05901         if (!AddExponent(m, 1)) return 0;
05902         m->Prec = m->frac[0] = 1;
05903     } else {
05904         VpNmlz(m);
05905     }
05906     return 1;
05907 }
05908 
05909 /*
05910  *  y = x - fix(x)
05911  */
05912 VP_EXPORT void
05913 VpFrac(Real *y, Real *x)
05914 {
05915     size_t my, ind_y, ind_x;
05916 
05917     if(!VpHasVal(x)) {
05918         VpAsgn(y,x,1);
05919         goto Exit;
05920     }
05921 
05922     if (x->exponent > 0 && (size_t)x->exponent >= x->Prec) {
05923         VpSetZero(y,VpGetSign(x));
05924         goto Exit;
05925     }
05926     else if(x->exponent <= 0) {
05927         VpAsgn(y, x, 1);
05928         goto Exit;
05929     }
05930 
05931     /* satisfy: x->exponent > 0 */
05932 
05933     y->Prec = x->Prec - (size_t)x->exponent;
05934     y->Prec = Min(y->Prec, y->MaxPrec);
05935     y->exponent = 0;
05936     VpSetSign(y,VpGetSign(x));
05937     ind_y = 0;
05938     my = y->Prec;
05939     ind_x = x->exponent;
05940     while(ind_y < my) {
05941         y->frac[ind_y] = x->frac[ind_x];
05942         ++ind_y;
05943         ++ind_x;
05944     }
05945     VpNmlz(y);
05946 
05947 Exit:
05948 #ifdef BIGDECIMAL_DEBUG
05949     if(gfDebug) {
05950         VPrint(stdout, "VpFrac y=%\n", y);
05951         VPrint(stdout, "    x=%\n", x);
05952     }
05953 #endif /* BIGDECIMAL_DEBUG */
05954     return;
05955 }
05956 
05957 /*
05958  *   y = x ** n
05959  */
05960 VP_EXPORT int
05961 VpPower(Real *y, Real *x, SIGNED_VALUE n)
05962 {
05963     size_t s, ss;
05964     ssize_t sign;
05965     Real *w1 = NULL;
05966     Real *w2 = NULL;
05967 
05968     if(VpIsZero(x)) {
05969         if(n==0) {
05970            VpSetOne(y);
05971            goto Exit;
05972         }
05973         sign = VpGetSign(x);
05974         if(n<0) {
05975            n = -n;
05976            if(sign<0) sign = (n%2)?(-1):(1);
05977            VpSetInf (y,sign);
05978         } else {
05979            if(sign<0) sign = (n%2)?(-1):(1);
05980            VpSetZero(y,sign);
05981         }
05982         goto Exit;
05983     }
05984     if(VpIsNaN(x)) {
05985         VpSetNaN(y);
05986         goto Exit;
05987     }
05988     if(VpIsInf(x)) {
05989         if(n==0) {
05990             VpSetOne(y);
05991             goto Exit;
05992         }
05993         if(n>0) {
05994             VpSetInf(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1);
05995             goto Exit;
05996         }
05997         VpSetZero(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1);
05998         goto Exit;
05999     }
06000 
06001     if((x->exponent == 1) &&(x->Prec == 1) &&(x->frac[0] == 1)) {
06002         /* abs(x) = 1 */
06003         VpSetOne(y);
06004         if(VpGetSign(x) > 0) goto Exit;
06005         if((n % 2) == 0) goto Exit;
06006         VpSetSign(y, -1);
06007         goto Exit;
06008     }
06009 
06010     if(n > 0) sign = 1;
06011     else if(n < 0) {
06012         sign = -1;
06013         n = -n;
06014     } else {
06015         VpSetOne(y);
06016         goto Exit;
06017     }
06018 
06019     /* Allocate working variables  */
06020 
06021     w1 = VpAlloc((y->MaxPrec + 2) * BASE_FIG, "#0");
06022     w2 = VpAlloc((w1->MaxPrec * 2 + 1) * BASE_FIG, "#0");
06023     /* calculation start */
06024 
06025     VpAsgn(y, x, 1);
06026     --n;
06027     while(n > 0) {
06028         VpAsgn(w1, x, 1);
06029         s = 1;
06030         while (ss = s, (s += s) <= (size_t)n) {
06031             VpMult(w2, w1, w1);
06032             VpAsgn(w1, w2, 1);
06033         }
06034         n -= (SIGNED_VALUE)ss;
06035         VpMult(w2, y, w1);
06036         VpAsgn(y, w2, 1);
06037     }
06038     if(sign < 0) {
06039         VpDivd(w1, w2, VpConstOne, y);
06040         VpAsgn(y, w1, 1);
06041     }
06042 
06043 Exit:
06044 #ifdef BIGDECIMAL_DEBUG
06045     if(gfDebug) {
06046         VPrint(stdout, "VpPower y=%\n", y);
06047         VPrint(stdout, "VpPower x=%\n", x);
06048         printf("  n=%d\n", n);
06049     }
06050 #endif /* BIGDECIMAL_DEBUG */
06051     VpFree(w2);
06052     VpFree(w1);
06053     return 1;
06054 }
06055 
06056 #ifdef BIGDECIMAL_DEBUG
06057 int
06058 VpVarCheck(Real * v)
06059 /*
06060  * Checks the validity of the Real variable v.
06061  * [Input]
06062  *   v ... Real *, variable to be checked.
06063  * [Returns]
06064  *   0  ... correct v.
06065  *   other ... error
06066  */
06067 {
06068     size_t i;
06069 
06070     if(v->MaxPrec <= 0) {
06071         printf("ERROR(VpVarCheck): Illegal Max. Precision(=%"PRIuSIZE")\n",
06072             v->MaxPrec);
06073         return 1;
06074     }
06075     if((v->Prec <= 0) ||((v->Prec) >(v->MaxPrec))) {
06076         printf("ERROR(VpVarCheck): Illegal Precision(=%"PRIuSIZE")\n", v->Prec);
06077         printf("       Max. Prec.=%"PRIuSIZE"\n", v->MaxPrec);
06078         return 2;
06079     }
06080     for(i = 0; i < v->Prec; ++i) {
06081         if((v->frac[i] >= BASE)) {
06082             printf("ERROR(VpVarCheck): Illegal fraction\n");
06083             printf("       Frac[%"PRIuSIZE"]=%lu\n", i, v->frac[i]);
06084             printf("       Prec.   =%"PRIuSIZE"\n", v->Prec);
06085             printf("       Exp. =%"PRIdVALUE"\n", v->exponent);
06086             printf("       BASE =%lu\n", BASE);
06087             return 3;
06088         }
06089     }
06090     return 0;
06091 }
06092 #endif /* BIGDECIMAL_DEBUG */
06093