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