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