File indexing completed on 2024-05-12 17:22:36
0001 /* floatnum.c: Arbitrary precision floating point numbers, based on bc. */ 0002 /* 0003 Copyright (C) 2007, 2008 Wolf Lammen. 0004 0005 This program is free software; you can redistribute it and/or modify 0006 it under the terms of the GNU General Public License as published by 0007 the Free Software Foundation; either version 2 of the License , or 0008 (at your option) any later version. 0009 0010 This program is distributed in the hope that it will be useful, 0011 but WITHOUT ANY WARRANTY; without even the implied warranty of 0012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0013 GNU General Public License for more details. 0014 0015 You should have received a copy of the GNU General Public License 0016 along with this program; see the file COPYING. If not, write to: 0017 0018 The Free Software Foundation, Inc. 0019 59 Temple Place, Suite 330 0020 Boston, MA 02111-1307 USA. 0021 0022 0023 You may contact the author by: 0024 e-mail: ookami1 <at> gmx <dot> de 0025 mail: Wolf Lammen 0026 Oertzweg 45 0027 22307 Hamburg 0028 Germany 0029 0030 *************************************************************************/ 0031 0032 /* a floating point engine based on bc's decimal 0033 fix point arithmetic. The speed is not overwhelming, 0034 but sufficient for calculators with limited demands. 0035 As bc's number.c is a portable engine, this should be 0036 portable as well. 0037 */ 0038 0039 #include "floatnum.h" 0040 #include "floatlong.h" 0041 #include <stdio.h> 0042 #include <string.h> 0043 0044 #define NOSPECIALVALUE 1 0045 0046 int maxdigits = MAXDIGITS; 0047 0048 Error float_error; 0049 int expmax = EXPMAX; 0050 int expmin = EXPMIN; 0051 0052 /* general helper routines */ 0053 0054 static int 0055 _max(int x, int y) 0056 { 0057 return x > y? x : y; 0058 } 0059 0060 static int 0061 _min(int x, int y) 0062 { 0063 return x < y? x : y; 0064 } 0065 0066 /* the return value points to the first character different 0067 from accept. */ 0068 const char* 0069 _memskip( 0070 const char* buf, 0071 const char* end, 0072 char accept) 0073 { 0074 for(--buf; ++buf != end && *buf == accept;); 0075 return buf; 0076 } 0077 0078 /* scans a value in a bc-num (or part of it) for the first 0079 occurence of a digit *different* from <digit>. The scan 0080 is limited to <count> bytes. Returns the offset of the 0081 matching byte, or <count>, if none was found */ 0082 static int 0083 _scan_digit( 0084 const char*p, 0085 int count, 0086 char digit) 0087 { 0088 const char* ps; 0089 0090 ps = p; 0091 for (; count-- > 0 && *p == digit; ++p); 0092 return p - ps; 0093 } 0094 0095 /* bc_num primitives */ 0096 0097 #define _scaleof(f) (f->significand->n_scale) 0098 #define _lenof(f) (f->significand->n_len) 0099 #define _valueof(f) (f->significand->n_value) 0100 #define _digit(f, offset) (*(_valueof(f) + offset)) 0101 #define _setscale(f, value) (_scaleof(f) = value) 0102 0103 /* converts floatnum's sign encodings (+1, -1) to bc_num 0104 sign encoding (PLUS, MINUS) */ 0105 static void 0106 _setsign( 0107 floatnum f, 0108 signed char value) 0109 { 0110 f->significand->n_sign = value < 0? MINUS : PLUS; 0111 #ifdef FLOATDEBUG 0112 f->value[0] = value < 0? '-' : '+'; 0113 #endif 0114 } 0115 0116 /* modifies a bc_num significand such that the last <count> 0117 (existing) digits of its value are not visible 0118 to bc_num operations any more. 0119 A negative <count> reverts this operation (unhide). 0120 pre: <count> <= n_scale*/ 0121 static void 0122 _hidelast( 0123 floatnum f, 0124 int count) 0125 { 0126 f->significand->n_scale -= count; 0127 } 0128 0129 /* modifies a bc_num significand such that the first <count> 0130 (existing) digits of its value are not visible 0131 to bc_num operations any more. 0132 A negative <count> reverts this operation (unhide). 0133 pre: <count> <= n_scale*/ 0134 static void 0135 _hidefirst( 0136 floatnum f, 0137 int count) 0138 { 0139 f->significand->n_value += count; 0140 f->significand->n_scale -= count; 0141 } 0142 0143 /* Usually, significands are normalized, which means they 0144 fulfill 1 <= x < 10. This operation moves the decimal 0145 point <count> places to the right, effectively multiplying 0146 the significand by a power of 10. A negative <count> 0147 reverts such an operation. 0148 pre: <count> < n_scale */ 0149 static void 0150 _movepoint( 0151 floatnum f, 0152 int digits) 0153 { 0154 f->significand->n_len += digits; 0155 f->significand->n_scale -= digits; 0156 } 0157 0158 /* floatstruct primitives */ 0159 0160 /* a quick check for NaN and 0 */ 0161 static char 0162 _is_special( 0163 cfloatnum f) 0164 { 0165 return f->significand == NULL; 0166 } 0167 0168 /* creates a shallow working copy of <source> in <dest>, 0169 using b as a container for the significand. 0170 On return, <dest> is equal in value to <source>. 0171 <dest> may then be modified freely in the course 0172 of an operation, without effecting source, 0173 *except* that the digits in *(dest->significand->n_value) 0174 have to be retained (or restored). 0175 <dest> must *never* be the destination of a 0176 float_xxx operation, in particiular, it must not be 0177 freed. Neither must b be modified (or freed) in a bc_num 0178 operation */ 0179 static void 0180 _copyfn( 0181 floatnum dest, 0182 cfloatnum source, 0183 bc_num b) 0184 { 0185 *dest = *source; 0186 if (!_is_special(source)) 0187 { 0188 *b = *(source->significand); 0189 b->n_ptr = 0; 0190 dest->significand = b; 0191 } 0192 } 0193 0194 /* If you want to execute a bc_num operation to a limited scale, 0195 it is a waste of computation time to pass operands 0196 with a longer scale, because bc lets the operand's 0197 scale override your limit. This function hides superfluous 0198 digits from bc, returning the original scale for restoring 0199 purposes */ 0200 static int 0201 _limit_scale( 0202 floatnum f, 0203 int newscale) 0204 { 0205 int oldscale; 0206 oldscale = _scaleof(f); 0207 _setscale(f, _min(oldscale, newscale)); 0208 return oldscale; 0209 } 0210 0211 /*============================ floatnum routines ===================*/ 0212 0213 int 0214 float_getrange() 0215 { 0216 return expmax; 0217 } 0218 0219 int 0220 float_getprecision() 0221 { 0222 return maxdigits; 0223 } 0224 0225 int 0226 float_setrange( 0227 int maxexp) 0228 { 0229 int result; 0230 0231 result = expmax; 0232 expmax = _max(_min(maxexp, MAXEXP), 1); 0233 expmin = -expmax - 1; 0234 return result; 0235 } 0236 0237 int 0238 float_setprecision( 0239 int digits) 0240 { 0241 int result; 0242 0243 result = maxdigits; 0244 maxdigits = _max(_min(digits, MAXDIGITS), 1); 0245 return result; 0246 } 0247 0248 /* checking the limits on exponents */ 0249 char 0250 float_isvalidexp( 0251 int exp) 0252 { 0253 return exp >= expmin && exp <= expmax; 0254 } 0255 0256 /* clears the error state as well */ 0257 Error 0258 float_geterror() 0259 { 0260 Error tmp; 0261 0262 tmp = float_error; 0263 float_error = Success; 0264 return tmp; 0265 } 0266 0267 /* the first error blocks all others as it may be the source 0268 of a cascade of dependent errors */ 0269 void 0270 float_seterror( 0271 Error code) 0272 { 0273 if (float_error == Success) 0274 float_error = code; 0275 } 0276 0277 void 0278 floatnum_init() 0279 { 0280 bc_init_numbers(); 0281 float_geterror(); 0282 } 0283 0284 void 0285 float_create( 0286 floatnum f) 0287 { 0288 f->significand = NULL; 0289 f->exponent = EXPNAN; 0290 #ifdef FLOATDEBUG 0291 memcpy(f->value, "NaN", 4); 0292 #endif 0293 } 0294 0295 void 0296 float_setnan ( 0297 floatnum f) 0298 { 0299 bc_free_num(&(f->significand)); 0300 float_create(f); 0301 } 0302 0303 char 0304 _setnan( 0305 floatnum result) 0306 { 0307 float_setnan(result); 0308 return FALSE; 0309 } 0310 0311 char 0312 _seterror( 0313 floatnum result, 0314 Error code) 0315 { 0316 float_seterror(code); 0317 return _setnan(result); 0318 } 0319 0320 void 0321 float_setzero ( 0322 floatnum f) 0323 { 0324 bc_free_num(&(f->significand)); 0325 f->exponent = EXPZERO; 0326 #ifdef FLOATDEBUG 0327 f->value[0] ='0'; 0328 f->value[1] = 0; 0329 #endif 0330 } 0331 0332 char 0333 _setzero( 0334 floatnum result) 0335 { 0336 float_setzero(result); 0337 return TRUE; 0338 } 0339 0340 char 0341 float_isnan( 0342 cfloatnum f) 0343 { 0344 return _is_special(f) && f->exponent != EXPZERO; 0345 } 0346 0347 char 0348 float_iszero( 0349 cfloatnum f) 0350 { 0351 return _is_special(f) && f->exponent == EXPZERO; 0352 } 0353 0354 int 0355 float_getlength( 0356 cfloatnum f) 0357 { 0358 return _is_special(f)? 0 : _scaleof(f) + 1; 0359 } 0360 0361 char 0362 float_getdigit( 0363 cfloatnum f, 0364 int ofs) 0365 { 0366 if (ofs >= float_getlength(f) || ofs < 0) 0367 return 0; 0368 return _digit(f, ofs); 0369 } 0370 0371 /* checks whether f is a NaN and sets the float_error 0372 variable accordingly. Used in parameter checks when 0373 float_xxx calls are executed. 0374 FALSE is returned if a NaN is encountered. */ 0375 char 0376 _checknan( 0377 cfloatnum f) 0378 { 0379 if (!float_isnan(f)) 0380 return TRUE; 0381 float_seterror(NoOperand); 0382 return FALSE; 0383 } 0384 0385 /* checks whether <digits> is positive and 0386 sets the float_error variable accordingly. Used in 0387 parameter checks when float_xxx calls are executed. 0388 Some operations accept a special value like EXACT, 0389 that has to pass this check, even though its numerical 0390 encoding violates the boundaries. 0391 If a function does not accept a special value, 0392 use NOSPECIALVALUE as a parameter for <specialval>. 0393 TRUE is returned if the check is passed. 0394 The check for the limit MAXDIGITS is not executed 0395 here, because some intermediate operations have to succeed 0396 on more than MAXDIGITS digits */ 0397 static char 0398 _checkdigits( 0399 int digits, 0400 int specialval) 0401 { 0402 if ((digits > 0 && digits <= maxdigits) || digits == specialval) 0403 return TRUE; 0404 float_seterror(InvalidPrecision); 0405 return FALSE; 0406 } 0407 0408 /* backward-scans the significand in <f>. Returns the number of 0409 digits equal to <digit> beginning with the <scale>+1-th digit of 0410 f->significand->n_value. */ 0411 static int 0412 _bscandigit( 0413 cfloatnum f, 0414 int scale, 0415 char digit) 0416 { 0417 char* p; 0418 char* ps; 0419 0420 ps = _valueof(f); 0421 for (p = ps + scale + 1; p-- != ps && *p == digit;); 0422 return scale - (p - ps); 0423 } 0424 0425 /* scans two significands for the first occurence 0426 of a pair of differnt digits. Returns the number 0427 of equal digits at the beginning */ 0428 static int 0429 _scan_equal( 0430 floatnum v1, 0431 floatnum v2) 0432 { 0433 int count, i; 0434 char* p1; 0435 char* p2; 0436 0437 count = _min(_scaleof(v1), _scaleof(v2)); 0438 p1 = _valueof(v1); 0439 p2 = _valueof(v2); 0440 i = 0; 0441 for (; *(p1++) == *(p2++) && ++i <= count;); 0442 return i; 0443 } 0444 0445 /* scans two significands until it finds a digit different 0446 from 0 in the first significand, or a digit different from 0447 9 in the second operand. The scan is limited by <count> 0448 compares and starts with the *second* digit in the 0449 significands. Returns the number of found (0,9) pairs. */ 0450 static int 0451 _scan_09pairs( 0452 floatnum f1, 0453 floatnum f2, 0454 int count) 0455 { 0456 char* p; 0457 char* p1; 0458 char* p2; 0459 0460 p1 = _valueof(f1) + 1; 0461 p2 = _valueof(f2) + 1; 0462 p = p1; 0463 for (; count-- > 0 && *p1 == 0 && *(p2++) == 9; ++p1); 0464 return p1 - p; 0465 } 0466 0467 signed char 0468 float_getsign( 0469 cfloatnum f) 0470 { 0471 if(_is_special(f)) 0472 return 0; 0473 return f->significand->n_sign == PLUS? 1 : -1; 0474 } 0475 0476 void 0477 float_setsign( 0478 floatnum f, 0479 signed char s) 0480 { 0481 if (s == 1 || s == -1) 0482 { 0483 if(!_is_special(f)) 0484 _setsign(f, s); 0485 } 0486 else if (s != 0) 0487 float_setnan(f); 0488 } 0489 0490 char 0491 float_neg( 0492 floatnum f) 0493 { 0494 float_setsign(f, -float_getsign(f)); 0495 return _checknan(f); 0496 } 0497 0498 char 0499 float_abs( 0500 floatnum f) 0501 { 0502 if(float_getsign(f) == -1) 0503 float_neg(f); 0504 return _checknan(f); 0505 } 0506 0507 signed char 0508 float_cmp( 0509 cfloatnum val1, 0510 cfloatnum val2) 0511 { 0512 signed char sgn1; 0513 0514 if (!_checknan(val1) || !_checknan(val2)) 0515 return UNORDERED; 0516 sgn1 = float_getsign(val1); 0517 if (float_getsign(val2) != sgn1) 0518 { 0519 if (sgn1 != 0) 0520 return sgn1; 0521 return -float_getsign(val2); 0522 } 0523 if (val1->exponent > val2->exponent) 0524 return sgn1; 0525 if (val1->exponent < val2->exponent) 0526 return -sgn1; 0527 if (_is_special(val1)) 0528 return 0; 0529 return (bc_compare(val1->significand, val2->significand)); 0530 } 0531 0532 /* normalizing process: 0533 hides leading zeros in a significand and corrects the 0534 exponent accordingly */ 0535 static void 0536 _corr_lead_zero( 0537 floatnum f) 0538 { 0539 int count; 0540 0541 count = _scan_digit(_valueof(f), float_getlength(f), 0); 0542 _hidefirst(f, count); 0543 f->exponent-=count; 0544 } 0545 0546 /* normalizing process: 0547 if the significand is > 10 in magnitude, this function 0548 corrects this */ 0549 static void 0550 _corr_overflow( 0551 floatnum f) 0552 { 0553 int shift; 0554 0555 shift = _lenof(f) - 1; 0556 _movepoint(f, -shift); 0557 f->exponent += shift; 0558 } 0559 0560 /* cuts off trailing zeros at the end of a significand */ 0561 static void 0562 _corr_trailing_zeros( 0563 floatnum f) 0564 { 0565 _hidelast(f, _bscandigit(f, _scaleof(f), 0)); 0566 } 0567 0568 static char hexdigits[] = "0123456789ABCDEF"; 0569 0570 int 0571 float_getsignificand( 0572 char* buf, 0573 int bufsz, 0574 cfloatnum f) 0575 { 0576 int idx, lg; 0577 0578 if (bufsz <= 0) 0579 return 0; 0580 if (float_isnan(f)) 0581 { 0582 *buf = 'N'; 0583 return 1; 0584 } 0585 if (float_iszero(f)) 0586 { 0587 *buf = '0'; 0588 return 1; 0589 } 0590 idx = -1; 0591 lg = _min(bufsz, float_getlength(f)); 0592 for(; ++idx < lg;) 0593 *(buf++) = hexdigits[(int)float_getdigit(f, idx)]; 0594 return lg; 0595 } 0596 0597 int 0598 float_getexponent( 0599 cfloatnum f) 0600 { 0601 if (_is_special(f)) 0602 return 0; 0603 return f->exponent; 0604 } 0605 0606 int float_getscientific( 0607 char* buf, 0608 int bufsz, 0609 cfloatnum f) 0610 { 0611 char b[42]; /* supports exponents encoded in up to 128 bits */ 0612 int sgnlg, explg, mlg; 0613 0614 /* handle special cases */ 0615 if(float_isnan(f)) 0616 { 0617 if (bufsz < 4) 0618 return -1; 0619 memcpy(buf, "NaN\0", 4); 0620 return 3; 0621 } 0622 if(float_iszero(f)) 0623 { 0624 if (bufsz < 2) 0625 return -1; 0626 *buf = '0'; 0627 *(buf+1) = '\0'; 0628 return 1; 0629 } 0630 0631 /* set one byte aside for sign? */ 0632 sgnlg = 0; 0633 if(float_getsign(f) < 0) 0634 sgnlg = 1; 0635 0636 /* convert the exponent */ 0637 sprintf(b, "%d", float_getexponent(f)); 0638 explg = strlen(b); 0639 0640 /* 3 extra bytes for dot, exp char and terminating \0 */ 0641 bufsz -= explg + sgnlg + 3; /* rest is for significand */ 0642 if (bufsz <= 0) 0643 /* buffer too small */ 0644 return-1; 0645 0646 if(sgnlg > 0) 0647 *(buf++) = '-'; 0648 0649 /* get the digit sequence of the significand, trailing zeros cut off */ 0650 mlg = float_getsignificand(++buf, bufsz, f) - 1; 0651 0652 /* move the first digit one byte to the front and fill 0653 the gap with a dot */ 0654 *(buf-1) = *buf; 0655 *(buf++) = '.'; 0656 0657 /* append the exponent */ 0658 *(buf+mlg) = 'e'; 0659 memcpy(buf+mlg+1, b, explg); 0660 0661 /* the trailing \0 */ 0662 *(buf+mlg+explg+1) = '\0'; 0663 return sgnlg + mlg + explg + 3; 0664 } 0665 0666 #ifdef FLOATDEBUG 0667 void _setvalue_(floatnum f) 0668 { 0669 f->value[float_getsignificand(f->value+2, sizeof(f->value)-3, f)+2] = 0; 0670 f->value[1] = f->value[2]; 0671 f->value[2] = '.'; 0672 f->value[0] = float_getsign(f) < 0? '-' : '+'; 0673 } 0674 #endif 0675 0676 int 0677 float_setsignificand( 0678 floatnum f, 0679 int* leadingzeros, 0680 const char* buf, 0681 int bufsz) 0682 { 0683 const char* p; 0684 const char* dot; 0685 const char* last; 0686 const char* b; 0687 char* bcp; 0688 int zeros; 0689 int lg; 0690 char c; 0691 0692 float_setnan(f); 0693 if (bufsz == NULLTERMINATED) 0694 bufsz = strlen(buf); 0695 0696 /* initialize the output parameters for all 0697 early out branches */ 0698 if (leadingzeros != NULL) 0699 *leadingzeros = 0; 0700 0701 if (bufsz <= 0) 0702 return -1; 0703 0704 dot = memchr(buf, '.', bufsz); 0705 /* do not accept more than 1 dots */ 0706 if (dot != NULL && memchr(dot + 1, '.', bufsz - (dot - buf)) != NULL) 0707 return -1; 0708 0709 last = buf + bufsz; /* points behind the input buffer */ 0710 0711 /* skip all leading zeros */ 0712 b = _memskip(buf, last, '0'); 0713 0714 /* is the first non-zero character found a dot? */ 0715 if (b == dot) 0716 /* then skip all zeros following the dot */ 0717 b = _memskip(b+1, last, '0'); 0718 0719 /* the 'leading zeros' */ 0720 zeros = b - buf - (dot == NULL || dot >= b? 0:1); 0721 0722 /* only zeros found? */ 0723 if (b == last) 0724 { 0725 /* indicate no dot, no leading zeros, because 0726 this does not matter in case of zero */ 0727 if (bufsz > (dot == NULL? 0:1)) 0728 float_setzero(f); 0729 /* do not accept a dot without any zero */ 0730 return -1; 0731 } 0732 0733 /* size of the rest buffer without leading zeros */ 0734 bufsz -= b - buf; 0735 0736 /* does the rest buffer contain a dot? */ 0737 lg = dot >= b && dot - b < maxdigits? 1 : 0; 0738 0739 /* points behind the last significant digit */ 0740 p = b + _min(maxdigits + lg, bufsz); 0741 0742 /* digits, limited by MAXDIGITS */ 0743 lg = _min(maxdigits, bufsz - lg); 0744 0745 /* reduce lg by the number of trailing zeros */ 0746 for (; *--p == '0'; --lg); 0747 if (*(p--) == '.') 0748 for (; *(p--) == '0'; --lg); 0749 0750 /* get a bc_num of sufficient size */ 0751 f->significand = bc_new_num(1, lg - 1); 0752 if (f->significand == NULL) 0753 return -1; 0754 0755 /* exponent is forced to 0 */ 0756 f->exponent = 0; 0757 0758 /* copy lg digits into bc_num buffer, 0759 scan the rest for invalid characters */ 0760 bcp = _valueof(f); 0761 for(; --bufsz >= 0;) 0762 { 0763 c = *(b++); 0764 if (c != '.') /* ignore a dot */ 0765 { 0766 if (c < '0' || c > '9') 0767 { 0768 /* invalid character */ 0769 float_setnan(f); 0770 return -1; 0771 } 0772 if (--lg >= 0) 0773 *(bcp++) = c - '0'; 0774 } 0775 } 0776 0777 if (leadingzeros != NULL) 0778 *leadingzeros = zeros; 0779 #ifdef FLOATDEBUG 0780 _setvalue_(f); 0781 #endif 0782 return dot == NULL? -1 : dot - buf; 0783 } 0784 0785 void 0786 float_setexponent( 0787 floatnum f, 0788 int exponent) 0789 { 0790 if (!_is_special(f)) 0791 { 0792 if (!float_isvalidexp(exponent)) 0793 float_setnan(f); 0794 else 0795 f->exponent = exponent; 0796 } 0797 } 0798 0799 void 0800 float_setscientific( 0801 floatnum f, 0802 const char* buf, 0803 int bufsz) 0804 { 0805 int exppos; 0806 int zeros; 0807 int dotpos; 0808 unsigned exp, ovfl; 0809 signed char expsign, msign; 0810 const char* expptr; 0811 const char* last; 0812 char c; 0813 0814 float_setnan(f); 0815 0816 if (bufsz == NULLTERMINATED) 0817 bufsz = strlen(buf); 0818 0819 /* find the offset of the exponent character, 0820 or -1, if not found */ 0821 for(exppos = bufsz; --exppos >= 0;) 0822 if ((c = *(buf+exppos)) == 'E' || c == 'e') 0823 break; 0824 0825 /* marks the end of the exponent string */ 0826 last = buf + bufsz; 0827 0828 /* pre-set exponent to +0, which is the right value, 0829 if there is no exponent. */ 0830 exp = 0; 0831 expsign = 1; 0832 ovfl = 0; 0833 if (exppos >= 0) 0834 { 0835 /* points behind the exponent character */ 0836 expptr = buf + (exppos + 1); 0837 if (expptr == last) 0838 /* do not accept an exponent char without an integer */ 0839 return; 0840 0841 /* get the exponent sign */ 0842 switch(*expptr) 0843 { 0844 case '-': 0845 expsign = -1; /* and fall through */ 0846 case '+': 0847 ++expptr; 0848 } 0849 if (expptr == last) 0850 /* do not accept a sign without a digit following */ 0851 return; 0852 0853 /* encode the sequence of digits into an unsignedeger */ 0854 for (;expptr != last ;) 0855 { 0856 if (*expptr < '0' || *expptr > '9') 0857 /* invalid char encountered */ 0858 return; 0859 ovfl = 10; 0860 if (_longmul(&exp, &ovfl)) 0861 { 0862 ovfl = *(expptr++) - '0'; 0863 _longadd(&exp, &ovfl); 0864 } 0865 if (ovfl != 0 || exp > EXPMAX+1) 0866 { 0867 /* do not return immediately, because the 0868 significand can be zero */ 0869 ovfl = 1; 0870 break; 0871 } 0872 } 0873 /* move the last pointer to the exponent char.*/ 0874 last = buf + exppos; 0875 } 0876 /* last points behind the significand part. 0877 exp is at most -EXPMIN */ 0878 0879 /* get the sign of the significand */ 0880 msign = 1; 0881 if (buf != last) 0882 switch(*buf) 0883 { 0884 case '-': 0885 msign = -1; /* fall through */ 0886 case '+': 0887 ++buf; 0888 } 0889 0890 /* let setsignificand convert the sequence of digits 0891 into a significand. If a dot is found, its position 0892 is given in dotpos, -1 otherwise. 0893 zeros are the count of leading '0' digits before 0894 the first non_zero digit. */ 0895 dotpos = float_setsignificand(f, &zeros, buf, last-buf); 0896 if (_is_special(f)) 0897 /* setsignificand either found a zero or encountered 0898 invalid characters */ 0899 return; 0900 0901 /* if we did not find a dot, we assume an integer, 0902 and put the dot after last digit */ 0903 if (dotpos == -1) 0904 dotpos = last - buf; 0905 0906 /* leading zeros shift the dot to the left. 0907 dotpos is now the exponent that results 0908 from the position of the dot in the significand. */ 0909 dotpos -= zeros+1; 0910 0911 /* combine the dot position with the explicit exponent */ 0912 if (ovfl != 0 || !_checkadd(&dotpos, expsign * (int)exp)) 0913 /* exponent overflow */ 0914 float_setnan(f); 0915 0916 float_setexponent(f, dotpos); 0917 float_setsign(f, msign); 0918 } 0919 0920 /* normalizes a significand such that 1 <= x < 10. 0921 If the exponent overflows during this operation 0922 this is notified. */ 0923 static char 0924 _normalize( 0925 floatnum f) 0926 { 0927 _corr_lead_zero(f); 0928 if (f->significand != NULL && _lenof(f) > 1) 0929 _corr_overflow(f); 0930 if (f->significand != NULL) 0931 _corr_trailing_zeros(f); 0932 if (f->significand != NULL && !float_isvalidexp(f->exponent)) 0933 { 0934 float_seterror(Underflow); 0935 if (f->exponent > 0) 0936 float_seterror(Overflow); 0937 float_setnan(f); 0938 } 0939 #ifdef FLOATDEBUG 0940 if (f->significand != NULL) 0941 _setvalue_(f); 0942 #endif 0943 return f->significand != NULL; 0944 } 0945 0946 void 0947 float_setinteger(floatnum dest, int value) 0948 { 0949 char buf[BITS_IN_UNSIGNED/3 + 3]; 0950 0951 sprintf(buf, "%d", value); 0952 float_setscientific(dest, buf, NULLTERMINATED); 0953 } 0954 0955 void 0956 float_move( 0957 floatnum dest, 0958 floatnum source) 0959 { 0960 if (dest != source) 0961 { 0962 float_setnan(dest); 0963 *dest = *source; 0964 float_create(source); 0965 } 0966 } 0967 0968 /* creates a copy of <source> and assigns it to 0969 <dest>. The significand is guaranteed to have 0970 <scale>+1 digits. The <dest> significand is 0971 truncated, or padded with zeros to the right, 0972 to achieve the desired length. 0973 <scale> may assume the special value EXACT, in 0974 which case a true copy is generated. 0975 This function allows an in-place copy 0976 (dest == source). */ 0977 static void 0978 _scaled_clone( 0979 floatnum dest, 0980 cfloatnum source, 0981 int scale) 0982 { 0983 /* dest == source allowed! */ 0984 0985 bc_num mant; 0986 unsigned exp; 0987 signed char sign; 0988 0989 mant = NULL; 0990 if(scale == EXACT) 0991 scale = _scaleof(source); 0992 if (dest == source && scale <= _scaleof(source)) 0993 { 0994 _setscale(dest, scale); 0995 return; 0996 } 0997 mant = bc_new_num(1, scale); 0998 scale = _min(scale, _scaleof(source)); 0999 memcpy(mant->n_value, _valueof(source), scale+1); 1000 sign = float_getsign(source); 1001 exp = source->exponent; 1002 float_setnan(dest); 1003 dest->exponent = exp; 1004 dest->significand = mant; 1005 #ifdef FLOATDEBUG 1006 _setvalue_(dest); 1007 #endif 1008 float_setsign(dest, sign); 1009 } 1010 1011 char 1012 float_copy( 1013 floatnum dest, 1014 cfloatnum source, 1015 int digits) 1016 { 1017 int scale, save; 1018 1019 if (digits == EXACT) 1020 digits = _max(1, float_getlength(source)); 1021 if (!_checkdigits(digits, NOSPECIALVALUE)) 1022 return _seterror(dest, InvalidPrecision); 1023 if (_is_special(source)) 1024 { 1025 if (dest != source) 1026 float_free(dest); 1027 *dest = *source; 1028 } 1029 else 1030 { 1031 // invariant: source has to be restored, if it is != dest 1032 scale = _min(digits - 1, _scaleof(source)); 1033 save = _limit_scale((floatnum)source, scale); 1034 _corr_trailing_zeros((floatnum)source); 1035 _scaled_clone(dest, source, EXACT); 1036 if (dest != source) 1037 _setscale(source, save); 1038 } 1039 return TRUE; 1040 } 1041 1042 /* rounding a value towards zero */ 1043 static void 1044 _trunc( 1045 floatnum dest, 1046 cfloatnum x, 1047 int scale) 1048 { 1049 scale -= _bscandigit(x, scale, 0); 1050 _scaled_clone(dest, x, scale); 1051 #ifdef FLOATDEBUG 1052 _setvalue_(dest); 1053 #endif 1054 } 1055 1056 /* rounding a value towards infinity */ 1057 static char 1058 _roundup( 1059 floatnum dest, 1060 cfloatnum x, 1061 int scale) 1062 { 1063 scale -= _bscandigit(x, scale, 9); 1064 _scaled_clone(dest, x, _max(0, scale)); 1065 if (scale < 0) 1066 { 1067 *_valueof(dest) = 1; 1068 if (!float_isvalidexp(++dest->exponent)) 1069 return FALSE; 1070 } 1071 else 1072 { 1073 ++*(_valueof(dest) + scale); 1074 } 1075 #ifdef FLOATDEBUG 1076 _setvalue_(dest); 1077 #endif 1078 return TRUE; 1079 } 1080 1081 char 1082 float_round( 1083 floatnum dest, 1084 cfloatnum src, 1085 int digits, 1086 roundmode mode) 1087 { 1088 int scalediff, scale; 1089 char digit; 1090 signed char sign, updown; 1091 1092 if (mode > TOMINUSINFINITY) 1093 return _seterror(dest, InvalidParam); 1094 if (!_checkdigits(digits, NOSPECIALVALUE)) 1095 return _setnan(dest); 1096 if (float_isnan(src)) 1097 return _seterror(dest, NoOperand); 1098 updown = 0; 1099 scale = digits - 1; 1100 if (float_getlength(src) > digits) 1101 { 1102 sign = float_getsign(src); 1103 switch(mode) 1104 { 1105 case TONEAREST: 1106 scalediff = _scaleof(src) - scale; 1107 if (scalediff > 0) 1108 { 1109 digit = _digit(src, digits); 1110 if (digit < 5 1111 || (digit == 5 1112 && scalediff == 1 1113 && (_digit(src, scale) & 1) == 0)) 1114 updown = -1; 1115 else 1116 updown = 1; 1117 } 1118 break; 1119 case TOZERO: 1120 updown = -1; 1121 break; 1122 case TOINFINITY: 1123 updown = 1; 1124 break; 1125 case TOPLUSINFINITY: 1126 updown = sign; 1127 break; 1128 case TOMINUSINFINITY: 1129 updown = -sign; 1130 break; 1131 } 1132 } 1133 switch (updown) 1134 { 1135 case 1: 1136 if (!_roundup(dest, src, scale)) 1137 return _seterror(dest, Overflow); 1138 break; 1139 case 0: 1140 float_copy(dest, src, digits); 1141 break; 1142 case -1: 1143 _trunc(dest, src, scale); 1144 break; 1145 } 1146 return TRUE; 1147 } 1148 1149 char 1150 float_int( 1151 floatnum f) 1152 { 1153 if (!_checknan(f)) 1154 return FALSE; 1155 if (f->exponent < 0) 1156 float_setzero(f); 1157 else if (!float_iszero(f)) 1158 float_round(f, f, f->exponent+1, TOZERO); 1159 return TRUE; 1160 } 1161 1162 char 1163 float_frac( 1164 floatnum f) 1165 { 1166 if (!_checknan(f) || float_iszero(f) || f->exponent < 0) 1167 return !float_isnan(f); 1168 if (_scaleof(f) <= f->exponent) 1169 float_setzero(f); 1170 else 1171 { 1172 int leadingzeros = 0; 1173 _hidefirst(f, f->exponent + 1); 1174 leadingzeros = _scan_digit(_valueof(f), float_getlength(f), 0); 1175 _hidefirst(f, leadingzeros); 1176 f->exponent = -leadingzeros - 1; 1177 #ifdef FLOATDEBUG 1178 _setvalue_(f); 1179 #endif 1180 } 1181 return TRUE; 1182 } 1183 1184 /* the general purpose add/subtract routine that deals with 1185 the ordinary case. */ 1186 static char 1187 _addsub_normal( 1188 floatnum dest, 1189 floatnum summand1, 1190 floatnum summand2, 1191 int digits) 1192 { 1193 floatstruct tmp; 1194 int expdiff; 1195 int scale; 1196 int fulllength; 1197 int extradigit; 1198 1199 /* the operands are ordered by their exponent */ 1200 1201 expdiff = (unsigned)(summand1->exponent - summand2->exponent); 1202 1203 /* the full length of the sum (without carry) */ 1204 fulllength = _max(expdiff + _scaleof(summand2), _scaleof(summand1)) + 1; 1205 1206 extradigit = 0; 1207 if (digits == EXACT || digits > fulllength) 1208 digits = fulllength; 1209 else 1210 { 1211 if (float_getsign(summand1) + float_getsign(summand2) == 0) 1212 extradigit = 1; /* a true subtraction needs no space for a carry */ 1213 if (expdiff > digits + extradigit) 1214 /* second operand underflows due to exponent diff */ 1215 return float_copy(dest, summand1, digits+extradigit); 1216 } 1217 1218 if (digits > maxdigits) 1219 return _seterror(dest, InvalidPrecision); 1220 1221 /* we cannot add the operands directly 1222 because of possibly different exponents. 1223 So we assume the second operand "to be OK" 1224 and shift the decimal point of the first 1225 appropriately to the right. 1226 There is a cheap way to do this: 1227 increment len by expdiff and decrement 1228 scale by the same amount. 1229 But: Check the operand is long enough 1230 to do this. */ 1231 1232 float_create(&tmp); 1233 if (_scaleof(summand1) < expdiff) 1234 { 1235 _scaled_clone(&tmp, summand1, expdiff); 1236 summand1 = &tmp; 1237 } 1238 scale = digits + extradigit - (int)expdiff - 1; 1239 _movepoint(summand1, expdiff); 1240 1241 /* truncate overly long operands */ 1242 _limit_scale(summand1, scale); 1243 _limit_scale(summand2, scale); 1244 1245 /* add */ 1246 dest->exponent = summand2->exponent; 1247 bc_add(summand1->significand, 1248 summand2->significand, 1249 &(dest->significand), 1250 scale); 1251 1252 float_free(&tmp); 1253 1254 return _normalize(dest); 1255 } 1256 1257 static char 1258 _sub_checkborrow( 1259 floatnum dest, 1260 floatnum summand1, 1261 floatnum summand2, 1262 int digits) 1263 { 1264 /* the operands have opposite signs, the same exponent, 1265 but their first digit of the significand differ. 1266 The operands are ordered by this digit. */ 1267 int result; 1268 int borrow; 1269 int scale1, scale2; 1270 char save; 1271 char* v1; 1272 char* v2; 1273 1274 /* Cancellation occurs, when the operands are of type 1275 p.000...yyy - q.999...xxx, p-q == 1, because a borrow 1276 propagates from the difference ..0yyy.. - ..9xxx.., 1277 leaving zeros in the first part. We check for this here. */ 1278 1279 borrow = 0; 1280 if (_digit(summand1, 0) - _digit(summand2, 0) == 1) 1281 { 1282 scale1 = _scaleof(summand1); 1283 scale2 = _scaleof(summand2); 1284 if (scale1 == 0) 1285 /* the special case of a one-digit first operand 1286 p. - q.999..xxx */ 1287 borrow = _scan_digit(_valueof(summand2) + 1, scale2, 9); 1288 else if (scale2 > 0) 1289 /* count the 0 - 9 pairs after the first digit, the area 1290 where a borrow can propagate */ 1291 borrow = _scan_09pairs(summand1, summand2, _min(scale1, scale2)); 1292 1293 /* In case of a one-digit second operand (p.yyy.. - q. == 1.yyy..), 1294 nothing is cancelled out. Borrow is already set to 0, and this is 1295 the correct value for this case, so nothing has to be done here */ 1296 1297 if (borrow > 0) 1298 { 1299 /* we have cancellation here. We skip all digits, that 1300 cancel out due to a propagating borrow. These 1301 include the first digit and all following 1302 (0,9) digit pairs, except the last one. The last 1303 pair may be subject to cancellation or not, we do not 1304 care, because the following digit pair either yields a 1305 non-zero difference, or creates no borrow. Our standard 1306 adder is good enough to deal with such a limited 1307 cancelling effect. We will replace the last (0,9) 1308 digit pair with a (9,8) pair. This prevents the 1309 creation of a borrow, and yet, will yield the correct 1310 result */ 1311 1312 /* hide all digits until the last found 0 - 9 pair */ 1313 summand2->exponent -= borrow; 1314 summand1->exponent -= borrow; 1315 /* in case of a one-digit significand, there is nothing to hide */ 1316 if (scale1 > 0) 1317 _hidefirst(summand1, borrow); 1318 _hidefirst(summand2, borrow); 1319 1320 /* we replace the last found 0 - 9 pair by a 9 - 8 pair, 1321 avoiding a carry, yet yielding the correct result */ 1322 save = *(v1 = _valueof(summand1)); 1323 *v1 = 9; 1324 *(v2 = _valueof(summand2)) = 8; 1325 } 1326 } 1327 result = _addsub_normal(dest, summand1, summand2, digits); 1328 1329 /* restore the modified digits */ 1330 if (borrow > 0) 1331 { 1332 if (summand1 != dest) 1333 *v1 = save; 1334 if (summand2 != dest) 1335 *v2 = 9; 1336 } 1337 return result; 1338 } 1339 1340 static char 1341 _sub_expdiff0( 1342 floatnum dest, 1343 floatnum summand1, 1344 floatnum summand2, 1345 int digits) 1346 { 1347 int eq; 1348 int result; 1349 1350 /* the operands are ordered by their significand length, 1351 and have the same exponent, and different sign */ 1352 1353 /* One type of cancellation is when both significands set out 1354 with the same digits. Since these digits cancel out 1355 during subtraction, we look out for the first pair 1356 of different digits. eq receives the number of 1357 equal digits, which may be 0 */ 1358 eq = _scan_equal(summand1, summand2); 1359 if (float_getlength(summand2) == eq) 1360 { 1361 /* the complete second operand is cancelled out */ 1362 if (float_getlength(summand1) == eq) 1363 /* op1 == -op2 */ 1364 return _setzero(dest); 1365 /* If xxx.. denotes the second operand, the (longer) 1366 first one is of form xxx..yyy.., since it has 1367 the same digits in the beginning. During 1368 subtraction the xxx... part is cancelled out, and 1369 this leaves yyy... as the subtraction result. 1370 By copying the yyy... to the result, we can shortcut the 1371 subtraction. 1372 But before doing so, we have to check yyy... for 1373 leading zeros, because the cancellation continues in 1374 this case. */ 1375 eq += _scan_digit(_valueof(summand1) + eq, 1376 float_getlength(summand1) - eq, 0); 1377 _hidefirst(summand1, eq); 1378 result = float_copy(dest, summand1, digits); 1379 dest->exponent -= eq; 1380 return result != 0 && _normalize(dest); 1381 } 1382 /* hide the identical digits, and do the 1383 subtraction without them. */ 1384 summand1->exponent -= eq; 1385 summand2->exponent -= eq; 1386 _hidefirst(summand1, eq); 1387 _hidefirst(summand2, eq); 1388 1389 /* order the operands by their first digit */ 1390 if (_digit(summand1, 0) >= _digit(summand2, 0)) 1391 return _sub_checkborrow(dest, summand1, summand2, digits); 1392 return _sub_checkborrow(dest, summand2, summand1, digits); 1393 } 1394 1395 static char 1396 _sub_expdiff1( 1397 floatnum dest, 1398 floatnum summand1, 1399 floatnum summand2, 1400 int digits) 1401 { 1402 /* Cancellation occurs when subtracting 0.9xxx from 1403 1.0yyy */ 1404 1405 int result; 1406 char singledigit; 1407 char* v1; 1408 char* v2; 1409 1410 /* the operands have different sign, are ordered by their 1411 exponent, and the difference of the exponents is 1 */ 1412 1413 /* 1.0yyy may be given as a single digit 1 or as a string of 1414 digits starting with 1.0 */ 1415 singledigit = _scaleof(summand1) == 0; 1416 if (_digit(summand1, 0) != 1 1417 || _digit(summand2, 0) != 9 1418 || (!singledigit && _digit(summand1, 1) != 0)) 1419 return _addsub_normal(dest, summand1, summand2, digits); 1420 1421 /* we have cancellation here. We transform this 1422 case into that of equal exponents. */ 1423 1424 /* we align both operands by hiding the first digit (1) of the 1425 greater operand. This leaves .0yyy which matches the 1426 second operand .9xxx. Unfortunately, if the first operand 1427 has only one digit, we cannot hide it, so we have to 1428 work around this then. */ 1429 if (!singledigit) 1430 _hidefirst(summand1, 1); 1431 /* we change the leading digits into a '9' and a '8' resp. 1432 So we finally subtract .8xxx from .9yyy, yielding 1433 the correct result. */ 1434 v1 = _valueof(summand1); 1435 v2 = _valueof(summand2); 1436 *v1 = 9; 1437 *v2 = 8; 1438 summand1->exponent--; 1439 result = _sub_checkborrow(dest, summand1, summand2, digits); 1440 1441 /* restore the original digits */ 1442 if (summand1 != dest) 1443 *v1 = singledigit? 1 : 0; 1444 if (summand2 != dest) 1445 *v2 = 9; 1446 return result; 1447 } 1448 1449 static char 1450 _sub_ordered( 1451 floatnum dest, 1452 floatnum summand1, 1453 floatnum summand2, 1454 int digits) 1455 { 1456 /* we have to check for cancellation when subtracting. 1457 Cancellation occurs when the operands are almost 1458 equal in magnitude. E.g. in 1.234 - 1.226 = 0.008, 1459 the result is on quite a different scale than the 1460 operands. Actually, this is a big problem, because 1461 it means that the (true) subtraction is numerically not 1462 stable. There is no way to get around this; you always 1463 have to take this into account, when subtracting. 1464 We make the best out of it, and check for cancellation 1465 in advance, so the result is at least valid to all digits, 1466 if the operands are known to be exact. 1467 Cancellation occurs only, if the difference between the 1468 exponents is 1 at most. We prepare the critical cases in 1469 specialized routines, and let the standard routine do the 1470 rest. */ 1471 1472 /* the operands are ordered by their exponent */ 1473 1474 unsigned expdiff; 1475 1476 expdiff = (unsigned)(summand1->exponent - summand2->exponent); 1477 switch (expdiff) 1478 { 1479 case 0: 1480 /* order the operands by their length of the significands */ 1481 if (float_getlength(summand1) >= float_getlength(summand2)) 1482 return _sub_expdiff0(dest, summand1, summand2, digits); 1483 return _sub_expdiff0(dest, summand2, summand1, digits); 1484 case 1: 1485 return _sub_expdiff1(dest, summand1, summand2, digits); 1486 } 1487 return _addsub_normal(dest, summand1, summand2, digits); 1488 } 1489 1490 static char 1491 _addsub_ordered( 1492 floatnum dest, 1493 floatnum summand1, 1494 floatnum summand2, 1495 int digits) 1496 { 1497 /* operands are ordered by their exponent */ 1498 1499 /* handle a bunch of special cases */ 1500 if (!_checkdigits(digits, EXACT) || !_checknan(summand1)) 1501 return _setnan(dest); 1502 if (float_iszero(summand2)) 1503 return float_copy(dest, summand1, digits); 1504 1505 /* separate true addition from true subtraction */ 1506 if (float_getsign(summand1) == float_getsign(summand2)) 1507 return _addsub_normal(dest, summand1, summand2, digits); 1508 return _sub_ordered(dest, summand1, summand2, digits); 1509 } 1510 1511 char 1512 float_add( 1513 floatnum dest, 1514 cfloatnum summand1, 1515 cfloatnum summand2, 1516 int digits) 1517 { 1518 bc_struct bc1, bc2; 1519 floatstruct tmp1, tmp2; 1520 floatnum s1, s2; 1521 1522 /* the adder may occasionally adjust operands to 1523 his needs. Hence we work on temporary structures */ 1524 s1 = dest; 1525 s2 = dest; 1526 if (dest != summand1) 1527 { 1528 _copyfn(&tmp1, summand1, &bc1); 1529 s1 = &tmp1; 1530 } 1531 if (dest != summand2) 1532 { 1533 _copyfn(&tmp2, summand2, &bc2); 1534 s2 = &tmp2; 1535 } 1536 1537 /* order the operands by their exponent. This should 1538 bring a NaN always to the front, and keeps zeros after any 1539 other number. */ 1540 if (s1->exponent >= s2->exponent) 1541 return _addsub_ordered(dest, s1, s2, digits); 1542 return _addsub_ordered(dest, s2, s1, digits); 1543 } 1544 1545 char 1546 float_sub( 1547 floatnum dest, 1548 cfloatnum minuend, 1549 cfloatnum subtrahend, 1550 int scale) 1551 { 1552 int result; 1553 if (minuend == subtrahend) 1554 { 1555 /* changing the sign of one operand would change that of 1556 the other as well. So this is a special case */ 1557 if(!_checknan(minuend)) 1558 return FALSE; 1559 _setzero(dest); 1560 } 1561 /* do not use float_neg, because it may change float_error */ 1562 float_setsign((floatnum)subtrahend, -float_getsign(subtrahend)); 1563 result = float_add(dest, minuend, subtrahend, scale); 1564 if (dest != subtrahend) 1565 float_setsign((floatnum)subtrahend, -float_getsign(subtrahend)); 1566 return result; 1567 } 1568 1569 char 1570 float_mul( 1571 floatnum dest, 1572 cfloatnum factor1, 1573 cfloatnum factor2, 1574 int digits) 1575 { 1576 int result; 1577 int fullscale; 1578 int savescale1, savescale2; 1579 int scale; 1580 1581 /* handle a bunch of special cases */ 1582 if (!_checkdigits(digits, EXACT) 1583 || !_checknan(factor1) 1584 || !_checknan(factor2)) 1585 /* invalid scale value or NaN operand */ 1586 return _setnan(dest); 1587 if (float_iszero(factor1) || float_iszero(factor2)) 1588 return _setzero(dest); 1589 1590 scale = digits - 1; 1591 fullscale = _scaleof(factor1) + _scaleof(factor2); 1592 if (digits == EXACT || scale > fullscale) 1593 scale = fullscale; 1594 1595 if (scale >= maxdigits) 1596 /* scale too large */ 1597 return _seterror(dest, InvalidPrecision); 1598 1599 /* limit the scale of the operands to sane sizes */ 1600 savescale1 = _limit_scale((floatnum)factor1, scale); 1601 savescale2 = _limit_scale((floatnum)factor2, scale); 1602 1603 /* multiply */ 1604 dest->exponent = factor1->exponent + factor2->exponent; 1605 bc_multiply(factor1->significand, factor2->significand, &(dest->significand), scale); 1606 result = _normalize(dest); 1607 1608 /* reverse order is necessary in case factor1 == factor2 */ 1609 if (dest != factor2) 1610 _setscale(factor2, savescale2); 1611 if (dest != factor1) 1612 _setscale(factor1, savescale1); 1613 return result; 1614 } 1615 1616 char 1617 float_div( 1618 floatnum dest, 1619 cfloatnum dividend, 1620 cfloatnum divisor, 1621 int digits) 1622 { 1623 int result; 1624 int savescale1, savescale2; 1625 int exp; 1626 1627 /* handle a bunch of special cases */ 1628 if (!_checkdigits(digits, INTQUOT) || !_checknan(dividend) 1629 || !_checknan(divisor)) 1630 return _setnan(dest); 1631 if (float_iszero(divisor)) 1632 return _seterror(dest, ZeroDivide); 1633 if (float_iszero(dividend)) 1634 /* 0/x == 0 */ 1635 return _setzero(dest); 1636 1637 exp = dividend->exponent - divisor->exponent; 1638 1639 /* check for integer quotient */ 1640 if(digits == INTQUOT) 1641 { 1642 if (exp < 0) 1643 return _setzero(dest); 1644 digits = exp; 1645 } 1646 1647 /* scale OK? */ 1648 if(digits > maxdigits) 1649 return _seterror(dest, InvalidPrecision); 1650 1651 /* limit the scale of the operands to sane sizes */ 1652 savescale1 = _limit_scale((floatnum)dividend, digits); 1653 savescale2 = _limit_scale((floatnum)divisor, digits); 1654 1655 /* divide */ 1656 result = TRUE; 1657 dest->exponent = exp; 1658 bc_divide(dividend->significand, 1659 divisor->significand, 1660 &(dest->significand), 1661 digits); 1662 if (bc_is_zero(dest->significand)) 1663 float_setzero(dest); 1664 else 1665 result = _normalize(dest); 1666 1667 /* reverse order is necessary in case divisor == dividend */ 1668 if (dest != divisor) 1669 _setscale(divisor, savescale2); 1670 if (dest != dividend) 1671 _setscale(dividend, savescale1); 1672 return result; 1673 } 1674 1675 char 1676 float_divmod( 1677 floatnum quotient, 1678 floatnum remainder, 1679 cfloatnum dividend, 1680 cfloatnum divisor, 1681 int digits) 1682 { 1683 int exp, exp1; 1684 1685 if (!_checkdigits(digits, INTQUOT) || !_checknan(dividend) 1686 || !_checknan(divisor) || quotient == remainder 1687 || float_iszero(divisor) || float_getlength(divisor) > maxdigits) 1688 { 1689 if (quotient == remainder) 1690 float_seterror(InvalidParam); 1691 if (float_getlength(divisor) > maxdigits) 1692 float_seterror(TooExpensive); 1693 if (float_iszero(divisor)) 1694 float_seterror(ZeroDivide); 1695 float_setnan(quotient); 1696 return _setnan(remainder); 1697 } 1698 if (float_iszero(dividend)) 1699 { 1700 float_setzero(quotient); 1701 return _setzero(remainder); 1702 } 1703 exp1 = dividend->exponent; 1704 exp = exp1 - divisor->exponent; 1705 if(digits-- == INTQUOT) 1706 { 1707 if (exp < 0) 1708 { 1709 if (float_copy(remainder, dividend, EXACT)) 1710 return _setzero(quotient); 1711 return _setnan(quotient); 1712 } 1713 digits = exp; 1714 } 1715 if (digits > maxdigits) 1716 { 1717 float_setnan(quotient); 1718 return _seterror(remainder, TooExpensive); 1719 } 1720 1721 /* divide */ 1722 quotient->exponent = exp; 1723 remainder->exponent = exp1; 1724 bc_divmod(dividend->significand, 1725 divisor->significand, 1726 &(quotient->significand), 1727 &(remainder->significand), 1728 digits); 1729 1730 /* if something goes wrong (one of the results overflows 1731 or underflows), always set both quotient and remainder 1732 to NaN */ 1733 if (bc_is_zero(remainder->significand)) 1734 float_setzero(remainder); 1735 else if (!_normalize(remainder)) 1736 return _setnan(quotient); 1737 if (bc_is_zero(quotient->significand)) 1738 float_setzero(quotient); 1739 else if (!_normalize(quotient)) 1740 return _setnan(remainder); 1741 return TRUE; 1742 } 1743 1744 char 1745 float_sqrt(floatnum value, int digits) 1746 { 1747 if (!_checkdigits(digits, NOSPECIALVALUE) || !_checknan(value)) 1748 return _setnan(value); 1749 switch (float_getsign(value)) 1750 { 1751 case -1: 1752 return _seterror(value, OutOfDomain); 1753 case 0: 1754 return TRUE; 1755 } 1756 if ((value->exponent & 1) != 0) 1757 { 1758 if (float_getlength(value) == 1) 1759 _scaled_clone(value, value, 1); 1760 _movepoint(value, 1); 1761 } 1762 bc_sqrt(&value->significand, digits - 1); 1763 #ifdef FLOATDEBUG 1764 _setvalue_(value); 1765 #endif 1766 if (value->exponent >= 0) 1767 value->exponent >>= 1; 1768 else 1769 value->exponent = -((1-value->exponent) >> 1); 1770 return TRUE; 1771 }