diff options
Diffstat (limited to 'sysdeps/powerpc/q_addsub.c')
| -rw-r--r-- | sysdeps/powerpc/q_addsub.c | 549 |
1 files changed, 549 insertions, 0 deletions
diff --git a/sysdeps/powerpc/q_addsub.c b/sysdeps/powerpc/q_addsub.c new file mode 100644 index 0000000000..e4ef6d8165 --- /dev/null +++ b/sysdeps/powerpc/q_addsub.c @@ -0,0 +1,549 @@ +/* Add or subtract two 128-bit floating point values. C prototype. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* Add 'a' to 'b' and put the result in 'result', but treat a[0]=axx, + b[0]=bxx. bxx differs from b[0] only in the high bit, similarly axx. */ +/* Exceptions to raise: + - Invalid (SNaN) + - Invalid (Inf-Inf) + - Overflow + - Underflow + - Inexact + */ + +/* Handle cases where exponent of a or b is maximum. */ +static void +handle_max_exponent(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + const unsigned ax, /* axx >> 16 & 0x7fff. */ + const unsigned bx) /* bxx >> 16 & 0x7fff. */ +{ + int ax_ismax, bx_ismax; + unsigned a1,a2,a3, b1,b2,b3; + int a_zeromant, b_zeromant; + + ax_ismax = ax == 0x7fff; + bx_ismax = bx == 0x7fff; + + assert(ax_ismax || bx_ismax); + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + a_zeromant = (axx & 0xffff | a1 | a2 | a3) == 0; + b_zeromant = (bxx & 0xffff | b1 | b2 | b3) == 0; + + /* Deal with SNaNs. */ + if ( ax_ismax && !a_zeromant && (axx & 0x8000) == 0 + || bx_ismax && !b_zeromant && (bxx & 0x8000) == 0) + { + set_fpscr_bit(FPSCR_VXSNAN); + axx |= 0x8000; /* Demote the SNaN to a QNaN (whichever of */ + bxx |= 0x8000; /* a or b it was). */ + } + /* Deal with Inf-Inf. */ + else if (a_zeromant && b_zeromant && (axx ^ bxx) == 0x80000000) + { + set_fpscr_bit(FPSCR_VXISI); + bxx |= 0x8000; /* Return an appropriate QNaN. */ + } + + /* Return the lexicographically larger of a or b, ignoring the sign + bits. */ + if ((axx & 0x7fffffff) > (bxx & 0x7fffffff)) goto return_a; + else if ((axx & 0x7fffffff) < (bxx & 0x7fffffff)) goto return_b; + else if (a1 > b1) goto return_a; + else if (a1 < b1) goto return_b; + else if (a2 > b2) goto return_a; + else if (a2 < b2) goto return_b; + else if (a3 > b3) goto return_a; /* I've clearly been writing too */ + else if (a3 < b3) goto return_b; /* much Fortran... */ + + /* If they are equal except for the sign bits, return 'b'. */ + +return_b: + result[0] = bxx; result[1] = b1; result[2] = b2; result[3] = b3; + return; + +return_a: + result[0] = axx; result[1] = a1; result[2] = a2; result[3] = a3; + return; +} + +/* Renormalise and output a FP number. */ +static void +renormalise_value(unsigned result[4], + const unsigned axx, + unsigned ax, + unsigned r0, + unsigned r1, + unsigned r2, + unsigned r3) +{ + int rshift; + if (r0 != 0 || cntlzw(a1) < 16 || 32 > ax-1) + { + rshift = cntlzw(r0)-15 + (-(cntlzw(r0) >> 5) & cntlzw(a1)); + assert(rshift < 32); + if (rshift > ax-1) + { + ax--; + rshift = ax; + } + + result[0] = (axx & 0x80000000 + | ax-rshift << 16 + | r0 << rshift & 0xffff + | a1 >> 32-rshift & 0xffff); + result[1] = a1 << rshift | a2 >> 32-rshift; + result[2] = a2 << rshift | a3 >> 32-rshift; + result[3] = a3 << rshift; + return; + } + result[3] = 0; + /* Special case for zero. */ + if (a1 == 0 && a2 == 0 && a3 == 0) + { + result[0] = axx & 0x80000000; + result[1] = result[2] = 0; + return; + } + while (a1 != 0 && cntlzw(a2) >= 16 && 64 <= ax-1) + { + ax -= 32; + a1 = a2; a2 = a3; a3 = 0; + } + rshift = cntlzw(a1)-15 + (-(cntlzw(a1) >> 5) & cntlzw(a2)); + assert(rshift < 32); + if (rshift > ax-1-32) + { + ax--; + rshift = ax-32; + } + + result[0] = (axx & 0x80000000 + | ax-rshift-32 << 16 + | a1 << rshift & 0xffff + | a2 >> 32-rshift & 0xffff); + result[1] = a2 << rshift | a3 >> 32-rshift; + result[2] = a3 << rshift; + return; +} + +/* Handle the case where one or both numbers are denormalised or zero. + This case almost never happens, so we don't slow the main code + with it. */ +static void +handle_min_exponent(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + const unsigned ax, /* axx >> 16 & 0x7fff. */ + const unsigned bx) /* bxx >> 16 & 0x7fff. */ +{ + int ax_denorm, bx_denorm; + unsigned a1,a2,a3, b1,b2,b3; + int a_zeromant, b_zeromant; + + ax_denorm = ax == 0; + bx_denorm = bx == 0; + + assert(ax_denorm || bx_denorm); + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + +} + +/* Add a+b+cin modulo 2^32, put result in 'r' and carry in 'cout'. */ +#define addc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)+(b)+(cin); + (cout) = addc_tmp >> 32; + (r) = addc_tmp; + } + +/* Calculate a+~b+cin modulo 2^32, put result in 'r' and carry in 'cout'. */ +#define subc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)-(b)+(cin)-1; + (cout) = addc_tmp >> 63; + (r) = addc_tmp; + } + +/* Handle the case where both exponents are the same. This requires quite + a different algorithm than the general case. */ +static void +handle_equal_exponents(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + unsigned ax) /* [ab]xx >> 16 & 0x7fff. */ +{ + unsigned a1,a2,a3, b1,b2,b3; + int roundmode; + unsigned carry, r0; + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + if ((int)(axx ^ bxx) >= 0) + { + int roundmode; + + /* Adding. */ + roundmode = fegetround(); + + /* What about overflow? */ + if (ax == 0x7ffe) + { + /* Oh no! Too big! */ + /* Result: + rounding result + -------- ------ + nearest return Inf with sign of a,b + zero return nearest possible non-Inf value with + sign of a,b + +Inf return +Inf if a,b>0, otherwise return + value just before -Inf. + -Inf return +Inf if a,b>0, otherwise return + value just before -Inf. + */ + set_fpscr_bit(FPSCR_OX); + /* Overflow always produces inexact result. */ + set_fpscr_bit(FPSCR_XX); + + if ( roundmode == FE_TONEAREST + || roundmode == FE_UPWARD && (int)axx >= 0 + || roundmode == FE_DOWNWARD && (int)axx < 0) + { + result[3] = result[2] = result[1] = 0; + result[0] = axx & 0xffff0000 | 0x7fff0000; + } + else + { + result[3] = result[2] = result[1] = 0xffffffff; + result[0] = axx & 0xfffe0000 | 0x7ffeffff; + } + return; + } + + /* We need to worry about rounding/inexact here. Do it like this: */ + if (a3 + b3 & 1) + { + /* Need to round. Upwards? */ + set_fpscr_bit(FPSCR_XX); + carry = ( roundmode == FE_NEAREST && (a3 + b3 & 2) != 0 + || roundmode == FE_UPWARD && (int)axx >= 0 + || roundmode == FE_DOWNWARD && (int)axx < 0); + } + else + carry = 0; /* Result will be exact. */ + + /* Perform the addition. */ + addc(a3,carry,a3,b3,carry); + addc(a2,carry,a2,b2,carry); + addc(a1,carry,a1,b1,carry); + r0 = (axx & 0xffff) + (bxx & 0xffff) + carry; + + /* Shift right by 1. */ + result[3] = a3 >> 1 | a2 << 31; + result[2] = a2 >> 1 | a1 << 31; + result[1] = a1 >> 1 | r0 << 31; + /* Exponent of result is exponent of inputs plus 1. + Sign of result is common sign of inputs. */ + result[0] = r0 >> 1 & 0xffff | axx + 0x10000 & 0xffff0000; + } + else + { + /* Subtracting. */ + + /* Perform the subtraction, a-b. */ + subc(a3,carry,a3,b3,0); + subc(a2,carry,a2,b2,carry); + subc(a1,carry,a1,b1,carry); + subc(r0,carry,a0&0xffff,b0&0xffff,carry); + + /* Maybe we should have calculated b-a... */ + if (carry) + { + subc(a3,carry,0,a3,0); + subc(a2,carry,0,a2,carry); + subc(a1,carry,0,a1,carry); + subc(r0,carry,0,r0,carry); + axx ^= 0x80000000; + } + + renormalise_value(result, axx, ax, r0, a1, a2, a3); + } +} + + +static void +add(unsigned result[4], const unsigned a[4], const unsigned b[4], + unsigned axx, unsigned bxx) +{ + int ax, bx, diff, carry; + unsigned a0,a1,a2,a3, b0,b1,b2,b3,b4, sdiff; + + ax = axx >> 16 & 0x7fff; + bx = bxx >> 16 & 0x7fff; + + /* Deal with NaNs and Inf. */ + if (ax == 0x7fff || bx == 0x7fff) + { + handle_max_exponent(result, a, b, axx, bxx, ax, bx); + return; + } + /* Deal with denorms and zero. */ + if (ax == 0 || bx == 0) + { + handle_min_exponent(result, a, b, axx, bxx, ax, bx); + return; + } + /* Finally, one special case, when both exponents are equal. */ + if (ax == bx) + { + handle_equal_exponents(result, a, b, axx, bxx, ax); + return; + } + + sdiff = axx ^ bxx; + /* Swap a and b if b has a larger magnitude than a, so that a will have + the larger magnitude. */ + if (ax < bx) + { + const unsigned *t; + t = b; b = a; a = t; + diff = bx - ax; + ax = bx; + axx = bxx; + } + else + diff = ax - bx; + + a0 = a[0] & 0xffff | 0x10000; a1 = a[1]; a2 = a[2]; a3 = a[3]; + b0 = b[0] & 0xffff | 0x10000; b1 = b[1]; b2 = b[2]; b3 = b[3]; + if (diff < 32) + { + b4 = b3 << 32-diff; + b3 = b3 >> diff | b2 << 32-biff; + b2 = b2 >> diff | b1 << 32-diff; + b1 = b1 >> diff | b0 << 32-diff; + b0 = b0 >> diff; + } + else if (diff < 64) + { + diff -= 32; + b4 = b3 & 1 | b3 >> (diff == 32) | b2 << 32-biff; + b3 = b2 >> diff | b1 << 32-diff; + b2 = b1 >> diff | b0 << 32-diff; + b1 = b0 >> diff; + b0 = 0; + } + else if (diff < 96) + { + b4 = b2 | b3 | b1 << 32-diff; + b3 = b1 >> diff | b0 << 32-diff; + b2 = b0 >> diff; + b1 = b0 = 0; + } + else if (diff < 128) + { + b4 = b1 | b2 | b3 | b0 << 32-diff; + b3 = b0 >> diff; + b2 = b1 = b0 = 0; + } + else + { + b4 = b0|b1|b2|b3; + b3 = b2 = b1 = b0 = 0; + } + + /* Now, two cases: one for addition, one for subtraction. */ + if ((int)sdiff >= 0) + { + /* Addition. */ + + /* + + /* Perform the addition. */ + addc(a3,carry,a3,b3,0); + addc(a2,carry,a2,b2,carry); + addc(a1,carry,a1,b1,carry); + addc(a0,carry,a0,b0,carry); + + + + if (a0 & 0x20000) + { + /* Need to renormalise by shifting right. */ + /* Shift right by 1. */ + b4 = b4 | a3 << 31; + a3 = a3 >> 1 | a2 << 31; + a2 = a2 >> 1 | a1 << 31; + result[1] = a1 >> 1 | r0 << 31; + /* Exponent of result is exponent of inputs plus 1. + Sign of result is common sign of inputs. */ + result[0] = r0 >> 1 & 0xffff | axx + 0x10000 & 0xffff0000; + } + + + } + else + { + /* Subtraction. */ + + } +} + +/* Add the absolute values of two 128-bit floating point values, + give the result the sign of one of them. The only exception this + can raise is for SNaN. */ +static void +aadd(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + unsigned ax, bx, xd; + const unsigned *sml; + unsigned t0,t1,t2,t3,tx, s0,s1,s2,s3,s4, carry; + int rmode, xdelta, shift; + + ax = a[0] >> 16 & 0x7fff; + bx = b[0] >> 16 & 0x7fff; + + /* Deal with . */ + if (ax == 0x7fff) + { + t0 = a[0]; t1 = a[1]; t2 = a[2]; t3 = a[3]; + /* Check for SNaN. */ + if ((t0 & 0x8000) == 0 + && (t0 & 0x7fff | t1 | t2 | t3) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Return b. */ + result[0] = t0; result[1] = t1; result[2] = t2; result[3] = t3; + return; + } + /* Deal with b==Inf or b==NaN. */ + if (bx == 0x7fff) + { + t0 = b[0]; t1 = b[1]; t2 = b[2]; t3 = b[3]; + /* Check for SNaN. */ + if ((t0 & 0x8000) == 0 + && (t0 & 0x7fff | t1 | t2 | t3) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Return b. */ + result[0] = t0; result[1] = t1; result[2] = t2; result[3] = t3; + return; + } + + /* Choose the larger of the two to be 't', and the smaller to be 's'. */ + if (ax > bx) + { + t0 = a[0] & 0xffff | (ax != 0) << 16; + t1 = a[1]; t2 = a[2]; t3 = a[3]; tx = ax; + s0 = b[0] & 0xffff | (bx != 0) << 16; + s1 = b[1]; s2 = b[2]; s3 = b[3]; + xd = ax-bx; + } + else + { + t0 = b[0] & 0xffff | (bx != 0) << 16; + t1 = b[1]; t2 = b[2]; t3 = b[3]; tx = bx; + s0 = a[0] & 0xffff | (ax != 0) << 16; + s1 = a[1]; s2 = a[2]; s3 = a[3]; + sml = a; + xd = bx-ax; + } + + /* Shift 's2' right by 'xd' bits. */ + switch (xd >> 5) + { + case 0: + s4 = 0; + break; + case 1: + s4 = s3; s3 = s2; s2 = s1; s1 = s0; s0 = 0; + break; + case 2: + s4 = s2 | s3 != 0; + s3 = s1; s2 = s0; s1 = s0 = 0; + break; + case 3: + s4 = s1 | (s3|s2) != 0; + s3 = s0; s2 = s1 = s0 = 0; + break; + default: + s4 = s0 | (s3|s2|s1) != 0; + s3 = s2 = s1 = s0 = 0; + } + xd = xd & 0x1f; + if (xd != 0) + { + s4 = s4 >> xd | (s4 << 32-xd) != 0 | s3 << 32-xd; + s3 = s3 >> xd | s2 << 32-xd; + s2 = s2 >> xd | s1 << 32-xd; + s1 = s1 >> xd | s0 << 32-xd; + s0 = s0 >> xd; + } + + /* Do the addition. */ +#define addc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)+(b)+(cin); + (cout) = addc_tmp >> 32; + (r) = addc_tmp; + } + addc(t3,carry,t3,s3,0); + addc(t2,carry,t2,s2,carry); + addc(t1,carry,t1,s1,carry); + t0 = t0 + s0 + carry; + + /* Renormalise. */ + xdelta = 15-cntlzw(t0); + if (tx + xdelta <= 0x7fff) + shift = xdelta; + else + { + } +} + +/* Add two 128-bit floating point values. */ +void +__q_add(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + if ((a[0] ^ b[0]) >= 0) + aadd(result, a, b); + else + asubtract(result, a, b); +} + +/* Subtract two 128-bit floating point values. */ +void +__q_sub(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + if ((a[0] ^ b[0]) < 0) + aadd(result, a, b); + else + asubtract(result, a, b); +} |