Rev 804 | Rev 829 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 804 | Rev 828 | ||
---|---|---|---|
Line 26... | Line 26... | ||
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
27 | */ |
27 | */ |
28 | 28 | ||
29 | #include<sftypes.h> |
29 | #include<sftypes.h> |
30 | #include<add.h> |
30 | #include<add.h> |
- | 31 | #include<div.h> |
|
31 | #include<comparison.h> |
32 | #include<comparison.h> |
- | 33 | #include<mul.h> |
|
32 | 34 | ||
33 | float32 divFloat32(float32 a, float32 b) |
35 | float32 divFloat32(float32 a, float32 b) |
34 | { |
36 | { |
35 | float32 result; |
37 | float32 result; |
36 | __s32 aexp, bexp, cexp; |
38 | __s32 aexp, bexp, cexp; |
Line 138... | Line 140... | ||
138 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
140 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
139 | } |
141 | } |
140 | 142 | ||
141 | /* pack and round */ |
143 | /* pack and round */ |
142 | 144 | ||
143 | /* TODO: find first nonzero digit and shift result and detect possibly underflow */ |
145 | /* find first nonzero digit and shift result and detect possibly underflow */ |
144 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
146 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
145 | cexp--; |
147 | cexp--; |
146 | cfrac <<= 1; |
148 | cfrac <<= 1; |
147 | /* TODO: fix underflow */ |
149 | /* TODO: fix underflow */ |
148 | }; |
150 | }; |
149 | 151 | ||
150 | - | ||
151 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
152 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
152 | 153 | ||
153 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
154 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
154 | ++cexp; |
155 | ++cexp; |
155 | cfrac >>= 1; |
156 | cfrac >>= 1; |
Line 183... | Line 184... | ||
183 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
184 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
184 | 185 | ||
185 | return result; |
186 | return result; |
186 | } |
187 | } |
187 | 188 | ||
- | 189 | float64 divFloat64(float64 a, float64 b) |
|
- | 190 | { |
|
- | 191 | float64 result; |
|
- | 192 | __s32 aexp, bexp, cexp; |
|
- | 193 | __u64 afrac, bfrac, cfrac; |
|
- | 194 | __u64 remlo, remhi; |
|
- | 195 | ||
- | 196 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
|
- | 197 | ||
- | 198 | if (isFloat64NaN(a)) { |
|
- | 199 | if (isFloat64SigNaN(a)) { |
|
- | 200 | /*FIXME: SigNaN*/ |
|
- | 201 | } |
|
- | 202 | /*NaN*/ |
|
- | 203 | return a; |
|
- | 204 | } |
|
- | 205 | ||
- | 206 | if (isFloat64NaN(b)) { |
|
- | 207 | if (isFloat64SigNaN(b)) { |
|
- | 208 | /*FIXME: SigNaN*/ |
|
- | 209 | } |
|
- | 210 | /*NaN*/ |
|
- | 211 | return b; |
|
- | 212 | } |
|
- | 213 | ||
- | 214 | if (isFloat64Infinity(a)) { |
|
- | 215 | if (isFloat64Infinity(b)) { |
|
- | 216 | /*FIXME: inf / inf */ |
|
- | 217 | result.binary = FLOAT64_NAN; |
|
- | 218 | return result; |
|
- | 219 | } |
|
- | 220 | /* inf / num */ |
|
- | 221 | result.parts.exp = a.parts.exp; |
|
- | 222 | result.parts.fraction = a.parts.fraction; |
|
- | 223 | return result; |
|
- | 224 | } |
|
- | 225 | ||
- | 226 | if (isFloat64Infinity(b)) { |
|
- | 227 | if (isFloat64Zero(a)) { |
|
- | 228 | /* FIXME 0 / inf */ |
|
- | 229 | result.parts.exp = 0; |
|
- | 230 | result.parts.fraction = 0; |
|
- | 231 | return result; |
|
- | 232 | } |
|
- | 233 | /* FIXME: num / inf*/ |
|
- | 234 | result.parts.exp = 0; |
|
- | 235 | result.parts.fraction = 0; |
|
- | 236 | return result; |
|
- | 237 | } |
|
- | 238 | ||
- | 239 | if (isFloat64Zero(b)) { |
|
- | 240 | if (isFloat64Zero(a)) { |
|
- | 241 | /*FIXME: 0 / 0*/ |
|
- | 242 | result.binary = FLOAT64_NAN; |
|
- | 243 | return result; |
|
- | 244 | } |
|
- | 245 | /* FIXME: division by zero */ |
|
- | 246 | result.parts.exp = 0; |
|
- | 247 | result.parts.fraction = 0; |
|
- | 248 | return result; |
|
- | 249 | } |
|
- | 250 | ||
- | 251 | ||
- | 252 | afrac = a.parts.fraction; |
|
- | 253 | aexp = a.parts.exp; |
|
- | 254 | bfrac = b.parts.fraction; |
|
- | 255 | bexp = b.parts.exp; |
|
- | 256 | ||
- | 257 | /* denormalized numbers */ |
|
- | 258 | if (aexp == 0) { |
|
- | 259 | if (afrac == 0) { |
|
- | 260 | result.parts.exp = 0; |
|
- | 261 | result.parts.fraction = 0; |
|
- | 262 | return result; |
|
- | 263 | } |
|
- | 264 | /* normalize it*/ |
|
- | 265 | ||
- | 266 | afrac <<= 1; |
|
- | 267 | /* afrac is nonzero => it must stop */ |
|
- | 268 | while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
|
- | 269 | afrac <<= 1; |
|
- | 270 | aexp--; |
|
- | 271 | } |
|
- | 272 | } |
|
- | 273 | ||
- | 274 | if (bexp == 0) { |
|
- | 275 | bfrac <<= 1; |
|
- | 276 | /* bfrac is nonzero => it must stop */ |
|
- | 277 | while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
|
- | 278 | bfrac <<= 1; |
|
- | 279 | bexp--; |
|
- | 280 | } |
|
- | 281 | } |
|
- | 282 | ||
- | 283 | afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); |
|
- | 284 | bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); |
|
- | 285 | ||
- | 286 | if ( bfrac <= (afrac << 1) ) { |
|
- | 287 | afrac >>= 1; |
|
- | 288 | aexp++; |
|
- | 289 | } |
|
- | 290 | ||
- | 291 | cexp = aexp - bexp + FLOAT64_BIAS - 2; |
|
- | 292 | ||
- | 293 | cfrac = divFloat64estim(afrac, bfrac); |
|
- | 294 | ||
- | 295 | if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ |
|
- | 296 | mul64integers( bfrac, cfrac, &remlo, &remhi); |
|
- | 297 | /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ |
|
- | 298 | remhi = afrac - remhi - ( remlo > 0); |
|
- | 299 | remlo = - remlo; |
|
- | 300 | ||
- | 301 | while ((__s64) remhi < 0) { |
|
- | 302 | cfrac--; |
|
- | 303 | remlo += bfrac; |
|
- | 304 | remhi += ( remlo < bfrac ); |
|
- | 305 | } |
|
- | 306 | cfrac |= ( remlo != 0 ); |
|
- | 307 | } |
|
- | 308 | ||
- | 309 | /* pack and round */ |
|
- | 310 | ||
- | 311 | /* find first nonzero digit and shift result and detect possibly underflow */ |
|
- | 312 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ) )))) { |
|
- | 313 | cexp--; |
|
- | 314 | cfrac <<= 1; |
|
- | 315 | /* TODO: fix underflow */ |
|
- | 316 | }; |
|
- | 317 | ||
- | 318 | ||
- | 319 | cfrac >>= 1; |
|
- | 320 | ++cexp; |
|
- | 321 | cfrac += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3)); |
|
- | 322 | ||
- | 323 | if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) { |
|
- | 324 | ++cexp; |
|
- | 325 | cfrac >>= 1; |
|
- | 326 | } |
|
- | 327 | ||
- | 328 | /* check overflow */ |
|
- | 329 | if (cexp >= FLOAT64_MAX_EXPONENT ) { |
|
- | 330 | /* FIXME: overflow, return infinity */ |
|
- | 331 | result.parts.exp = FLOAT64_MAX_EXPONENT; |
|
- | 332 | result.parts.fraction = 0; |
|
- | 333 | return result; |
|
- | 334 | } |
|
- | 335 | ||
- | 336 | if (cexp < 0) { |
|
- | 337 | /* FIXME: underflow */ |
|
- | 338 | result.parts.exp = 0; |
|
- | 339 | if ((cexp + FLOAT64_FRACTION_SIZE) < 0) { |
|
- | 340 | result.parts.fraction = 0; |
|
- | 341 | return result; |
|
- | 342 | } |
|
- | 343 | cfrac >>= 1; |
|
- | 344 | while (cexp < 0) { |
|
- | 345 | cexp ++; |
|
- | 346 | cfrac >>= 1; |
|
- | 347 | } |
|
- | 348 | return result; |
|
- | 349 | ||
- | 350 | } else { |
|
- | 351 | cexp ++; /*normalized*/ |
|
- | 352 | result.parts.exp = (__u32)cexp; |
|
- | 353 | } |
|
- | 354 | ||
- | 355 | result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK)); |
|
- | 356 | ||
- | 357 | return result; |
|
- | 358 | } |
|
- | 359 | ||
- | 360 | __u64 divFloat64estim(__u64 a, __u64 b) |
|
- | 361 | { |
|
- | 362 | __u64 bhi; |
|
- | 363 | __u64 remhi, remlo; |
|
- | 364 | __u64 result; |
|
- | 365 | ||
- | 366 | if ( b <= a ) { |
|
- | 367 | return 0xFFFFFFFFFFFFFFFFull; |
|
- | 368 | } |
|
- | 369 | ||
- | 370 | bhi = b >> 32; |
|
- | 371 | result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; |
|
- | 372 | mul64integers(b, result, &remlo, &remhi); |
|
- | 373 | ||
- | 374 | remhi = a - remhi - (remlo > 0); |
|
- | 375 | remlo = - remlo; |
|
- | 376 | ||
- | 377 | b <<= 32; |
|
- | 378 | while ( (__s64) remhi < 0 ) { |
|
- | 379 | result -= 0x1ll << 32; |
|
- | 380 | remlo += b; |
|
- | 381 | remhi += bhi + ( remlo < b ); |
|
- | 382 | } |
|
- | 383 | remhi = (remhi << 32) | (remlo >> 32); |
|
- | 384 | if (( bhi << 32) <= remhi) { |
|
- | 385 | result |= 0xFFFFFFFF; |
|
- | 386 | } else { |
|
- | 387 | result |= remhi / bhi; |
|
- | 388 | } |
|
- | 389 | ||
- | 390 | ||
- | 391 | return result; |
|
- | 392 | } |
|
- | 393 |