File indexing completed on 2024-05-12 17:22:35
0001 /* floatpower.c: power operation, based on floatnum. */ 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 #include "floatipower.h" 0033 #include "floatconst.h" 0034 #include "floatcommon.h" 0035 #include "floatlong.h" 0036 0037 /* for radix conversion, we need to get a result, 0038 even though it might slightly overflow or underflow. 0039 That is why we keep the exponent separate. 0040 For limited exponents (< 1024) the relative 0041 error for a 100 digit calculation is < 1e-99 0042 x != 0 */ 0043 char 0044 _raiseposi( 0045 floatnum x, 0046 int* expx, 0047 unsigned exponent, 0048 int digits) 0049 { 0050 int exppwr, extra; 0051 floatstruct pwr; 0052 0053 float_create(&pwr); 0054 float_move(&pwr, x); 0055 exppwr = float_getexponent(&pwr); 0056 float_setexponent(&pwr, 0); 0057 if ((exponent & 1) != 0) 0058 { 0059 float_copy(x, &pwr, digits+1); 0060 *expx = exppwr; 0061 } 0062 else 0063 { 0064 float_copy(x, &c1, EXACT); 0065 *expx = 0; 0066 } 0067 extra = _findfirstbit(exponent)/3+1; 0068 while((exponent >>= 1) != 0) 0069 { 0070 float_mul(&pwr, &pwr, &pwr, digits+extra); 0071 if (!_checkadd(&exppwr, exppwr)) 0072 break; 0073 exppwr += float_getexponent(&pwr); 0074 float_setexponent(&pwr, 0); 0075 if((exponent & 1) != 0) 0076 { 0077 float_mul(x, x, &pwr, digits+extra); 0078 *expx += exppwr + float_getexponent(x); 0079 float_setexponent(x, 0); 0080 } 0081 } 0082 float_free(&pwr); 0083 return exponent == 0; 0084 } 0085 0086 char 0087 _raiseposi_( 0088 floatnum x, 0089 unsigned exponent, 0090 int digits) 0091 { 0092 int exp; 0093 char result; 0094 0095 result = _raiseposi(x, &exp, exponent, digits); 0096 float_setexponent(x, exp); 0097 return result; 0098 } 0099 0100 /* raises a non-zero x to the exponent-th power */ 0101 char 0102 _raisei( 0103 floatnum x, 0104 int exponent, 0105 int digits) 0106 { 0107 int expx; 0108 signed char sgn; 0109 char negexp; 0110 0111 switch (exponent) 0112 { 0113 case 0: 0114 float_copy(x, &c1, EXACT); /* fall through */ 0115 case 1: 0116 return 1; 0117 } 0118 if (float_getlength(x) == 1 && float_getdigit(x, 0) == 1) 0119 { 0120 /* power of ten */ 0121 if (!_checkmul(&exponent, float_getexponent(x)) 0122 || exponent < EXPMIN || exponent > EXPMAX) 0123 return 0; 0124 sgn = float_getsign(x) < 0? -(exponent & 1) : 1; 0125 float_setexponent(x, exponent); 0126 float_setsign(x, sgn); 0127 return 1; 0128 } 0129 negexp = exponent < 0; 0130 if (negexp) 0131 exponent = -exponent; 0132 if (!_raiseposi(x, &expx, exponent, digits) 0133 || expx < EXPMIN || expx > EXPMAX) 0134 return 0; 0135 float_setexponent(x, expx); 0136 return negexp? float_reciprocal(x, digits) : 1; 0137 }