Ruby 1.9.2p136(2010-12-25revision30365)

- Main Page
- Modules
- Data Structures
- Files

- File List
- Globals

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