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