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