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