File indexing completed on 2024-05-12 17:22:35
0001 /* floatlong.c: portable double size integer arithmetic. */ 0002 /* 0003 Copyright (C) 2007 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 "floatlong.h" 0033 0034 #define HALFSIZE (sizeof(unsigned) * 4) 0035 #define LOWMASK ((1 << HALFSIZE) - 1) 0036 #define MSB (1 << (2*HALFSIZE-1)) 0037 0038 /*************** functions handling a single unsigned ********************/ 0039 0040 int 0041 _findfirstbit( 0042 unsigned value) 0043 { 0044 int result; 0045 0046 result = -1; 0047 while (value != 0) 0048 { 0049 value >>= 1; 0050 ++result; 0051 } 0052 return result; 0053 } 0054 0055 /* pre: value != 0 */ 0056 static int 0057 _revfindfirstbit( 0058 unsigned value) 0059 { 0060 int result; 0061 0062 result = 0; 0063 while ((value & 1) == 0) 0064 { 0065 value >>= 1; 0066 ++result; 0067 } 0068 return result; 0069 } 0070 0071 /******************* double unsigned functions ********************/ 0072 0073 static void 0074 _longsplit( 0075 unsigned value, 0076 unsigned* low, 0077 unsigned* high) 0078 { 0079 *high = value >> HALFSIZE; 0080 *low = value & LOWMASK; 0081 } 0082 0083 static unsigned 0084 _longcat( 0085 unsigned low, 0086 unsigned high) 0087 { 0088 return (high << HALFSIZE) + low; 0089 } 0090 0091 char 0092 _longadd( 0093 unsigned* s1, 0094 unsigned* s2) 0095 { 0096 unsigned s1h, s1l, s2h, s2l; 0097 0098 _longsplit(*s1, &s1l, &s1h); 0099 _longsplit(*s2, &s2l, &s2h); 0100 s1l += s2l; 0101 _longsplit(s1l, &s1l, &s2l); 0102 s1h += s2h + s2l; 0103 _longsplit(s1h, &s1h, s2); 0104 *s1 = _longcat(s1l, s1h); 0105 return *s2 == 0; 0106 } 0107 0108 char 0109 _longmul( 0110 unsigned* f1, 0111 unsigned* f2) 0112 { 0113 unsigned f1h, f1l, f2h, f2l; 0114 0115 _longsplit(*f1, &f1l, &f1h); 0116 _longsplit(*f2, &f2l, &f2h); 0117 *f1 = f1l * f2l; 0118 *f2 = f1h * f2h; 0119 f1l *= f2h; 0120 f2l *= f1h; 0121 _longadd(&f1l, &f2l); 0122 _longsplit(f1l, &f1l, &f1h); 0123 f1l <<= HALFSIZE; 0124 _longadd(f1, &f1l); 0125 *f2 += f1l + (f2l << HALFSIZE) + f1h; 0126 return *f2 == 0; 0127 } 0128 0129 unsigned 0130 _longshr( 0131 unsigned low, 0132 unsigned high, 0133 char shift) 0134 { 0135 if (shift == 0) 0136 return low; 0137 return (low >> shift) | (high << ((2*HALFSIZE)-shift)); 0138 } 0139 0140 unsigned 0141 _longshl( 0142 unsigned low, 0143 unsigned high, 0144 char shift) 0145 { 0146 if (shift == 0) 0147 return high; 0148 return (low >> ((2*HALFSIZE)-shift)) | (high << shift); 0149 } 0150 0151 /***************** unsigned array functions *****************/ 0152 0153 unsigned 0154 _longarrayadd( 0155 unsigned* uarray, 0156 int lg, 0157 unsigned incr) 0158 { 0159 for(; lg-- > 0 && incr != 0;) 0160 _longadd(uarray++, &incr); 0161 return incr; 0162 } 0163 0164 unsigned 0165 _longarraymul( 0166 unsigned* uarray, 0167 int lg, 0168 unsigned factor) 0169 { 0170 unsigned ovfl, carry; 0171 0172 carry = 0; 0173 ovfl = 0; 0174 for (; lg-- > 0;) 0175 { 0176 carry += ovfl; 0177 ovfl = factor; 0178 _longmul(uarray, &ovfl); 0179 _longadd(uarray++, &carry); 0180 } 0181 return carry + ovfl; 0182 } 0183 0184 unsigned 0185 _bitsubstr( 0186 unsigned* uarray, 0187 int ofs) 0188 { 0189 int idx; 0190 0191 if (ofs <= 0) 0192 return *uarray << -ofs; 0193 idx = ofs / BITS_IN_UNSIGNED; 0194 return _longshr(*(uarray+idx), *(uarray+idx+1), 0195 ofs - idx * BITS_IN_UNSIGNED); 0196 } 0197 0198 void 0199 _orsubstr( 0200 unsigned* uarray, 0201 int bitofs, 0202 unsigned value) 0203 { 0204 int idx; 0205 int ofs; 0206 0207 idx = bitofs / BITS_IN_UNSIGNED; 0208 ofs = bitofs - idx * BITS_IN_UNSIGNED; 0209 if (ofs == 0) 0210 *(uarray + idx) |= value; 0211 else 0212 { 0213 *(uarray + idx) |= value << ofs; 0214 *(uarray + idx + 1) |= value >> (BITS_IN_UNSIGNED - ofs); 0215 } 0216 } 0217 0218 /**************** longint functions ****************/ 0219 0220 char 0221 _isfull( 0222 t_longint* l) 0223 { 0224 return l->length >= (int)UARRAYLG - 1; 0225 } 0226 0227 unsigned 0228 _bitlength( 0229 t_longint* l) 0230 { 0231 if (l->length == 0) 0232 return 0; 0233 return (l->length - 1) * BITS_IN_UNSIGNED 0234 + _findfirstbit(l->value[l->length-1]) + 1; 0235 } 0236 0237 unsigned 0238 _lastnonzerobit( 0239 t_longint* l) 0240 { 0241 int i; 0242 0243 if (l->length == 0) 0244 return -1; 0245 i = -1; 0246 for (; ++i < l->length && l->value[i] == 0;); 0247 return i * BITS_IN_UNSIGNED + _revfindfirstbit(l->value[i]); 0248 } 0249 0250 unsigned 0251 _longintadd( 0252 t_longint* l, 0253 unsigned summand) 0254 { 0255 unsigned ovfl; 0256 0257 ovfl = _longarrayadd(l->value, l->length, summand); 0258 if (ovfl != 0 && !_isfull(l)) 0259 { 0260 l->value[l->length] = ovfl; 0261 ++l->length; 0262 ovfl = 0; 0263 } 0264 return ovfl; 0265 } 0266 0267 unsigned 0268 _longintmul( 0269 t_longint* l, 0270 unsigned factor) 0271 { 0272 unsigned ovfl; 0273 0274 ovfl = _longarraymul(l->value, l->length, factor); 0275 if (ovfl != 0 && !_isfull(l)) 0276 { 0277 l->value[l->length] = ovfl; 0278 ++l->length; 0279 ovfl = 0; 0280 } 0281 return ovfl; 0282 } 0283 0284 char 0285 _longintsetsize( 0286 t_longint* l, 0287 unsigned bitlength) 0288 { 0289 int i; 0290 0291 if (bitlength == 0) 0292 l->length = 0; 0293 else 0294 l->length = (bitlength - 1) / BITS_IN_UNSIGNED + 1; 0295 if (l->length >= (int)UARRAYLG) 0296 return 0; 0297 for(i = l->length; i >= 0; --i) 0298 l->value[i] = 0; 0299 return 1; 0300 } 0301 0302 /**************** overflow checking integer functions *************/ 0303 0304 char 0305 _checkadd( 0306 int* s1, 0307 int s2) 0308 { 0309 unsigned u1, u2; 0310 0311 if ((*s1 < 0) ^ (s2 < 0)) 0312 { 0313 /* you can always add summands of opposite sign */ 0314 *s1 += s2; 0315 return 1; 0316 } 0317 if (s2 < 0) 0318 { 0319 u1 = -*s1; 0320 u2 = -s2; 0321 _longadd(&u1, &u2); 0322 *s1 = -(int)u1; 0323 return u2 == 0 && *s1 < 0; 0324 } 0325 u1 = *s1; 0326 u2 = s2; 0327 _longadd(&u1, &u2); 0328 *s1 = u1; 0329 return *s1 >= 0; 0330 } 0331 0332 char 0333 _checkmul( 0334 int* f1, 0335 int f2) 0336 { 0337 unsigned x1, x2; 0338 signed char sgn = 1; 0339 const unsigned msb = MSB; 0340 0341 if (*f1 >= 0) 0342 x1 = *f1; 0343 else 0344 { 0345 sgn = -1; 0346 x1 = -*f1; 0347 } 0348 if (f2 >= 0) 0349 x2 = f2; 0350 else 0351 { 0352 sgn = -sgn; 0353 x2 = -f2; 0354 } 0355 _longmul(&x1, &x2); 0356 if (sgn < 0) 0357 { 0358 *f1 = -(int)x1; 0359 return (x2 == 0 && x1 <= msb); 0360 } 0361 *f1 = (int)x1; 0362 return (x2 == 0 && x1 < msb); 0363 }