1 /* multi_arith.h: multi-precision integer arithmetic functions, needed 2 to do extended-precision floating point. 3 4 (c) 1998 David Huggins-Daines. 5 6 Somewhat based on arch/alpha/math-emu/ieee-math.c, which is (c) 7 David Mosberger-Tang. 8 9 You may copy, modify, and redistribute this file under the terms of 10 the GNU General Public License, version 2, or any later version, at 11 your convenience. */ 12 13 /* Note: 14 15 These are not general multi-precision math routines. Rather, they 16 implement the subset of integer arithmetic that we need in order to 17 multiply, divide, and normalize 128-bit unsigned mantissae. */ 18 19 #ifndef MULTI_ARITH_H 20 #define MULTI_ARITH_H 21 22 #if 0 /* old code... */ 23 24 /* Unsigned only, because we don't need signs to multiply and divide. */ 25 typedef unsigned int int128[4]; 26 27 /* Word order */ 28 enum { 29 MSW128, 30 NMSW128, 31 NLSW128, 32 LSW128 33 }; 34 35 /* big-endian */ 36 #define LO_WORD(ll) (((unsigned int *) &ll)[1]) 37 #define HI_WORD(ll) (((unsigned int *) &ll)[0]) 38 39 /* Convenience functions to stuff various integer values into int128s */ 40 41 static inline void zero128(int128 a) 42 { 43 a[LSW128] = a[NLSW128] = a[NMSW128] = a[MSW128] = 0; 44 } 45 46 /* Human-readable word order in the arguments */ 47 static inline void set128(unsigned int i3, unsigned int i2, unsigned int i1, 48 unsigned int i0, int128 a) 49 { 50 a[LSW128] = i0; 51 a[NLSW128] = i1; 52 a[NMSW128] = i2; 53 a[MSW128] = i3; 54 } 55 56 /* Convenience functions (for testing as well) */ 57 static inline void int64_to_128(unsigned long long src, int128 dest) 58 { 59 dest[LSW128] = (unsigned int) src; 60 dest[NLSW128] = src >> 32; 61 dest[NMSW128] = dest[MSW128] = 0; 62 } 63 64 static inline void int128_to_64(const int128 src, unsigned long long *dest) 65 { 66 *dest = src[LSW128] | (long long) src[NLSW128] << 32; 67 } 68 69 static inline void put_i128(const int128 a) 70 { 71 printk("%08x %08x %08x %08x\n", a[MSW128], a[NMSW128], 72 a[NLSW128], a[LSW128]); 73 } 74 75 /* Internal shifters: 76 77 Note that these are only good for 0 < count < 32. 78 */ 79 80 static inline void _lsl128(unsigned int count, int128 a) 81 { 82 a[MSW128] = (a[MSW128] << count) | (a[NMSW128] >> (32 - count)); 83 a[NMSW128] = (a[NMSW128] << count) | (a[NLSW128] >> (32 - count)); 84 a[NLSW128] = (a[NLSW128] << count) | (a[LSW128] >> (32 - count)); 85 a[LSW128] <<= count; 86 } 87 88 static inline void _lsr128(unsigned int count, int128 a) 89 { 90 a[LSW128] = (a[LSW128] >> count) | (a[NLSW128] << (32 - count)); 91 a[NLSW128] = (a[NLSW128] >> count) | (a[NMSW128] << (32 - count)); 92 a[NMSW128] = (a[NMSW128] >> count) | (a[MSW128] << (32 - count)); 93 a[MSW128] >>= count; 94 } 95 96 /* Should be faster, one would hope */ 97 98 static inline void lslone128(int128 a) 99 { 100 asm volatile ("lsl.l #1,%0\n" 101 "roxl.l #1,%1\n" 102 "roxl.l #1,%2\n" 103 "roxl.l #1,%3\n" 104 : 105 "=d" (a[LSW128]), 106 "=d"(a[NLSW128]), 107 "=d"(a[NMSW128]), 108 "=d"(a[MSW128]) 109 : 110 "0"(a[LSW128]), 111 "1"(a[NLSW128]), 112 "2"(a[NMSW128]), 113 "3"(a[MSW128])); 114 } 115 116 static inline void lsrone128(int128 a) 117 { 118 asm volatile ("lsr.l #1,%0\n" 119 "roxr.l #1,%1\n" 120 "roxr.l #1,%2\n" 121 "roxr.l #1,%3\n" 122 : 123 "=d" (a[MSW128]), 124 "=d"(a[NMSW128]), 125 "=d"(a[NLSW128]), 126 "=d"(a[LSW128]) 127 : 128 "0"(a[MSW128]), 129 "1"(a[NMSW128]), 130 "2"(a[NLSW128]), 131 "3"(a[LSW128])); 132 } 133 134 /* Generalized 128-bit shifters: 135 136 These bit-shift to a multiple of 32, then move whole longwords. */ 137 138 static inline void lsl128(unsigned int count, int128 a) 139 { 140 int wordcount, i; 141 142 if (count % 32) 143 _lsl128(count % 32, a); 144 145 if (0 == (wordcount = count / 32)) 146 return; 147 148 /* argh, gak, endian-sensitive */ 149 for (i = 0; i < 4 - wordcount; i++) { 150 a[i] = a[i + wordcount]; 151 } 152 for (i = 3; i >= 4 - wordcount; --i) { 153 a[i] = 0; 154 } 155 } 156 157 static inline void lsr128(unsigned int count, int128 a) 158 { 159 int wordcount, i; 160 161 if (count % 32) 162 _lsr128(count % 32, a); 163 164 if (0 == (wordcount = count / 32)) 165 return; 166 167 for (i = 3; i >= wordcount; --i) { 168 a[i] = a[i - wordcount]; 169 } 170 for (i = 0; i < wordcount; i++) { 171 a[i] = 0; 172 } 173 } 174 175 static inline int orl128(int a, int128 b) 176 { 177 b[LSW128] |= a; 178 } 179 180 static inline int btsthi128(const int128 a) 181 { 182 return a[MSW128] & 0x80000000; 183 } 184 185 /* test bits (numbered from 0 = LSB) up to and including "top" */ 186 static inline int bftestlo128(int top, const int128 a) 187 { 188 int r = 0; 189 190 if (top > 31) 191 r |= a[LSW128]; 192 if (top > 63) 193 r |= a[NLSW128]; 194 if (top > 95) 195 r |= a[NMSW128]; 196 197 r |= a[3 - (top / 32)] & ((1 << (top % 32 + 1)) - 1); 198 199 return (r != 0); 200 } 201 202 /* Aargh. We need these because GCC is broken */ 203 /* FIXME: do them in assembly, for goodness' sake! */ 204 static inline void mask64(int pos, unsigned long long *mask) 205 { 206 *mask = 0; 207 208 if (pos < 32) { 209 LO_WORD(*mask) = (1 << pos) - 1; 210 return; 211 } 212 LO_WORD(*mask) = -1; 213 HI_WORD(*mask) = (1 << (pos - 32)) - 1; 214 } 215 216 static inline void bset64(int pos, unsigned long long *dest) 217 { 218 /* This conditional will be optimized away. Thanks, GCC! */ 219 if (pos < 32) 220 asm volatile ("bset %1,%0":"=m" 221 (LO_WORD(*dest)):"id"(pos)); 222 else 223 asm volatile ("bset %1,%0":"=m" 224 (HI_WORD(*dest)):"id"(pos - 32)); 225 } 226 227 static inline int btst64(int pos, unsigned long long dest) 228 { 229 if (pos < 32) 230 return (0 != (LO_WORD(dest) & (1 << pos))); 231 else 232 return (0 != (HI_WORD(dest) & (1 << (pos - 32)))); 233 } 234 235 static inline void lsl64(int count, unsigned long long *dest) 236 { 237 if (count < 32) { 238 HI_WORD(*dest) = (HI_WORD(*dest) << count) 239 | (LO_WORD(*dest) >> count); 240 LO_WORD(*dest) <<= count; 241 return; 242 } 243 count -= 32; 244 HI_WORD(*dest) = LO_WORD(*dest) << count; 245 LO_WORD(*dest) = 0; 246 } 247 248 static inline void lsr64(int count, unsigned long long *dest) 249 { 250 if (count < 32) { 251 LO_WORD(*dest) = (LO_WORD(*dest) >> count) 252 | (HI_WORD(*dest) << (32 - count)); 253 HI_WORD(*dest) >>= count; 254 return; 255 } 256 count -= 32; 257 LO_WORD(*dest) = HI_WORD(*dest) >> count; 258 HI_WORD(*dest) = 0; 259 } 260 #endif 261 262 static inline void fp_denormalize(struct fp_ext *reg, unsigned int cnt) 263 { 264 reg->exp += cnt; 265 266 switch (cnt) { 267 case 0 ... 8: 268 reg->lowmant = reg->mant.m32[1] << (8 - cnt); 269 reg->mant.m32[1] = (reg->mant.m32[1] >> cnt) | 270 (reg->mant.m32[0] << (32 - cnt)); 271 reg->mant.m32[0] = reg->mant.m32[0] >> cnt; 272 break; 273 case 9 ... 32: 274 reg->lowmant = reg->mant.m32[1] >> (cnt - 8); 275 if (reg->mant.m32[1] << (40 - cnt)) 276 reg->lowmant |= 1; 277 reg->mant.m32[1] = (reg->mant.m32[1] >> cnt) | 278 (reg->mant.m32[0] << (32 - cnt)); 279 reg->mant.m32[0] = reg->mant.m32[0] >> cnt; 280 break; 281 case 33 ... 39: 282 asm volatile ("bfextu %1{%2,#8},%0" : "=d" (reg->lowmant) 283 : "m" (reg->mant.m32[0]), "d" (64 - cnt)); 284 if (reg->mant.m32[1] << (40 - cnt)) 285 reg->lowmant |= 1; 286 reg->mant.m32[1] = reg->mant.m32[0] >> (cnt - 32); 287 reg->mant.m32[0] = 0; 288 break; 289 case 40 ... 71: 290 reg->lowmant = reg->mant.m32[0] >> (cnt - 40); 291 if ((reg->mant.m32[0] << (72 - cnt)) || reg->mant.m32[1]) 292 reg->lowmant |= 1; 293 reg->mant.m32[1] = reg->mant.m32[0] >> (cnt - 32); 294 reg->mant.m32[0] = 0; 295 break; 296 default: 297 reg->lowmant = reg->mant.m32[0] || reg->mant.m32[1]; 298 reg->mant.m32[0] = 0; 299 reg->mant.m32[1] = 0; 300 break; 301 } 302 } 303 304 static inline int fp_overnormalize(struct fp_ext *reg) 305 { 306 int shift; 307 308 if (reg->mant.m32[0]) { 309 asm ("bfffo %1{#0,#32},%0" : "=d" (shift) : "dm" (reg->mant.m32[0])); 310 reg->mant.m32[0] = (reg->mant.m32[0] << shift) | (reg->mant.m32[1] >> (32 - shift)); 311 reg->mant.m32[1] = (reg->mant.m32[1] << shift); 312 } else { 313 asm ("bfffo %1{#0,#32},%0" : "=d" (shift) : "dm" (reg->mant.m32[1])); 314 reg->mant.m32[0] = (reg->mant.m32[1] << shift); 315 reg->mant.m32[1] = 0; 316 shift += 32; 317 } 318 319 return shift; 320 } 321 322 static inline int fp_addmant(struct fp_ext *dest, struct fp_ext *src) 323 { 324 int carry; 325 326 /* we assume here, gcc only insert move and a clr instr */ 327 asm volatile ("add.b %1,%0" : "=d,g" (dest->lowmant) 328 : "g,d" (src->lowmant), "0,0" (dest->lowmant)); 329 asm volatile ("addx.l %1,%0" : "=d" (dest->mant.m32[1]) 330 : "d" (src->mant.m32[1]), "0" (dest->mant.m32[1])); 331 asm volatile ("addx.l %1,%0" : "=d" (dest->mant.m32[0]) 332 : "d" (src->mant.m32[0]), "0" (dest->mant.m32[0])); 333 asm volatile ("addx.l %0,%0" : "=d" (carry) : "0" (0)); 334 335 return carry; 336 } 337 338 static inline int fp_addcarry(struct fp_ext *reg) 339 { 340 if (++reg->exp == 0x7fff) { 341 if (reg->mant.m64) 342 fp_set_sr(FPSR_EXC_INEX2); 343 reg->mant.m64 = 0; 344 fp_set_sr(FPSR_EXC_OVFL); 345 return 0; 346 } 347 reg->lowmant = (reg->mant.m32[1] << 7) | (reg->lowmant ? 1 : 0); 348 reg->mant.m32[1] = (reg->mant.m32[1] >> 1) | 349 (reg->mant.m32[0] << 31); 350 reg->mant.m32[0] = (reg->mant.m32[0] >> 1) | 0x80000000; 351 352 return 1; 353 } 354 355 static inline void fp_submant(struct fp_ext *dest, struct fp_ext *src1, 356 struct fp_ext *src2) 357 { 358 /* we assume here, gcc only insert move and a clr instr */ 359 asm volatile ("sub.b %1,%0" : "=d,g" (dest->lowmant) 360 : "g,d" (src2->lowmant), "0,0" (src1->lowmant)); 361 asm volatile ("subx.l %1,%0" : "=d" (dest->mant.m32[1]) 362 : "d" (src2->mant.m32[1]), "0" (src1->mant.m32[1])); 363 asm volatile ("subx.l %1,%0" : "=d" (dest->mant.m32[0]) 364 : "d" (src2->mant.m32[0]), "0" (src1->mant.m32[0])); 365 } 366 367 #define fp_mul64(desth, destl, src1, src2) ({ \ 368 asm ("mulu.l %2,%1:%0" : "=d" (destl), "=d" (desth) \ 369 : "dm" (src1), "0" (src2)); \ 370 }) 371 #define fp_div64(quot, rem, srch, srcl, div) \ 372 asm ("divu.l %2,%1:%0" : "=d" (quot), "=d" (rem) \ 373 : "dm" (div), "1" (srch), "0" (srcl)) 374 #define fp_add64(dest1, dest2, src1, src2) ({ \ 375 asm ("add.l %1,%0" : "=d,dm" (dest2) \ 376 : "dm,d" (src2), "0,0" (dest2)); \ 377 asm ("addx.l %1,%0" : "=d" (dest1) \ 378 : "d" (src1), "0" (dest1)); \ 379 }) 380 #define fp_addx96(dest, src) ({ \ 381 /* we assume here, gcc only insert move and a clr instr */ \ 382 asm volatile ("add.l %1,%0" : "=d,g" (dest->m32[2]) \ 383 : "g,d" (temp.m32[1]), "0,0" (dest->m32[2])); \ 384 asm volatile ("addx.l %1,%0" : "=d" (dest->m32[1]) \ 385 : "d" (temp.m32[0]), "0" (dest->m32[1])); \ 386 asm volatile ("addx.l %1,%0" : "=d" (dest->m32[0]) \ 387 : "d" (0), "0" (dest->m32[0])); \ 388 }) 389 #define fp_sub64(dest, src) ({ \ 390 asm ("sub.l %1,%0" : "=d,dm" (dest.m32[1]) \ 391 : "dm,d" (src.m32[1]), "0,0" (dest.m32[1])); \ 392 asm ("subx.l %1,%0" : "=d" (dest.m32[0]) \ 393 : "d" (src.m32[0]), "0" (dest.m32[0])); \ 394 }) 395 #define fp_sub96c(dest, srch, srcm, srcl) ({ \ 396 char carry; \ 397 asm ("sub.l %1,%0" : "=d,dm" (dest.m32[2]) \ 398 : "dm,d" (srcl), "0,0" (dest.m32[2])); \ 399 asm ("subx.l %1,%0" : "=d" (dest.m32[1]) \ 400 : "d" (srcm), "0" (dest.m32[1])); \ 401 asm ("subx.l %2,%1; scs %0" : "=d" (carry), "=d" (dest.m32[0]) \ 402 : "d" (srch), "1" (dest.m32[0])); \ 403 carry; \ 404 }) 405 406 static inline void fp_multiplymant(union fp_mant128 *dest, struct fp_ext *src1, 407 struct fp_ext *src2) 408 { 409 union fp_mant64 temp; 410 411 fp_mul64(dest->m32[0], dest->m32[1], src1->mant.m32[0], src2->mant.m32[0]); 412 fp_mul64(dest->m32[2], dest->m32[3], src1->mant.m32[1], src2->mant.m32[1]); 413 414 fp_mul64(temp.m32[0], temp.m32[1], src1->mant.m32[0], src2->mant.m32[1]); 415 fp_addx96(dest, temp); 416 417 fp_mul64(temp.m32[0], temp.m32[1], src1->mant.m32[1], src2->mant.m32[0]); 418 fp_addx96(dest, temp); 419 } 420 421 static inline void fp_dividemant(union fp_mant128 *dest, struct fp_ext *src, 422 struct fp_ext *div) 423 { 424 union fp_mant128 tmp; 425 union fp_mant64 tmp64; 426 unsigned long *mantp = dest->m32; 427 unsigned long fix, rem, first, dummy; 428 int i; 429 430 /* the algorithm below requires dest to be smaller than div, 431 but both have the high bit set */ 432 if (src->mant.m64 >= div->mant.m64) { 433 fp_sub64(src->mant, div->mant); 434 *mantp = 1; 435 } else 436 *mantp = 0; 437 mantp++; 438 439 /* basic idea behind this algorithm: we can't divide two 64bit numbers 440 (AB/CD) directly, but we can calculate AB/C0, but this means this 441 quotient is off by C0/CD, so we have to multiply the first result 442 to fix the result, after that we have nearly the correct result 443 and only a few corrections are needed. */ 444 445 /* C0/CD can be precalculated, but it's an 64bit division again, but 446 we can make it a bit easier, by dividing first through C so we get 447 10/1D and now only a single shift and the value fits into 32bit. */ 448 fix = 0x80000000; 449 dummy = div->mant.m32[1] / div->mant.m32[0] + 1; 450 dummy = (dummy >> 1) | fix; 451 fp_div64(fix, dummy, fix, 0, dummy); 452 fix--; 453 454 for (i = 0; i < 3; i++, mantp++) { 455 if (src->mant.m32[0] == div->mant.m32[0]) { 456 fp_div64(first, rem, 0, src->mant.m32[1], div->mant.m32[0]); 457 458 fp_mul64(*mantp, dummy, first, fix); 459 *mantp += fix; 460 } else { 461 fp_div64(first, rem, src->mant.m32[0], src->mant.m32[1], div->mant.m32[0]); 462 463 fp_mul64(*mantp, dummy, first, fix); 464 } 465 466 fp_mul64(tmp.m32[0], tmp.m32[1], div->mant.m32[0], first - *mantp); 467 fp_add64(tmp.m32[0], tmp.m32[1], 0, rem); 468 tmp.m32[2] = 0; 469 470 fp_mul64(tmp64.m32[0], tmp64.m32[1], *mantp, div->mant.m32[1]); 471 fp_sub96c(tmp, 0, tmp64.m32[0], tmp64.m32[1]); 472 473 src->mant.m32[0] = tmp.m32[1]; 474 src->mant.m32[1] = tmp.m32[2]; 475 476 while (!fp_sub96c(tmp, 0, div->mant.m32[0], div->mant.m32[1])) { 477 src->mant.m32[0] = tmp.m32[1]; 478 src->mant.m32[1] = tmp.m32[2]; 479 *mantp += 1; 480 } 481 } 482 } 483 484 #if 0 485 static inline unsigned int fp_fls128(union fp_mant128 *src) 486 { 487 unsigned long data; 488 unsigned int res, off; 489 490 if ((data = src->m32[0])) 491 off = 0; 492 else if ((data = src->m32[1])) 493 off = 32; 494 else if ((data = src->m32[2])) 495 off = 64; 496 else if ((data = src->m32[3])) 497 off = 96; 498 else 499 return 128; 500 501 asm ("bfffo %1{#0,#32},%0" : "=d" (res) : "dm" (data)); 502 return res + off; 503 } 504 505 static inline void fp_shiftmant128(union fp_mant128 *src, int shift) 506 { 507 unsigned long sticky; 508 509 switch (shift) { 510 case 0: 511 return; 512 case 1: 513 asm volatile ("lsl.l #1,%0" 514 : "=d" (src->m32[3]) : "0" (src->m32[3])); 515 asm volatile ("roxl.l #1,%0" 516 : "=d" (src->m32[2]) : "0" (src->m32[2])); 517 asm volatile ("roxl.l #1,%0" 518 : "=d" (src->m32[1]) : "0" (src->m32[1])); 519 asm volatile ("roxl.l #1,%0" 520 : "=d" (src->m32[0]) : "0" (src->m32[0])); 521 return; 522 case 2 ... 31: 523 src->m32[0] = (src->m32[0] << shift) | (src->m32[1] >> (32 - shift)); 524 src->m32[1] = (src->m32[1] << shift) | (src->m32[2] >> (32 - shift)); 525 src->m32[2] = (src->m32[2] << shift) | (src->m32[3] >> (32 - shift)); 526 src->m32[3] = (src->m32[3] << shift); 527 return; 528 case 32 ... 63: 529 shift -= 32; 530 src->m32[0] = (src->m32[1] << shift) | (src->m32[2] >> (32 - shift)); 531 src->m32[1] = (src->m32[2] << shift) | (src->m32[3] >> (32 - shift)); 532 src->m32[2] = (src->m32[3] << shift); 533 src->m32[3] = 0; 534 return; 535 case 64 ... 95: 536 shift -= 64; 537 src->m32[0] = (src->m32[2] << shift) | (src->m32[3] >> (32 - shift)); 538 src->m32[1] = (src->m32[3] << shift); 539 src->m32[2] = src->m32[3] = 0; 540 return; 541 case 96 ... 127: 542 shift -= 96; 543 src->m32[0] = (src->m32[3] << shift); 544 src->m32[1] = src->m32[2] = src->m32[3] = 0; 545 return; 546 case -31 ... -1: 547 shift = -shift; 548 sticky = 0; 549 if (src->m32[3] << (32 - shift)) 550 sticky = 1; 551 src->m32[3] = (src->m32[3] >> shift) | (src->m32[2] << (32 - shift)) | sticky; 552 src->m32[2] = (src->m32[2] >> shift) | (src->m32[1] << (32 - shift)); 553 src->m32[1] = (src->m32[1] >> shift) | (src->m32[0] << (32 - shift)); 554 src->m32[0] = (src->m32[0] >> shift); 555 return; 556 case -63 ... -32: 557 shift = -shift - 32; 558 sticky = 0; 559 if ((src->m32[2] << (32 - shift)) || src->m32[3]) 560 sticky = 1; 561 src->m32[3] = (src->m32[2] >> shift) | (src->m32[1] << (32 - shift)) | sticky; 562 src->m32[2] = (src->m32[1] >> shift) | (src->m32[0] << (32 - shift)); 563 src->m32[1] = (src->m32[0] >> shift); 564 src->m32[0] = 0; 565 return; 566 case -95 ... -64: 567 shift = -shift - 64; 568 sticky = 0; 569 if ((src->m32[1] << (32 - shift)) || src->m32[2] || src->m32[3]) 570 sticky = 1; 571 src->m32[3] = (src->m32[1] >> shift) | (src->m32[0] << (32 - shift)) | sticky; 572 src->m32[2] = (src->m32[0] >> shift); 573 src->m32[1] = src->m32[0] = 0; 574 return; 575 case -127 ... -96: 576 shift = -shift - 96; 577 sticky = 0; 578 if ((src->m32[0] << (32 - shift)) || src->m32[1] || src->m32[2] || src->m32[3]) 579 sticky = 1; 580 src->m32[3] = (src->m32[0] >> shift) | sticky; 581 src->m32[2] = src->m32[1] = src->m32[0] = 0; 582 return; 583 } 584 585 if (shift < 0 && (src->m32[0] || src->m32[1] || src->m32[2] || src->m32[3])) 586 src->m32[3] = 1; 587 else 588 src->m32[3] = 0; 589 src->m32[2] = 0; 590 src->m32[1] = 0; 591 src->m32[0] = 0; 592 } 593 #endif 594 595 static inline void fp_putmant128(struct fp_ext *dest, union fp_mant128 *src, 596 int shift) 597 { 598 unsigned long tmp; 599 600 switch (shift) { 601 case 0: 602 dest->mant.m64 = src->m64[0]; 603 dest->lowmant = src->m32[2] >> 24; 604 if (src->m32[3] || (src->m32[2] << 8)) 605 dest->lowmant |= 1; 606 break; 607 case 1: 608 asm volatile ("lsl.l #1,%0" 609 : "=d" (tmp) : "0" (src->m32[2])); 610 asm volatile ("roxl.l #1,%0" 611 : "=d" (dest->mant.m32[1]) : "0" (src->m32[1])); 612 asm volatile ("roxl.l #1,%0" 613 : "=d" (dest->mant.m32[0]) : "0" (src->m32[0])); 614 dest->lowmant = tmp >> 24; 615 if (src->m32[3] || (tmp << 8)) 616 dest->lowmant |= 1; 617 break; 618 case 31: 619 asm volatile ("lsr.l #1,%1; roxr.l #1,%0" 620 : "=d" (dest->mant.m32[0]) 621 : "d" (src->m32[0]), "0" (src->m32[1])); 622 asm volatile ("roxr.l #1,%0" 623 : "=d" (dest->mant.m32[1]) : "0" (src->m32[2])); 624 asm volatile ("roxr.l #1,%0" 625 : "=d" (tmp) : "0" (src->m32[3])); 626 dest->lowmant = tmp >> 24; 627 if (src->m32[3] << 7) 628 dest->lowmant |= 1; 629 break; 630 case 32: 631 dest->mant.m32[0] = src->m32[1]; 632 dest->mant.m32[1] = src->m32[2]; 633 dest->lowmant = src->m32[3] >> 24; 634 if (src->m32[3] << 8) 635 dest->lowmant |= 1; 636 break; 637 } 638 } 639 640 #if 0 /* old code... */ 641 static inline int fls(unsigned int a) 642 { 643 int r; 644 645 asm volatile ("bfffo %1{#0,#32},%0" 646 : "=d" (r) : "md" (a)); 647 return r; 648 } 649 650 /* fls = "find last set" (cf. ffs(3)) */ 651 static inline int fls128(const int128 a) 652 { 653 if (a[MSW128]) 654 return fls(a[MSW128]); 655 if (a[NMSW128]) 656 return fls(a[NMSW128]) + 32; 657 /* XXX: it probably never gets beyond this point in actual 658 use, but that's indicative of a more general problem in the 659 algorithm (i.e. as per the actual 68881 implementation, we 660 really only need at most 67 bits of precision [plus 661 overflow]) so I'm not going to fix it. */ 662 if (a[NLSW128]) 663 return fls(a[NLSW128]) + 64; 664 if (a[LSW128]) 665 return fls(a[LSW128]) + 96; 666 else 667 return -1; 668 } 669 670 static inline int zerop128(const int128 a) 671 { 672 return !(a[LSW128] | a[NLSW128] | a[NMSW128] | a[MSW128]); 673 } 674 675 static inline int nonzerop128(const int128 a) 676 { 677 return (a[LSW128] | a[NLSW128] | a[NMSW128] | a[MSW128]); 678 } 679 680 /* Addition and subtraction */ 681 /* Do these in "pure" assembly, because "extended" asm is unmanageable 682 here */ 683 static inline void add128(const int128 a, int128 b) 684 { 685 /* rotating carry flags */ 686 unsigned int carry[2]; 687 688 carry[0] = a[LSW128] > (0xffffffff - b[LSW128]); 689 b[LSW128] += a[LSW128]; 690 691 carry[1] = a[NLSW128] > (0xffffffff - b[NLSW128] - carry[0]); 692 b[NLSW128] = a[NLSW128] + b[NLSW128] + carry[0]; 693 694 carry[0] = a[NMSW128] > (0xffffffff - b[NMSW128] - carry[1]); 695 b[NMSW128] = a[NMSW128] + b[NMSW128] + carry[1]; 696 697 b[MSW128] = a[MSW128] + b[MSW128] + carry[0]; 698 } 699 700 /* Note: assembler semantics: "b -= a" */ 701 static inline void sub128(const int128 a, int128 b) 702 { 703 /* rotating borrow flags */ 704 unsigned int borrow[2]; 705 706 borrow[0] = b[LSW128] < a[LSW128]; 707 b[LSW128] -= a[LSW128]; 708 709 borrow[1] = b[NLSW128] < a[NLSW128] + borrow[0]; 710 b[NLSW128] = b[NLSW128] - a[NLSW128] - borrow[0]; 711 712 borrow[0] = b[NMSW128] < a[NMSW128] + borrow[1]; 713 b[NMSW128] = b[NMSW128] - a[NMSW128] - borrow[1]; 714 715 b[MSW128] = b[MSW128] - a[MSW128] - borrow[0]; 716 } 717 718 /* Poor man's 64-bit expanding multiply */ 719 static inline void mul64(unsigned long long a, unsigned long long b, int128 c) 720 { 721 unsigned long long acc; 722 int128 acc128; 723 724 zero128(acc128); 725 zero128(c); 726 727 /* first the low words */ 728 if (LO_WORD(a) && LO_WORD(b)) { 729 acc = (long long) LO_WORD(a) * LO_WORD(b); 730 c[NLSW128] = HI_WORD(acc); 731 c[LSW128] = LO_WORD(acc); 732 } 733 /* Next the high words */ 734 if (HI_WORD(a) && HI_WORD(b)) { 735 acc = (long long) HI_WORD(a) * HI_WORD(b); 736 c[MSW128] = HI_WORD(acc); 737 c[NMSW128] = LO_WORD(acc); 738 } 739 /* The middle words */ 740 if (LO_WORD(a) && HI_WORD(b)) { 741 acc = (long long) LO_WORD(a) * HI_WORD(b); 742 acc128[NMSW128] = HI_WORD(acc); 743 acc128[NLSW128] = LO_WORD(acc); 744 add128(acc128, c); 745 } 746 /* The first and last words */ 747 if (HI_WORD(a) && LO_WORD(b)) { 748 acc = (long long) HI_WORD(a) * LO_WORD(b); 749 acc128[NMSW128] = HI_WORD(acc); 750 acc128[NLSW128] = LO_WORD(acc); 751 add128(acc128, c); 752 } 753 } 754 755 /* Note: unsigned */ 756 static inline int cmp128(int128 a, int128 b) 757 { 758 if (a[MSW128] < b[MSW128]) 759 return -1; 760 if (a[MSW128] > b[MSW128]) 761 return 1; 762 if (a[NMSW128] < b[NMSW128]) 763 return -1; 764 if (a[NMSW128] > b[NMSW128]) 765 return 1; 766 if (a[NLSW128] < b[NLSW128]) 767 return -1; 768 if (a[NLSW128] > b[NLSW128]) 769 return 1; 770 771 return (signed) a[LSW128] - b[LSW128]; 772 } 773 774 inline void div128(int128 a, int128 b, int128 c) 775 { 776 int128 mask; 777 778 /* Algorithm: 779 780 Shift the divisor until it's at least as big as the 781 dividend, keeping track of the position to which we've 782 shifted it, i.e. the power of 2 which we've multiplied it 783 by. 784 785 Then, for this power of 2 (the mask), and every one smaller 786 than it, subtract the mask from the dividend and add it to 787 the quotient until the dividend is smaller than the raised 788 divisor. At this point, divide the dividend and the mask 789 by 2 (i.e. shift one place to the right). Lather, rinse, 790 and repeat, until there are no more powers of 2 left. */ 791 792 /* FIXME: needless to say, there's room for improvement here too. */ 793 794 /* Shift up */ 795 /* XXX: since it just has to be "at least as big", we can 796 probably eliminate this horribly wasteful loop. I will 797 have to prove this first, though */ 798 set128(0, 0, 0, 1, mask); 799 while (cmp128(b, a) < 0 && !btsthi128(b)) { 800 lslone128(b); 801 lslone128(mask); 802 } 803 804 /* Shift down */ 805 zero128(c); 806 do { 807 if (cmp128(a, b) >= 0) { 808 sub128(b, a); 809 add128(mask, c); 810 } 811 lsrone128(mask); 812 lsrone128(b); 813 } while (nonzerop128(mask)); 814 815 /* The remainder is in a... */ 816 } 817 #endif 818 819 #endif /* MULTI_ARITH_H */ 820