Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 1 | /* This file is distributed under the University of Illinois Open Source |
| 2 | * License. See LICENSE.TXT for details. |
| 3 | */ |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 4 | |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 5 | /* long double __gcc_qdiv(long double x, long double y); |
| 6 | * This file implements the PowerPC 128-bit double-double division operation. |
| 7 | * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) |
| 8 | */ |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 9 | |
| 10 | #include "DD.h" |
| 11 | |
| 12 | long double __gcc_qdiv(long double a, long double b) |
| 13 | { |
| 14 | static const uint32_t infinityHi = UINT32_C(0x7ff00000); |
| 15 | DD dst = { .ld = a }, src = { .ld = b }; |
| 16 | |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 17 | register double x = dst.s.hi, x1 = dst.s.lo, |
| 18 | y = src.s.hi, y1 = src.s.lo; |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 19 | |
| 20 | double yHi, yLo, qHi, qLo; |
| 21 | double yq, tmp, q; |
| 22 | |
| 23 | q = x / y; |
| 24 | |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 25 | /* Detect special cases */ |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 26 | if (q == 0.0) { |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 27 | dst.s.hi = q; |
| 28 | dst.s.lo = 0.0; |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 29 | return dst.ld; |
| 30 | } |
| 31 | |
| 32 | const doublebits qBits = { .d = q }; |
| 33 | if (((uint32_t)(qBits.x >> 32) & infinityHi) == infinityHi) { |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 34 | dst.s.hi = q; |
| 35 | dst.s.lo = 0.0; |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 36 | return dst.ld; |
| 37 | } |
| 38 | |
| 39 | yHi = high26bits(y); |
| 40 | qHi = high26bits(q); |
| 41 | |
| 42 | yq = y * q; |
| 43 | yLo = y - yHi; |
| 44 | qLo = q - qHi; |
| 45 | |
| 46 | tmp = LOWORDER(yq, yHi, yLo, qHi, qLo); |
| 47 | tmp = (x - yq) - tmp; |
| 48 | tmp = ((tmp + x1) - y1 * q) / y; |
| 49 | x = q + tmp; |
| 50 | |
Edward O'Callaghan | ccf4813 | 2009-08-09 18:41:02 +0000 | [diff] [blame] | 51 | dst.s.lo = (q - x) + tmp; |
| 52 | dst.s.hi = x; |
Daniel Dunbar | fd08999 | 2009-06-26 16:47:03 +0000 | [diff] [blame] | 53 | |
| 54 | return dst.ld; |
| 55 | } |