1/* 2 * 3 * Ruby BigDecimal(Variable decimal precision) extension library. 4 * 5 * Copyright(C) 2002 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp) 6 * 7 * You may distribute under the terms of either the GNU General Public 8 * License or the Artistic License, as specified in the README file 9 * of this BigDecimal distribution. 10 * 11 * NOTE: Change log in this source removed to reduce source code size. 12 * See rev. 1.25 if needed. 13 * 14 */ 15 16/* #define BIGDECIMAL_DEBUG 1 */ 17#ifdef BIGDECIMAL_DEBUG 18# define BIGDECIMAL_ENABLE_VPRINT 1 19#endif 20#include "bigdecimal.h" 21 22#ifndef BIGDECIMAL_DEBUG 23# define NDEBUG 24#endif 25#include <assert.h> 26 27#include <ctype.h> 28#include <stdio.h> 29#include <stdlib.h> 30#include <string.h> 31#include <errno.h> 32#include <math.h> 33#include "math.h" 34 35#ifdef HAVE_IEEEFP_H 36#include <ieeefp.h> 37#endif 38 39/* #define ENABLE_NUMERIC_STRING */ 40 41VALUE rb_cBigDecimal; 42VALUE rb_mBigMath; 43 44static ID id_BigDecimal_exception_mode; 45static ID id_BigDecimal_rounding_mode; 46static ID id_BigDecimal_precision_limit; 47 48static ID id_up; 49static ID id_down; 50static ID id_truncate; 51static ID id_half_up; 52static ID id_default; 53static ID id_half_down; 54static ID id_half_even; 55static ID id_banker; 56static ID id_ceiling; 57static ID id_ceil; 58static ID id_floor; 59static ID id_to_r; 60static ID id_eq; 61 62/* MACRO's to guard objects from GC by keeping them in stack */ 63#define ENTER(n) volatile VALUE RB_UNUSED_VAR(vStack[n]);int iStack=0 64#define PUSH(x) vStack[iStack++] = (VALUE)(x); 65#define SAVE(p) PUSH(p->obj); 66#define GUARD_OBJ(p,y) {p=y;SAVE(p);} 67 68#define BASE_FIG RMPD_COMPONENT_FIGURES 69#define BASE RMPD_BASE 70 71#define HALF_BASE (BASE/2) 72#define BASE1 (BASE/10) 73 74#ifndef DBLE_FIG 75#define DBLE_FIG (DBL_DIG+1) /* figure of double */ 76#endif 77 78#ifndef RBIGNUM_ZERO_P 79# define RBIGNUM_ZERO_P(x) (RBIGNUM_LEN(x) == 0 || \ 80 (RBIGNUM_DIGITS(x)[0] == 0 && \ 81 (RBIGNUM_LEN(x) == 1 || bigzero_p(x)))) 82#endif 83 84static inline int 85bigzero_p(VALUE x) 86{ 87 long i; 88 BDIGIT *ds = RBIGNUM_DIGITS(x); 89 90 for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) { 91 if (ds[i]) return 0; 92 } 93 return 1; 94} 95 96#ifndef RRATIONAL_ZERO_P 97# define RRATIONAL_ZERO_P(x) (FIXNUM_P(RRATIONAL(x)->num) && \ 98 FIX2LONG(RRATIONAL(x)->num) == 0) 99#endif 100 101#ifndef RRATIONAL_NEGATIVE_P 102# define RRATIONAL_NEGATIVE_P(x) RTEST(rb_funcall((x), '<', 1, INT2FIX(0))) 103#endif 104 105#ifndef DECIMAL_SIZE_OF_BITS 106#define DECIMAL_SIZE_OF_BITS(n) (((n) * 3010 + 9998) / 9999) 107/* an approximation of ceil(n * log10(2)), upto 65536 at least */ 108#endif 109 110#ifdef PRIsVALUE 111# define RB_OBJ_CLASSNAME(obj) rb_obj_class(obj) 112# define RB_OBJ_STRING(obj) (obj) 113#else 114# define PRIsVALUE "s" 115# define RB_OBJ_CLASSNAME(obj) rb_obj_classname(obj) 116# define RB_OBJ_STRING(obj) StringValueCStr(obj) 117#endif 118 119/* 120 * ================== Ruby Interface part ========================== 121 */ 122#define DoSomeOne(x,y,f) rb_num_coerce_bin(x,y,f) 123 124/* 125 * Returns the BigDecimal version number. 126 */ 127static VALUE 128BigDecimal_version(VALUE self) 129{ 130 /* 131 * 1.0.0: Ruby 1.8.0 132 * 1.0.1: Ruby 1.8.1 133 * 1.1.0: Ruby 1.9.3 134 */ 135 return rb_str_new2("1.1.0"); 136} 137 138/* 139 * VP routines used in BigDecimal part 140 */ 141static unsigned short VpGetException(void); 142static void VpSetException(unsigned short f); 143static void VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v); 144static int VpLimitRound(Real *c, size_t ixDigit); 145static Real *VpCopy(Real *pv, Real const* const x); 146 147/* 148 * **** BigDecimal part **** 149 */ 150 151static void 152BigDecimal_delete(void *pv) 153{ 154 VpFree(pv); 155} 156 157static size_t 158BigDecimal_memsize(const void *ptr) 159{ 160 const Real *pv = ptr; 161 return pv ? (sizeof(*pv) + pv->MaxPrec * sizeof(BDIGIT)) : 0; 162} 163 164static const rb_data_type_t BigDecimal_data_type = { 165 "BigDecimal", 166 {0, BigDecimal_delete, BigDecimal_memsize,}, 167}; 168 169static inline int 170is_kind_of_BigDecimal(VALUE const v) 171{ 172 return rb_typeddata_is_kind_of(v, &BigDecimal_data_type); 173} 174 175static VALUE 176ToValue(Real *p) 177{ 178 if(VpIsNaN(p)) { 179 VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'(Not a Number)",0); 180 } else if(VpIsPosInf(p)) { 181 VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0); 182 } else if(VpIsNegInf(p)) { 183 VpException(VP_EXCEPTION_INFINITY,"Computation results to '-Infinity'",0); 184 } 185 return p->obj; 186} 187 188NORETURN(static void cannot_be_coerced_into_BigDecimal(VALUE, VALUE)); 189 190static void 191cannot_be_coerced_into_BigDecimal(VALUE exc_class, VALUE v) 192{ 193 VALUE str; 194 195 if (rb_special_const_p(v)) { 196 str = rb_inspect(v); 197 } 198 else { 199 str = rb_class_name(rb_obj_class(v)); 200 } 201 202 str = rb_str_cat2(rb_str_dup(str), " can't be coerced into BigDecimal"); 203 rb_exc_raise(rb_exc_new3(exc_class, str)); 204} 205 206static VALUE BigDecimal_div2(int, VALUE*, VALUE); 207 208static Real* 209GetVpValueWithPrec(VALUE v, long prec, int must) 210{ 211 Real *pv; 212 VALUE num, bg, args[2]; 213 char szD[128]; 214 VALUE orig = Qundef; 215 216again: 217 switch(TYPE(v)) 218 { 219 case T_FLOAT: 220 if (prec < 0) goto unable_to_coerce_without_prec; 221 if (prec > DBL_DIG+1)goto SomeOneMayDoIt; 222 v = rb_funcall(v, id_to_r, 0); 223 goto again; 224 case T_RATIONAL: 225 if (prec < 0) goto unable_to_coerce_without_prec; 226 227 if (orig == Qundef ? (orig = v, 1) : orig != v) { 228 num = RRATIONAL(v)->num; 229 pv = GetVpValueWithPrec(num, -1, must); 230 if (pv == NULL) goto SomeOneMayDoIt; 231 232 args[0] = RRATIONAL(v)->den; 233 args[1] = LONG2NUM(prec); 234 v = BigDecimal_div2(2, args, ToValue(pv)); 235 goto again; 236 } 237 238 v = orig; 239 goto SomeOneMayDoIt; 240 241 case T_DATA: 242 if (is_kind_of_BigDecimal(v)) { 243 pv = DATA_PTR(v); 244 return pv; 245 } 246 else { 247 goto SomeOneMayDoIt; 248 } 249 break; 250 251 case T_FIXNUM: 252 sprintf(szD, "%ld", FIX2LONG(v)); 253 return VpCreateRbObject(VpBaseFig() * 2 + 1, szD); 254 255#ifdef ENABLE_NUMERIC_STRING 256 case T_STRING: 257 SafeStringValue(v); 258 return VpCreateRbObject(strlen(RSTRING_PTR(v)) + VpBaseFig() + 1, 259 RSTRING_PTR(v)); 260#endif /* ENABLE_NUMERIC_STRING */ 261 262 case T_BIGNUM: 263 bg = rb_big2str(v, 10); 264 return VpCreateRbObject(strlen(RSTRING_PTR(bg)) + VpBaseFig() + 1, 265 RSTRING_PTR(bg)); 266 default: 267 goto SomeOneMayDoIt; 268 } 269 270SomeOneMayDoIt: 271 if (must) { 272 cannot_be_coerced_into_BigDecimal(rb_eTypeError, v); 273 } 274 return NULL; /* NULL means to coerce */ 275 276unable_to_coerce_without_prec: 277 if (must) { 278 rb_raise(rb_eArgError, 279 "%"PRIsVALUE" can't be coerced into BigDecimal without a precision", 280 RB_OBJ_CLASSNAME(v)); 281 } 282 return NULL; 283} 284 285static Real* 286GetVpValue(VALUE v, int must) 287{ 288 return GetVpValueWithPrec(v, -1, must); 289} 290 291/* call-seq: 292 * BigDecimal.double_fig 293 * 294 * The BigDecimal.double_fig class method returns the number of digits a 295 * Float number is allowed to have. The result depends upon the CPU and OS 296 * in use. 297 */ 298static VALUE 299BigDecimal_double_fig(VALUE self) 300{ 301 return INT2FIX(VpDblFig()); 302} 303 304/* call-seq: 305 * precs 306 * 307 * Returns an Array of two Integer values. 308 * 309 * The first value is the current number of significant digits in the 310 * BigDecimal. The second value is the maximum number of significant digits 311 * for the BigDecimal. 312 */ 313static VALUE 314BigDecimal_prec(VALUE self) 315{ 316 ENTER(1); 317 Real *p; 318 VALUE obj; 319 320 GUARD_OBJ(p,GetVpValue(self,1)); 321 obj = rb_assoc_new(INT2NUM(p->Prec*VpBaseFig()), 322 INT2NUM(p->MaxPrec*VpBaseFig())); 323 return obj; 324} 325 326/* 327 * call-seq: hash 328 * 329 * Creates a hash for this BigDecimal. 330 * 331 * Two BigDecimals with equal sign, 332 * fractional part and exponent have the same hash. 333 */ 334static VALUE 335BigDecimal_hash(VALUE self) 336{ 337 ENTER(1); 338 Real *p; 339 st_index_t hash; 340 341 GUARD_OBJ(p,GetVpValue(self,1)); 342 hash = (st_index_t)p->sign; 343 /* hash!=2: the case for 0(1),NaN(0) or +-Infinity(3) is sign itself */ 344 if(hash == 2 || hash == (st_index_t)-2) { 345 hash ^= rb_memhash(p->frac, sizeof(BDIGIT)*p->Prec); 346 hash += p->exponent; 347 } 348 return INT2FIX(hash); 349} 350 351/* 352 * call-seq: _dump 353 * 354 * Method used to provide marshalling support. 355 * 356 * inf = BigDecimal.new('Infinity') 357 * => #<BigDecimal:1e16fa8,'Infinity',9(9)> 358 * BigDecimal._load(inf._dump) 359 * => #<BigDecimal:1df8dc8,'Infinity',9(9)> 360 * 361 * See the Marshal module. 362 */ 363static VALUE 364BigDecimal_dump(int argc, VALUE *argv, VALUE self) 365{ 366 ENTER(5); 367 Real *vp; 368 char *psz; 369 VALUE dummy; 370 volatile VALUE dump; 371 372 rb_scan_args(argc, argv, "01", &dummy); 373 GUARD_OBJ(vp,GetVpValue(self,1)); 374 dump = rb_str_new(0,VpNumOfChars(vp,"E")+50); 375 psz = RSTRING_PTR(dump); 376 sprintf(psz, "%"PRIuSIZE":", VpMaxPrec(vp)*VpBaseFig()); 377 VpToString(vp, psz+strlen(psz), 0, 0); 378 rb_str_resize(dump, strlen(psz)); 379 return dump; 380} 381 382/* 383 * Internal method used to provide marshalling support. See the Marshal module. 384 */ 385static VALUE 386BigDecimal_load(VALUE self, VALUE str) 387{ 388 ENTER(2); 389 Real *pv; 390 unsigned char *pch; 391 unsigned char ch; 392 unsigned long m=0; 393 394 SafeStringValue(str); 395 pch = (unsigned char *)RSTRING_PTR(str); 396 /* First get max prec */ 397 while((*pch)!=(unsigned char)'\0' && (ch=*pch++)!=(unsigned char)':') { 398 if(!ISDIGIT(ch)) { 399 rb_raise(rb_eTypeError, "load failed: invalid character in the marshaled string"); 400 } 401 m = m*10 + (unsigned long)(ch-'0'); 402 } 403 if(m>VpBaseFig()) m -= VpBaseFig(); 404 GUARD_OBJ(pv,VpNewRbClass(m,(char *)pch,self)); 405 m /= VpBaseFig(); 406 if(m && pv->MaxPrec>m) pv->MaxPrec = m+1; 407 return ToValue(pv); 408} 409 410static unsigned short 411check_rounding_mode(VALUE const v) 412{ 413 unsigned short sw; 414 ID id; 415 switch (TYPE(v)) { 416 case T_SYMBOL: 417 id = SYM2ID(v); 418 if (id == id_up) 419 return VP_ROUND_UP; 420 if (id == id_down || id == id_truncate) 421 return VP_ROUND_DOWN; 422 if (id == id_half_up || id == id_default) 423 return VP_ROUND_HALF_UP; 424 if (id == id_half_down) 425 return VP_ROUND_HALF_DOWN; 426 if (id == id_half_even || id == id_banker) 427 return VP_ROUND_HALF_EVEN; 428 if (id == id_ceiling || id == id_ceil) 429 return VP_ROUND_CEIL; 430 if (id == id_floor) 431 return VP_ROUND_FLOOR; 432 rb_raise(rb_eArgError, "invalid rounding mode"); 433 434 default: 435 break; 436 } 437 438 Check_Type(v, T_FIXNUM); 439 sw = (unsigned short)FIX2UINT(v); 440 if (!VpIsRoundMode(sw)) { 441 rb_raise(rb_eArgError, "invalid rounding mode"); 442 } 443 return sw; 444} 445 446/* call-seq: 447 * BigDecimal.mode(mode, value) 448 * 449 * Controls handling of arithmetic exceptions and rounding. If no value 450 * is supplied, the current value is returned. 451 * 452 * Six values of the mode parameter control the handling of arithmetic 453 * exceptions: 454 * 455 * BigDecimal::EXCEPTION_NaN 456 * BigDecimal::EXCEPTION_INFINITY 457 * BigDecimal::EXCEPTION_UNDERFLOW 458 * BigDecimal::EXCEPTION_OVERFLOW 459 * BigDecimal::EXCEPTION_ZERODIVIDE 460 * BigDecimal::EXCEPTION_ALL 461 * 462 * For each mode parameter above, if the value set is false, computation 463 * continues after an arithmetic exception of the appropriate type. 464 * When computation continues, results are as follows: 465 * 466 * EXCEPTION_NaN:: NaN 467 * EXCEPTION_INFINITY:: +infinity or -infinity 468 * EXCEPTION_UNDERFLOW:: 0 469 * EXCEPTION_OVERFLOW:: +infinity or -infinity 470 * EXCEPTION_ZERODIVIDE:: +infinity or -infinity 471 * 472 * One value of the mode parameter controls the rounding of numeric values: 473 * BigDecimal::ROUND_MODE. The values it can take are: 474 * 475 * ROUND_UP, :up:: round away from zero 476 * ROUND_DOWN, :down, :truncate:: round towards zero (truncate) 477 * ROUND_HALF_UP, :half_up, :default:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round away from zero. (default) 478 * ROUND_HALF_DOWN, :half_down:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards zero. 479 * 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) 480 * ROUND_CEILING, :ceiling, :ceil:: round towards positive infinity (ceil) 481 * ROUND_FLOOR, :floor:: round towards negative infinity (floor) 482 * 483 */ 484static VALUE 485BigDecimal_mode(int argc, VALUE *argv, VALUE self) 486{ 487 VALUE which; 488 VALUE val; 489 unsigned long f,fo; 490 491 if(rb_scan_args(argc,argv,"11",&which,&val)==1) val = Qnil; 492 493 Check_Type(which, T_FIXNUM); 494 f = (unsigned long)FIX2INT(which); 495 496 if(f&VP_EXCEPTION_ALL) { 497 /* Exception mode setting */ 498 fo = VpGetException(); 499 if(val==Qnil) return INT2FIX(fo); 500 if(val!=Qfalse && val!=Qtrue) { 501 rb_raise(rb_eArgError, "second argument must be true or false"); 502 return Qnil; /* Not reached */ 503 } 504 if(f&VP_EXCEPTION_INFINITY) { 505 VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_INFINITY): 506 (fo&(~VP_EXCEPTION_INFINITY)))); 507 } 508 fo = VpGetException(); 509 if(f&VP_EXCEPTION_NaN) { 510 VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_NaN): 511 (fo&(~VP_EXCEPTION_NaN)))); 512 } 513 fo = VpGetException(); 514 if(f&VP_EXCEPTION_UNDERFLOW) { 515 VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_UNDERFLOW): 516 (fo&(~VP_EXCEPTION_UNDERFLOW)))); 517 } 518 fo = VpGetException(); 519 if(f&VP_EXCEPTION_ZERODIVIDE) { 520 VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_ZERODIVIDE): 521 (fo&(~VP_EXCEPTION_ZERODIVIDE)))); 522 } 523 fo = VpGetException(); 524 return INT2FIX(fo); 525 } 526 if (VP_ROUND_MODE == f) { 527 /* Rounding mode setting */ 528 unsigned short sw; 529 fo = VpGetRoundMode(); 530 if (NIL_P(val)) return INT2FIX(fo); 531 sw = check_rounding_mode(val); 532 fo = VpSetRoundMode(sw); 533 return INT2FIX(fo); 534 } 535 rb_raise(rb_eTypeError, "first argument for BigDecimal#mode invalid"); 536 return Qnil; 537} 538 539static size_t 540GetAddSubPrec(Real *a, Real *b) 541{ 542 size_t mxs; 543 size_t mx = a->Prec; 544 SIGNED_VALUE d; 545 546 if(!VpIsDef(a) || !VpIsDef(b)) return (size_t)-1L; 547 if(mx < b->Prec) mx = b->Prec; 548 if(a->exponent!=b->exponent) { 549 mxs = mx; 550 d = a->exponent - b->exponent; 551 if (d < 0) d = -d; 552 mx = mx + (size_t)d; 553 if (mx<mxs) { 554 return VpException(VP_EXCEPTION_INFINITY,"Exponent overflow",0); 555 } 556 } 557 return mx; 558} 559 560static SIGNED_VALUE 561GetPositiveInt(VALUE v) 562{ 563 SIGNED_VALUE n; 564 Check_Type(v, T_FIXNUM); 565 n = FIX2INT(v); 566 if (n < 0) { 567 rb_raise(rb_eArgError, "argument must be positive"); 568 } 569 return n; 570} 571 572VP_EXPORT Real * 573VpNewRbClass(size_t mx, const char *str, VALUE klass) 574{ 575 Real *pv = VpAlloc(mx,str); 576 pv->obj = TypedData_Wrap_Struct(klass, &BigDecimal_data_type, pv); 577 return pv; 578} 579 580VP_EXPORT Real * 581VpCreateRbObject(size_t mx, const char *str) 582{ 583 Real *pv = VpAlloc(mx,str); 584 pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv); 585 return pv; 586} 587 588#define VpAllocReal(prec) (Real *)VpMemAlloc(offsetof(Real, frac) + (prec) * sizeof(BDIGIT)) 589#define VpReallocReal(ptr, prec) (Real *)VpMemRealloc((ptr), offsetof(Real, frac) + (prec) * sizeof(BDIGIT)) 590 591static Real * 592VpCopy(Real *pv, Real const* const x) 593{ 594 assert(x != NULL); 595 596 pv = VpReallocReal(pv, x->MaxPrec); 597 pv->MaxPrec = x->MaxPrec; 598 pv->Prec = x->Prec; 599 pv->exponent = x->exponent; 600 pv->sign = x->sign; 601 pv->flag = x->flag; 602 MEMCPY(pv->frac, x->frac, BDIGIT, pv->MaxPrec); 603 604 return pv; 605} 606 607/* Returns True if the value is Not a Number */ 608static VALUE 609BigDecimal_IsNaN(VALUE self) 610{ 611 Real *p = GetVpValue(self,1); 612 if(VpIsNaN(p)) return Qtrue; 613 return Qfalse; 614} 615 616/* Returns nil, -1, or +1 depending on whether the value is finite, 617 * -infinity, or +infinity. 618 */ 619static VALUE 620BigDecimal_IsInfinite(VALUE self) 621{ 622 Real *p = GetVpValue(self,1); 623 if(VpIsPosInf(p)) return INT2FIX(1); 624 if(VpIsNegInf(p)) return INT2FIX(-1); 625 return Qnil; 626} 627 628/* Returns True if the value is finite (not NaN or infinite) */ 629static VALUE 630BigDecimal_IsFinite(VALUE self) 631{ 632 Real *p = GetVpValue(self,1); 633 if(VpIsNaN(p)) return Qfalse; 634 if(VpIsInf(p)) return Qfalse; 635 return Qtrue; 636} 637 638static void 639BigDecimal_check_num(Real *p) 640{ 641 if(VpIsNaN(p)) { 642 VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'(Not a Number)",1); 643 } else if(VpIsPosInf(p)) { 644 VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",1); 645 } else if(VpIsNegInf(p)) { 646 VpException(VP_EXCEPTION_INFINITY,"Computation results to '-Infinity'",1); 647 } 648} 649 650static VALUE BigDecimal_split(VALUE self); 651 652/* Returns the value as an integer (Fixnum or Bignum). 653 * 654 * If the BigNumber is infinity or NaN, raises FloatDomainError. 655 */ 656static VALUE 657BigDecimal_to_i(VALUE self) 658{ 659 ENTER(5); 660 ssize_t e, nf; 661 Real *p; 662 663 GUARD_OBJ(p,GetVpValue(self,1)); 664 BigDecimal_check_num(p); 665 666 e = VpExponent10(p); 667 if(e<=0) return INT2FIX(0); 668 nf = VpBaseFig(); 669 if(e<=nf) { 670 return LONG2NUM((long)(VpGetSign(p)*(BDIGIT_DBL_SIGNED)p->frac[0])); 671 } 672 else { 673 VALUE a = BigDecimal_split(self); 674 VALUE digits = RARRAY_PTR(a)[1]; 675 VALUE numerator = rb_funcall(digits, rb_intern("to_i"), 0); 676 VALUE ret; 677 ssize_t dpower = e - (ssize_t)RSTRING_LEN(digits); 678 679 if (VpGetSign(p) < 0) { 680 numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1)); 681 } 682 if (dpower < 0) { 683 ret = rb_funcall(numerator, rb_intern("div"), 1, 684 rb_funcall(INT2FIX(10), rb_intern("**"), 1, 685 INT2FIX(-dpower))); 686 } 687 else 688 ret = rb_funcall(numerator, '*', 1, 689 rb_funcall(INT2FIX(10), rb_intern("**"), 1, 690 INT2FIX(dpower))); 691 if (RB_TYPE_P(ret, T_FLOAT)) 692 rb_raise(rb_eFloatDomainError, "Infinity"); 693 return ret; 694 } 695} 696 697/* Returns a new Float object having approximately the same value as the 698 * BigDecimal number. Normal accuracy limits and built-in errors of binary 699 * Float arithmetic apply. 700 */ 701static VALUE 702BigDecimal_to_f(VALUE self) 703{ 704 ENTER(1); 705 Real *p; 706 double d; 707 SIGNED_VALUE e; 708 char *buf; 709 volatile VALUE str; 710 711 GUARD_OBJ(p, GetVpValue(self, 1)); 712 if (VpVtoD(&d, &e, p) != 1) 713 return rb_float_new(d); 714 if (e > (SIGNED_VALUE)(DBL_MAX_10_EXP+BASE_FIG)) 715 goto overflow; 716 if (e < (SIGNED_VALUE)(DBL_MIN_10_EXP-BASE_FIG)) 717 goto underflow; 718 719 str = rb_str_new(0, VpNumOfChars(p,"E")); 720 buf = RSTRING_PTR(str); 721 VpToString(p, buf, 0, 0); 722 errno = 0; 723 d = strtod(buf, 0); 724 if (errno == ERANGE) { 725 if (d == 0.0) goto underflow; 726 if (fabs(d) >= HUGE_VAL) goto overflow; 727 } 728 return rb_float_new(d); 729 730overflow: 731 VpException(VP_EXCEPTION_OVERFLOW, "BigDecimal to Float conversion", 0); 732 if (p->sign >= 0) 733 return rb_float_new(VpGetDoublePosInf()); 734 else 735 return rb_float_new(VpGetDoubleNegInf()); 736 737underflow: 738 VpException(VP_EXCEPTION_UNDERFLOW, "BigDecimal to Float conversion", 0); 739 if (p->sign >= 0) 740 return rb_float_new(0.0); 741 else 742 return rb_float_new(-0.0); 743} 744 745 746/* Converts a BigDecimal to a Rational. 747 */ 748static VALUE 749BigDecimal_to_r(VALUE self) 750{ 751 Real *p; 752 ssize_t sign, power, denomi_power; 753 VALUE a, digits, numerator; 754 755 p = GetVpValue(self,1); 756 BigDecimal_check_num(p); 757 758 sign = VpGetSign(p); 759 power = VpExponent10(p); 760 a = BigDecimal_split(self); 761 digits = RARRAY_PTR(a)[1]; 762 denomi_power = power - RSTRING_LEN(digits); 763 numerator = rb_funcall(digits, rb_intern("to_i"), 0); 764 765 if (sign < 0) { 766 numerator = rb_funcall(numerator, '*', 1, INT2FIX(-1)); 767 } 768 if (denomi_power < 0) { 769 return rb_Rational(numerator, 770 rb_funcall(INT2FIX(10), rb_intern("**"), 1, 771 INT2FIX(-denomi_power))); 772 } 773 else { 774 return rb_Rational1(rb_funcall(numerator, '*', 1, 775 rb_funcall(INT2FIX(10), rb_intern("**"), 1, 776 INT2FIX(denomi_power)))); 777 } 778} 779 780/* The coerce method provides support for Ruby type coercion. It is not 781 * enabled by default. 782 * 783 * This means that binary operations like + * / or - can often be performed 784 * on a BigDecimal and an object of another type, if the other object can 785 * be coerced into a BigDecimal value. 786 * 787 * e.g. 788 * a = BigDecimal.new("1.0") 789 * b = a / 2.0 -> 0.5 790 * 791 * Note that coercing a String to a BigDecimal is not supported by default; 792 * it requires a special compile-time option when building Ruby. 793 */ 794static VALUE 795BigDecimal_coerce(VALUE self, VALUE other) 796{ 797 ENTER(2); 798 VALUE obj; 799 Real *b; 800 801 if (RB_TYPE_P(other, T_FLOAT)) { 802 obj = rb_assoc_new(other, BigDecimal_to_f(self)); 803 } 804 else { 805 if (RB_TYPE_P(other, T_RATIONAL)) { 806 Real* pv = DATA_PTR(self); 807 GUARD_OBJ(b, GetVpValueWithPrec(other, pv->Prec*VpBaseFig(), 1)); 808 } 809 else { 810 GUARD_OBJ(b, GetVpValue(other, 1)); 811 } 812 obj = rb_assoc_new(b->obj, self); 813 } 814 815 return obj; 816} 817 818/* 819 * call-seq: +@ 820 * 821 * Return self. 822 * 823 * e.g. 824 * b = +a # b == a 825 */ 826static VALUE 827BigDecimal_uplus(VALUE self) 828{ 829 return self; 830} 831 832 /* 833 * Document-method: BigDecimal#add 834 * Document-method: BigDecimal#+ 835 * 836 * call-seq: 837 * add(value, digits) 838 * 839 * Add the specified value. 840 * 841 * e.g. 842 * c = a.add(b,n) 843 * c = a + b 844 * 845 * digits:: If specified and less than the number of significant digits of the 846 * result, the result is rounded to that number of digits, according to 847 * BigDecimal.mode. 848 */ 849static VALUE 850BigDecimal_add(VALUE self, VALUE r) 851{ 852 ENTER(5); 853 Real *c, *a, *b; 854 size_t mx; 855 856 GUARD_OBJ(a, GetVpValue(self, 1)); 857 if (RB_TYPE_P(r, T_FLOAT)) { 858 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 859 } 860 else if (RB_TYPE_P(r, T_RATIONAL)) { 861 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 862 } 863 else { 864 b = GetVpValue(r,0); 865 } 866 867 if (!b) return DoSomeOne(self,r,'+'); 868 SAVE(b); 869 870 if (VpIsNaN(b)) return b->obj; 871 if (VpIsNaN(a)) return a->obj; 872 873 mx = GetAddSubPrec(a, b); 874 if (mx == (size_t)-1L) { 875 GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0")); 876 VpAddSub(c, a, b, 1); 877 } 878 else { 879 GUARD_OBJ(c, VpCreateRbObject(mx * (VpBaseFig() + 1), "0")); 880 if(!mx) { 881 VpSetInf(c, VpGetSign(a)); 882 } 883 else { 884 VpAddSub(c, a, b, 1); 885 } 886 } 887 return ToValue(c); 888} 889 890 /* call-seq: 891 * sub(value, digits) 892 * 893 * Subtract the specified value. 894 * 895 * e.g. 896 * c = a.sub(b,n) 897 * c = a - b 898 * 899 * digits:: If specified and less than the number of significant digits of the 900 * result, the result is rounded to that number of digits, according to 901 * BigDecimal.mode. 902 */ 903static VALUE 904BigDecimal_sub(VALUE self, VALUE r) 905{ 906 ENTER(5); 907 Real *c, *a, *b; 908 size_t mx; 909 910 GUARD_OBJ(a,GetVpValue(self,1)); 911 if (RB_TYPE_P(r, T_FLOAT)) { 912 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 913 } 914 else if (RB_TYPE_P(r, T_RATIONAL)) { 915 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 916 } 917 else { 918 b = GetVpValue(r,0); 919 } 920 921 if(!b) return DoSomeOne(self,r,'-'); 922 SAVE(b); 923 924 if(VpIsNaN(b)) return b->obj; 925 if(VpIsNaN(a)) return a->obj; 926 927 mx = GetAddSubPrec(a,b); 928 if (mx == (size_t)-1L) { 929 GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0")); 930 VpAddSub(c, a, b, -1); 931 } else { 932 GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0")); 933 if(!mx) { 934 VpSetInf(c,VpGetSign(a)); 935 } else { 936 VpAddSub(c, a, b, -1); 937 } 938 } 939 return ToValue(c); 940} 941 942static VALUE 943BigDecimalCmp(VALUE self, VALUE r,char op) 944{ 945 ENTER(5); 946 SIGNED_VALUE e; 947 Real *a, *b=0; 948 GUARD_OBJ(a,GetVpValue(self,1)); 949 switch (TYPE(r)) { 950 case T_DATA: 951 if (!is_kind_of_BigDecimal(r)) break; 952 /* fall through */ 953 case T_FIXNUM: 954 /* fall through */ 955 case T_BIGNUM: 956 GUARD_OBJ(b, GetVpValue(r,0)); 957 break; 958 959 case T_FLOAT: 960 GUARD_OBJ(b, GetVpValueWithPrec(r, DBL_DIG+1, 0)); 961 break; 962 963 case T_RATIONAL: 964 GUARD_OBJ(b, GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 0)); 965 break; 966 967 default: 968 break; 969 } 970 if (b == NULL) { 971 ID f = 0; 972 973 switch (op) { 974 case '*': 975 return rb_num_coerce_cmp(self, r, rb_intern("<=>")); 976 977 case '=': 978 return RTEST(rb_num_coerce_cmp(self, r, rb_intern("=="))) ? Qtrue : Qfalse; 979 980 case 'G': 981 f = rb_intern(">="); 982 break; 983 984 case 'L': 985 f = rb_intern("<="); 986 break; 987 988 case '>': 989 /* fall through */ 990 case '<': 991 f = (ID)op; 992 break; 993 994 default: 995 break; 996 } 997 return rb_num_coerce_relop(self, r, f); 998 } 999 SAVE(b); 1000 e = VpComp(a, b); 1001 if (e == 999) 1002 return (op == '*') ? Qnil : Qfalse; 1003 switch (op) { 1004 case '*': 1005 return INT2FIX(e); /* any op */ 1006 1007 case '=': 1008 if(e==0) return Qtrue; 1009 return Qfalse; 1010 1011 case 'G': 1012 if(e>=0) return Qtrue; 1013 return Qfalse; 1014 1015 case '>': 1016 if(e> 0) return Qtrue; 1017 return Qfalse; 1018 1019 case 'L': 1020 if(e<=0) return Qtrue; 1021 return Qfalse; 1022 1023 case '<': 1024 if(e< 0) return Qtrue; 1025 return Qfalse; 1026 1027 default: 1028 break; 1029 } 1030 1031 rb_bug("Undefined operation in BigDecimalCmp()"); 1032 1033 UNREACHABLE; 1034} 1035 1036/* Returns True if the value is zero. */ 1037static VALUE 1038BigDecimal_zero(VALUE self) 1039{ 1040 Real *a = GetVpValue(self,1); 1041 return VpIsZero(a) ? Qtrue : Qfalse; 1042} 1043 1044/* Returns self if the value is non-zero, nil otherwise. */ 1045static VALUE 1046BigDecimal_nonzero(VALUE self) 1047{ 1048 Real *a = GetVpValue(self,1); 1049 return VpIsZero(a) ? Qnil : self; 1050} 1051 1052/* The comparison operator. 1053 * a <=> b is 0 if a == b, 1 if a > b, -1 if a < b. 1054 */ 1055static VALUE 1056BigDecimal_comp(VALUE self, VALUE r) 1057{ 1058 return BigDecimalCmp(self, r, '*'); 1059} 1060 1061/* 1062 * Tests for value equality; returns true if the values are equal. 1063 * 1064 * The == and === operators and the eql? method have the same implementation 1065 * for BigDecimal. 1066 * 1067 * Values may be coerced to perform the comparison: 1068 * 1069 * BigDecimal.new('1.0') == 1.0 -> true 1070 */ 1071static VALUE 1072BigDecimal_eq(VALUE self, VALUE r) 1073{ 1074 return BigDecimalCmp(self, r, '='); 1075} 1076 1077/* call-seq: 1078 * a < b 1079 * 1080 * Returns true if a is less than b. 1081 * 1082 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce). 1083 */ 1084static VALUE 1085BigDecimal_lt(VALUE self, VALUE r) 1086{ 1087 return BigDecimalCmp(self, r, '<'); 1088} 1089 1090/* call-seq: 1091 * a <= b 1092 * 1093 * Returns true if a is less than or equal to b. 1094 * 1095 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce). 1096 */ 1097static VALUE 1098BigDecimal_le(VALUE self, VALUE r) 1099{ 1100 return BigDecimalCmp(self, r, 'L'); 1101} 1102 1103/* call-seq: 1104 * a > b 1105 * 1106 * Returns true if a is greater than b. 1107 * 1108 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce). 1109 */ 1110static VALUE 1111BigDecimal_gt(VALUE self, VALUE r) 1112{ 1113 return BigDecimalCmp(self, r, '>'); 1114} 1115 1116/* call-seq: 1117 * a >= b 1118 * 1119 * Returns true if a is greater than or equal to b. 1120 * 1121 * Values may be coerced to perform the comparison (see ==, BigDecimal#coerce) 1122 */ 1123static VALUE 1124BigDecimal_ge(VALUE self, VALUE r) 1125{ 1126 return BigDecimalCmp(self, r, 'G'); 1127} 1128 1129/* 1130 * call-seq: -@ 1131 * 1132 * Return the negation of self. 1133 * 1134 * e.g. 1135 * b = -a 1136 * b == a * -1 1137 */ 1138static VALUE 1139BigDecimal_neg(VALUE self) 1140{ 1141 ENTER(5); 1142 Real *c, *a; 1143 GUARD_OBJ(a,GetVpValue(self,1)); 1144 GUARD_OBJ(c,VpCreateRbObject(a->Prec *(VpBaseFig() + 1), "0")); 1145 VpAsgn(c, a, -1); 1146 return ToValue(c); 1147} 1148 1149 /* 1150 * Document-method: BigDecimal#mult 1151 * 1152 * call-seq: mult(value, digits) 1153 * 1154 * Multiply by the specified value. 1155 * 1156 * e.g. 1157 * c = a.mult(b,n) 1158 * c = a * b 1159 * 1160 * digits:: If specified and less than the number of significant digits of the 1161 * result, the result is rounded to that number of digits, according to 1162 * BigDecimal.mode. 1163 */ 1164static VALUE 1165BigDecimal_mult(VALUE self, VALUE r) 1166{ 1167 ENTER(5); 1168 Real *c, *a, *b; 1169 size_t mx; 1170 1171 GUARD_OBJ(a,GetVpValue(self,1)); 1172 if (RB_TYPE_P(r, T_FLOAT)) { 1173 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 1174 } 1175 else if (RB_TYPE_P(r, T_RATIONAL)) { 1176 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 1177 } 1178 else { 1179 b = GetVpValue(r,0); 1180 } 1181 1182 if(!b) return DoSomeOne(self,r,'*'); 1183 SAVE(b); 1184 1185 mx = a->Prec + b->Prec; 1186 GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0")); 1187 VpMult(c, a, b); 1188 return ToValue(c); 1189} 1190 1191static VALUE 1192BigDecimal_divide(Real **c, Real **res, Real **div, VALUE self, VALUE r) 1193/* For c = self.div(r): with round operation */ 1194{ 1195 ENTER(5); 1196 Real *a, *b; 1197 size_t mx; 1198 1199 GUARD_OBJ(a,GetVpValue(self,1)); 1200 if (RB_TYPE_P(r, T_FLOAT)) { 1201 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 1202 } 1203 else if (RB_TYPE_P(r, T_RATIONAL)) { 1204 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 1205 } 1206 else { 1207 b = GetVpValue(r,0); 1208 } 1209 1210 if(!b) return DoSomeOne(self,r,'/'); 1211 SAVE(b); 1212 1213 *div = b; 1214 mx = a->Prec + vabs(a->exponent); 1215 if(mx<b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent); 1216 mx =(mx + 1) * VpBaseFig(); 1217 GUARD_OBJ((*c),VpCreateRbObject(mx, "#0")); 1218 GUARD_OBJ((*res),VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0")); 1219 VpDivd(*c, *res, a, b); 1220 return (VALUE)0; 1221} 1222 1223 /* call-seq: 1224 * div(value, digits) 1225 * quo(value) 1226 * 1227 * Divide by the specified value. 1228 * 1229 * e.g. 1230 * c = a.div(b,n) 1231 * 1232 * digits:: If specified and less than the number of significant digits of the 1233 * result, the result is rounded to that number of digits, according to 1234 * BigDecimal.mode. 1235 * 1236 * If digits is 0, the result is the same as the / operator. If not, the 1237 * result is an integer BigDecimal, by analogy with Float#div. 1238 * 1239 * The alias quo is provided since <code>div(value, 0)</code> is the same as 1240 * computing the quotient; see BigDecimal#divmod. 1241 */ 1242static VALUE 1243BigDecimal_div(VALUE self, VALUE r) 1244/* For c = self/r: with round operation */ 1245{ 1246 ENTER(5); 1247 Real *c=NULL, *res=NULL, *div = NULL; 1248 r = BigDecimal_divide(&c, &res, &div, self, r); 1249 if(r!=(VALUE)0) return r; /* coerced by other */ 1250 SAVE(c);SAVE(res);SAVE(div); 1251 /* a/b = c + r/b */ 1252 /* c xxxxx 1253 r 00000yyyyy ==> (y/b)*BASE >= HALF_BASE 1254 */ 1255 /* Round */ 1256 if(VpHasVal(div)) { /* frac[0] must be zero for NaN,INF,Zero */ 1257 VpInternalRound(c, 0, c->frac[c->Prec-1], (BDIGIT)(VpBaseVal()*(BDIGIT_DBL)res->frac[0]/div->frac[0])); 1258 } 1259 return ToValue(c); 1260} 1261 1262/* 1263 * %: mod = a%b = a - (a.to_f/b).floor * b 1264 * div = (a.to_f/b).floor 1265 */ 1266static VALUE 1267BigDecimal_DoDivmod(VALUE self, VALUE r, Real **div, Real **mod) 1268{ 1269 ENTER(8); 1270 Real *c=NULL, *d=NULL, *res=NULL; 1271 Real *a, *b; 1272 size_t mx; 1273 1274 GUARD_OBJ(a,GetVpValue(self,1)); 1275 if (RB_TYPE_P(r, T_FLOAT)) { 1276 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 1277 } 1278 else if (RB_TYPE_P(r, T_RATIONAL)) { 1279 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 1280 } 1281 else { 1282 b = GetVpValue(r,0); 1283 } 1284 1285 if(!b) return Qfalse; 1286 SAVE(b); 1287 1288 if(VpIsNaN(a) || VpIsNaN(b)) goto NaN; 1289 if(VpIsInf(a) && VpIsInf(b)) goto NaN; 1290 if(VpIsZero(b)) { 1291 rb_raise(rb_eZeroDivError, "divided by 0"); 1292 } 1293 if(VpIsInf(a)) { 1294 GUARD_OBJ(d,VpCreateRbObject(1, "0")); 1295 VpSetInf(d, (SIGNED_VALUE)(VpGetSign(a) == VpGetSign(b) ? 1 : -1)); 1296 GUARD_OBJ(c,VpCreateRbObject(1, "NaN")); 1297 *div = d; 1298 *mod = c; 1299 return Qtrue; 1300 } 1301 if(VpIsInf(b)) { 1302 GUARD_OBJ(d,VpCreateRbObject(1, "0")); 1303 *div = d; 1304 *mod = a; 1305 return Qtrue; 1306 } 1307 if(VpIsZero(a)) { 1308 GUARD_OBJ(c,VpCreateRbObject(1, "0")); 1309 GUARD_OBJ(d,VpCreateRbObject(1, "0")); 1310 *div = d; 1311 *mod = c; 1312 return Qtrue; 1313 } 1314 1315 mx = a->Prec + vabs(a->exponent); 1316 if(mx<b->Prec + vabs(b->exponent)) mx = b->Prec + vabs(b->exponent); 1317 mx =(mx + 1) * VpBaseFig(); 1318 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1319 GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0")); 1320 VpDivd(c, res, a, b); 1321 mx = c->Prec *(VpBaseFig() + 1); 1322 GUARD_OBJ(d,VpCreateRbObject(mx, "0")); 1323 VpActiveRound(d,c,VP_ROUND_DOWN,0); 1324 VpMult(res,d,b); 1325 VpAddSub(c,a,res,-1); 1326 if(!VpIsZero(c) && (VpGetSign(a)*VpGetSign(b)<0)) { 1327 VpAddSub(res,d,VpOne(),-1); 1328 GUARD_OBJ(d,VpCreateRbObject(GetAddSubPrec(c, b)*(VpBaseFig() + 1), "0")); 1329 VpAddSub(d ,c,b, 1); 1330 *div = res; 1331 *mod = d; 1332 } else { 1333 *div = d; 1334 *mod = c; 1335 } 1336 return Qtrue; 1337 1338NaN: 1339 GUARD_OBJ(c,VpCreateRbObject(1, "NaN")); 1340 GUARD_OBJ(d,VpCreateRbObject(1, "NaN")); 1341 *div = d; 1342 *mod = c; 1343 return Qtrue; 1344} 1345 1346/* call-seq: 1347 * a % b 1348 * a.modulo(b) 1349 * 1350 * Returns the modulus from dividing by b. 1351 * 1352 * See BigDecimal#divmod. 1353 */ 1354static VALUE 1355BigDecimal_mod(VALUE self, VALUE r) /* %: a%b = a - (a.to_f/b).floor * b */ 1356{ 1357 ENTER(3); 1358 Real *div=NULL, *mod=NULL; 1359 1360 if(BigDecimal_DoDivmod(self,r,&div,&mod)) { 1361 SAVE(div); SAVE(mod); 1362 return ToValue(mod); 1363 } 1364 return DoSomeOne(self,r,'%'); 1365} 1366 1367static VALUE 1368BigDecimal_divremain(VALUE self, VALUE r, Real **dv, Real **rv) 1369{ 1370 ENTER(10); 1371 size_t mx; 1372 Real *a=NULL, *b=NULL, *c=NULL, *res=NULL, *d=NULL, *rr=NULL, *ff=NULL; 1373 Real *f=NULL; 1374 1375 GUARD_OBJ(a,GetVpValue(self,1)); 1376 if (RB_TYPE_P(r, T_FLOAT)) { 1377 b = GetVpValueWithPrec(r, DBL_DIG+1, 1); 1378 } 1379 else if (RB_TYPE_P(r, T_RATIONAL)) { 1380 b = GetVpValueWithPrec(r, a->Prec*VpBaseFig(), 1); 1381 } 1382 else { 1383 b = GetVpValue(r,0); 1384 } 1385 1386 if(!b) return DoSomeOne(self,r,rb_intern("remainder")); 1387 SAVE(b); 1388 1389 mx =(a->MaxPrec + b->MaxPrec) *VpBaseFig(); 1390 GUARD_OBJ(c ,VpCreateRbObject(mx, "0")); 1391 GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0")); 1392 GUARD_OBJ(rr ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0")); 1393 GUARD_OBJ(ff ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0")); 1394 1395 VpDivd(c, res, a, b); 1396 1397 mx = c->Prec *(VpBaseFig() + 1); 1398 1399 GUARD_OBJ(d,VpCreateRbObject(mx, "0")); 1400 GUARD_OBJ(f,VpCreateRbObject(mx, "0")); 1401 1402 VpActiveRound(d,c,VP_ROUND_DOWN,0); /* 0: round off */ 1403 1404 VpFrac(f, c); 1405 VpMult(rr,f,b); 1406 VpAddSub(ff,res,rr,1); 1407 1408 *dv = d; 1409 *rv = ff; 1410 return (VALUE)0; 1411} 1412 1413/* Returns the remainder from dividing by the value. 1414 * 1415 * x.remainder(y) means x-y*(x/y).truncate 1416 */ 1417static VALUE 1418BigDecimal_remainder(VALUE self, VALUE r) /* remainder */ 1419{ 1420 VALUE f; 1421 Real *d,*rv=0; 1422 f = BigDecimal_divremain(self,r,&d,&rv); 1423 if(f!=(VALUE)0) return f; 1424 return ToValue(rv); 1425} 1426 1427/* Divides by the specified value, and returns the quotient and modulus 1428 * as BigDecimal numbers. The quotient is rounded towards negative infinity. 1429 * 1430 * For example: 1431 * 1432 * require 'bigdecimal' 1433 * 1434 * a = BigDecimal.new("42") 1435 * b = BigDecimal.new("9") 1436 * 1437 * q,m = a.divmod(b) 1438 * 1439 * c = q * b + m 1440 * 1441 * a == c -> true 1442 * 1443 * The quotient q is (a/b).floor, and the modulus is the amount that must be 1444 * added to q * b to get a. 1445 */ 1446static VALUE 1447BigDecimal_divmod(VALUE self, VALUE r) 1448{ 1449 ENTER(5); 1450 Real *div=NULL, *mod=NULL; 1451 1452 if(BigDecimal_DoDivmod(self,r,&div,&mod)) { 1453 SAVE(div); SAVE(mod); 1454 return rb_assoc_new(ToValue(div), ToValue(mod)); 1455 } 1456 return DoSomeOne(self,r,rb_intern("divmod")); 1457} 1458 1459/* 1460 * See BigDecimal#quo 1461 */ 1462static VALUE 1463BigDecimal_div2(int argc, VALUE *argv, VALUE self) 1464{ 1465 ENTER(5); 1466 VALUE b,n; 1467 int na = rb_scan_args(argc,argv,"11",&b,&n); 1468 if(na==1) { /* div in Float sense */ 1469 Real *div=NULL; 1470 Real *mod; 1471 if(BigDecimal_DoDivmod(self,b,&div,&mod)) { 1472 return BigDecimal_to_i(ToValue(div)); 1473 } 1474 return DoSomeOne(self,b,rb_intern("div")); 1475 } else { /* div in BigDecimal sense */ 1476 SIGNED_VALUE ix = GetPositiveInt(n); 1477 if (ix == 0) return BigDecimal_div(self, b); 1478 else { 1479 Real *res=NULL; 1480 Real *av=NULL, *bv=NULL, *cv=NULL; 1481 size_t mx = (ix+VpBaseFig()*2); 1482 size_t pl = VpSetPrecLimit(0); 1483 1484 GUARD_OBJ(cv,VpCreateRbObject(mx,"0")); 1485 GUARD_OBJ(av,GetVpValue(self,1)); 1486 GUARD_OBJ(bv,GetVpValue(b,1)); 1487 mx = av->Prec + bv->Prec + 2; 1488 if(mx <= cv->MaxPrec) mx = cv->MaxPrec+1; 1489 GUARD_OBJ(res,VpCreateRbObject((mx * 2 + 2)*VpBaseFig(), "#0")); 1490 VpDivd(cv,res,av,bv); 1491 VpSetPrecLimit(pl); 1492 VpLeftRound(cv,VpGetRoundMode(),ix); 1493 return ToValue(cv); 1494 } 1495 } 1496} 1497 1498static VALUE 1499BigDecimal_add2(VALUE self, VALUE b, VALUE n) 1500{ 1501 ENTER(2); 1502 Real *cv; 1503 SIGNED_VALUE mx = GetPositiveInt(n); 1504 if (mx == 0) return BigDecimal_add(self, b); 1505 else { 1506 size_t pl = VpSetPrecLimit(0); 1507 VALUE c = BigDecimal_add(self,b); 1508 VpSetPrecLimit(pl); 1509 GUARD_OBJ(cv,GetVpValue(c,1)); 1510 VpLeftRound(cv,VpGetRoundMode(),mx); 1511 return ToValue(cv); 1512 } 1513} 1514 1515static VALUE 1516BigDecimal_sub2(VALUE self, VALUE b, VALUE n) 1517{ 1518 ENTER(2); 1519 Real *cv; 1520 SIGNED_VALUE mx = GetPositiveInt(n); 1521 if (mx == 0) return BigDecimal_sub(self, b); 1522 else { 1523 size_t pl = VpSetPrecLimit(0); 1524 VALUE c = BigDecimal_sub(self,b); 1525 VpSetPrecLimit(pl); 1526 GUARD_OBJ(cv,GetVpValue(c,1)); 1527 VpLeftRound(cv,VpGetRoundMode(),mx); 1528 return ToValue(cv); 1529 } 1530} 1531 1532static VALUE 1533BigDecimal_mult2(VALUE self, VALUE b, VALUE n) 1534{ 1535 ENTER(2); 1536 Real *cv; 1537 SIGNED_VALUE mx = GetPositiveInt(n); 1538 if (mx == 0) return BigDecimal_mult(self, b); 1539 else { 1540 size_t pl = VpSetPrecLimit(0); 1541 VALUE c = BigDecimal_mult(self,b); 1542 VpSetPrecLimit(pl); 1543 GUARD_OBJ(cv,GetVpValue(c,1)); 1544 VpLeftRound(cv,VpGetRoundMode(),mx); 1545 return ToValue(cv); 1546 } 1547} 1548 1549/* Returns the absolute value. 1550 * 1551 * BigDecimal('5').abs -> 5 1552 * 1553 * BigDecimal('-3').abs -> 3 1554 */ 1555static VALUE 1556BigDecimal_abs(VALUE self) 1557{ 1558 ENTER(5); 1559 Real *c, *a; 1560 size_t mx; 1561 1562 GUARD_OBJ(a,GetVpValue(self,1)); 1563 mx = a->Prec *(VpBaseFig() + 1); 1564 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1565 VpAsgn(c, a, 1); 1566 VpChangeSign(c, 1); 1567 return ToValue(c); 1568} 1569 1570/* call-seq: 1571 * sqrt(n) 1572 * 1573 * Returns the square root of the value. 1574 * 1575 * Result has at least n significant digits. 1576 */ 1577static VALUE 1578BigDecimal_sqrt(VALUE self, VALUE nFig) 1579{ 1580 ENTER(5); 1581 Real *c, *a; 1582 size_t mx, n; 1583 1584 GUARD_OBJ(a,GetVpValue(self,1)); 1585 mx = a->Prec *(VpBaseFig() + 1); 1586 1587 n = GetPositiveInt(nFig) + VpDblFig() + 1; 1588 if(mx <= n) mx = n; 1589 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1590 VpSqrt(c, a); 1591 return ToValue(c); 1592} 1593 1594/* Return the integer part of the number. 1595 */ 1596static VALUE 1597BigDecimal_fix(VALUE self) 1598{ 1599 ENTER(5); 1600 Real *c, *a; 1601 size_t mx; 1602 1603 GUARD_OBJ(a,GetVpValue(self,1)); 1604 mx = a->Prec *(VpBaseFig() + 1); 1605 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1606 VpActiveRound(c,a,VP_ROUND_DOWN,0); /* 0: round off */ 1607 return ToValue(c); 1608} 1609 1610/* call-seq: 1611 * round(n, mode) 1612 * 1613 * Round to the nearest 1 (by default), returning the result as a BigDecimal. 1614 * 1615 * BigDecimal('3.14159').round #=> 3 1616 * BigDecimal('8.7').round #=> 9 1617 * 1618 * If n is specified and positive, the fractional part of the result has no 1619 * more than that many digits. 1620 * 1621 * If n is specified and negative, at least that many digits to the left of the 1622 * decimal point will be 0 in the result. 1623 * 1624 * BigDecimal('3.14159').round(3) #=> 3.142 1625 * BigDecimal('13345.234').round(-2) #=> 13300.0 1626 * 1627 * The value of the optional mode argument can be used to determine how 1628 * rounding is performed; see BigDecimal.mode. 1629 */ 1630static VALUE 1631BigDecimal_round(int argc, VALUE *argv, VALUE self) 1632{ 1633 ENTER(5); 1634 Real *c, *a; 1635 int iLoc = 0; 1636 VALUE vLoc; 1637 VALUE vRound; 1638 size_t mx, pl; 1639 1640 unsigned short sw = VpGetRoundMode(); 1641 1642 switch (rb_scan_args(argc, argv, "02", &vLoc, &vRound)) { 1643 case 0: 1644 iLoc = 0; 1645 break; 1646 case 1: 1647 Check_Type(vLoc, T_FIXNUM); 1648 iLoc = FIX2INT(vLoc); 1649 break; 1650 case 2: 1651 Check_Type(vLoc, T_FIXNUM); 1652 iLoc = FIX2INT(vLoc); 1653 sw = check_rounding_mode(vRound); 1654 break; 1655 } 1656 1657 pl = VpSetPrecLimit(0); 1658 GUARD_OBJ(a,GetVpValue(self,1)); 1659 mx = a->Prec *(VpBaseFig() + 1); 1660 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1661 VpSetPrecLimit(pl); 1662 VpActiveRound(c,a,sw,iLoc); 1663 if (argc == 0) { 1664 return BigDecimal_to_i(ToValue(c)); 1665 } 1666 return ToValue(c); 1667} 1668 1669/* call-seq: 1670 * truncate(n) 1671 * 1672 * Truncate to the nearest 1, returning the result as a BigDecimal. 1673 * 1674 * BigDecimal('3.14159').truncate #=> 3 1675 * BigDecimal('8.7').truncate #=> 8 1676 * 1677 * If n is specified and positive, the fractional part of the result has no 1678 * more than that many digits. 1679 * 1680 * If n is specified and negative, at least that many digits to the left of the 1681 * decimal point will be 0 in the result. 1682 * 1683 * BigDecimal('3.14159').truncate(3) #=> 3.141 1684 * BigDecimal('13345.234').truncate(-2) #=> 13300.0 1685 */ 1686static VALUE 1687BigDecimal_truncate(int argc, VALUE *argv, VALUE self) 1688{ 1689 ENTER(5); 1690 Real *c, *a; 1691 int iLoc; 1692 VALUE vLoc; 1693 size_t mx, pl = VpSetPrecLimit(0); 1694 1695 if(rb_scan_args(argc,argv,"01",&vLoc)==0) { 1696 iLoc = 0; 1697 } else { 1698 Check_Type(vLoc, T_FIXNUM); 1699 iLoc = FIX2INT(vLoc); 1700 } 1701 1702 GUARD_OBJ(a,GetVpValue(self,1)); 1703 mx = a->Prec *(VpBaseFig() + 1); 1704 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1705 VpSetPrecLimit(pl); 1706 VpActiveRound(c,a,VP_ROUND_DOWN,iLoc); /* 0: truncate */ 1707 if (argc == 0) { 1708 return BigDecimal_to_i(ToValue(c)); 1709 } 1710 return ToValue(c); 1711} 1712 1713/* Return the fractional part of the number. 1714 */ 1715static VALUE 1716BigDecimal_frac(VALUE self) 1717{ 1718 ENTER(5); 1719 Real *c, *a; 1720 size_t mx; 1721 1722 GUARD_OBJ(a,GetVpValue(self,1)); 1723 mx = a->Prec *(VpBaseFig() + 1); 1724 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1725 VpFrac(c, a); 1726 return ToValue(c); 1727} 1728 1729/* call-seq: 1730 * floor(n) 1731 * 1732 * Return the largest integer less than or equal to the value, as a BigDecimal. 1733 * 1734 * BigDecimal('3.14159').floor #=> 3 1735 * BigDecimal('-9.1').floor #=> -10 1736 * 1737 * If n is specified and positive, the fractional part of the result has no 1738 * more than that many digits. 1739 * 1740 * If n is specified and negative, at least that 1741 * many digits to the left of the decimal point will be 0 in the result. 1742 * 1743 * BigDecimal('3.14159').floor(3) #=> 3.141 1744 * BigDecimal('13345.234').floor(-2) #=> 13300.0 1745 */ 1746static VALUE 1747BigDecimal_floor(int argc, VALUE *argv, VALUE self) 1748{ 1749 ENTER(5); 1750 Real *c, *a; 1751 int iLoc; 1752 VALUE vLoc; 1753 size_t mx, pl = VpSetPrecLimit(0); 1754 1755 if(rb_scan_args(argc,argv,"01",&vLoc)==0) { 1756 iLoc = 0; 1757 } else { 1758 Check_Type(vLoc, T_FIXNUM); 1759 iLoc = FIX2INT(vLoc); 1760 } 1761 1762 GUARD_OBJ(a,GetVpValue(self,1)); 1763 mx = a->Prec *(VpBaseFig() + 1); 1764 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1765 VpSetPrecLimit(pl); 1766 VpActiveRound(c,a,VP_ROUND_FLOOR,iLoc); 1767#ifdef BIGDECIMAL_DEBUG 1768 VPrint(stderr, "floor: c=%\n", c); 1769#endif 1770 if (argc == 0) { 1771 return BigDecimal_to_i(ToValue(c)); 1772 } 1773 return ToValue(c); 1774} 1775 1776/* call-seq: 1777 * ceil(n) 1778 * 1779 * Return the smallest integer greater than or equal to the value, as a BigDecimal. 1780 * 1781 * BigDecimal('3.14159').ceil #=> 4 1782 * BigDecimal('-9.1').ceil #=> -9 1783 * 1784 * If n is specified and positive, the fractional part of the result has no 1785 * more than that many digits. 1786 * 1787 * If n is specified and negative, at least that 1788 * many digits to the left of the decimal point will be 0 in the result. 1789 * 1790 * BigDecimal('3.14159').ceil(3) #=> 3.142 1791 * BigDecimal('13345.234').ceil(-2) #=> 13400.0 1792 */ 1793static VALUE 1794BigDecimal_ceil(int argc, VALUE *argv, VALUE self) 1795{ 1796 ENTER(5); 1797 Real *c, *a; 1798 int iLoc; 1799 VALUE vLoc; 1800 size_t mx, pl = VpSetPrecLimit(0); 1801 1802 if(rb_scan_args(argc,argv,"01",&vLoc)==0) { 1803 iLoc = 0; 1804 } else { 1805 Check_Type(vLoc, T_FIXNUM); 1806 iLoc = FIX2INT(vLoc); 1807 } 1808 1809 GUARD_OBJ(a,GetVpValue(self,1)); 1810 mx = a->Prec *(VpBaseFig() + 1); 1811 GUARD_OBJ(c,VpCreateRbObject(mx, "0")); 1812 VpSetPrecLimit(pl); 1813 VpActiveRound(c,a,VP_ROUND_CEIL,iLoc); 1814 if (argc == 0) { 1815 return BigDecimal_to_i(ToValue(c)); 1816 } 1817 return ToValue(c); 1818} 1819 1820/* call-seq: 1821 * to_s(s) 1822 * 1823 * Converts the value to a string. 1824 * 1825 * The default format looks like 0.xxxxEnn. 1826 * 1827 * The optional parameter s consists of either an integer; or an optional '+' 1828 * or ' ', followed by an optional number, followed by an optional 'E' or 'F'. 1829 * 1830 * If there is a '+' at the start of s, positive values are returned with 1831 * a leading '+'. 1832 * 1833 * A space at the start of s returns positive values with a leading space. 1834 * 1835 * If s contains a number, a space is inserted after each group of that many 1836 * fractional digits. 1837 * 1838 * If s ends with an 'E', engineering notation (0.xxxxEnn) is used. 1839 * 1840 * If s ends with an 'F', conventional floating point notation is used. 1841 * 1842 * Examples: 1843 * 1844 * BigDecimal.new('-123.45678901234567890').to_s('5F') 1845 * #=> '-123.45678 90123 45678 9' 1846 * 1847 * BigDecimal.new('123.45678901234567890').to_s('+8F') 1848 * #=> '+123.45678901 23456789' 1849 * 1850 * BigDecimal.new('123.45678901234567890').to_s(' F') 1851 * #=> ' 123.4567890123456789' 1852 */ 1853static VALUE 1854BigDecimal_to_s(int argc, VALUE *argv, VALUE self) 1855{ 1856 ENTER(5); 1857 int fmt = 0; /* 0:E format */ 1858 int fPlus = 0; /* =0:default,=1: set ' ' before digits ,set '+' before digits. */ 1859 Real *vp; 1860 volatile VALUE str; 1861 char *psz; 1862 char ch; 1863 size_t nc, mc = 0; 1864 VALUE f; 1865 1866 GUARD_OBJ(vp,GetVpValue(self,1)); 1867 1868 if (rb_scan_args(argc,argv,"01",&f)==1) { 1869 if (RB_TYPE_P(f, T_STRING)) { 1870 SafeStringValue(f); 1871 psz = RSTRING_PTR(f); 1872 if (*psz == ' ') { 1873 fPlus = 1; 1874 psz++; 1875 } 1876 else if (*psz == '+') { 1877 fPlus = 2; 1878 psz++; 1879 } 1880 while ((ch = *psz++) != 0) { 1881 if (ISSPACE(ch)) { 1882 continue; 1883 } 1884 if (!ISDIGIT(ch)) { 1885 if (ch == 'F' || ch == 'f') { 1886 fmt = 1; /* F format */ 1887 } 1888 break; 1889 } 1890 mc = mc * 10 + ch - '0'; 1891 } 1892 } 1893 else { 1894 mc = (size_t)GetPositiveInt(f); 1895 } 1896 } 1897 if (fmt) { 1898 nc = VpNumOfChars(vp, "F"); 1899 } 1900 else { 1901 nc = VpNumOfChars(vp, "E"); 1902 } 1903 if (mc > 0) { 1904 nc += (nc + mc - 1) / mc + 1; 1905 } 1906 1907 str = rb_str_new(0, nc); 1908 psz = RSTRING_PTR(str); 1909 1910 if (fmt) { 1911 VpToFString(vp, psz, mc, fPlus); 1912 } 1913 else { 1914 VpToString (vp, psz, mc, fPlus); 1915 } 1916 rb_str_resize(str, strlen(psz)); 1917 return str; 1918} 1919 1920/* Splits a BigDecimal number into four parts, returned as an array of values. 1921 * 1922 * The first value represents the sign of the BigDecimal, and is -1 or 1, or 0 1923 * if the BigDecimal is Not a Number. 1924 * 1925 * The second value is a string representing the significant digits of the 1926 * BigDecimal, with no leading zeros. 1927 * 1928 * The third value is the base used for arithmetic (currently always 10) as an 1929 * Integer. 1930 * 1931 * The fourth value is an Integer exponent. 1932 * 1933 * If the BigDecimal can be represented as 0.xxxxxx*10**n, then xxxxxx is the 1934 * string of significant digits with no leading zeros, and n is the exponent. 1935 * 1936 * From these values, you can translate a BigDecimal to a float as follows: 1937 * 1938 * sign, significant_digits, base, exponent = a.split 1939 * f = sign * "0.#{significant_digits}".to_f * (base ** exponent) 1940 * 1941 * (Note that the to_f method is provided as a more convenient way to translate 1942 * a BigDecimal to a Float.) 1943 */ 1944static VALUE 1945BigDecimal_split(VALUE self) 1946{ 1947 ENTER(5); 1948 Real *vp; 1949 VALUE obj,str; 1950 ssize_t e, s; 1951 char *psz1; 1952 1953 GUARD_OBJ(vp,GetVpValue(self,1)); 1954 str = rb_str_new(0, VpNumOfChars(vp,"E")); 1955 psz1 = RSTRING_PTR(str); 1956 VpSzMantissa(vp,psz1); 1957 s = 1; 1958 if(psz1[0]=='-') { 1959 size_t len = strlen(psz1+1); 1960 1961 memmove(psz1, psz1+1, len); 1962 psz1[len] = '\0'; 1963 s = -1; 1964 } 1965 if(psz1[0]=='N') s=0; /* NaN */ 1966 e = VpExponent10(vp); 1967 obj = rb_ary_new2(4); 1968 rb_ary_push(obj, INT2FIX(s)); 1969 rb_ary_push(obj, str); 1970 rb_str_resize(str, strlen(psz1)); 1971 rb_ary_push(obj, INT2FIX(10)); 1972 rb_ary_push(obj, INT2NUM(e)); 1973 return obj; 1974} 1975 1976/* Returns the exponent of the BigDecimal number, as an Integer. 1977 * 1978 * If the number can be represented as 0.xxxxxx*10**n where xxxxxx is a string 1979 * of digits with no leading zeros, then n is the exponent. 1980 */ 1981static VALUE 1982BigDecimal_exponent(VALUE self) 1983{ 1984 ssize_t e = VpExponent10(GetVpValue(self, 1)); 1985 return INT2NUM(e); 1986} 1987 1988/* Returns debugging information about the value as a string of comma-separated 1989 * values in angle brackets with a leading #: 1990 * 1991 * BigDecimal.new("1234.5678").inspect -> 1992 * "#<BigDecimal:b7ea1130,'0.12345678E4',8(12)>" 1993 * 1994 * The first part is the address, the second is the value as a string, and 1995 * the final part ss(mm) is the current number of significant digits and the 1996 * maximum number of significant digits, respectively. 1997 */ 1998static VALUE 1999BigDecimal_inspect(VALUE self) 2000{ 2001 ENTER(5); 2002 Real *vp; 2003 volatile VALUE obj; 2004 size_t nc; 2005 char *psz, *tmp; 2006 2007 GUARD_OBJ(vp,GetVpValue(self,1)); 2008 nc = VpNumOfChars(vp,"E"); 2009 nc +=(nc + 9) / 10; 2010 2011 obj = rb_str_new(0, nc+256); 2012 psz = RSTRING_PTR(obj); 2013 sprintf(psz,"#<BigDecimal:%"PRIxVALUE",'",self); 2014 tmp = psz + strlen(psz); 2015 VpToString(vp, tmp, 10, 0); 2016 tmp += strlen(tmp); 2017 sprintf(tmp, "',%"PRIuSIZE"(%"PRIuSIZE")>", VpPrec(vp)*VpBaseFig(), VpMaxPrec(vp)*VpBaseFig()); 2018 rb_str_resize(obj, strlen(psz)); 2019 return obj; 2020} 2021 2022static VALUE BigMath_s_exp(VALUE, VALUE, VALUE); 2023static VALUE BigMath_s_log(VALUE, VALUE, VALUE); 2024 2025#define BigMath_exp(x, n) BigMath_s_exp(rb_mBigMath, (x), (n)) 2026#define BigMath_log(x, n) BigMath_s_log(rb_mBigMath, (x), (n)) 2027 2028inline static int 2029is_integer(VALUE x) 2030{ 2031 return (RB_TYPE_P(x, T_FIXNUM) || RB_TYPE_P(x, T_BIGNUM)); 2032} 2033 2034inline static int 2035is_negative(VALUE x) 2036{ 2037 if (FIXNUM_P(x)) { 2038 return FIX2LONG(x) < 0; 2039 } 2040 else if (RB_TYPE_P(x, T_BIGNUM)) { 2041 return RBIGNUM_NEGATIVE_P(x); 2042 } 2043 else if (RB_TYPE_P(x, T_FLOAT)) { 2044 return RFLOAT_VALUE(x) < 0.0; 2045 } 2046 return RTEST(rb_funcall(x, '<', 1, INT2FIX(0))); 2047} 2048 2049#define is_positive(x) (!is_negative(x)) 2050 2051inline static int 2052is_zero(VALUE x) 2053{ 2054 VALUE num; 2055 2056 switch (TYPE(x)) { 2057 case T_FIXNUM: 2058 return FIX2LONG(x) == 0; 2059 2060 case T_BIGNUM: 2061 return Qfalse; 2062 2063 case T_RATIONAL: 2064 num = RRATIONAL(x)->num; 2065 return FIXNUM_P(num) && FIX2LONG(num) == 0; 2066 2067 default: 2068 break; 2069 } 2070 2071 return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(0))); 2072} 2073 2074inline static int 2075is_one(VALUE x) 2076{ 2077 VALUE num, den; 2078 2079 switch (TYPE(x)) { 2080 case T_FIXNUM: 2081 return FIX2LONG(x) == 1; 2082 2083 case T_BIGNUM: 2084 return Qfalse; 2085 2086 case T_RATIONAL: 2087 num = RRATIONAL(x)->num; 2088 den = RRATIONAL(x)->den; 2089 return FIXNUM_P(den) && FIX2LONG(den) == 1 && 2090 FIXNUM_P(num) && FIX2LONG(num) == 1; 2091 2092 default: 2093 break; 2094 } 2095 2096 return RTEST(rb_funcall(x, id_eq, 1, INT2FIX(1))); 2097} 2098 2099inline static int 2100is_even(VALUE x) 2101{ 2102 switch (TYPE(x)) { 2103 case T_FIXNUM: 2104 return (FIX2LONG(x) % 2) == 0; 2105 2106 case T_BIGNUM: 2107 return (RBIGNUM_DIGITS(x)[0] % 2) == 0; 2108 2109 default: 2110 break; 2111 } 2112 2113 return 0; 2114} 2115 2116static VALUE 2117rmpd_power_by_big_decimal(Real const* x, Real const* exp, ssize_t const n) 2118{ 2119 VALUE log_x, multiplied, y; 2120 volatile VALUE obj = exp->obj; 2121 2122 if (VpIsZero(exp)) { 2123 return ToValue(VpCreateRbObject(n, "1")); 2124 } 2125 2126 log_x = BigMath_log(x->obj, SSIZET2NUM(n+1)); 2127 multiplied = BigDecimal_mult2(exp->obj, log_x, SSIZET2NUM(n+1)); 2128 y = BigMath_exp(multiplied, SSIZET2NUM(n)); 2129 RB_GC_GUARD(obj); 2130 2131 return y; 2132} 2133 2134/* call-seq: 2135 * power(n) 2136 * power(n, prec) 2137 * 2138 * Returns the value raised to the power of n. 2139 * 2140 * Note that n must be an Integer. 2141 * 2142 * Also available as the operator ** 2143 */ 2144static VALUE 2145BigDecimal_power(int argc, VALUE*argv, VALUE self) 2146{ 2147 ENTER(5); 2148 VALUE vexp, prec; 2149 Real* exp = NULL; 2150 Real *x, *y; 2151 ssize_t mp, ma, n; 2152 SIGNED_VALUE int_exp; 2153 double d; 2154 2155 rb_scan_args(argc, argv, "11", &vexp, &prec); 2156 2157 GUARD_OBJ(x, GetVpValue(self, 1)); 2158 n = NIL_P(prec) ? (ssize_t)(x->Prec*VpBaseFig()) : NUM2SSIZET(prec); 2159 2160 if (VpIsNaN(x)) { 2161 y = VpCreateRbObject(n, "0#"); 2162 RB_GC_GUARD(y->obj); 2163 VpSetNaN(y); 2164 return ToValue(y); 2165 } 2166 2167retry: 2168 switch (TYPE(vexp)) { 2169 case T_FIXNUM: 2170 break; 2171 2172 case T_BIGNUM: 2173 break; 2174 2175 case T_FLOAT: 2176 d = RFLOAT_VALUE(vexp); 2177 if (d == round(d)) { 2178 vexp = LL2NUM((LONG_LONG)round(d)); 2179 goto retry; 2180 } 2181 exp = GetVpValueWithPrec(vexp, DBL_DIG+1, 1); 2182 break; 2183 2184 case T_RATIONAL: 2185 if (is_zero(RRATIONAL(vexp)->num)) { 2186 if (is_positive(vexp)) { 2187 vexp = INT2FIX(0); 2188 goto retry; 2189 } 2190 } 2191 else if (is_one(RRATIONAL(vexp)->den)) { 2192 vexp = RRATIONAL(vexp)->num; 2193 goto retry; 2194 } 2195 exp = GetVpValueWithPrec(vexp, n, 1); 2196 break; 2197 2198 case T_DATA: 2199 if (is_kind_of_BigDecimal(vexp)) { 2200 VALUE zero = INT2FIX(0); 2201 VALUE rounded = BigDecimal_round(1, &zero, vexp); 2202 if (RTEST(BigDecimal_eq(vexp, rounded))) { 2203 vexp = BigDecimal_to_i(vexp); 2204 goto retry; 2205 } 2206 exp = DATA_PTR(vexp); 2207 break; 2208 } 2209 /* fall through */ 2210 default: 2211 rb_raise(rb_eTypeError, 2212 "wrong argument type %"PRIsVALUE" (expected scalar Numeric)", 2213 RB_OBJ_CLASSNAME(vexp)); 2214 } 2215 2216 if (VpIsZero(x)) { 2217 if (is_negative(vexp)) { 2218 y = VpCreateRbObject(n, "#0"); 2219 RB_GC_GUARD(y->obj); 2220 if (VpGetSign(x) < 0) { 2221 if (is_integer(vexp)) { 2222 if (is_even(vexp)) { 2223 /* (-0) ** (-even_integer) -> Infinity */ 2224 VpSetPosInf(y); 2225 } 2226 else { 2227 /* (-0) ** (-odd_integer) -> -Infinity */ 2228 VpSetNegInf(y); 2229 } 2230 } 2231 else { 2232 /* (-0) ** (-non_integer) -> Infinity */ 2233 VpSetPosInf(y); 2234 } 2235 } 2236 else { 2237 /* (+0) ** (-num) -> Infinity */ 2238 VpSetPosInf(y); 2239 } 2240 return ToValue(y); 2241 } 2242 else if (is_zero(vexp)) { 2243 return ToValue(VpCreateRbObject(n, "1")); 2244 } 2245 else { 2246 return ToValue(VpCreateRbObject(n, "0")); 2247 } 2248 } 2249 2250 if (is_zero(vexp)) { 2251 return ToValue(VpCreateRbObject(n, "1")); 2252 } 2253 else if (is_one(vexp)) { 2254 return self; 2255 } 2256 2257 if (VpIsInf(x)) { 2258 if (is_negative(vexp)) { 2259 if (VpGetSign(x) < 0) { 2260 if (is_integer(vexp)) { 2261 if (is_even(vexp)) { 2262 /* (-Infinity) ** (-even_integer) -> +0 */ 2263 return ToValue(VpCreateRbObject(n, "0")); 2264 } 2265 else { 2266 /* (-Infinity) ** (-odd_integer) -> -0 */ 2267 return ToValue(VpCreateRbObject(n, "-0")); 2268 } 2269 } 2270 else { 2271 /* (-Infinity) ** (-non_integer) -> -0 */ 2272 return ToValue(VpCreateRbObject(n, "-0")); 2273 } 2274 } 2275 else { 2276 return ToValue(VpCreateRbObject(n, "0")); 2277 } 2278 } 2279 else { 2280 y = VpCreateRbObject(n, "0#"); 2281 if (VpGetSign(x) < 0) { 2282 if (is_integer(vexp)) { 2283 if (is_even(vexp)) { 2284 VpSetPosInf(y); 2285 } 2286 else { 2287 VpSetNegInf(y); 2288 } 2289 } 2290 else { 2291 /* TODO: support complex */ 2292 rb_raise(rb_eMathDomainError, 2293 "a non-integral exponent for a negative base"); 2294 } 2295 } 2296 else { 2297 VpSetPosInf(y); 2298 } 2299 return ToValue(y); 2300 } 2301 } 2302 2303 if (exp != NULL) { 2304 return rmpd_power_by_big_decimal(x, exp, n); 2305 } 2306 else if (RB_TYPE_P(vexp, T_BIGNUM)) { 2307 VALUE abs_value = BigDecimal_abs(self); 2308 if (is_one(abs_value)) { 2309 return ToValue(VpCreateRbObject(n, "1")); 2310 } 2311 else if (RTEST(rb_funcall(abs_value, '<', 1, INT2FIX(1)))) { 2312 if (is_negative(vexp)) { 2313 y = VpCreateRbObject(n, "0#"); 2314 if (is_even(vexp)) { 2315 VpSetInf(y, VpGetSign(x)); 2316 } 2317 else { 2318 VpSetInf(y, -VpGetSign(x)); 2319 } 2320 return ToValue(y); 2321 } 2322 else if (VpGetSign(x) < 0 && is_even(vexp)) { 2323 return ToValue(VpCreateRbObject(n, "-0")); 2324 } 2325 else { 2326 return ToValue(VpCreateRbObject(n, "0")); 2327 } 2328 } 2329 else { 2330 if (is_positive(vexp)) { 2331 y = VpCreateRbObject(n, "0#"); 2332 if (is_even(vexp)) { 2333 VpSetInf(y, VpGetSign(x)); 2334 } 2335 else { 2336 VpSetInf(y, -VpGetSign(x)); 2337 } 2338 return ToValue(y); 2339 } 2340 else if (VpGetSign(x) < 0 && is_even(vexp)) { 2341 return ToValue(VpCreateRbObject(n, "-0")); 2342 } 2343 else { 2344 return ToValue(VpCreateRbObject(n, "0")); 2345 } 2346 } 2347 } 2348 2349 int_exp = FIX2INT(vexp); 2350 ma = int_exp; 2351 if (ma < 0) ma = -ma; 2352 if (ma == 0) ma = 1; 2353 2354 if (VpIsDef(x)) { 2355 mp = x->Prec * (VpBaseFig() + 1); 2356 GUARD_OBJ(y, VpCreateRbObject(mp * (ma + 1), "0")); 2357 } 2358 else { 2359 GUARD_OBJ(y, VpCreateRbObject(1, "0")); 2360 } 2361 VpPower(y, x, int_exp); 2362 return ToValue(y); 2363} 2364 2365/* call-seq: 2366 * big_decimal ** exp -> big_decimal 2367 * 2368 * It is a synonym of BigDecimal#power(exp). 2369 */ 2370static VALUE 2371BigDecimal_power_op(VALUE self, VALUE exp) 2372{ 2373 return BigDecimal_power(1, &exp, self); 2374} 2375 2376static VALUE 2377BigDecimal_s_allocate(VALUE klass) 2378{ 2379 return VpNewRbClass(0, NULL, klass)->obj; 2380} 2381 2382static Real *BigDecimal_new(int argc, VALUE *argv); 2383 2384/* call-seq: 2385 * new(initial, digits) 2386 * 2387 * Create a new BigDecimal object. 2388 * 2389 * initial:: The initial value, as an Integer, a Float, a Rational, 2390 * a BigDecimal, or a String. 2391 * 2392 * If it is a String, spaces are ignored and unrecognized characters 2393 * terminate the value. 2394 * 2395 * digits:: The number of significant digits, as a Fixnum. If omitted or 0, 2396 * the number of significant digits is determined from the initial 2397 * value. 2398 * 2399 * The actual number of significant digits used in computation is usually 2400 * larger than the specified number. 2401 */ 2402static VALUE 2403BigDecimal_initialize(int argc, VALUE *argv, VALUE self) 2404{ 2405 ENTER(1); 2406 Real *pv = rb_check_typeddata(self, &BigDecimal_data_type); 2407 Real *x; 2408 2409 GUARD_OBJ(x, BigDecimal_new(argc, argv)); 2410 if (ToValue(x)) { 2411 pv = VpCopy(pv, x); 2412 } 2413 else { 2414 VpFree(pv); 2415 pv = x; 2416 } 2417 DATA_PTR(self) = pv; 2418 pv->obj = self; 2419 return self; 2420} 2421 2422/* :nodoc: 2423 * 2424 * private method to dup and clone the provided BigDecimal +other+ 2425 */ 2426static VALUE 2427BigDecimal_initialize_copy(VALUE self, VALUE other) 2428{ 2429 Real *pv = rb_check_typeddata(self, &BigDecimal_data_type); 2430 Real *x = rb_check_typeddata(other, &BigDecimal_data_type); 2431 2432 if (self != other) { 2433 DATA_PTR(self) = VpCopy(pv, x); 2434 } 2435 return self; 2436} 2437 2438static Real * 2439BigDecimal_new(int argc, VALUE *argv) 2440{ 2441 size_t mf; 2442 VALUE nFig; 2443 VALUE iniValue; 2444 2445 if (rb_scan_args(argc, argv, "11", &iniValue, &nFig) == 1) { 2446 mf = 0; 2447 } 2448 else { 2449 mf = GetPositiveInt(nFig); 2450 } 2451 2452 switch (TYPE(iniValue)) { 2453 case T_DATA: 2454 if (is_kind_of_BigDecimal(iniValue)) { 2455 return DATA_PTR(iniValue); 2456 } 2457 break; 2458 2459 case T_FIXNUM: 2460 /* fall through */ 2461 case T_BIGNUM: 2462 return GetVpValue(iniValue, 1); 2463 2464 case T_FLOAT: 2465 if (mf > DBL_DIG+1) { 2466 rb_raise(rb_eArgError, "precision too large."); 2467 } 2468 /* fall through */ 2469 case T_RATIONAL: 2470 if (NIL_P(nFig)) { 2471 rb_raise(rb_eArgError, 2472 "can't omit precision for a %"PRIsVALUE".", 2473 RB_OBJ_CLASSNAME(iniValue)); 2474 } 2475 return GetVpValueWithPrec(iniValue, mf, 1); 2476 2477 case T_STRING: 2478 /* fall through */ 2479 default: 2480 break; 2481 } 2482 StringValueCStr(iniValue); 2483 return VpAlloc(mf, RSTRING_PTR(iniValue)); 2484} 2485 2486static VALUE 2487BigDecimal_global_new(int argc, VALUE *argv, VALUE self) 2488{ 2489 ENTER(1); 2490 Real *pv; 2491 2492 GUARD_OBJ(pv, BigDecimal_new(argc, argv)); 2493 if (ToValue(pv)) pv = VpCopy(NULL, pv); 2494 pv->obj = TypedData_Wrap_Struct(rb_cBigDecimal, &BigDecimal_data_type, pv); 2495 return pv->obj; 2496} 2497 2498 /* call-seq: 2499 * BigDecimal.limit(digits) 2500 * 2501 * Limit the number of significant digits in newly created BigDecimal 2502 * numbers to the specified value. Rounding is performed as necessary, 2503 * as specified by BigDecimal.mode. 2504 * 2505 * A limit of 0, the default, means no upper limit. 2506 * 2507 * The limit specified by this method takes less priority over any limit 2508 * specified to instance methods such as ceil, floor, truncate, or round. 2509 */ 2510static VALUE 2511BigDecimal_limit(int argc, VALUE *argv, VALUE self) 2512{ 2513 VALUE nFig; 2514 VALUE nCur = INT2NUM(VpGetPrecLimit()); 2515 2516 if(rb_scan_args(argc,argv,"01",&nFig)==1) { 2517 int nf; 2518 if(nFig==Qnil) return nCur; 2519 Check_Type(nFig, T_FIXNUM); 2520 nf = FIX2INT(nFig); 2521 if(nf<0) { 2522 rb_raise(rb_eArgError, "argument must be positive"); 2523 } 2524 VpSetPrecLimit(nf); 2525 } 2526 return nCur; 2527} 2528 2529/* Returns the sign of the value. 2530 * 2531 * Returns a positive value if > 0, a negative value if < 0, and a 2532 * zero if == 0. 2533 * 2534 * The specific value returned indicates the type and sign of the BigDecimal, 2535 * as follows: 2536 * 2537 * BigDecimal::SIGN_NaN:: value is Not a Number 2538 * BigDecimal::SIGN_POSITIVE_ZERO:: value is +0 2539 * BigDecimal::SIGN_NEGATIVE_ZERO:: value is -0 2540 * BigDecimal::SIGN_POSITIVE_INFINITE:: value is +infinity 2541 * BigDecimal::SIGN_NEGATIVE_INFINITE:: value is -infinity 2542 * BigDecimal::SIGN_POSITIVE_FINITE:: value is positive 2543 * BigDecimal::SIGN_NEGATIVE_FINITE:: value is negative 2544 */ 2545static VALUE 2546BigDecimal_sign(VALUE self) 2547{ /* sign */ 2548 int s = GetVpValue(self,1)->sign; 2549 return INT2FIX(s); 2550} 2551 2552/* 2553 * call-seq: BigDecimal.save_exception_mode { ... } 2554 * 2555 * Excecute the provided block, but preserve the exception mode 2556 * 2557 * BigDecimal.save_exception_mode do 2558 * BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false) 2559 * BigDecimal.mode(BigDecimal::EXCEPTION_NaN, false) 2560 * 2561 * BigDecimal.new(BigDecimal('Infinity')) 2562 * BigDecimal.new(BigDecimal('-Infinity')) 2563 * BigDecimal(BigDecimal.new('NaN')) 2564 * end 2565 * 2566 * For use with the BigDecimal::EXCEPTION_* 2567 * 2568 * See BigDecimal.mode 2569 */ 2570static VALUE 2571BigDecimal_save_exception_mode(VALUE self) 2572{ 2573 unsigned short const exception_mode = VpGetException(); 2574 int state; 2575 VALUE ret = rb_protect(rb_yield, Qnil, &state); 2576 VpSetException(exception_mode); 2577 if (state) rb_jump_tag(state); 2578 return ret; 2579} 2580 2581/* 2582 * call-seq: BigDecimal.save_rounding_mode { ... } 2583 * 2584 * Excecute the provided block, but preserve the rounding mode 2585 * 2586 * BigDecimal.save_exception_mode do 2587 * BigDecimal.mode(BigDecimal::ROUND_MODE, :up) 2588 * puts BigDecimal.mode(BigDecimal::ROUND_MODE) 2589 * end 2590 * 2591 * For use with the BigDecimal::ROUND_* 2592 * 2593 * See BigDecimal.mode 2594 */ 2595static VALUE 2596BigDecimal_save_rounding_mode(VALUE self) 2597{ 2598 unsigned short const round_mode = VpGetRoundMode(); 2599 int state; 2600 VALUE ret = rb_protect(rb_yield, Qnil, &state); 2601 VpSetRoundMode(round_mode); 2602 if (state) rb_jump_tag(state); 2603 return ret; 2604} 2605 2606/* 2607 * call-seq: BigDecimal.save_limit { ... } 2608 * 2609 * Excecute the provided block, but preserve the precision limit 2610 * 2611 * BigDecimal.limit(100) 2612 * puts BigDecimal.limit 2613 * BigDecimal.save_limit do 2614 * BigDecimal.limit(200) 2615 * puts BigDecimal.limit 2616 * end 2617 * puts BigDecimal.limit 2618 * 2619 */ 2620static VALUE 2621BigDecimal_save_limit(VALUE self) 2622{ 2623 size_t const limit = VpGetPrecLimit(); 2624 int state; 2625 VALUE ret = rb_protect(rb_yield, Qnil, &state); 2626 VpSetPrecLimit(limit); 2627 if (state) rb_jump_tag(state); 2628 return ret; 2629} 2630 2631/* call-seq: 2632 * BigMath.exp(x, prec) 2633 * 2634 * Computes the value of e (the base of natural logarithms) raised to the 2635 * power of x, to the specified number of digits of precision. 2636 * 2637 * If x is infinity, returns Infinity. 2638 * 2639 * If x is NaN, returns NaN. 2640 */ 2641static VALUE 2642BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec) 2643{ 2644 ssize_t prec, n, i; 2645 Real* vx = NULL; 2646 VALUE one, d, x1, y, z; 2647 int negative = 0; 2648 int infinite = 0; 2649 int nan = 0; 2650 double flo; 2651 2652 prec = NUM2SSIZET(vprec); 2653 if (prec <= 0) { 2654 rb_raise(rb_eArgError, "Zero or negative precision for exp"); 2655 } 2656 2657 /* TODO: the following switch statement is almostly the same as one in the 2658 * BigDecimalCmp function. */ 2659 switch (TYPE(x)) { 2660 case T_DATA: 2661 if (!is_kind_of_BigDecimal(x)) break; 2662 vx = DATA_PTR(x); 2663 negative = VpGetSign(vx) < 0; 2664 infinite = VpIsPosInf(vx) || VpIsNegInf(vx); 2665 nan = VpIsNaN(vx); 2666 break; 2667 2668 case T_FIXNUM: 2669 /* fall through */ 2670 case T_BIGNUM: 2671 vx = GetVpValue(x, 0); 2672 break; 2673 2674 case T_FLOAT: 2675 flo = RFLOAT_VALUE(x); 2676 negative = flo < 0; 2677 infinite = isinf(flo); 2678 nan = isnan(flo); 2679 if (!infinite && !nan) { 2680 vx = GetVpValueWithPrec(x, DBL_DIG+1, 0); 2681 } 2682 break; 2683 2684 case T_RATIONAL: 2685 vx = GetVpValueWithPrec(x, prec, 0); 2686 break; 2687 2688 default: 2689 break; 2690 } 2691 if (infinite) { 2692 if (negative) { 2693 return ToValue(GetVpValueWithPrec(INT2NUM(0), prec, 1)); 2694 } 2695 else { 2696 Real* vy; 2697 vy = VpCreateRbObject(prec, "#0"); 2698 RB_GC_GUARD(vy->obj); 2699 VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE); 2700 return ToValue(vy); 2701 } 2702 } 2703 else if (nan) { 2704 Real* vy; 2705 vy = VpCreateRbObject(prec, "#0"); 2706 RB_GC_GUARD(vy->obj); 2707 VpSetNaN(vy); 2708 return ToValue(vy); 2709 } 2710 else if (vx == NULL) { 2711 cannot_be_coerced_into_BigDecimal(rb_eArgError, x); 2712 } 2713 x = RB_GC_GUARD(vx->obj); 2714 2715 n = prec + rmpd_double_figures(); 2716 negative = VpGetSign(vx) < 0; 2717 if (negative) { 2718 VpSetSign(vx, 1); 2719 } 2720 2721 RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1")); 2722 RB_GC_GUARD(x1) = one; 2723 RB_GC_GUARD(y) = one; 2724 RB_GC_GUARD(d) = y; 2725 RB_GC_GUARD(z) = one; 2726 i = 0; 2727 2728 while (!VpIsZero((Real*)DATA_PTR(d))) { 2729 VALUE argv[2]; 2730 SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y)); 2731 SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d)); 2732 ssize_t m = n - vabs(ey - ed); 2733 if (m <= 0) { 2734 break; 2735 } 2736 else if ((size_t)m < rmpd_double_figures()) { 2737 m = rmpd_double_figures(); 2738 } 2739 2740 x1 = BigDecimal_mult2(x1, x, SSIZET2NUM(n)); 2741 ++i; 2742 z = BigDecimal_mult(z, SSIZET2NUM(i)); 2743 argv[0] = z; 2744 argv[1] = SSIZET2NUM(m); 2745 d = BigDecimal_div2(2, argv, x1); 2746 y = BigDecimal_add(y, d); 2747 } 2748 2749 if (negative) { 2750 VALUE argv[2]; 2751 argv[0] = y; 2752 argv[1] = vprec; 2753 return BigDecimal_div2(2, argv, one); 2754 } 2755 else { 2756 vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y))); 2757 return BigDecimal_round(1, &vprec, y); 2758 } 2759} 2760 2761/* call-seq: 2762 * BigMath.log(x, prec) 2763 * 2764 * Computes the natural logarithm of x to the specified number of digits of 2765 * precision. 2766 * 2767 * If x is zero or negative, raises Math::DomainError. 2768 * 2769 * If x is positive infinity, returns Infinity. 2770 * 2771 * If x is NaN, returns NaN. 2772 */ 2773static VALUE 2774BigMath_s_log(VALUE klass, VALUE x, VALUE vprec) 2775{ 2776 ssize_t prec, n, i; 2777 SIGNED_VALUE expo; 2778 Real* vx = NULL; 2779 VALUE argv[2], vn, one, two, w, x2, y, d; 2780 int zero = 0; 2781 int negative = 0; 2782 int infinite = 0; 2783 int nan = 0; 2784 double flo; 2785 long fix; 2786 2787 if (!is_integer(vprec)) { 2788 rb_raise(rb_eArgError, "precision must be an Integer"); 2789 } 2790 2791 prec = NUM2SSIZET(vprec); 2792 if (prec <= 0) { 2793 rb_raise(rb_eArgError, "Zero or negative precision for exp"); 2794 } 2795 2796 /* TODO: the following switch statement is almostly the same as one in the 2797 * BigDecimalCmp function. */ 2798 switch (TYPE(x)) { 2799 case T_DATA: 2800 if (!is_kind_of_BigDecimal(x)) break; 2801 vx = DATA_PTR(x); 2802 zero = VpIsZero(vx); 2803 negative = VpGetSign(vx) < 0; 2804 infinite = VpIsPosInf(vx) || VpIsNegInf(vx); 2805 nan = VpIsNaN(vx); 2806 break; 2807 2808 case T_FIXNUM: 2809 fix = FIX2LONG(x); 2810 zero = fix == 0; 2811 negative = fix < 0; 2812 goto get_vp_value; 2813 2814 case T_BIGNUM: 2815 zero = RBIGNUM_ZERO_P(x); 2816 negative = RBIGNUM_NEGATIVE_P(x); 2817get_vp_value: 2818 if (zero || negative) break; 2819 vx = GetVpValue(x, 0); 2820 break; 2821 2822 case T_FLOAT: 2823 flo = RFLOAT_VALUE(x); 2824 zero = flo == 0; 2825 negative = flo < 0; 2826 infinite = isinf(flo); 2827 nan = isnan(flo); 2828 if (!zero && !negative && !infinite && !nan) { 2829 vx = GetVpValueWithPrec(x, DBL_DIG+1, 1); 2830 } 2831 break; 2832 2833 case T_RATIONAL: 2834 zero = RRATIONAL_ZERO_P(x); 2835 negative = RRATIONAL_NEGATIVE_P(x); 2836 if (zero || negative) break; 2837 vx = GetVpValueWithPrec(x, prec, 1); 2838 break; 2839 2840 case T_COMPLEX: 2841 rb_raise(rb_eMathDomainError, 2842 "Complex argument for BigMath.log"); 2843 2844 default: 2845 break; 2846 } 2847 if (infinite && !negative) { 2848 Real* vy; 2849 vy = VpCreateRbObject(prec, "#0"); 2850 RB_GC_GUARD(vy->obj); 2851 VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE); 2852 return ToValue(vy); 2853 } 2854 else if (nan) { 2855 Real* vy; 2856 vy = VpCreateRbObject(prec, "#0"); 2857 RB_GC_GUARD(vy->obj); 2858 VpSetNaN(vy); 2859 return ToValue(vy); 2860 } 2861 else if (zero || negative) { 2862 rb_raise(rb_eMathDomainError, 2863 "Zero or negative argument for log"); 2864 } 2865 else if (vx == NULL) { 2866 cannot_be_coerced_into_BigDecimal(rb_eArgError, x); 2867 } 2868 x = ToValue(vx); 2869 2870 RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1")); 2871 RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2")); 2872 2873 n = prec + rmpd_double_figures(); 2874 RB_GC_GUARD(vn) = SSIZET2NUM(n); 2875 expo = VpExponent10(vx); 2876 if (expo < 0 || expo >= 3) { 2877 char buf[16]; 2878 snprintf(buf, 16, "1E%"PRIdVALUE, -expo); 2879 x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn); 2880 } 2881 else { 2882 expo = 0; 2883 } 2884 w = BigDecimal_sub(x, one); 2885 argv[0] = BigDecimal_add(x, one); 2886 argv[1] = vn; 2887 x = BigDecimal_div2(2, argv, w); 2888 RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn); 2889 RB_GC_GUARD(y) = x; 2890 RB_GC_GUARD(d) = y; 2891 i = 1; 2892 while (!VpIsZero((Real*)DATA_PTR(d))) { 2893 SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y)); 2894 SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d)); 2895 ssize_t m = n - vabs(ey - ed); 2896 if (m <= 0) { 2897 break; 2898 } 2899 else if ((size_t)m < rmpd_double_figures()) { 2900 m = rmpd_double_figures(); 2901 } 2902 2903 x = BigDecimal_mult2(x2, x, vn); 2904 i += 2; 2905 argv[0] = SSIZET2NUM(i); 2906 argv[1] = SSIZET2NUM(m); 2907 d = BigDecimal_div2(2, argv, x); 2908 y = BigDecimal_add(y, d); 2909 } 2910 2911 y = BigDecimal_mult(y, two); 2912 if (expo != 0) { 2913 VALUE log10, vexpo, dy; 2914 log10 = BigMath_s_log(klass, INT2FIX(10), vprec); 2915 vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1)); 2916 dy = BigDecimal_mult(log10, vexpo); 2917 y = BigDecimal_add(y, dy); 2918 } 2919 2920 return y; 2921} 2922 2923/* Document-class: BigDecimal 2924 * BigDecimal provides arbitrary-precision floating point decimal arithmetic. 2925 * 2926 * Copyright (C) 2002 by Shigeo Kobayashi <shigeo@tinyforest.gr.jp>. 2927 * 2928 * You may distribute under the terms of either the GNU General Public 2929 * License or the Artistic License, as specified in the README file 2930 * of the BigDecimal distribution. 2931 * 2932 * Documented by mathew <meta@pobox.com>. 2933 * 2934 * = Introduction 2935 * 2936 * Ruby provides built-in support for arbitrary precision integer arithmetic. 2937 * 2938 * For example: 2939 * 2940 * 42**13 #=> 1265437718438866624512 2941 * 2942 * BigDecimal provides similar support for very large or very accurate floating 2943 * point numbers. 2944 * 2945 * Decimal arithmetic is also useful for general calculation, because it 2946 * provides the correct answers people expect--whereas normal binary floating 2947 * point arithmetic often introduces subtle errors because of the conversion 2948 * between base 10 and base 2. 2949 * 2950 * For example, try: 2951 * 2952 * sum = 0 2953 * for i in (1..10000) 2954 * sum = sum + 0.0001 2955 * end 2956 * print sum #=> 0.9999999999999062 2957 * 2958 * and contrast with the output from: 2959 * 2960 * require 'bigdecimal' 2961 * 2962 * sum = BigDecimal.new("0") 2963 * for i in (1..10000) 2964 * sum = sum + BigDecimal.new("0.0001") 2965 * end 2966 * print sum #=> 0.1E1 2967 * 2968 * Similarly: 2969 * 2970 * (BigDecimal.new("1.2") - BigDecimal("1.0")) == BigDecimal("0.2") #=> true 2971 * 2972 * (1.2 - 1.0) == 0.2 #=> false 2973 * 2974 * = Special features of accurate decimal arithmetic 2975 * 2976 * Because BigDecimal is more accurate than normal binary floating point 2977 * arithmetic, it requires some special values. 2978 * 2979 * == Infinity 2980 * 2981 * BigDecimal sometimes needs to return infinity, for example if you divide 2982 * a value by zero. 2983 * 2984 * BigDecimal.new("1.0") / BigDecimal.new("0.0") #=> infinity 2985 * BigDecimal.new("-1.0") / BigDecimal.new("0.0") #=> -infinity 2986 * 2987 * You can represent infinite numbers to BigDecimal using the strings 2988 * <code>'Infinity'</code>, <code>'+Infinity'</code> and 2989 * <code>'-Infinity'</code> (case-sensitive) 2990 * 2991 * == Not a Number 2992 * 2993 * When a computation results in an undefined value, the special value +NaN+ 2994 * (for 'not a number') is returned. 2995 * 2996 * Example: 2997 * 2998 * BigDecimal.new("0.0") / BigDecimal.new("0.0") #=> NaN 2999 * 3000 * You can also create undefined values. 3001 * 3002 * NaN is never considered to be the same as any other value, even NaN itself: 3003 * 3004 * n = BigDecimal.new('NaN') 3005 * n == 0.0 #=> nil 3006 * n == n #=> nil 3007 * 3008 * == Positive and negative zero 3009 * 3010 * If a computation results in a value which is too small to be represented as 3011 * a BigDecimal within the currently specified limits of precision, zero must 3012 * be returned. 3013 * 3014 * If the value which is too small to be represented is negative, a BigDecimal 3015 * value of negative zero is returned. 3016 * 3017 * BigDecimal.new("1.0") / BigDecimal.new("-Infinity") #=> -0.0 3018 * 3019 * If the value is positive, a value of positive zero is returned. 3020 * 3021 * BigDecimal.new("1.0") / BigDecimal.new("Infinity") #=> 0.0 3022 * 3023 * (See BigDecimal.mode for how to specify limits of precision.) 3024 * 3025 * Note that +-0.0+ and +0.0+ are considered to be the same for the purposes of 3026 * comparison. 3027 * 3028 * Note also that in mathematics, there is no particular concept of negative 3029 * or positive zero; true mathematical zero has no sign. 3030 */ 3031void 3032Init_bigdecimal(void) 3033{ 3034 VALUE arg; 3035 3036 id_BigDecimal_exception_mode = rb_intern_const("BigDecimal.exception_mode"); 3037 id_BigDecimal_rounding_mode = rb_intern_const("BigDecimal.rounding_mode"); 3038 id_BigDecimal_precision_limit = rb_intern_const("BigDecimal.precision_limit"); 3039 3040 /* Initialize VP routines */ 3041 VpInit(0UL); 3042 3043 /* Class and method registration */ 3044 rb_cBigDecimal = rb_define_class("BigDecimal",rb_cNumeric); 3045 rb_define_alloc_func(rb_cBigDecimal, BigDecimal_s_allocate); 3046 3047 /* Global function */ 3048 rb_define_global_function("BigDecimal", BigDecimal_global_new, -1); 3049 3050 /* Class methods */ 3051 rb_define_singleton_method(rb_cBigDecimal, "mode", BigDecimal_mode, -1); 3052 rb_define_singleton_method(rb_cBigDecimal, "limit", BigDecimal_limit, -1); 3053 rb_define_singleton_method(rb_cBigDecimal, "double_fig", BigDecimal_double_fig, 0); 3054 rb_define_singleton_method(rb_cBigDecimal, "_load", BigDecimal_load, 1); 3055 rb_define_singleton_method(rb_cBigDecimal, "ver", BigDecimal_version, 0); 3056 3057 rb_define_singleton_method(rb_cBigDecimal, "save_exception_mode", BigDecimal_save_exception_mode, 0); 3058 rb_define_singleton_method(rb_cBigDecimal, "save_rounding_mode", BigDecimal_save_rounding_mode, 0); 3059 rb_define_singleton_method(rb_cBigDecimal, "save_limit", BigDecimal_save_limit, 0); 3060 3061 /* Constants definition */ 3062 3063 /* 3064 * Base value used in internal calculations. On a 32 bit system, BASE 3065 * is 10000, indicating that calculation is done in groups of 4 digits. 3066 * (If it were larger, BASE**2 wouldn't fit in 32 bits, so you couldn't 3067 * guarantee that two groups could always be multiplied together without 3068 * overflow.) 3069 */ 3070 rb_define_const(rb_cBigDecimal, "BASE", INT2FIX((SIGNED_VALUE)VpBaseVal())); 3071 3072 /* Exceptions */ 3073 3074 /* 3075 * 0xff: Determines whether overflow, underflow or zero divide result in 3076 * an exception being thrown. See BigDecimal.mode. 3077 */ 3078 rb_define_const(rb_cBigDecimal, "EXCEPTION_ALL",INT2FIX(VP_EXCEPTION_ALL)); 3079 3080 /* 3081 * 0x02: Determines what happens when the result of a computation is not a 3082 * number (NaN). See BigDecimal.mode. 3083 */ 3084 rb_define_const(rb_cBigDecimal, "EXCEPTION_NaN",INT2FIX(VP_EXCEPTION_NaN)); 3085 3086 /* 3087 * 0x01: Determines what happens when the result of a computation is 3088 * infinity. See BigDecimal.mode. 3089 */ 3090 rb_define_const(rb_cBigDecimal, "EXCEPTION_INFINITY",INT2FIX(VP_EXCEPTION_INFINITY)); 3091 3092 /* 3093 * 0x04: Determines what happens when the result of a computation is an 3094 * underflow (a result too small to be represented). See BigDecimal.mode. 3095 */ 3096 rb_define_const(rb_cBigDecimal, "EXCEPTION_UNDERFLOW",INT2FIX(VP_EXCEPTION_UNDERFLOW)); 3097 3098 /* 3099 * 0x01: Determines what happens when the result of a computation is an 3100 * overflow (a result too large to be represented). See BigDecimal.mode. 3101 */ 3102 rb_define_const(rb_cBigDecimal, "EXCEPTION_OVERFLOW",INT2FIX(VP_EXCEPTION_OVERFLOW)); 3103 3104 /* 3105 * 0x01: Determines what happens when a division by zero is performed. 3106 * See BigDecimal.mode. 3107 */ 3108 rb_define_const(rb_cBigDecimal, "EXCEPTION_ZERODIVIDE",INT2FIX(VP_EXCEPTION_ZERODIVIDE)); 3109 3110 /* 3111 * 0x100: Determines what happens when a result must be rounded in order to 3112 * fit in the appropriate number of significant digits. See 3113 * BigDecimal.mode. 3114 */ 3115 rb_define_const(rb_cBigDecimal, "ROUND_MODE",INT2FIX(VP_ROUND_MODE)); 3116 3117 /* 1: Indicates that values should be rounded away from zero. See 3118 * BigDecimal.mode. 3119 */ 3120 rb_define_const(rb_cBigDecimal, "ROUND_UP",INT2FIX(VP_ROUND_UP)); 3121 3122 /* 2: Indicates that values should be rounded towards zero. See 3123 * BigDecimal.mode. 3124 */ 3125 rb_define_const(rb_cBigDecimal, "ROUND_DOWN",INT2FIX(VP_ROUND_DOWN)); 3126 3127 /* 3: Indicates that digits >= 5 should be rounded up, others rounded down. 3128 * See BigDecimal.mode. */ 3129 rb_define_const(rb_cBigDecimal, "ROUND_HALF_UP",INT2FIX(VP_ROUND_HALF_UP)); 3130 3131 /* 4: Indicates that digits >= 6 should be rounded up, others rounded down. 3132 * See BigDecimal.mode. 3133 */ 3134 rb_define_const(rb_cBigDecimal, "ROUND_HALF_DOWN",INT2FIX(VP_ROUND_HALF_DOWN)); 3135 /* 5: Round towards +infinity. See BigDecimal.mode. */ 3136 rb_define_const(rb_cBigDecimal, "ROUND_CEILING",INT2FIX(VP_ROUND_CEIL)); 3137 3138 /* 6: Round towards -infinity. See BigDecimal.mode. */ 3139 rb_define_const(rb_cBigDecimal, "ROUND_FLOOR",INT2FIX(VP_ROUND_FLOOR)); 3140 3141 /* 7: Round towards the even neighbor. See BigDecimal.mode. */ 3142 rb_define_const(rb_cBigDecimal, "ROUND_HALF_EVEN",INT2FIX(VP_ROUND_HALF_EVEN)); 3143 3144 /* 0: Indicates that a value is not a number. See BigDecimal.sign. */ 3145 rb_define_const(rb_cBigDecimal, "SIGN_NaN",INT2FIX(VP_SIGN_NaN)); 3146 3147 /* 1: Indicates that a value is +0. See BigDecimal.sign. */ 3148 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_ZERO",INT2FIX(VP_SIGN_POSITIVE_ZERO)); 3149 3150 /* -1: Indicates that a value is -0. See BigDecimal.sign. */ 3151 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_ZERO",INT2FIX(VP_SIGN_NEGATIVE_ZERO)); 3152 3153 /* 2: Indicates that a value is positive and finite. See BigDecimal.sign. */ 3154 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_FINITE",INT2FIX(VP_SIGN_POSITIVE_FINITE)); 3155 3156 /* -2: Indicates that a value is negative and finite. See BigDecimal.sign. */ 3157 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_FINITE",INT2FIX(VP_SIGN_NEGATIVE_FINITE)); 3158 3159 /* 3: Indicates that a value is positive and infinite. See BigDecimal.sign. */ 3160 rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_INFINITE",INT2FIX(VP_SIGN_POSITIVE_INFINITE)); 3161 3162 /* -3: Indicates that a value is negative and infinite. See BigDecimal.sign. */ 3163 rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_INFINITE",INT2FIX(VP_SIGN_NEGATIVE_INFINITE)); 3164 3165 arg = rb_str_new2("+Infinity"); 3166 /* Positive infinity value. */ 3167 rb_define_const(rb_cBigDecimal, "INFINITY", BigDecimal_global_new(1, &arg, rb_cBigDecimal)); 3168 arg = rb_str_new2("NaN"); 3169 /* 'Not a Number' value. */ 3170 rb_define_const(rb_cBigDecimal, "NAN", BigDecimal_global_new(1, &arg, rb_cBigDecimal)); 3171 3172 3173 /* instance methods */ 3174 rb_define_method(rb_cBigDecimal, "initialize", BigDecimal_initialize, -1); 3175 rb_define_method(rb_cBigDecimal, "initialize_copy", BigDecimal_initialize_copy, 1); 3176 rb_define_method(rb_cBigDecimal, "precs", BigDecimal_prec, 0); 3177 3178 rb_define_method(rb_cBigDecimal, "add", BigDecimal_add2, 2); 3179 rb_define_method(rb_cBigDecimal, "sub", BigDecimal_sub2, 2); 3180 rb_define_method(rb_cBigDecimal, "mult", BigDecimal_mult2, 2); 3181 rb_define_method(rb_cBigDecimal, "div", BigDecimal_div2, -1); 3182 rb_define_method(rb_cBigDecimal, "hash", BigDecimal_hash, 0); 3183 rb_define_method(rb_cBigDecimal, "to_s", BigDecimal_to_s, -1); 3184 rb_define_method(rb_cBigDecimal, "to_i", BigDecimal_to_i, 0); 3185 rb_define_method(rb_cBigDecimal, "to_int", BigDecimal_to_i, 0); 3186 rb_define_method(rb_cBigDecimal, "to_r", BigDecimal_to_r, 0); 3187 rb_define_method(rb_cBigDecimal, "split", BigDecimal_split, 0); 3188 rb_define_method(rb_cBigDecimal, "+", BigDecimal_add, 1); 3189 rb_define_method(rb_cBigDecimal, "-", BigDecimal_sub, 1); 3190 rb_define_method(rb_cBigDecimal, "+@", BigDecimal_uplus, 0); 3191 rb_define_method(rb_cBigDecimal, "-@", BigDecimal_neg, 0); 3192 rb_define_method(rb_cBigDecimal, "*", BigDecimal_mult, 1); 3193 rb_define_method(rb_cBigDecimal, "/", BigDecimal_div, 1); 3194 rb_define_method(rb_cBigDecimal, "quo", BigDecimal_div, 1); 3195 rb_define_method(rb_cBigDecimal, "%", BigDecimal_mod, 1); 3196 rb_define_method(rb_cBigDecimal, "modulo", BigDecimal_mod, 1); 3197 rb_define_method(rb_cBigDecimal, "remainder", BigDecimal_remainder, 1); 3198 rb_define_method(rb_cBigDecimal, "divmod", BigDecimal_divmod, 1); 3199 /* rb_define_method(rb_cBigDecimal, "dup", BigDecimal_dup, 0); */ 3200 rb_define_method(rb_cBigDecimal, "to_f", BigDecimal_to_f, 0); 3201 rb_define_method(rb_cBigDecimal, "abs", BigDecimal_abs, 0); 3202 rb_define_method(rb_cBigDecimal, "sqrt", BigDecimal_sqrt, 1); 3203 rb_define_method(rb_cBigDecimal, "fix", BigDecimal_fix, 0); 3204 rb_define_method(rb_cBigDecimal, "round", BigDecimal_round, -1); 3205 rb_define_method(rb_cBigDecimal, "frac", BigDecimal_frac, 0); 3206 rb_define_method(rb_cBigDecimal, "floor", BigDecimal_floor, -1); 3207 rb_define_method(rb_cBigDecimal, "ceil", BigDecimal_ceil, -1); 3208 rb_define_method(rb_cBigDecimal, "power", BigDecimal_power, -1); 3209 rb_define_method(rb_cBigDecimal, "**", BigDecimal_power_op, 1); 3210 rb_define_method(rb_cBigDecimal, "<=>", BigDecimal_comp, 1); 3211 rb_define_method(rb_cBigDecimal, "==", BigDecimal_eq, 1); 3212 rb_define_method(rb_cBigDecimal, "===", BigDecimal_eq, 1); 3213 rb_define_method(rb_cBigDecimal, "eql?", BigDecimal_eq, 1); 3214 rb_define_method(rb_cBigDecimal, "<", BigDecimal_lt, 1); 3215 rb_define_method(rb_cBigDecimal, "<=", BigDecimal_le, 1); 3216 rb_define_method(rb_cBigDecimal, ">", BigDecimal_gt, 1); 3217 rb_define_method(rb_cBigDecimal, ">=", BigDecimal_ge, 1); 3218 rb_define_method(rb_cBigDecimal, "zero?", BigDecimal_zero, 0); 3219 rb_define_method(rb_cBigDecimal, "nonzero?", BigDecimal_nonzero, 0); 3220 rb_define_method(rb_cBigDecimal, "coerce", BigDecimal_coerce, 1); 3221 rb_define_method(rb_cBigDecimal, "inspect", BigDecimal_inspect, 0); 3222 rb_define_method(rb_cBigDecimal, "exponent", BigDecimal_exponent, 0); 3223 rb_define_method(rb_cBigDecimal, "sign", BigDecimal_sign, 0); 3224 rb_define_method(rb_cBigDecimal, "nan?", BigDecimal_IsNaN, 0); 3225 rb_define_method(rb_cBigDecimal, "infinite?", BigDecimal_IsInfinite, 0); 3226 rb_define_method(rb_cBigDecimal, "finite?", BigDecimal_IsFinite, 0); 3227 rb_define_method(rb_cBigDecimal, "truncate", BigDecimal_truncate, -1); 3228 rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1); 3229 3230 /* mathematical functions */ 3231 rb_mBigMath = rb_define_module("BigMath"); 3232 rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2); 3233 rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2); 3234 3235 id_up = rb_intern_const("up"); 3236 id_down = rb_intern_const("down"); 3237 id_truncate = rb_intern_const("truncate"); 3238 id_half_up = rb_intern_const("half_up"); 3239 id_default = rb_intern_const("default"); 3240 id_half_down = rb_intern_const("half_down"); 3241 id_half_even = rb_intern_const("half_even"); 3242 id_banker = rb_intern_const("banker"); 3243 id_ceiling = rb_intern_const("ceiling"); 3244 id_ceil = rb_intern_const("ceil"); 3245 id_floor = rb_intern_const("floor"); 3246 id_to_r = rb_intern_const("to_r"); 3247 id_eq = rb_intern_const("=="); 3248} 3249 3250/* 3251 * 3252 * ============================================================================ 3253 * 3254 * vp_ routines begin from here. 3255 * 3256 * ============================================================================ 3257 * 3258 */ 3259#ifdef BIGDECIMAL_DEBUG 3260static int gfDebug = 1; /* Debug switch */ 3261#if 0 3262static int gfCheckVal = 1; /* Value checking flag in VpNmlz() */ 3263#endif 3264#endif /* BIGDECIMAL_DEBUG */ 3265 3266static Real *VpConstOne; /* constant 1.0 */ 3267static Real *VpPt5; /* constant 0.5 */ 3268#define maxnr 100UL /* Maximum iterations for calcurating sqrt. */ 3269 /* used in VpSqrt() */ 3270 3271/* ETC */ 3272#define MemCmp(x,y,z) memcmp(x,y,z) 3273#define StrCmp(x,y) strcmp(x,y) 3274 3275static int VpIsDefOP(Real *c,Real *a,Real *b,int sw); 3276static int AddExponent(Real *a, SIGNED_VALUE n); 3277static BDIGIT VpAddAbs(Real *a,Real *b,Real *c); 3278static BDIGIT VpSubAbs(Real *a,Real *b,Real *c); 3279static 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); 3280static int VpNmlz(Real *a); 3281static void VpFormatSt(char *psz, size_t fFmt); 3282static int VpRdup(Real *m, size_t ind_m); 3283 3284#ifdef BIGDECIMAL_DEBUG 3285static int gnAlloc=0; /* Memory allocation counter */ 3286#endif /* BIGDECIMAL_DEBUG */ 3287 3288VP_EXPORT void * 3289VpMemAlloc(size_t mb) 3290{ 3291 void *p = xmalloc(mb); 3292 if (!p) { 3293 VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1); 3294 } 3295 memset(p, 0, mb); 3296#ifdef BIGDECIMAL_DEBUG 3297 gnAlloc++; /* Count allocation call */ 3298#endif /* BIGDECIMAL_DEBUG */ 3299 return p; 3300} 3301 3302VP_EXPORT void * 3303VpMemRealloc(void *ptr, size_t mb) 3304{ 3305 void *p = xrealloc(ptr, mb); 3306 if (!p) { 3307 VpException(VP_EXCEPTION_MEMORY, "failed to allocate memory", 1); 3308 } 3309 return p; 3310} 3311 3312VP_EXPORT void 3313VpFree(Real *pv) 3314{ 3315 if(pv != NULL) { 3316 xfree(pv); 3317#ifdef BIGDECIMAL_DEBUG 3318 gnAlloc--; /* Decrement allocation count */ 3319 if(gnAlloc==0) { 3320 printf(" *************** All memories allocated freed ****************"); 3321 getchar(); 3322 } 3323 if(gnAlloc<0) { 3324 printf(" ??????????? Too many memory free calls(%d) ?????????????\n",gnAlloc); 3325 getchar(); 3326 } 3327#endif /* BIGDECIMAL_DEBUG */ 3328 } 3329} 3330 3331/* 3332 * EXCEPTION Handling. 3333 */ 3334 3335#define rmpd_set_thread_local_exception_mode(mode) \ 3336 rb_thread_local_aset( \ 3337 rb_thread_current(), \ 3338 id_BigDecimal_exception_mode, \ 3339 INT2FIX((int)(mode)) \ 3340 ) 3341 3342static unsigned short 3343VpGetException (void) 3344{ 3345 VALUE const vmode = rb_thread_local_aref( 3346 rb_thread_current(), 3347 id_BigDecimal_exception_mode 3348 ); 3349 3350 if (NIL_P(vmode)) { 3351 rmpd_set_thread_local_exception_mode(RMPD_EXCEPTION_MODE_DEFAULT); 3352 return RMPD_EXCEPTION_MODE_DEFAULT; 3353 } 3354 3355 return (unsigned short)FIX2UINT(vmode); 3356} 3357 3358static void 3359VpSetException(unsigned short f) 3360{ 3361 rmpd_set_thread_local_exception_mode(f); 3362} 3363 3364/* 3365 * Precision limit. 3366 */ 3367 3368#define rmpd_set_thread_local_precision_limit(limit) \ 3369 rb_thread_local_aset( \ 3370 rb_thread_current(), \ 3371 id_BigDecimal_precision_limit, \ 3372 SIZET2NUM(limit) \ 3373 ) 3374#define RMPD_PRECISION_LIMIT_DEFAULT ((size_t)0) 3375 3376/* These 2 functions added at v1.1.7 */ 3377VP_EXPORT size_t 3378VpGetPrecLimit(void) 3379{ 3380 VALUE const vlimit = rb_thread_local_aref( 3381 rb_thread_current(), 3382 id_BigDecimal_precision_limit 3383 ); 3384 3385 if (NIL_P(vlimit)) { 3386 rmpd_set_thread_local_precision_limit(RMPD_PRECISION_LIMIT_DEFAULT); 3387 return RMPD_PRECISION_LIMIT_DEFAULT; 3388 } 3389 3390 return NUM2SIZET(vlimit); 3391} 3392 3393VP_EXPORT size_t 3394VpSetPrecLimit(size_t n) 3395{ 3396 size_t const s = VpGetPrecLimit(); 3397 rmpd_set_thread_local_precision_limit(n); 3398 return s; 3399} 3400 3401/* 3402 * Rounding mode. 3403 */ 3404 3405#define rmpd_set_thread_local_rounding_mode(mode) \ 3406 rb_thread_local_aset( \ 3407 rb_thread_current(), \ 3408 id_BigDecimal_rounding_mode, \ 3409 INT2FIX((int)(mode)) \ 3410 ) 3411 3412VP_EXPORT unsigned short 3413VpGetRoundMode(void) 3414{ 3415 VALUE const vmode = rb_thread_local_aref( 3416 rb_thread_current(), 3417 id_BigDecimal_rounding_mode 3418 ); 3419 3420 if (NIL_P(vmode)) { 3421 rmpd_set_thread_local_rounding_mode(RMPD_ROUNDING_MODE_DEFAULT); 3422 return RMPD_ROUNDING_MODE_DEFAULT; 3423 } 3424 3425 return (unsigned short)FIX2INT(vmode); 3426} 3427 3428VP_EXPORT int 3429VpIsRoundMode(unsigned short n) 3430{ 3431 switch (n) { 3432 case VP_ROUND_UP: 3433 case VP_ROUND_DOWN: 3434 case VP_ROUND_HALF_UP: 3435 case VP_ROUND_HALF_DOWN: 3436 case VP_ROUND_CEIL: 3437 case VP_ROUND_FLOOR: 3438 case VP_ROUND_HALF_EVEN: 3439 return 1; 3440 3441 default: 3442 return 0; 3443 } 3444} 3445 3446VP_EXPORT unsigned short 3447VpSetRoundMode(unsigned short n) 3448{ 3449 if (VpIsRoundMode(n)) { 3450 rmpd_set_thread_local_rounding_mode(n); 3451 return n; 3452 } 3453 3454 return VpGetRoundMode(); 3455} 3456 3457/* 3458 * 0.0 & 1.0 generator 3459 * These gZero_..... and gOne_..... can be any name 3460 * referenced from nowhere except Zero() and One(). 3461 * gZero_..... and gOne_..... must have global scope 3462 * (to let the compiler know they may be changed in outside 3463 * (... but not actually..)). 3464 */ 3465volatile const double gZero_ABCED9B1_CE73__00400511F31D = 0.0; 3466volatile const double gOne_ABCED9B4_CE73__00400511F31D = 1.0; 3467static double 3468Zero(void) 3469{ 3470 return gZero_ABCED9B1_CE73__00400511F31D; 3471} 3472 3473static double 3474One(void) 3475{ 3476 return gOne_ABCED9B4_CE73__00400511F31D; 3477} 3478 3479/* 3480 ---------------------------------------------------------------- 3481 Value of sign in Real structure is reserved for future use. 3482 short sign; 3483 ==0 : NaN 3484 1 : Positive zero 3485 -1 : Negative zero 3486 2 : Positive number 3487 -2 : Negative number 3488 3 : Positive infinite number 3489 -3 : Negative infinite number 3490 ---------------------------------------------------------------- 3491*/ 3492 3493VP_EXPORT double 3494VpGetDoubleNaN(void) /* Returns the value of NaN */ 3495{ 3496 static double fNaN = 0.0; 3497 if(fNaN==0.0) fNaN = Zero()/Zero(); 3498 return fNaN; 3499} 3500 3501VP_EXPORT double 3502VpGetDoublePosInf(void) /* Returns the value of +Infinity */ 3503{ 3504 static double fInf = 0.0; 3505 if(fInf==0.0) fInf = One()/Zero(); 3506 return fInf; 3507} 3508 3509VP_EXPORT double 3510VpGetDoubleNegInf(void) /* Returns the value of -Infinity */ 3511{ 3512 static double fInf = 0.0; 3513 if(fInf==0.0) fInf = -(One()/Zero()); 3514 return fInf; 3515} 3516 3517VP_EXPORT double 3518VpGetDoubleNegZero(void) /* Returns the value of -0 */ 3519{ 3520 static double nzero = 1000.0; 3521 if(nzero!=0.0) nzero = (One()/VpGetDoubleNegInf()); 3522 return nzero; 3523} 3524 3525#if 0 /* unused */ 3526VP_EXPORT int 3527VpIsNegDoubleZero(double v) 3528{ 3529 double z = VpGetDoubleNegZero(); 3530 return MemCmp(&v,&z,sizeof(v))==0; 3531} 3532#endif 3533 3534VP_EXPORT int 3535VpException(unsigned short f, const char *str,int always) 3536{ 3537 VALUE exc; 3538 int fatal=0; 3539 unsigned short const exception_mode = VpGetException(); 3540 3541 if(f==VP_EXCEPTION_OP || f==VP_EXCEPTION_MEMORY) always = 1; 3542 3543 if (always || (exception_mode & f)) { 3544 switch(f) 3545 { 3546 /* 3547 case VP_EXCEPTION_OVERFLOW: 3548 */ 3549 case VP_EXCEPTION_ZERODIVIDE: 3550 case VP_EXCEPTION_INFINITY: 3551 case VP_EXCEPTION_NaN: 3552 case VP_EXCEPTION_UNDERFLOW: 3553 case VP_EXCEPTION_OP: 3554 exc = rb_eFloatDomainError; 3555 goto raise; 3556 case VP_EXCEPTION_MEMORY: 3557 fatal = 1; 3558 goto raise; 3559 default: 3560 fatal = 1; 3561 goto raise; 3562 } 3563 } 3564 return 0; /* 0 Means VpException() raised no exception */ 3565 3566raise: 3567 if(fatal) rb_fatal("%s", str); 3568 else rb_raise(exc, "%s", str); 3569 return 0; 3570} 3571 3572/* Throw exception or returns 0,when resulting c is Inf or NaN */ 3573/* sw=1:+ 2:- 3:* 4:/ */ 3574static int 3575VpIsDefOP(Real *c,Real *a,Real *b,int sw) 3576{ 3577 if(VpIsNaN(a) || VpIsNaN(b)) { 3578 /* at least a or b is NaN */ 3579 VpSetNaN(c); 3580 goto NaN; 3581 } 3582 3583 if(VpIsInf(a)) { 3584 if(VpIsInf(b)) { 3585 switch(sw) 3586 { 3587 case 1: /* + */ 3588 if(VpGetSign(a)==VpGetSign(b)) { 3589 VpSetInf(c,VpGetSign(a)); 3590 goto Inf; 3591 } else { 3592 VpSetNaN(c); 3593 goto NaN; 3594 } 3595 case 2: /* - */ 3596 if(VpGetSign(a)!=VpGetSign(b)) { 3597 VpSetInf(c,VpGetSign(a)); 3598 goto Inf; 3599 } else { 3600 VpSetNaN(c); 3601 goto NaN; 3602 } 3603 break; 3604 case 3: /* * */ 3605 VpSetInf(c,VpGetSign(a)*VpGetSign(b)); 3606 goto Inf; 3607 break; 3608 case 4: /* / */ 3609 VpSetNaN(c); 3610 goto NaN; 3611 } 3612 VpSetNaN(c); 3613 goto NaN; 3614 } 3615 /* Inf op Finite */ 3616 switch(sw) 3617 { 3618 case 1: /* + */ 3619 case 2: /* - */ 3620 VpSetInf(c,VpGetSign(a)); 3621 break; 3622 case 3: /* * */ 3623 if(VpIsZero(b)) { 3624 VpSetNaN(c); 3625 goto NaN; 3626 } 3627 VpSetInf(c,VpGetSign(a)*VpGetSign(b)); 3628 break; 3629 case 4: /* / */ 3630 VpSetInf(c,VpGetSign(a)*VpGetSign(b)); 3631 } 3632 goto Inf; 3633 } 3634 3635 if(VpIsInf(b)) { 3636 switch(sw) 3637 { 3638 case 1: /* + */ 3639 VpSetInf(c,VpGetSign(b)); 3640 break; 3641 case 2: /* - */ 3642 VpSetInf(c,-VpGetSign(b)); 3643 break; 3644 case 3: /* * */ 3645 if(VpIsZero(a)) { 3646 VpSetNaN(c); 3647 goto NaN; 3648 } 3649 VpSetInf(c,VpGetSign(a)*VpGetSign(b)); 3650 break; 3651 case 4: /* / */ 3652 VpSetZero(c,VpGetSign(a)*VpGetSign(b)); 3653 } 3654 goto Inf; 3655 } 3656 return 1; /* Results OK */ 3657 3658Inf: 3659 return VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0); 3660NaN: 3661 return VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'",0); 3662} 3663 3664/* 3665 ---------------------------------------------------------------- 3666*/ 3667 3668/* 3669 * returns number of chars needed to represent vp in specified format. 3670 */ 3671VP_EXPORT size_t 3672VpNumOfChars(Real *vp,const char *pszFmt) 3673{ 3674 SIGNED_VALUE ex; 3675 size_t nc; 3676 3677 if(vp == NULL) return BASE_FIG*2+6; 3678 if(!VpIsDef(vp)) return 32; /* not sure,may be OK */ 3679 3680 switch(*pszFmt) 3681 { 3682 case 'F': 3683 nc = BASE_FIG*(vp->Prec + 1)+2; 3684 ex = vp->exponent; 3685 if(ex < 0) { 3686 nc += BASE_FIG*(size_t)(-ex); 3687 } 3688 else { 3689 if((size_t)ex > vp->Prec) { 3690 nc += BASE_FIG*((size_t)ex - vp->Prec); 3691 } 3692 } 3693 break; 3694 case 'E': 3695 default: 3696 nc = BASE_FIG*(vp->Prec + 2)+6; /* 3: sign + exponent chars */ 3697 } 3698 return nc; 3699} 3700 3701/* 3702 * Initializer for Vp routines and constants used. 3703 * [Input] 3704 * BaseVal: Base value(assigned to BASE) for Vp calculation. 3705 * It must be the form BaseVal=10**n.(n=1,2,3,...) 3706 * If Base <= 0L,then the BASE will be calcurated so 3707 * that BASE is as large as possible satisfying the 3708 * relation MaxVal <= BASE*(BASE+1). Where the value 3709 * MaxVal is the largest value which can be represented 3710 * by one BDIGIT word in the computer used. 3711 * 3712 * [Returns] 3713 * 1+DBL_DIG ... OK 3714 */ 3715VP_EXPORT size_t 3716VpInit(BDIGIT BaseVal) 3717{ 3718 /* Setup +/- Inf NaN -0 */ 3719 VpGetDoubleNaN(); 3720 VpGetDoublePosInf(); 3721 VpGetDoubleNegInf(); 3722 VpGetDoubleNegZero(); 3723 3724 /* Allocates Vp constants. */ 3725 VpConstOne = VpAlloc(1UL, "1"); 3726 VpPt5 = VpAlloc(1UL, ".5"); 3727 3728#ifdef BIGDECIMAL_DEBUG 3729 gnAlloc = 0; 3730#endif /* BIGDECIMAL_DEBUG */ 3731 3732#ifdef BIGDECIMAL_DEBUG 3733 if(gfDebug) { 3734 printf("VpInit: BaseVal = %lu\n", BaseVal); 3735 printf(" BASE = %lu\n", BASE); 3736 printf(" HALF_BASE = %lu\n", HALF_BASE); 3737 printf(" BASE1 = %lu\n", BASE1); 3738 printf(" BASE_FIG = %u\n", BASE_FIG); 3739 printf(" DBLE_FIG = %d\n", DBLE_FIG); 3740 } 3741#endif /* BIGDECIMAL_DEBUG */ 3742 3743 return rmpd_double_figures(); 3744} 3745 3746VP_EXPORT Real * 3747VpOne(void) 3748{ 3749 return VpConstOne; 3750} 3751 3752/* If exponent overflows,then raise exception or returns 0 */ 3753static int 3754AddExponent(Real *a, SIGNED_VALUE n) 3755{ 3756 SIGNED_VALUE e = a->exponent; 3757 SIGNED_VALUE m = e+n; 3758 SIGNED_VALUE eb, mb; 3759 if(e>0) { 3760 if(n>0) { 3761 mb = m*(SIGNED_VALUE)BASE_FIG; 3762 eb = e*(SIGNED_VALUE)BASE_FIG; 3763 if(mb<eb) goto overflow; 3764 } 3765 } else if(n<0) { 3766 mb = m*(SIGNED_VALUE)BASE_FIG; 3767 eb = e*(SIGNED_VALUE)BASE_FIG; 3768 if(mb>eb) goto underflow; 3769 } 3770 a->exponent = m; 3771 return 1; 3772 3773/* Overflow/Underflow ==> Raise exception or returns 0 */ 3774underflow: 3775 VpSetZero(a,VpGetSign(a)); 3776 return VpException(VP_EXCEPTION_UNDERFLOW,"Exponent underflow",0); 3777 3778overflow: 3779 VpSetInf(a,VpGetSign(a)); 3780 return VpException(VP_EXCEPTION_OVERFLOW,"Exponent overflow",0); 3781} 3782 3783/* 3784 * Allocates variable. 3785 * [Input] 3786 * mx ... allocation unit, if zero then mx is determined by szVal. 3787 * The mx is the number of effective digits can to be stored. 3788 * szVal ... value assigned(char). If szVal==NULL,then zero is assumed. 3789 * If szVal[0]=='#' then Max. Prec. will not be considered(1.1.7), 3790 * full precision specified by szVal is allocated. 3791 * 3792 * [Returns] 3793 * Pointer to the newly allocated variable, or 3794 * NULL be returned if memory allocation is failed,or any error. 3795 */ 3796VP_EXPORT Real * 3797VpAlloc(size_t mx, const char *szVal) 3798{ 3799 size_t i, ni, ipn, ipf, nf, ipe, ne, nalloc; 3800 char v,*psz; 3801 int sign=1; 3802 Real *vp = NULL; 3803 size_t mf = VpGetPrecLimit(); 3804 VALUE buf; 3805 3806 mx = (mx + BASE_FIG - 1) / BASE_FIG + 1; /* Determine allocation unit. */ 3807 if (szVal) { 3808 while (ISSPACE(*szVal)) szVal++; 3809 if (*szVal != '#') { 3810 if (mf) { 3811 mf = (mf + BASE_FIG - 1) / BASE_FIG + 2; /* Needs 1 more for div */ 3812 if (mx > mf) { 3813 mx = mf; 3814 } 3815 } 3816 } 3817 else { 3818 ++szVal; 3819 } 3820 } 3821 else { 3822 /* necessary to be able to store */ 3823 /* at least mx digits. */ 3824 /* szVal==NULL ==> allocate zero value. */ 3825 vp = VpAllocReal(mx); 3826 /* xmalloc() alway returns(or throw interruption) */ 3827 vp->MaxPrec = mx; /* set max precision */ 3828 VpSetZero(vp,1); /* initialize vp to zero. */ 3829 return vp; 3830 } 3831 3832 /* Skip all '_' after digit: 2006-6-30 */ 3833 ni = 0; 3834 buf = rb_str_tmp_new(strlen(szVal)+1); 3835 psz = RSTRING_PTR(buf); 3836 i = 0; 3837 ipn = 0; 3838 while ((psz[i]=szVal[ipn]) != 0) { 3839 if (ISDIGIT(psz[i])) ++ni; 3840 if (psz[i] == '_') { 3841 if (ni > 0) { ipn++; continue; } 3842 psz[i] = 0; 3843 break; 3844 } 3845 ++i; 3846 ++ipn; 3847 } 3848 /* Skip trailing spaces */ 3849 while (--i > 0) { 3850 if (ISSPACE(psz[i])) psz[i] = 0; 3851 else break; 3852 } 3853 szVal = psz; 3854 3855 /* Check on Inf & NaN */ 3856 if (StrCmp(szVal, SZ_PINF) == 0 || 3857 StrCmp(szVal, SZ_INF) == 0 ) { 3858 vp = VpAllocReal(1); 3859 vp->MaxPrec = 1; /* set max precision */ 3860 VpSetPosInf(vp); 3861 return vp; 3862 } 3863 if (StrCmp(szVal, SZ_NINF) == 0) { 3864 vp = VpAllocReal(1); 3865 vp->MaxPrec = 1; /* set max precision */ 3866 VpSetNegInf(vp); 3867 return vp; 3868 } 3869 if (StrCmp(szVal, SZ_NaN) == 0) { 3870 vp = VpAllocReal(1); 3871 vp->MaxPrec = 1; /* set max precision */ 3872 VpSetNaN(vp); 3873 return vp; 3874 } 3875 3876 /* check on number szVal[] */ 3877 ipn = i = 0; 3878 if (szVal[i] == '-') { sign=-1; ++i; } 3879 else if (szVal[i] == '+') ++i; 3880 /* Skip digits */ 3881 ni = 0; /* digits in mantissa */ 3882 while ((v = szVal[i]) != 0) { 3883 if (!ISDIGIT(v)) break; 3884 ++i; 3885 ++ni; 3886 } 3887 nf = 0; 3888 ipf = 0; 3889 ipe = 0; 3890 ne = 0; 3891 if (v) { 3892 /* other than digit nor \0 */ 3893 if (szVal[i] == '.') { /* xxx. */ 3894 ++i; 3895 ipf = i; 3896 while ((v = szVal[i]) != 0) { /* get fraction part. */ 3897 if (!ISDIGIT(v)) break; 3898 ++i; 3899 ++nf; 3900 } 3901 } 3902 ipe = 0; /* Exponent */ 3903 3904 switch (szVal[i]) { 3905 case '\0': 3906 break; 3907 case 'e': case 'E': 3908 case 'd': case 'D': 3909 ++i; 3910 ipe = i; 3911 v = szVal[i]; 3912 if ((v == '-') || (v == '+')) ++i; 3913 while ((v=szVal[i]) != 0) { 3914 if (!ISDIGIT(v)) break; 3915 ++i; 3916 ++ne; 3917 } 3918 break; 3919 default: 3920 break; 3921 } 3922 } 3923 nalloc = (ni + nf + BASE_FIG - 1) / BASE_FIG + 1; /* set effective allocation */ 3924 /* units for szVal[] */ 3925 if (mx <= 0) mx = 1; 3926 nalloc = Max(nalloc, mx); 3927 mx = nalloc; 3928 vp = VpAllocReal(mx); 3929 /* xmalloc() alway returns(or throw interruption) */ 3930 vp->MaxPrec = mx; /* set max precision */ 3931 VpSetZero(vp, sign); 3932 VpCtoV(vp, &szVal[ipn], ni, &szVal[ipf], nf, &szVal[ipe], ne); 3933 rb_str_resize(buf, 0); 3934 return vp; 3935} 3936 3937/* 3938 * Assignment(c=a). 3939 * [Input] 3940 * a ... RHSV 3941 * isw ... switch for assignment. 3942 * c = a when isw > 0 3943 * c = -a when isw < 0 3944 * if c->MaxPrec < a->Prec,then round operation 3945 * will be performed. 3946 * [Output] 3947 * c ... LHSV 3948 */ 3949VP_EXPORT size_t 3950VpAsgn(Real *c, Real *a, int isw) 3951{ 3952 size_t n; 3953 if(VpIsNaN(a)) { 3954 VpSetNaN(c); 3955 return 0; 3956 } 3957 if(VpIsInf(a)) { 3958 VpSetInf(c,isw*VpGetSign(a)); 3959 return 0; 3960 } 3961 3962 /* check if the RHS is zero */ 3963 if(!VpIsZero(a)) { 3964 c->exponent = a->exponent; /* store exponent */ 3965 VpSetSign(c,(isw*VpGetSign(a))); /* set sign */ 3966 n =(a->Prec < c->MaxPrec) ?(a->Prec) :(c->MaxPrec); 3967 c->Prec = n; 3968 memcpy(c->frac, a->frac, n * sizeof(BDIGIT)); 3969 /* Needs round ? */ 3970 if(isw!=10) { 3971 /* Not in ActiveRound */ 3972 if(c->Prec < a->Prec) { 3973 VpInternalRound(c,n,(n>0)?a->frac[n-1]:0,a->frac[n]); 3974 } else { 3975 VpLimitRound(c,0); 3976 } 3977 } 3978 } else { 3979 /* The value of 'a' is zero. */ 3980 VpSetZero(c,isw*VpGetSign(a)); 3981 return 1; 3982 } 3983 return c->Prec*BASE_FIG; 3984} 3985 3986/* 3987 * c = a + b when operation = 1 or 2 3988 * = a - b when operation = -1 or -2. 3989 * Returns number of significant digits of c 3990 */ 3991VP_EXPORT size_t 3992VpAddSub(Real *c, Real *a, Real *b, int operation) 3993{ 3994 short sw, isw; 3995 Real *a_ptr, *b_ptr; 3996 size_t n, na, nb, i; 3997 BDIGIT mrv; 3998 3999#ifdef BIGDECIMAL_DEBUG 4000 if(gfDebug) { 4001 VPrint(stdout, "VpAddSub(enter) a=% \n", a); 4002 VPrint(stdout, " b=% \n", b); 4003 printf(" operation=%d\n", operation); 4004 } 4005#endif /* BIGDECIMAL_DEBUG */ 4006 4007 if(!VpIsDefOP(c,a,b,(operation>0)?1:2)) return 0; /* No significant digits */ 4008 4009 /* check if a or b is zero */ 4010 if(VpIsZero(a)) { 4011 /* a is zero,then assign b to c */ 4012 if(!VpIsZero(b)) { 4013 VpAsgn(c, b, operation); 4014 } else { 4015 /* Both a and b are zero. */ 4016 if(VpGetSign(a)<0 && operation*VpGetSign(b)<0) { 4017 /* -0 -0 */ 4018 VpSetZero(c,-1); 4019 } else { 4020 VpSetZero(c,1); 4021 } 4022 return 1; /* 0: 1 significant digits */ 4023 } 4024 return c->Prec*BASE_FIG; 4025 } 4026 if(VpIsZero(b)) { 4027 /* b is zero,then assign a to c. */ 4028 VpAsgn(c, a, 1); 4029 return c->Prec*BASE_FIG; 4030 } 4031 4032 if(operation < 0) sw = -1; 4033 else sw = 1; 4034 4035 /* compare absolute value. As a result,|a_ptr|>=|b_ptr| */ 4036 if(a->exponent > b->exponent) { 4037 a_ptr = a; 4038 b_ptr = b; 4039 } /* |a|>|b| */ 4040 else if(a->exponent < b->exponent) { 4041 a_ptr = b; 4042 b_ptr = a; 4043 } /* |a|<|b| */ 4044 else { 4045 /* Exponent part of a and b is the same,then compare fraction */ 4046 /* part */ 4047 na = a->Prec; 4048 nb = b->Prec; 4049 n = Min(na, nb); 4050 for(i=0;i < n; ++i) { 4051 if(a->frac[i] > b->frac[i]) { 4052 a_ptr = a; 4053 b_ptr = b; 4054 goto end_if; 4055 } else if(a->frac[i] < b->frac[i]) { 4056 a_ptr = b; 4057 b_ptr = a; 4058 goto end_if; 4059 } 4060 } 4061 if(na > nb) { 4062 a_ptr = a; 4063 b_ptr = b; 4064 goto end_if; 4065 } else if(na < nb) { 4066 a_ptr = b; 4067 b_ptr = a; 4068 goto end_if; 4069 } 4070 /* |a| == |b| */ 4071 if(VpGetSign(a) + sw *VpGetSign(b) == 0) { 4072 VpSetZero(c,1); /* abs(a)=abs(b) and operation = '-' */ 4073 return c->Prec*BASE_FIG; 4074 } 4075 a_ptr = a; 4076 b_ptr = b; 4077 } 4078 4079end_if: 4080 isw = VpGetSign(a) + sw *VpGetSign(b); 4081 /* 4082 * isw = 0 ...( 1)+(-1),( 1)-( 1),(-1)+(1),(-1)-(-1) 4083 * = 2 ...( 1)+( 1),( 1)-(-1) 4084 * =-2 ...(-1)+(-1),(-1)-( 1) 4085 * If isw==0, then c =(Sign a_ptr)(|a_ptr|-|b_ptr|) 4086 * else c =(Sign ofisw)(|a_ptr|+|b_ptr|) 4087 */ 4088 if(isw) { /* addition */ 4089 VpSetSign(c, 1); 4090 mrv = VpAddAbs(a_ptr, b_ptr, c); 4091 VpSetSign(c, isw / 2); 4092 } else { /* subtraction */ 4093 VpSetSign(c, 1); 4094 mrv = VpSubAbs(a_ptr, b_ptr, c); 4095 if(a_ptr == a) { 4096 VpSetSign(c,VpGetSign(a)); 4097 } else { 4098 VpSetSign(c,VpGetSign(a_ptr) * sw); 4099 } 4100 } 4101 VpInternalRound(c,0,(c->Prec>0)?c->frac[c->Prec-1]:0,mrv); 4102 4103#ifdef BIGDECIMAL_DEBUG 4104 if(gfDebug) { 4105 VPrint(stdout, "VpAddSub(result) c=% \n", c); 4106 VPrint(stdout, " a=% \n", a); 4107 VPrint(stdout, " b=% \n", b); 4108 printf(" operation=%d\n", operation); 4109 } 4110#endif /* BIGDECIMAL_DEBUG */ 4111 return c->Prec*BASE_FIG; 4112} 4113 4114/* 4115 * Addition of two variable precisional variables 4116 * a and b assuming abs(a)>abs(b). 4117 * c = abs(a) + abs(b) ; where |a|>=|b| 4118 */ 4119static BDIGIT 4120VpAddAbs(Real *a, Real *b, Real *c) 4121{ 4122 size_t word_shift; 4123 size_t ap; 4124 size_t bp; 4125 size_t cp; 4126 size_t a_pos; 4127 size_t b_pos, b_pos_with_word_shift; 4128 size_t c_pos; 4129 BDIGIT av, bv, carry, mrv; 4130 4131#ifdef BIGDECIMAL_DEBUG 4132 if(gfDebug) { 4133 VPrint(stdout, "VpAddAbs called: a = %\n", a); 4134 VPrint(stdout, " b = %\n", b); 4135 } 4136#endif /* BIGDECIMAL_DEBUG */ 4137 4138 word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv); 4139 a_pos = ap; 4140 b_pos = bp; 4141 c_pos = cp; 4142 if(word_shift==(size_t)-1L) return 0; /* Overflow */ 4143 if(b_pos == (size_t)-1L) goto Assign_a; 4144 4145 mrv = av + bv; /* Most right val. Used for round. */ 4146 4147 /* Just assign the last few digits of b to c because a has no */ 4148 /* corresponding digits to be added. */ 4149 while(b_pos + word_shift > a_pos) { 4150 --c_pos; 4151 if(b_pos > 0) { 4152 c->frac[c_pos] = b->frac[--b_pos]; 4153 } else { 4154 --word_shift; 4155 c->frac[c_pos] = 0; 4156 } 4157 } 4158 4159 /* Just assign the last few digits of a to c because b has no */ 4160 /* corresponding digits to be added. */ 4161 b_pos_with_word_shift = b_pos + word_shift; 4162 while(a_pos > b_pos_with_word_shift) { 4163 c->frac[--c_pos] = a->frac[--a_pos]; 4164 } 4165 carry = 0; /* set first carry be zero */ 4166 4167 /* Now perform addition until every digits of b will be */ 4168 /* exhausted. */ 4169 while(b_pos > 0) { 4170 c->frac[--c_pos] = a->frac[--a_pos] + b->frac[--b_pos] + carry; 4171 if(c->frac[c_pos] >= BASE) { 4172 c->frac[c_pos] -= BASE; 4173 carry = 1; 4174 } else { 4175 carry = 0; 4176 } 4177 } 4178 4179 /* Just assign the first few digits of a with considering */ 4180 /* the carry obtained so far because b has been exhausted. */ 4181 while(a_pos > 0) { 4182 c->frac[--c_pos] = a->frac[--a_pos] + carry; 4183 if(c->frac[c_pos] >= BASE) { 4184 c->frac[c_pos] -= BASE; 4185 carry = 1; 4186 } else { 4187 carry = 0; 4188 } 4189 } 4190 if(c_pos) c->frac[c_pos - 1] += carry; 4191 goto Exit; 4192 4193Assign_a: 4194 VpAsgn(c, a, 1); 4195 mrv = 0; 4196 4197Exit: 4198 4199#ifdef BIGDECIMAL_DEBUG 4200 if(gfDebug) { 4201 VPrint(stdout, "VpAddAbs exit: c=% \n", c); 4202 } 4203#endif /* BIGDECIMAL_DEBUG */ 4204 return mrv; 4205} 4206 4207/* 4208 * c = abs(a) - abs(b) 4209 */ 4210static BDIGIT 4211VpSubAbs(Real *a, Real *b, Real *c) 4212{ 4213 size_t word_shift; 4214 size_t ap; 4215 size_t bp; 4216 size_t cp; 4217 size_t a_pos; 4218 size_t b_pos, b_pos_with_word_shift; 4219 size_t c_pos; 4220 BDIGIT av, bv, borrow, mrv; 4221 4222#ifdef BIGDECIMAL_DEBUG 4223 if(gfDebug) { 4224 VPrint(stdout, "VpSubAbs called: a = %\n", a); 4225 VPrint(stdout, " b = %\n", b); 4226 } 4227#endif /* BIGDECIMAL_DEBUG */ 4228 4229 word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv); 4230 a_pos = ap; 4231 b_pos = bp; 4232 c_pos = cp; 4233 if(word_shift==(size_t)-1L) return 0; /* Overflow */ 4234 if(b_pos == (size_t)-1L) goto Assign_a; 4235 4236 if(av >= bv) { 4237 mrv = av - bv; 4238 borrow = 0; 4239 } else { 4240 mrv = 0; 4241 borrow = 1; 4242 } 4243 4244 /* Just assign the values which are the BASE subtracted by */ 4245 /* each of the last few digits of the b because the a has no */ 4246 /* corresponding digits to be subtracted. */ 4247 if(b_pos + word_shift > a_pos) { 4248 while(b_pos + word_shift > a_pos) { 4249 --c_pos; 4250 if(b_pos > 0) { 4251 c->frac[c_pos] = BASE - b->frac[--b_pos] - borrow; 4252 } else { 4253 --word_shift; 4254 c->frac[c_pos] = BASE - borrow; 4255 } 4256 borrow = 1; 4257 } 4258 } 4259 /* Just assign the last few digits of a to c because b has no */ 4260 /* corresponding digits to subtract. */ 4261 4262 b_pos_with_word_shift = b_pos + word_shift; 4263 while(a_pos > b_pos_with_word_shift) { 4264 c->frac[--c_pos] = a->frac[--a_pos]; 4265 } 4266 4267 /* Now perform subtraction until every digits of b will be */ 4268 /* exhausted. */ 4269 while(b_pos > 0) { 4270 --c_pos; 4271 if(a->frac[--a_pos] < b->frac[--b_pos] + borrow) { 4272 c->frac[c_pos] = BASE + a->frac[a_pos] - b->frac[b_pos] - borrow; 4273 borrow = 1; 4274 } else { 4275 c->frac[c_pos] = a->frac[a_pos] - b->frac[b_pos] - borrow; 4276 borrow = 0; 4277 } 4278 } 4279 4280 /* Just assign the first few digits of a with considering */ 4281 /* the borrow obtained so far because b has been exhausted. */ 4282 while(a_pos > 0) { 4283 --c_pos; 4284 if(a->frac[--a_pos] < borrow) { 4285 c->frac[c_pos] = BASE + a->frac[a_pos] - borrow; 4286 borrow = 1; 4287 } else { 4288 c->frac[c_pos] = a->frac[a_pos] - borrow; 4289 borrow = 0; 4290 } 4291 } 4292 if(c_pos) c->frac[c_pos - 1] -= borrow; 4293 goto Exit; 4294 4295Assign_a: 4296 VpAsgn(c, a, 1); 4297 mrv = 0; 4298 4299Exit: 4300#ifdef BIGDECIMAL_DEBUG 4301 if(gfDebug) { 4302 VPrint(stdout, "VpSubAbs exit: c=% \n", c); 4303 } 4304#endif /* BIGDECIMAL_DEBUG */ 4305 return mrv; 4306} 4307 4308/* 4309 * Note: If(av+bv)>= HALF_BASE,then 1 will be added to the least significant 4310 * digit of c(In case of addition). 4311 * ------------------------- figure of output ----------------------------------- 4312 * a = xxxxxxxxxxx 4313 * b = xxxxxxxxxx 4314 * c =xxxxxxxxxxxxxxx 4315 * word_shift = | | 4316 * right_word = | | (Total digits in RHSV) 4317 * left_word = | | (Total digits in LHSV) 4318 * a_pos = | 4319 * b_pos = | 4320 * c_pos = | 4321 */ 4322static size_t 4323VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, BDIGIT *av, BDIGIT *bv) 4324{ 4325 size_t left_word, right_word, word_shift; 4326 c->frac[0] = 0; 4327 *av = *bv = 0; 4328 word_shift =((a->exponent) -(b->exponent)); 4329 left_word = b->Prec + word_shift; 4330 right_word = Max((a->Prec),left_word); 4331 left_word =(c->MaxPrec) - 1; /* -1 ... prepare for round up */ 4332 /* 4333 * check if 'round' is needed. 4334 */ 4335 if(right_word > left_word) { /* round ? */ 4336 /*--------------------------------- 4337 * Actual size of a = xxxxxxAxx 4338 * Actual size of b = xxxBxxxxx 4339 * Max. size of c = xxxxxx 4340 * Round off = |-----| 4341 * c_pos = | 4342 * right_word = | 4343 * a_pos = | 4344 */ 4345 *c_pos = right_word = left_word + 1; /* Set resulting precision */ 4346 /* be equal to that of c */ 4347 if((a->Prec) >=(c->MaxPrec)) { 4348 /* 4349 * a = xxxxxxAxxx 4350 * c = xxxxxx 4351 * a_pos = | 4352 */ 4353 *a_pos = left_word; 4354 *av = a->frac[*a_pos]; /* av is 'A' shown in above. */ 4355 } else { 4356 /* 4357 * a = xxxxxxx 4358 * c = xxxxxxxxxx 4359 * a_pos = | 4360 */ 4361 *a_pos = a->Prec; 4362 } 4363 if((b->Prec + word_shift) >= c->MaxPrec) { 4364 /* 4365 * a = xxxxxxxxx 4366 * b = xxxxxxxBxxx 4367 * c = xxxxxxxxxxx 4368 * b_pos = | 4369 */ 4370 if(c->MaxPrec >=(word_shift + 1)) { 4371 *b_pos = c->MaxPrec - word_shift - 1; 4372 *bv = b->frac[*b_pos]; 4373 } else { 4374 *b_pos = -1L; 4375 } 4376 } else { 4377 /* 4378 * a = xxxxxxxxxxxxxxxx 4379 * b = xxxxxx 4380 * c = xxxxxxxxxxxxx 4381 * b_pos = | 4382 */ 4383 *b_pos = b->Prec; 4384 } 4385 } else { /* The MaxPrec of c - 1 > The Prec of a + b */ 4386 /* 4387 * a = xxxxxxx 4388 * b = xxxxxx 4389 * c = xxxxxxxxxxx 4390 * c_pos = | 4391 */ 4392 *b_pos = b->Prec; 4393 *a_pos = a->Prec; 4394 *c_pos = right_word + 1; 4395 } 4396 c->Prec = *c_pos; 4397 c->exponent = a->exponent; 4398 if(!AddExponent(c,1)) return (size_t)-1L; 4399 return word_shift; 4400} 4401 4402/* 4403 * Return number og significant digits 4404 * c = a * b , Where a = a0a1a2 ... an 4405 * b = b0b1b2 ... bm 4406 * c = c0c1c2 ... cl 4407 * a0 a1 ... an * bm 4408 * a0 a1 ... an * bm-1 4409 * . . . 4410 * . . . 4411 * a0 a1 .... an * b0 4412 * +_____________________________ 4413 * c0 c1 c2 ...... cl 4414 * nc <---| 4415 * MaxAB |--------------------| 4416 */ 4417VP_EXPORT size_t 4418VpMult(Real *c, Real *a, Real *b) 4419{ 4420 size_t MxIndA, MxIndB, MxIndAB, MxIndC; 4421 size_t ind_c, i, ii, nc; 4422 size_t ind_as, ind_ae, ind_bs; 4423 BDIGIT carry; 4424 BDIGIT_DBL s; 4425 Real *w; 4426 4427#ifdef BIGDECIMAL_DEBUG 4428 if(gfDebug) { 4429 VPrint(stdout, "VpMult(Enter): a=% \n", a); 4430 VPrint(stdout, " b=% \n", b); 4431 } 4432#endif /* BIGDECIMAL_DEBUG */ 4433 4434 if(!VpIsDefOP(c,a,b,3)) return 0; /* No significant digit */ 4435 4436 if(VpIsZero(a) || VpIsZero(b)) { 4437 /* at least a or b is zero */ 4438 VpSetZero(c,VpGetSign(a)*VpGetSign(b)); 4439 return 1; /* 0: 1 significant digit */ 4440 } 4441 4442 if(VpIsOne(a)) { 4443 VpAsgn(c, b, VpGetSign(a)); 4444 goto Exit; 4445 } 4446 if(VpIsOne(b)) { 4447 VpAsgn(c, a, VpGetSign(b)); 4448 goto Exit; 4449 } 4450 if((b->Prec) >(a->Prec)) { 4451 /* Adjust so that digits(a)>digits(b) */ 4452 w = a; 4453 a = b; 4454 b = w; 4455 } 4456 w = NULL; 4457 MxIndA = a->Prec - 1; 4458 MxIndB = b->Prec - 1; 4459 MxIndC = c->MaxPrec - 1; 4460 MxIndAB = a->Prec + b->Prec - 1; 4461 4462 if(MxIndC < MxIndAB) { /* The Max. prec. of c < Prec(a)+Prec(b) */ 4463 w = c; 4464 c = VpAlloc((size_t)((MxIndAB + 1) * BASE_FIG), "#0"); 4465 MxIndC = MxIndAB; 4466 } 4467 4468 /* set LHSV c info */ 4469 4470 c->exponent = a->exponent; /* set exponent */ 4471 if(!AddExponent(c,b->exponent)) { 4472 if(w) VpFree(c); 4473 return 0; 4474 } 4475 VpSetSign(c,VpGetSign(a)*VpGetSign(b)); /* set sign */ 4476 carry = 0; 4477 nc = ind_c = MxIndAB; 4478 memset(c->frac, 0, (nc + 1) * sizeof(BDIGIT)); /* Initialize c */ 4479 c->Prec = nc + 1; /* set precision */ 4480 for(nc = 0; nc < MxIndAB; ++nc, --ind_c) { 4481 if(nc < MxIndB) { /* The left triangle of the Fig. */ 4482 ind_as = MxIndA - nc; 4483 ind_ae = MxIndA; 4484 ind_bs = MxIndB; 4485 } else if(nc <= MxIndA) { /* The middle rectangular of the Fig. */ 4486 ind_as = MxIndA - nc; 4487 ind_ae = MxIndA -(nc - MxIndB); 4488 ind_bs = MxIndB; 4489 } else if(nc > MxIndA) { /* The right triangle of the Fig. */ 4490 ind_as = 0; 4491 ind_ae = MxIndAB - nc - 1; 4492 ind_bs = MxIndB -(nc - MxIndA); 4493 } 4494 4495 for(i = ind_as; i <= ind_ae; ++i) { 4496 s = (BDIGIT_DBL)a->frac[i] * b->frac[ind_bs--]; 4497 carry = (BDIGIT)(s / BASE); 4498 s -= (BDIGIT_DBL)carry * BASE; 4499 c->frac[ind_c] += (BDIGIT)s; 4500 if(c->frac[ind_c] >= BASE) { 4501 s = c->frac[ind_c] / BASE; 4502 carry += (BDIGIT)s; 4503 c->frac[ind_c] -= (BDIGIT)(s * BASE); 4504 } 4505 if(carry) { 4506 ii = ind_c; 4507 while(ii-- > 0) { 4508 c->frac[ii] += carry; 4509 if(c->frac[ii] >= BASE) { 4510 carry = c->frac[ii] / BASE; 4511 c->frac[ii] -= (carry * BASE); 4512 } else { 4513 break; 4514 } 4515 } 4516 } 4517 } 4518 } 4519 if(w != NULL) { /* free work variable */ 4520 VpNmlz(c); 4521 VpAsgn(w, c, 1); 4522 VpFree(c); 4523 c = w; 4524 } else { 4525 VpLimitRound(c,0); 4526 } 4527 4528Exit: 4529#ifdef BIGDECIMAL_DEBUG 4530 if(gfDebug) { 4531 VPrint(stdout, "VpMult(c=a*b): c=% \n", c); 4532 VPrint(stdout, " a=% \n", a); 4533 VPrint(stdout, " b=% \n", b); 4534 } 4535#endif /*BIGDECIMAL_DEBUG */ 4536 return c->Prec*BASE_FIG; 4537} 4538 4539/* 4540 * c = a / b, remainder = r 4541 */ 4542VP_EXPORT size_t 4543VpDivd(Real *c, Real *r, Real *a, Real *b) 4544{ 4545 size_t word_a, word_b, word_c, word_r; 4546 size_t i, n, ind_a, ind_b, ind_c, ind_r; 4547 size_t nLoop; 4548 BDIGIT_DBL q, b1, b1p1, b1b2, b1b2p1, r1r2; 4549 BDIGIT borrow, borrow1, borrow2; 4550 BDIGIT_DBL qb; 4551 4552#ifdef BIGDECIMAL_DEBUG 4553 if(gfDebug) { 4554 VPrint(stdout, " VpDivd(c=a/b) a=% \n", a); 4555 VPrint(stdout, " b=% \n", b); 4556 } 4557#endif /*BIGDECIMAL_DEBUG */ 4558 4559 VpSetNaN(r); 4560 if(!VpIsDefOP(c,a,b,4)) goto Exit; 4561 if(VpIsZero(a)&&VpIsZero(b)) { 4562 VpSetNaN(c); 4563 return VpException(VP_EXCEPTION_NaN,"(VpDivd) 0/0 not defined(NaN)",0); 4564 } 4565 if(VpIsZero(b)) { 4566 VpSetInf(c,VpGetSign(a)*VpGetSign(b)); 4567 return VpException(VP_EXCEPTION_ZERODIVIDE,"(VpDivd) Divide by zero",0); 4568 } 4569 if(VpIsZero(a)) { 4570 /* numerator a is zero */ 4571 VpSetZero(c,VpGetSign(a)*VpGetSign(b)); 4572 VpSetZero(r,VpGetSign(a)*VpGetSign(b)); 4573 goto Exit; 4574 } 4575 if(VpIsOne(b)) { 4576 /* divide by one */ 4577 VpAsgn(c, a, VpGetSign(b)); 4578 VpSetZero(r,VpGetSign(a)); 4579 goto Exit; 4580 } 4581 4582 word_a = a->Prec; 4583 word_b = b->Prec; 4584 word_c = c->MaxPrec; 4585 word_r = r->MaxPrec; 4586 4587 ind_c = 0; 4588 ind_r = 1; 4589 4590 if(word_a >= word_r) goto space_error; 4591 4592 r->frac[0] = 0; 4593 while(ind_r <= word_a) { 4594 r->frac[ind_r] = a->frac[ind_r - 1]; 4595 ++ind_r; 4596 } 4597 4598 while(ind_r < word_r) r->frac[ind_r++] = 0; 4599 while(ind_c < word_c) c->frac[ind_c++] = 0; 4600 4601 /* initial procedure */ 4602 b1 = b1p1 = b->frac[0]; 4603 if(b->Prec <= 1) { 4604 b1b2p1 = b1b2 = b1p1 * BASE; 4605 } else { 4606 b1p1 = b1 + 1; 4607 b1b2p1 = b1b2 = b1 * BASE + b->frac[1]; 4608 if(b->Prec > 2) ++b1b2p1; 4609 } 4610 4611 /* */ 4612 /* loop start */ 4613 ind_c = word_r - 1; 4614 nLoop = Min(word_c,ind_c); 4615 ind_c = 1; 4616 while(ind_c < nLoop) { 4617 if(r->frac[ind_c] == 0) { 4618 ++ind_c; 4619 continue; 4620 } 4621 r1r2 = (BDIGIT_DBL)r->frac[ind_c] * BASE + r->frac[ind_c + 1]; 4622 if(r1r2 == b1b2) { 4623 /* The first two word digits is the same */ 4624 ind_b = 2; 4625 ind_a = ind_c + 2; 4626 while(ind_b < word_b) { 4627 if(r->frac[ind_a] < b->frac[ind_b]) goto div_b1p1; 4628 if(r->frac[ind_a] > b->frac[ind_b]) break; 4629 ++ind_a; 4630 ++ind_b; 4631 } 4632 /* The first few word digits of r and b is the same and */ 4633 /* the first different word digit of w is greater than that */ 4634 /* of b, so quotinet is 1 and just subtract b from r. */ 4635 borrow = 0; /* quotient=1, then just r-b */ 4636 ind_b = b->Prec - 1; 4637 ind_r = ind_c + ind_b; 4638 if(ind_r >= word_r) goto space_error; 4639 n = ind_b; 4640 for(i = 0; i <= n; ++i) { 4641 if(r->frac[ind_r] < b->frac[ind_b] + borrow) { 4642 r->frac[ind_r] += (BASE - (b->frac[ind_b] + borrow)); 4643 borrow = 1; 4644 } else { 4645 r->frac[ind_r] = r->frac[ind_r] - b->frac[ind_b] - borrow; 4646 borrow = 0; 4647 } 4648 --ind_r; 4649 --ind_b; 4650 } 4651 ++c->frac[ind_c]; 4652 goto carry; 4653 } 4654 /* The first two word digits is not the same, */ 4655 /* then compare magnitude, and divide actually. */ 4656 if(r1r2 >= b1b2p1) { 4657 q = r1r2 / b1b2p1; /* q == (BDIGIT)q */ 4658 c->frac[ind_c] += (BDIGIT)q; 4659 ind_r = b->Prec + ind_c - 1; 4660 goto sub_mult; 4661 } 4662 4663div_b1p1: 4664 if(ind_c + 1 >= word_c) goto out_side; 4665 q = r1r2 / b1p1; /* q == (BDIGIT)q */ 4666 c->frac[ind_c + 1] += (BDIGIT)q; 4667 ind_r = b->Prec + ind_c; 4668 4669sub_mult: 4670 borrow1 = borrow2 = 0; 4671 ind_b = word_b - 1; 4672 if(ind_r >= word_r) goto space_error; 4673 n = ind_b; 4674 for(i = 0; i <= n; ++i) { 4675 /* now, perform r = r - q * b */ 4676 qb = q * b->frac[ind_b]; 4677 if (qb < BASE) borrow1 = 0; 4678 else { 4679 borrow1 = (BDIGIT)(qb / BASE); 4680 qb -= (BDIGIT_DBL)borrow1 * BASE; /* get qb < BASE */ 4681 } 4682 if(r->frac[ind_r] < qb) { 4683 r->frac[ind_r] += (BDIGIT)(BASE - qb); 4684 borrow2 = borrow2 + borrow1 + 1; 4685 } else { 4686 r->frac[ind_r] -= (BDIGIT)qb; 4687 borrow2 += borrow1; 4688 } 4689 if(borrow2) { 4690 if(r->frac[ind_r - 1] < borrow2) { 4691 r->frac[ind_r - 1] += (BASE - borrow2); 4692 borrow2 = 1; 4693 } else { 4694 r->frac[ind_r - 1] -= borrow2; 4695 borrow2 = 0; 4696 } 4697 } 4698 --ind_r; 4699 --ind_b; 4700 } 4701 4702 r->frac[ind_r] -= borrow2; 4703carry: 4704 ind_r = ind_c; 4705 while(c->frac[ind_r] >= BASE) { 4706 c->frac[ind_r] -= BASE; 4707 --ind_r; 4708 ++c->frac[ind_r]; 4709 } 4710 } 4711 /* End of operation, now final arrangement */ 4712out_side: 4713 c->Prec = word_c; 4714 c->exponent = a->exponent; 4715 if(!AddExponent(c,2)) return 0; 4716 if(!AddExponent(c,-(b->exponent))) return 0; 4717 4718 VpSetSign(c,VpGetSign(a)*VpGetSign(b)); 4719 VpNmlz(c); /* normalize c */ 4720 r->Prec = word_r; 4721 r->exponent = a->exponent; 4722 if(!AddExponent(r,1)) return 0; 4723 VpSetSign(r,VpGetSign(a)); 4724 VpNmlz(r); /* normalize r(remainder) */ 4725 goto Exit; 4726 4727space_error: 4728#ifdef BIGDECIMAL_DEBUG 4729 if(gfDebug) { 4730 printf(" word_a=%lu\n", word_a); 4731 printf(" word_b=%lu\n", word_b); 4732 printf(" word_c=%lu\n", word_c); 4733 printf(" word_r=%lu\n", word_r); 4734 printf(" ind_r =%lu\n", ind_r); 4735 } 4736#endif /* BIGDECIMAL_DEBUG */ 4737 rb_bug("ERROR(VpDivd): space for remainder too small."); 4738 4739Exit: 4740#ifdef BIGDECIMAL_DEBUG 4741 if(gfDebug) { 4742 VPrint(stdout, " VpDivd(c=a/b), c=% \n", c); 4743 VPrint(stdout, " r=% \n", r); 4744 } 4745#endif /* BIGDECIMAL_DEBUG */ 4746 return c->Prec*BASE_FIG; 4747} 4748 4749/* 4750 * Input a = 00000xxxxxxxx En(5 preceeding zeros) 4751 * Output a = xxxxxxxx En-5 4752 */ 4753static int 4754VpNmlz(Real *a) 4755{ 4756 size_t ind_a, i; 4757 4758 if (!VpIsDef(a)) goto NoVal; 4759 if (VpIsZero(a)) goto NoVal; 4760 4761 ind_a = a->Prec; 4762 while (ind_a--) { 4763 if (a->frac[ind_a]) { 4764 a->Prec = ind_a + 1; 4765 i = 0; 4766 while (a->frac[i] == 0) ++i; /* skip the first few zeros */ 4767 if (i) { 4768 a->Prec -= i; 4769 if (!AddExponent(a, -(SIGNED_VALUE)i)) return 0; 4770 memmove(&a->frac[0], &a->frac[i], a->Prec*sizeof(BDIGIT)); 4771 } 4772 return 1; 4773 } 4774 } 4775 /* a is zero(no non-zero digit) */ 4776 VpSetZero(a, VpGetSign(a)); 4777 return 0; 4778 4779NoVal: 4780 a->frac[0] = 0; 4781 a->Prec = 1; 4782 return 0; 4783} 4784 4785/* 4786 * VpComp = 0 ... if a=b, 4787 * Pos ... a>b, 4788 * Neg ... a<b. 4789 * 999 ... result undefined(NaN) 4790 */ 4791VP_EXPORT int 4792VpComp(Real *a, Real *b) 4793{ 4794 int val; 4795 size_t mx, ind; 4796 int e; 4797 val = 0; 4798 if(VpIsNaN(a)||VpIsNaN(b)) return 999; 4799 if(!VpIsDef(a)) { 4800 if(!VpIsDef(b)) e = a->sign - b->sign; 4801 else e = a->sign; 4802 if(e>0) return 1; 4803 else if(e<0) return -1; 4804 else return 0; 4805 } 4806 if(!VpIsDef(b)) { 4807 e = -b->sign; 4808 if(e>0) return 1; 4809 else return -1; 4810 } 4811 /* Zero check */ 4812 if(VpIsZero(a)) { 4813 if(VpIsZero(b)) return 0; /* both zero */ 4814 val = -VpGetSign(b); 4815 goto Exit; 4816 } 4817 if(VpIsZero(b)) { 4818 val = VpGetSign(a); 4819 goto Exit; 4820 } 4821 4822 /* compare sign */ 4823 if(VpGetSign(a) > VpGetSign(b)) { 4824 val = 1; /* a>b */ 4825 goto Exit; 4826 } 4827 if(VpGetSign(a) < VpGetSign(b)) { 4828 val = -1; /* a<b */ 4829 goto Exit; 4830 } 4831 4832 /* a and b have same sign, && signe!=0,then compare exponent */ 4833 if((a->exponent) >(b->exponent)) { 4834 val = VpGetSign(a); 4835 goto Exit; 4836 } 4837 if((a->exponent) <(b->exponent)) { 4838 val = -VpGetSign(b); 4839 goto Exit; 4840 } 4841 4842 /* a and b have same exponent, then compare significand. */ 4843 mx =((a->Prec) <(b->Prec)) ?(a->Prec) :(b->Prec); 4844 ind = 0; 4845 while(ind < mx) { 4846 if((a->frac[ind]) >(b->frac[ind])) { 4847 val = VpGetSign(a); 4848 goto Exit; 4849 } 4850 if((a->frac[ind]) <(b->frac[ind])) { 4851 val = -VpGetSign(b); 4852 goto Exit; 4853 } 4854 ++ind; 4855 } 4856 if((a->Prec) >(b->Prec)) { 4857 val = VpGetSign(a); 4858 } else if((a->Prec) <(b->Prec)) { 4859 val = -VpGetSign(b); 4860 } 4861 4862Exit: 4863 if (val> 1) val = 1; 4864 else if(val<-1) val = -1; 4865 4866#ifdef BIGDECIMAL_DEBUG 4867 if(gfDebug) { 4868 VPrint(stdout, " VpComp a=%\n", a); 4869 VPrint(stdout, " b=%\n", b); 4870 printf(" ans=%d\n", val); 4871 } 4872#endif /* BIGDECIMAL_DEBUG */ 4873 return (int)val; 4874} 4875 4876#ifdef BIGDECIMAL_ENABLE_VPRINT 4877/* 4878 * cntl_chr ... ASCIIZ Character, print control characters 4879 * Available control codes: 4880 * % ... VP variable. To print '%', use '%%'. 4881 * \n ... new line 4882 * \b ... backspace 4883 * ... tab 4884 * Note: % must must not appear more than once 4885 * a ... VP variable to be printed 4886 */ 4887VP_EXPORT int 4888VPrint(FILE *fp, const char *cntl_chr, Real *a) 4889{ 4890 size_t i, j, nc, nd, ZeroSup; 4891 BDIGIT m, e, nn; 4892 4893 /* Check if NaN & Inf. */ 4894 if(VpIsNaN(a)) { 4895 fprintf(fp,SZ_NaN); 4896 return 8; 4897 } 4898 if(VpIsPosInf(a)) { 4899 fprintf(fp,SZ_INF); 4900 return 8; 4901 } 4902 if(VpIsNegInf(a)) { 4903 fprintf(fp,SZ_NINF); 4904 return 9; 4905 } 4906 if(VpIsZero(a)) { 4907 fprintf(fp,"0.0"); 4908 return 3; 4909 } 4910 4911 j = 0; 4912 nd = nc = 0; /* nd : number of digits in fraction part(every 10 digits, */ 4913 /* nd<=10). */ 4914 /* nc : number of caracters printed */ 4915 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */ 4916 while(*(cntl_chr + j)) { 4917 if((*(cntl_chr + j) == '%') &&(*(cntl_chr + j + 1) != '%')) { 4918 nc = 0; 4919 if(!VpIsZero(a)) { 4920 if(VpGetSign(a) < 0) { 4921 fprintf(fp, "-"); 4922 ++nc; 4923 } 4924 nc += fprintf(fp, "0."); 4925 for(i=0; i < a->Prec; ++i) { 4926 m = BASE1; 4927 e = a->frac[i]; 4928 while(m) { 4929 nn = e / m; 4930 if((!ZeroSup) || nn) { 4931 nc += fprintf(fp, "%lu", (unsigned long)nn); /* The leading zero(s) */ 4932 /* as 0.00xx will not */ 4933 /* be printed. */ 4934 ++nd; 4935 ZeroSup = 0; /* Set to print succeeding zeros */ 4936 } 4937 if(nd >= 10) { /* print ' ' after every 10 digits */ 4938 nd = 0; 4939 nc += fprintf(fp, " "); 4940 } 4941 e = e - nn * m; 4942 m /= 10; 4943 } 4944 } 4945 nc += fprintf(fp, "E%"PRIdSIZE, VpExponent10(a)); 4946 } else { 4947 nc += fprintf(fp, "0.0"); 4948 } 4949 } else { 4950 ++nc; 4951 if(*(cntl_chr + j) == '\\') { 4952 switch(*(cntl_chr + j + 1)) { 4953 case 'n': 4954 fprintf(fp, "\n"); 4955 ++j; 4956 break; 4957 case 't': 4958 fprintf(fp, "\t"); 4959 ++j; 4960 break; 4961 case 'b': 4962 fprintf(fp, "\n"); 4963 ++j; 4964 break; 4965 default: 4966 fprintf(fp, "%c", *(cntl_chr + j)); 4967 break; 4968 } 4969 } else { 4970 fprintf(fp, "%c", *(cntl_chr + j)); 4971 if(*(cntl_chr + j) == '%') ++j; 4972 } 4973 } 4974 j++; 4975 } 4976 return (int)nc; 4977} 4978#endif /* BIGDECIMAL_ENABLE_VPRINT */ 4979 4980static void 4981VpFormatSt(char *psz, size_t fFmt) 4982{ 4983 size_t ie, i, nf = 0; 4984 char ch; 4985 4986 if(fFmt<=0) return; 4987 4988 ie = strlen(psz); 4989 for(i = 0; i < ie; ++i) { 4990 ch = psz[i]; 4991 if(!ch) break; 4992 if(ISSPACE(ch) || ch=='-' || ch=='+') continue; 4993 if(ch == '.') { nf = 0;continue;} 4994 if(ch == 'E') break; 4995 nf++; 4996 if(nf > fFmt) { 4997 memmove(psz + i + 1, psz + i, ie - i + 1); 4998 ++ie; 4999 nf = 0; 5000 psz[i] = ' '; 5001 } 5002 } 5003} 5004 5005VP_EXPORT ssize_t 5006VpExponent10(Real *a) 5007{ 5008 ssize_t ex; 5009 size_t n; 5010 5011 if (!VpHasVal(a)) return 0; 5012 5013 ex = a->exponent * (ssize_t)BASE_FIG; 5014 n = BASE1; 5015 while ((a->frac[0] / n) == 0) { 5016 --ex; 5017 n /= 10; 5018 } 5019 return ex; 5020} 5021 5022VP_EXPORT void 5023VpSzMantissa(Real *a,char *psz) 5024{ 5025 size_t i, n, ZeroSup; 5026 BDIGIT_DBL m, e, nn; 5027 5028 if(VpIsNaN(a)) { 5029 sprintf(psz,SZ_NaN); 5030 return; 5031 } 5032 if(VpIsPosInf(a)) { 5033 sprintf(psz,SZ_INF); 5034 return; 5035 } 5036 if(VpIsNegInf(a)) { 5037 sprintf(psz,SZ_NINF); 5038 return; 5039 } 5040 5041 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */ 5042 if(!VpIsZero(a)) { 5043 if(VpGetSign(a) < 0) *psz++ = '-'; 5044 n = a->Prec; 5045 for (i=0; i < n; ++i) { 5046 m = BASE1; 5047 e = a->frac[i]; 5048 while (m) { 5049 nn = e / m; 5050 if((!ZeroSup) || nn) { 5051 sprintf(psz, "%lu", (unsigned long)nn); /* The leading zero(s) */ 5052 psz += strlen(psz); 5053 /* as 0.00xx will be ignored. */ 5054 ZeroSup = 0; /* Set to print succeeding zeros */ 5055 } 5056 e = e - nn * m; 5057 m /= 10; 5058 } 5059 } 5060 *psz = 0; 5061 while(psz[-1]=='0') *(--psz) = 0; 5062 } else { 5063 if(VpIsPosZero(a)) sprintf(psz, "0"); 5064 else sprintf(psz, "-0"); 5065 } 5066} 5067 5068VP_EXPORT int 5069VpToSpecialString(Real *a,char *psz,int fPlus) 5070/* fPlus =0:default, =1: set ' ' before digits , =2: set '+' before digits. */ 5071{ 5072 if(VpIsNaN(a)) { 5073 sprintf(psz,SZ_NaN); 5074 return 1; 5075 } 5076 5077 if(VpIsPosInf(a)) { 5078 if(fPlus==1) { 5079 *psz++ = ' '; 5080 } else if(fPlus==2) { 5081 *psz++ = '+'; 5082 } 5083 sprintf(psz,SZ_INF); 5084 return 1; 5085 } 5086 if(VpIsNegInf(a)) { 5087 sprintf(psz,SZ_NINF); 5088 return 1; 5089 } 5090 if(VpIsZero(a)) { 5091 if(VpIsPosZero(a)) { 5092 if(fPlus==1) sprintf(psz, " 0.0"); 5093 else if(fPlus==2) sprintf(psz, "+0.0"); 5094 else sprintf(psz, "0.0"); 5095 } else sprintf(psz, "-0.0"); 5096 return 1; 5097 } 5098 return 0; 5099} 5100 5101VP_EXPORT void 5102VpToString(Real *a, char *psz, size_t fFmt, int fPlus) 5103/* fPlus =0:default, =1: set ' ' before digits , =2:set '+' before digits. */ 5104{ 5105 size_t i, n, ZeroSup; 5106 BDIGIT shift, m, e, nn; 5107 char *pszSav = psz; 5108 ssize_t ex; 5109 5110 if (VpToSpecialString(a, psz, fPlus)) return; 5111 5112 ZeroSup = 1; /* Flag not to print the leading zeros as 0.00xxxxEnn */ 5113 5114 if (VpGetSign(a) < 0) *psz++ = '-'; 5115 else if (fPlus == 1) *psz++ = ' '; 5116 else if (fPlus == 2) *psz++ = '+'; 5117 5118 *psz++ = '0'; 5119 *psz++ = '.'; 5120 n = a->Prec; 5121 for(i=0;i < n;++i) { 5122 m = BASE1; 5123 e = a->frac[i]; 5124 while(m) { 5125 nn = e / m; 5126 if((!ZeroSup) || nn) { 5127 sprintf(psz, "%lu", (unsigned long)nn); /* The reading zero(s) */ 5128 psz += strlen(psz); 5129 /* as 0.00xx will be ignored. */ 5130 ZeroSup = 0; /* Set to print succeeding zeros */ 5131 } 5132 e = e - nn * m; 5133 m /= 10; 5134 } 5135 } 5136 ex = a->exponent * (ssize_t)BASE_FIG; 5137 shift = BASE1; 5138 while(a->frac[0] / shift == 0) { 5139 --ex; 5140 shift /= 10; 5141 } 5142 while(psz[-1]=='0') *(--psz) = 0; 5143 sprintf(psz, "E%"PRIdSIZE, ex); 5144 if(fFmt) VpFormatSt(pszSav, fFmt); 5145} 5146 5147VP_EXPORT void 5148VpToFString(Real *a, char *psz, size_t fFmt, int fPlus) 5149/* fPlus =0:default,=1: set ' ' before digits ,set '+' before digits. */ 5150{ 5151 size_t i, n; 5152 BDIGIT m, e, nn; 5153 char *pszSav = psz; 5154 ssize_t ex; 5155 5156 if(VpToSpecialString(a,psz,fPlus)) return; 5157 5158 if(VpGetSign(a) < 0) *psz++ = '-'; 5159 else if(fPlus==1) *psz++ = ' '; 5160 else if(fPlus==2) *psz++ = '+'; 5161 5162 n = a->Prec; 5163 ex = a->exponent; 5164 if(ex<=0) { 5165 *psz++ = '0';*psz++ = '.'; 5166 while(ex<0) { 5167 for(i=0;i<BASE_FIG;++i) *psz++ = '0'; 5168 ++ex; 5169 } 5170 ex = -1; 5171 } 5172 5173 for(i=0;i < n;++i) { 5174 --ex; 5175 if(i==0 && ex >= 0) { 5176 sprintf(psz, "%lu", (unsigned long)a->frac[i]); 5177 psz += strlen(psz); 5178 } else { 5179 m = BASE1; 5180 e = a->frac[i]; 5181 while(m) { 5182 nn = e / m; 5183 *psz++ = (char)(nn + '0'); 5184 e = e - nn * m; 5185 m /= 10; 5186 } 5187 } 5188 if(ex == 0) *psz++ = '.'; 5189 } 5190 while(--ex>=0) { 5191 m = BASE; 5192 while(m/=10) *psz++ = '0'; 5193 if(ex == 0) *psz++ = '.'; 5194 } 5195 *psz = 0; 5196 while(psz[-1]=='0') *(--psz) = 0; 5197 if(psz[-1]=='.') sprintf(psz, "0"); 5198 if(fFmt) VpFormatSt(pszSav, fFmt); 5199} 5200 5201/* 5202 * [Output] 5203 * a[] ... variable to be assigned the value. 5204 * [Input] 5205 * int_chr[] ... integer part(may include '+/-'). 5206 * ni ... number of characters in int_chr[],not including '+/-'. 5207 * frac[] ... fraction part. 5208 * nf ... number of characters in frac[]. 5209 * exp_chr[] ... exponent part(including '+/-'). 5210 * ne ... number of characters in exp_chr[],not including '+/-'. 5211 */ 5212VP_EXPORT int 5213VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne) 5214{ 5215 size_t i, j, ind_a, ma, mi, me; 5216 SIGNED_VALUE e, es, eb, ef; 5217 int sign, signe, exponent_overflow; 5218 5219 /* get exponent part */ 5220 e = 0; 5221 ma = a->MaxPrec; 5222 mi = ni; 5223 me = ne; 5224 signe = 1; 5225 exponent_overflow = 0; 5226 memset(a->frac, 0, ma * sizeof(BDIGIT)); 5227 if (ne > 0) { 5228 i = 0; 5229 if (exp_chr[0] == '-') { 5230 signe = -1; 5231 ++i; 5232 ++me; 5233 } 5234 else if (exp_chr[0] == '+') { 5235 ++i; 5236 ++me; 5237 } 5238 while (i < me) { 5239 es = e * (SIGNED_VALUE)BASE_FIG; 5240 e = e * 10 + exp_chr[i] - '0'; 5241 if (es > (SIGNED_VALUE)(e*BASE_FIG)) { 5242 exponent_overflow = 1; 5243 e = es; /* keep sign */ 5244 break; 5245 } 5246 ++i; 5247 } 5248 } 5249 5250 /* get integer part */ 5251 i = 0; 5252 sign = 1; 5253 if(1 /*ni >= 0*/) { 5254 if(int_chr[0] == '-') { 5255 sign = -1; 5256 ++i; 5257 ++mi; 5258 } else if(int_chr[0] == '+') { 5259 ++i; 5260 ++mi; 5261 } 5262 } 5263 5264 e = signe * e; /* e: The value of exponent part. */ 5265 e = e + ni; /* set actual exponent size. */ 5266 5267 if (e > 0) signe = 1; 5268 else signe = -1; 5269 5270 /* Adjust the exponent so that it is the multiple of BASE_FIG. */ 5271 j = 0; 5272 ef = 1; 5273 while (ef) { 5274 if (e >= 0) eb = e; 5275 else eb = -e; 5276 ef = eb / (SIGNED_VALUE)BASE_FIG; 5277 ef = eb - ef * (SIGNED_VALUE)BASE_FIG; 5278 if (ef) { 5279 ++j; /* Means to add one more preceeding zero */ 5280 ++e; 5281 } 5282 } 5283 5284 eb = e / (SIGNED_VALUE)BASE_FIG; 5285 5286 if (exponent_overflow) { 5287 int zero = 1; 5288 for ( ; i < mi && zero; i++) zero = int_chr[i] == '0'; 5289 for (i = 0; i < nf && zero; i++) zero = frac[i] == '0'; 5290 if (!zero && signe > 0) { 5291 VpSetInf(a, sign); 5292 VpException(VP_EXCEPTION_INFINITY, "exponent overflow",0); 5293 } 5294 else VpSetZero(a, sign); 5295 return 1; 5296 } 5297 5298 ind_a = 0; 5299 while (i < mi) { 5300 a->frac[ind_a] = 0; 5301 while ((j < BASE_FIG) && (i < mi)) { 5302 a->frac[ind_a] = a->frac[ind_a] * 10 + int_chr[i] - '0'; 5303 ++j; 5304 ++i; 5305 } 5306 if (i < mi) { 5307 ++ind_a; 5308 if (ind_a >= ma) goto over_flow; 5309 j = 0; 5310 } 5311 } 5312 5313 /* get fraction part */ 5314 5315 i = 0; 5316 while(i < nf) { 5317 while((j < BASE_FIG) && (i < nf)) { 5318 a->frac[ind_a] = a->frac[ind_a] * 10 + frac[i] - '0'; 5319 ++j; 5320 ++i; 5321 } 5322 if(i < nf) { 5323 ++ind_a; 5324 if(ind_a >= ma) goto over_flow; 5325 j = 0; 5326 } 5327 } 5328 goto Final; 5329 5330over_flow: 5331 rb_warn("Conversion from String to BigDecimal overflow (last few digits discarded)."); 5332 5333Final: 5334 if (ind_a >= ma) ind_a = ma - 1; 5335 while (j < BASE_FIG) { 5336 a->frac[ind_a] = a->frac[ind_a] * 10; 5337 ++j; 5338 } 5339 a->Prec = ind_a + 1; 5340 a->exponent = eb; 5341 VpSetSign(a,sign); 5342 VpNmlz(a); 5343 return 1; 5344} 5345 5346/* 5347 * [Input] 5348 * *m ... Real 5349 * [Output] 5350 * *d ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig. 5351 * *e ... exponent of m. 5352 * DBLE_FIG ... Number of digits in a double variable. 5353 * 5354 * m -> d*10**e, 0<d<BASE 5355 * [Returns] 5356 * 0 ... Zero 5357 * 1 ... Normal 5358 * 2 ... Infinity 5359 * -1 ... NaN 5360 */ 5361VP_EXPORT int 5362VpVtoD(double *d, SIGNED_VALUE *e, Real *m) 5363{ 5364 size_t ind_m, mm, fig; 5365 double div; 5366 int f = 1; 5367 5368 if(VpIsNaN(m)) { 5369 *d = VpGetDoubleNaN(); 5370 *e = 0; 5371 f = -1; /* NaN */ 5372 goto Exit; 5373 } else 5374 if(VpIsPosZero(m)) { 5375 *d = 0.0; 5376 *e = 0; 5377 f = 0; 5378 goto Exit; 5379 } else 5380 if(VpIsNegZero(m)) { 5381 *d = VpGetDoubleNegZero(); 5382 *e = 0; 5383 f = 0; 5384 goto Exit; 5385 } else 5386 if(VpIsPosInf(m)) { 5387 *d = VpGetDoublePosInf(); 5388 *e = 0; 5389 f = 2; 5390 goto Exit; 5391 } else 5392 if(VpIsNegInf(m)) { 5393 *d = VpGetDoubleNegInf(); 5394 *e = 0; 5395 f = 2; 5396 goto Exit; 5397 } 5398 /* Normal number */ 5399 fig =(DBLE_FIG + BASE_FIG - 1) / BASE_FIG; 5400 ind_m = 0; 5401 mm = Min(fig,(m->Prec)); 5402 *d = 0.0; 5403 div = 1.; 5404 while(ind_m < mm) { 5405 div /= (double)BASE; 5406 *d = *d + (double)m->frac[ind_m++] * div; 5407 } 5408 *e = m->exponent * (SIGNED_VALUE)BASE_FIG; 5409 *d *= VpGetSign(m); 5410 5411Exit: 5412#ifdef BIGDECIMAL_DEBUG 5413 if(gfDebug) { 5414 VPrint(stdout, " VpVtoD: m=%\n", m); 5415 printf(" d=%e * 10 **%ld\n", *d, *e); 5416 printf(" DBLE_FIG = %d\n", DBLE_FIG); 5417 } 5418#endif /*BIGDECIMAL_DEBUG */ 5419 return f; 5420} 5421 5422/* 5423 * m <- d 5424 */ 5425VP_EXPORT void 5426VpDtoV(Real *m, double d) 5427{ 5428 size_t ind_m, mm; 5429 SIGNED_VALUE ne; 5430 BDIGIT i; 5431 double val, val2; 5432 5433 if(isnan(d)) { 5434 VpSetNaN(m); 5435 goto Exit; 5436 } 5437 if(isinf(d)) { 5438 if(d>0.0) VpSetPosInf(m); 5439 else VpSetNegInf(m); 5440 goto Exit; 5441 } 5442 5443 if(d == 0.0) { 5444 VpSetZero(m,1); 5445 goto Exit; 5446 } 5447 val =(d > 0.) ? d :(-d); 5448 ne = 0; 5449 if(val >= 1.0) { 5450 while(val >= 1.0) { 5451 val /= (double)BASE; 5452 ++ne; 5453 } 5454 } else { 5455 val2 = 1.0 /(double)BASE; 5456 while(val < val2) { 5457 val *= (double)BASE; 5458 --ne; 5459 } 5460 } 5461 /* Now val = 0.xxxxx*BASE**ne */ 5462 5463 mm = m->MaxPrec; 5464 memset(m->frac, 0, mm * sizeof(BDIGIT)); 5465 for(ind_m = 0;val > 0.0 && ind_m < mm;ind_m++) { 5466 val *= (double)BASE; 5467 i = (BDIGIT)val; 5468 val -= (double)i; 5469 m->frac[ind_m] = i; 5470 } 5471 if(ind_m >= mm) ind_m = mm - 1; 5472 VpSetSign(m, (d > 0.0) ? 1 : -1); 5473 m->Prec = ind_m + 1; 5474 m->exponent = ne; 5475 5476 VpInternalRound(m, 0, (m->Prec > 0) ? m->frac[m->Prec-1] : 0, 5477 (BDIGIT)(val*(double)BASE)); 5478 5479Exit: 5480#ifdef BIGDECIMAL_DEBUG 5481 if(gfDebug) { 5482 printf("VpDtoV d=%30.30e\n", d); 5483 VPrint(stdout, " m=%\n", m); 5484 } 5485#endif /* BIGDECIMAL_DEBUG */ 5486 return; 5487} 5488 5489/* 5490 * m <- ival 5491 */ 5492#if 0 /* unused */ 5493VP_EXPORT void 5494VpItoV(Real *m, SIGNED_VALUE ival) 5495{ 5496 size_t mm, ind_m; 5497 size_t val, v1, v2, v; 5498 int isign; 5499 SIGNED_VALUE ne; 5500 5501 if(ival == 0) { 5502 VpSetZero(m,1); 5503 goto Exit; 5504 } 5505 isign = 1; 5506 val = ival; 5507 if(ival < 0) { 5508 isign = -1; 5509 val =(size_t)(-ival); 5510 } 5511 ne = 0; 5512 ind_m = 0; 5513 mm = m->MaxPrec; 5514 while(ind_m < mm) { 5515 m->frac[ind_m] = 0; 5516 ++ind_m; 5517 } 5518 ind_m = 0; 5519 while(val > 0) { 5520 if(val) { 5521 v1 = val; 5522 v2 = 1; 5523 while(v1 >= BASE) { 5524 v1 /= BASE; 5525 v2 *= BASE; 5526 } 5527 val = val - v2 * v1; 5528 v = v1; 5529 } else { 5530 v = 0; 5531 } 5532 m->frac[ind_m] = v; 5533 ++ind_m; 5534 ++ne; 5535 } 5536 m->Prec = ind_m - 1; 5537 m->exponent = ne; 5538 VpSetSign(m,isign); 5539 VpNmlz(m); 5540 5541Exit: 5542#ifdef BIGDECIMAL_DEBUG 5543 if(gfDebug) { 5544 printf(" VpItoV i=%d\n", ival); 5545 VPrint(stdout, " m=%\n", m); 5546 } 5547#endif /* BIGDECIMAL_DEBUG */ 5548 return; 5549} 5550#endif 5551 5552/* 5553 * y = SQRT(x), y*y - x =>0 5554 */ 5555VP_EXPORT int 5556VpSqrt(Real *y, Real *x) 5557{ 5558 Real *f = NULL; 5559 Real *r = NULL; 5560 size_t y_prec; 5561 SIGNED_VALUE n, e; 5562 SIGNED_VALUE prec; 5563 ssize_t nr; 5564 double val; 5565 5566 /* Zero, NaN or Infinity ? */ 5567 if(!VpHasVal(x)) { 5568 if(VpIsZero(x)||VpGetSign(x)>0) { 5569 VpAsgn(y,x,1); 5570 goto Exit; 5571 } 5572 VpSetNaN(y); 5573 return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(NaN or negative value)",0); 5574 goto Exit; 5575 } 5576 5577 /* Negative ? */ 5578 if(VpGetSign(x) < 0) { 5579 VpSetNaN(y); 5580 return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(negative value)",0); 5581 } 5582 5583 /* One ? */ 5584 if(VpIsOne(x)) { 5585 VpSetOne(y); 5586 goto Exit; 5587 } 5588 5589 n = (SIGNED_VALUE)y->MaxPrec; 5590 if (x->MaxPrec > (size_t)n) n = (ssize_t)x->MaxPrec; 5591 /* allocate temporally variables */ 5592 f = VpAlloc(y->MaxPrec * (BASE_FIG + 2), "#1"); 5593 r = VpAlloc((n + n) * (BASE_FIG + 2), "#1"); 5594 5595 nr = 0; 5596 y_prec = y->MaxPrec; 5597 5598 prec = x->exponent - (ssize_t)y_prec; 5599 if (x->exponent > 0) 5600 ++prec; 5601 else 5602 --prec; 5603 5604 VpVtoD(&val, &e, x); /* val <- x */ 5605 e /= (SIGNED_VALUE)BASE_FIG; 5606 n = e / 2; 5607 if (e - n * 2 != 0) { 5608 val /= BASE; 5609 n = (e + 1) / 2; 5610 } 5611 VpDtoV(y, sqrt(val)); /* y <- sqrt(val) */ 5612 y->exponent += n; 5613 n = (SIGNED_VALUE)((DBLE_FIG + BASE_FIG - 1) / BASE_FIG); 5614 y->MaxPrec = Min((size_t)n , y_prec); 5615 f->MaxPrec = y->MaxPrec + 1; 5616 n = (SIGNED_VALUE)(y_prec * BASE_FIG); 5617 if (n < (SIGNED_VALUE)maxnr) n = (SIGNED_VALUE)maxnr; 5618 do { 5619 y->MaxPrec *= 2; 5620 if (y->MaxPrec > y_prec) y->MaxPrec = y_prec; 5621 f->MaxPrec = y->MaxPrec; 5622 VpDivd(f, r, x, y); /* f = x/y */ 5623 VpAddSub(r, f, y, -1); /* r = f - y */ 5624 VpMult(f, VpPt5, r); /* f = 0.5*r */ 5625 if(VpIsZero(f)) goto converge; 5626 VpAddSub(r, f, y, 1); /* r = y + f */ 5627 VpAsgn(y, r, 1); /* y = r */ 5628 if(f->exponent <= prec) goto converge; 5629 } while(++nr < n); 5630 /* */ 5631#ifdef BIGDECIMAL_DEBUG 5632 if(gfDebug) { 5633 printf("ERROR(VpSqrt): did not converge within %ld iterations.\n", 5634 nr); 5635 } 5636#endif /* BIGDECIMAL_DEBUG */ 5637 y->MaxPrec = y_prec; 5638 5639converge: 5640 VpChangeSign(y, 1); 5641#ifdef BIGDECIMAL_DEBUG 5642 if(gfDebug) { 5643 VpMult(r, y, y); 5644 VpAddSub(f, x, r, -1); 5645 printf("VpSqrt: iterations = %"PRIdSIZE"\n", nr); 5646 VPrint(stdout, " y =% \n", y); 5647 VPrint(stdout, " x =% \n", x); 5648 VPrint(stdout, " x-y*y = % \n", f); 5649 } 5650#endif /* BIGDECIMAL_DEBUG */ 5651 y->MaxPrec = y_prec; 5652 5653Exit: 5654 VpFree(f); 5655 VpFree(r); 5656 return 1; 5657} 5658 5659/* 5660 * 5661 * nf: digit position for operation. 5662 * 5663 */ 5664VP_EXPORT int 5665VpMidRound(Real *y, unsigned short f, ssize_t nf) 5666/* 5667 * Round reletively from the decimal point. 5668 * f: rounding mode 5669 * nf: digit location to round from the the decimal point. 5670 */ 5671{ 5672 /* fracf: any positive digit under rounding position? */ 5673 /* fracf_1further: any positive digits under one further than the rounding position? */ 5674 /* exptoadd: number of digits needed to compensate negative nf */ 5675 int fracf, fracf_1further; 5676 ssize_t n,i,ix,ioffset, exptoadd; 5677 BDIGIT v, shifter; 5678 BDIGIT div; 5679 5680 nf += y->exponent * (ssize_t)BASE_FIG; 5681 exptoadd=0; 5682 if (nf < 0) { 5683 /* rounding position too left(large). */ 5684 if((f!=VP_ROUND_CEIL) && (f!=VP_ROUND_FLOOR)) { 5685 VpSetZero(y,VpGetSign(y)); /* truncate everything */ 5686 return 0; 5687 } 5688 exptoadd = -nf; 5689 nf = 0; 5690 } 5691 5692 ix = nf / (ssize_t)BASE_FIG; 5693 if ((size_t)ix >= y->Prec) return 0; /* rounding position too right(small). */ 5694 v = y->frac[ix]; 5695 5696 ioffset = nf - ix*(ssize_t)BASE_FIG; 5697 n = (ssize_t)BASE_FIG - ioffset - 1; 5698 for (shifter=1,i=0; i<n; ++i) shifter *= 10; 5699 5700 /* so the representation used (in y->frac) is an array of BDIGIT, where 5701 each BDIGIT contains a value between 0 and BASE-1, consisting of BASE_FIG 5702 decimal places. 5703 5704 (that numbers of decimal places are typed as ssize_t is somewhat confusing) 5705 5706 nf is now position (in decimal places) of the digit from the start of 5707 the array. 5708 ix is the position (in BDIGITS) of the BDIGIT containing the decimal digit, 5709 from the start of the array. 5710 v is the value of this BDIGIT 5711 ioffset is the number of extra decimal places along of this decimal digit 5712 within v. 5713 n is the number of decimal digits remaining within v after this decimal digit 5714 shifter is 10**n, 5715 v % shifter are the remaining digits within v 5716 v % (shifter * 10) are the digit together with the remaining digits within v 5717 v / shifter are the digit's predecessors together with the digit 5718 div = v / shifter / 10 is just the digit's precessors 5719 (v / shifter) - div*10 is just the digit, which is what v ends up being reassigned to. 5720 */ 5721 5722 fracf = (v % (shifter * 10) > 0); 5723 fracf_1further = ((v % shifter) > 0); 5724 5725 v /= shifter; 5726 div = v / 10; 5727 v = v - div*10; 5728 /* now v is just the digit required. 5729 now fracf is whether the digit or any of the remaining digits within v are non-zero 5730 now fracf_1further is whether any of the remaining digits within v are non-zero 5731 */ 5732 5733 /* now check all the remaining BDIGITS for zero-ness a whole BDIGIT at a time. 5734 if we spot any non-zeroness, that means that we foudn a positive digit under 5735 rounding position, and we also found a positive digit under one further than 5736 the rounding position, so both searches (to see if any such non-zero digit exists) 5737 can stop */ 5738 5739 for (i=ix+1; (size_t)i < y->Prec; i++) { 5740 if (y->frac[i] % BASE) { 5741 fracf = fracf_1further = 1; 5742 break; 5743 } 5744 } 5745 5746 /* now fracf = does any positive digit exist under the rounding position? 5747 now fracf_1further = does any positive digit exist under one further than the 5748 rounding position? 5749 now v = the first digit under the rounding position */ 5750 5751 /* drop digits after pointed digit */ 5752 memset(y->frac+ix+1, 0, (y->Prec - (ix+1)) * sizeof(BDIGIT)); 5753 5754 switch(f) { 5755 case VP_ROUND_DOWN: /* Truncate */ 5756 break; 5757 case VP_ROUND_UP: /* Roundup */ 5758 if (fracf) ++div; 5759 break; 5760 case VP_ROUND_HALF_UP: 5761 if (v>=5) ++div; 5762 break; 5763 case VP_ROUND_HALF_DOWN: 5764 if (v > 5 || (v == 5 && fracf_1further)) ++div; 5765 break; 5766 case VP_ROUND_CEIL: 5767 if (fracf && (VpGetSign(y)>0)) ++div; 5768 break; 5769 case VP_ROUND_FLOOR: 5770 if (fracf && (VpGetSign(y)<0)) ++div; 5771 break; 5772 case VP_ROUND_HALF_EVEN: /* Banker's rounding */ 5773 if (v > 5) ++div; 5774 else if (v == 5) { 5775 if (fracf_1further) { 5776 ++div; 5777 } 5778 else { 5779 if (ioffset == 0) { 5780 /* v is the first decimal digit of its BDIGIT; 5781 need to grab the previous BDIGIT if present 5782 to check for evenness of the previous decimal 5783 digit (which is same as that of the BDIGIT since 5784 base 10 has a factor of 2) */ 5785 if (ix && (y->frac[ix-1] % 2)) ++div; 5786 } 5787 else { 5788 if (div % 2) ++div; 5789 } 5790 } 5791 } 5792 break; 5793 } 5794 for (i=0; i<=n; ++i) div *= 10; 5795 if (div>=BASE) { 5796 if(ix) { 5797 y->frac[ix] = 0; 5798 VpRdup(y,ix); 5799 } else { 5800 short s = VpGetSign(y); 5801 SIGNED_VALUE e = y->exponent; 5802 VpSetOne(y); 5803 VpSetSign(y, s); 5804 y->exponent = e+1; 5805 } 5806 } else { 5807 y->frac[ix] = div; 5808 VpNmlz(y); 5809 } 5810 if (exptoadd > 0) { 5811 y->exponent += (SIGNED_VALUE)(exptoadd/BASE_FIG); 5812 exptoadd %= (ssize_t)BASE_FIG; 5813 for(i=0;i<exptoadd;i++) { 5814 y->frac[0] *= 10; 5815 if (y->frac[0] >= BASE) { 5816 y->frac[0] /= BASE; 5817 y->exponent++; 5818 } 5819 } 5820 } 5821 return 1; 5822} 5823 5824VP_EXPORT int 5825VpLeftRound(Real *y, unsigned short f, ssize_t nf) 5826/* 5827 * Round from the left hand side of the digits. 5828 */ 5829{ 5830 BDIGIT v; 5831 if (!VpHasVal(y)) return 0; /* Unable to round */ 5832 v = y->frac[0]; 5833 nf -= VpExponent(y)*(ssize_t)BASE_FIG; 5834 while ((v /= 10) != 0) nf--; 5835 nf += (ssize_t)BASE_FIG-1; 5836 return VpMidRound(y,f,nf); 5837} 5838 5839VP_EXPORT int 5840VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t nf) 5841{ 5842 /* First,assign whole value in truncation mode */ 5843 if (VpAsgn(y, x, 10) <= 1) return 0; /* Zero,NaN,or Infinity */ 5844 return VpMidRound(y,f,nf); 5845} 5846 5847static int 5848VpLimitRound(Real *c, size_t ixDigit) 5849{ 5850 size_t ix = VpGetPrecLimit(); 5851 if(!VpNmlz(c)) return -1; 5852 if(!ix) return 0; 5853 if(!ixDigit) ixDigit = c->Prec-1; 5854 if((ix+BASE_FIG-1)/BASE_FIG > ixDigit+1) return 0; 5855 return VpLeftRound(c, VpGetRoundMode(), (ssize_t)ix); 5856} 5857 5858/* If I understand correctly, this is only ever used to round off the final decimal 5859 digit of precision */ 5860static void 5861VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v) 5862{ 5863 int f = 0; 5864 5865 unsigned short const rounding_mode = VpGetRoundMode(); 5866 5867 if (VpLimitRound(c, ixDigit)) return; 5868 if (!v) return; 5869 5870 v /= BASE1; 5871 switch (rounding_mode) { 5872 case VP_ROUND_DOWN: 5873 break; 5874 case VP_ROUND_UP: 5875 if (v) f = 1; 5876 break; 5877 case VP_ROUND_HALF_UP: 5878 if (v >= 5) f = 1; 5879 break; 5880 case VP_ROUND_HALF_DOWN: 5881 /* this is ok - because this is the last digit of precision, 5882 the case where v == 5 and some further digits are nonzero 5883 will never occur */ 5884 if (v >= 6) f = 1; 5885 break; 5886 case VP_ROUND_CEIL: 5887 if (v && (VpGetSign(c) > 0)) f = 1; 5888 break; 5889 case VP_ROUND_FLOOR: 5890 if (v && (VpGetSign(c) < 0)) f = 1; 5891 break; 5892 case VP_ROUND_HALF_EVEN: /* Banker's rounding */ 5893 /* as per VP_ROUND_HALF_DOWN, because this is the last digit of precision, 5894 there is no case to worry about where v == 5 and some further digits are nonzero */ 5895 if (v > 5) f = 1; 5896 else if (v == 5 && vPrev % 2) f = 1; 5897 break; 5898 } 5899 if (f) { 5900 VpRdup(c, ixDigit); 5901 VpNmlz(c); 5902 } 5903} 5904 5905/* 5906 * Rounds up m(plus one to final digit of m). 5907 */ 5908static int 5909VpRdup(Real *m, size_t ind_m) 5910{ 5911 BDIGIT carry; 5912 5913 if (!ind_m) ind_m = m->Prec; 5914 5915 carry = 1; 5916 while (carry > 0 && (ind_m--)) { 5917 m->frac[ind_m] += carry; 5918 if (m->frac[ind_m] >= BASE) m->frac[ind_m] -= BASE; 5919 else carry = 0; 5920 } 5921 if(carry > 0) { /* Overflow,count exponent and set fraction part be 1 */ 5922 if (!AddExponent(m, 1)) return 0; 5923 m->Prec = m->frac[0] = 1; 5924 } else { 5925 VpNmlz(m); 5926 } 5927 return 1; 5928} 5929 5930/* 5931 * y = x - fix(x) 5932 */ 5933VP_EXPORT void 5934VpFrac(Real *y, Real *x) 5935{ 5936 size_t my, ind_y, ind_x; 5937 5938 if(!VpHasVal(x)) { 5939 VpAsgn(y,x,1); 5940 goto Exit; 5941 } 5942 5943 if (x->exponent > 0 && (size_t)x->exponent >= x->Prec) { 5944 VpSetZero(y,VpGetSign(x)); 5945 goto Exit; 5946 } 5947 else if(x->exponent <= 0) { 5948 VpAsgn(y, x, 1); 5949 goto Exit; 5950 } 5951 5952 /* satisfy: x->exponent > 0 */ 5953 5954 y->Prec = x->Prec - (size_t)x->exponent; 5955 y->Prec = Min(y->Prec, y->MaxPrec); 5956 y->exponent = 0; 5957 VpSetSign(y,VpGetSign(x)); 5958 ind_y = 0; 5959 my = y->Prec; 5960 ind_x = x->exponent; 5961 while(ind_y < my) { 5962 y->frac[ind_y] = x->frac[ind_x]; 5963 ++ind_y; 5964 ++ind_x; 5965 } 5966 VpNmlz(y); 5967 5968Exit: 5969#ifdef BIGDECIMAL_DEBUG 5970 if(gfDebug) { 5971 VPrint(stdout, "VpFrac y=%\n", y); 5972 VPrint(stdout, " x=%\n", x); 5973 } 5974#endif /* BIGDECIMAL_DEBUG */ 5975 return; 5976} 5977 5978/* 5979 * y = x ** n 5980 */ 5981VP_EXPORT int 5982VpPower(Real *y, Real *x, SIGNED_VALUE n) 5983{ 5984 size_t s, ss; 5985 ssize_t sign; 5986 Real *w1 = NULL; 5987 Real *w2 = NULL; 5988 5989 if(VpIsZero(x)) { 5990 if(n==0) { 5991 VpSetOne(y); 5992 goto Exit; 5993 } 5994 sign = VpGetSign(x); 5995 if(n<0) { 5996 n = -n; 5997 if(sign<0) sign = (n%2)?(-1):(1); 5998 VpSetInf (y,sign); 5999 } else { 6000 if(sign<0) sign = (n%2)?(-1):(1); 6001 VpSetZero(y,sign); 6002 } 6003 goto Exit; 6004 } 6005 if(VpIsNaN(x)) { 6006 VpSetNaN(y); 6007 goto Exit; 6008 } 6009 if(VpIsInf(x)) { 6010 if(n==0) { 6011 VpSetOne(y); 6012 goto Exit; 6013 } 6014 if(n>0) { 6015 VpSetInf(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1); 6016 goto Exit; 6017 } 6018 VpSetZero(y, (n%2==0 || VpIsPosInf(x)) ? 1 : -1); 6019 goto Exit; 6020 } 6021 6022 if((x->exponent == 1) &&(x->Prec == 1) &&(x->frac[0] == 1)) { 6023 /* abs(x) = 1 */ 6024 VpSetOne(y); 6025 if(VpGetSign(x) > 0) goto Exit; 6026 if((n % 2) == 0) goto Exit; 6027 VpSetSign(y, -1); 6028 goto Exit; 6029 } 6030 6031 if(n > 0) sign = 1; 6032 else if(n < 0) { 6033 sign = -1; 6034 n = -n; 6035 } else { 6036 VpSetOne(y); 6037 goto Exit; 6038 } 6039 6040 /* Allocate working variables */ 6041 6042 w1 = VpAlloc((y->MaxPrec + 2) * BASE_FIG, "#0"); 6043 w2 = VpAlloc((w1->MaxPrec * 2 + 1) * BASE_FIG, "#0"); 6044 /* calculation start */ 6045 6046 VpAsgn(y, x, 1); 6047 --n; 6048 while(n > 0) { 6049 VpAsgn(w1, x, 1); 6050 s = 1; 6051 while (ss = s, (s += s) <= (size_t)n) { 6052 VpMult(w2, w1, w1); 6053 VpAsgn(w1, w2, 1); 6054 } 6055 n -= (SIGNED_VALUE)ss; 6056 VpMult(w2, y, w1); 6057 VpAsgn(y, w2, 1); 6058 } 6059 if(sign < 0) { 6060 VpDivd(w1, w2, VpConstOne, y); 6061 VpAsgn(y, w1, 1); 6062 } 6063 6064Exit: 6065#ifdef BIGDECIMAL_DEBUG 6066 if(gfDebug) { 6067 VPrint(stdout, "VpPower y=%\n", y); 6068 VPrint(stdout, "VpPower x=%\n", x); 6069 printf(" n=%d\n", n); 6070 } 6071#endif /* BIGDECIMAL_DEBUG */ 6072 VpFree(w2); 6073 VpFree(w1); 6074 return 1; 6075} 6076 6077#ifdef BIGDECIMAL_DEBUG 6078int 6079VpVarCheck(Real * v) 6080/* 6081 * Checks the validity of the Real variable v. 6082 * [Input] 6083 * v ... Real *, variable to be checked. 6084 * [Returns] 6085 * 0 ... correct v. 6086 * other ... error 6087 */ 6088{ 6089 size_t i; 6090 6091 if(v->MaxPrec <= 0) { 6092 printf("ERROR(VpVarCheck): Illegal Max. Precision(=%"PRIuSIZE")\n", 6093 v->MaxPrec); 6094 return 1; 6095 } 6096 if((v->Prec <= 0) ||((v->Prec) >(v->MaxPrec))) { 6097 printf("ERROR(VpVarCheck): Illegal Precision(=%"PRIuSIZE")\n", v->Prec); 6098 printf(" Max. Prec.=%"PRIuSIZE"\n", v->MaxPrec); 6099 return 2; 6100 } 6101 for(i = 0; i < v->Prec; ++i) { 6102 if((v->frac[i] >= BASE)) { 6103 printf("ERROR(VpVarCheck): Illegal fraction\n"); 6104 printf(" Frac[%"PRIuSIZE"]=%lu\n", i, v->frac[i]); 6105 printf(" Prec. =%"PRIuSIZE"\n", v->Prec); 6106 printf(" Exp. =%"PRIdVALUE"\n", v->exponent); 6107 printf(" BASE =%lu\n", BASE); 6108 return 3; 6109 } 6110 } 6111 return 0; 6112} 6113#endif /* BIGDECIMAL_DEBUG */ 6114