URL
https://opencores.org/ocsvn/openrisc_me/openrisc_me/trunk
Subversion Repositories openrisc_me
[/] [openrisc/] [trunk/] [gnu-src/] [gcc-4.5.1/] [libgcc/] [config/] [libbid/] [bid_sqrt_macros.h] - Rev 300
Go to most recent revision | Compare with Previous | Blame | View Log
/* Copyright (C) 2007, 2009 Free Software Foundation, Inc. This file is part of GCC. GCC is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version. GCC 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 General Public License for more details. Under Section 7 of GPL version 3, you are granted additional permissions described in the GCC Runtime Library Exception, version 3.1, as published by the Free Software Foundation. You should have received a copy of the GNU General Public License and a copy of the GCC Runtime Library Exception along with this program; see the files COPYING3 and COPYING.RUNTIME respectively. If not, see <http://www.gnu.org/licenses/>. */ #ifndef _SQRT_MACROS_H_ #define _SQRT_MACROS_H_ #define FENCE __fence #if DOUBLE_EXTENDED_ON extern BINARY80 SQRT80 (BINARY80); __BID_INLINE__ UINT64 short_sqrt128 (UINT128 A10) { BINARY80 lx, ly, l64; int_float f64; // 2^64 f64.i = 0x5f800000; l64 = (BINARY80) f64.d; lx = (BINARY80) A10.w[1] * l64 + (BINARY80) A10.w[0]; ly = SQRT80 (lx); return (UINT64) ly; } __BID_INLINE__ void long_sqrt128 (UINT128 * pCS, UINT256 C256) { UINT256 C4; UINT128 CS; UINT64 X; SINT64 SE; BINARY80 l64, lm64, l128, lxL, lx, ly, lS, lSH, lSL, lE, l3, l2, l1, l0, lp, lCl; int_float fx, f64, fm64; int *ple = (int *) &lx; // 2^64 f64.i = 0x5f800000; l64 = (BINARY80) f64.d; l128 = l64 * l64; lx = l3 = (BINARY80) C256.w[3] * l64 * l128; l2 = (BINARY80) C256.w[2] * l128; lx = FENCE (lx + l2); l1 = (BINARY80) C256.w[1] * l64; lx = FENCE (lx + l1); l0 = (BINARY80) C256.w[0]; lx = FENCE (lx + l0); // sqrt(C256) lS = SQRT80 (lx); // get coefficient // 2^(-64) fm64.i = 0x1f800000; lm64 = (BINARY80) fm64.d; CS.w[1] = (UINT64) (lS * lm64); CS.w[0] = (UINT64) (lS - (BINARY80) CS.w[1] * l64); /////////////////////////////////////// // CAUTION! // little endian code only // add solution for big endian ////////////////////////////////////// lSH = lS; *((UINT64 *) & lSH) &= 0xffffffff00000000ull; // correction for C256 rounding lCl = FENCE (l3 - lx); lCl = FENCE (lCl + l2); lCl = FENCE (lCl + l1); lCl = FENCE (lCl + l0); lSL = lS - lSH; ////////////////////////////////////////// // Watch for compiler re-ordering // ///////////////////////////////////////// // C256-S^2 lxL = FENCE (lx - lSH * lSH); lp = lSH * lSL; lp += lp; lxL = FENCE (lxL - lp); lSL *= lSL; lxL = FENCE (lxL - lSL); lCl += lxL; // correction term lE = lCl / (lS + lS); // get low part of coefficient X = CS.w[0]; if (lCl >= 0) { SE = (SINT64) (lE); CS.w[0] += SE; if (CS.w[0] < X) CS.w[1]++; } else { SE = (SINT64) (-lE); CS.w[0] -= SE; if (CS.w[0] > X) CS.w[1]--; } pCS->w[0] = CS.w[0]; pCS->w[1] = CS.w[1]; } #else extern double sqrt (double); __BID_INLINE__ UINT64 short_sqrt128 (UINT128 A10) { UINT256 ARS, ARS0, AE0, AE, S; UINT64 MY, ES, CY; double lx, l64; int_double f64, ly; int ey, k; // 2^64 f64.i = 0x43f0000000000000ull; l64 = f64.d; lx = (double) A10.w[1] * l64 + (double) A10.w[0]; ly.d = 1.0 / sqrt (lx); MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull; ey = 0x3ff - (ly.i >> 52); // A10*RS^2 __mul_64x128_to_192 (ARS0, MY, A10); __mul_64x192_to_256 (ARS, MY, ARS0); // shr by 2*ey+40, to get a 64-bit value k = (ey << 1) + 104 - 64; if (k >= 128) { if (k > 128) ES = (ARS.w[2] >> (k - 128)) | (ARS.w[3] << (192 - k)); else ES = ARS.w[2]; } else { if (k >= 64) { ARS.w[0] = ARS.w[1]; ARS.w[1] = ARS.w[2]; k -= 64; } if (k) { __shr_128 (ARS, ARS, k); } ES = ARS.w[0]; } ES = ((SINT64) ES) >> 1; if (((SINT64) ES) < 0) { ES = -ES; // A*RS*eps (scaled by 2^64) __mul_64x192_to_256 (AE0, ES, ARS0); AE.w[0] = AE0.w[1]; AE.w[1] = AE0.w[2]; AE.w[2] = AE0.w[3]; __add_carry_out (S.w[0], CY, ARS0.w[0], AE.w[0]); __add_carry_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY); S.w[2] = ARS0.w[2] + AE.w[2] + CY; } else { // A*RS*eps (scaled by 2^64) __mul_64x192_to_256 (AE0, ES, ARS0); AE.w[0] = AE0.w[1]; AE.w[1] = AE0.w[2]; AE.w[2] = AE0.w[3]; __sub_borrow_out (S.w[0], CY, ARS0.w[0], AE.w[0]); __sub_borrow_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY); S.w[2] = ARS0.w[2] - AE.w[2] - CY; } k = ey + 51; if (k >= 64) { if (k >= 128) { S.w[0] = S.w[2]; S.w[1] = 0; k -= 128; } else { S.w[0] = S.w[1]; S.w[1] = S.w[2]; } k -= 64; } if (k) { __shr_128 (S, S, k); } return (UINT64) ((S.w[0] + 1) >> 1); } __BID_INLINE__ void long_sqrt128 (UINT128 * pCS, UINT256 C256) { UINT512 ARS0, ARS; UINT256 ARS00, AE, AE2, S; UINT128 ES, ES2, ARS1; UINT64 ES32, CY, MY; double l64, l128, lx, l2, l1, l0; int_double f64, ly; int ey, k, k2; // 2^64 f64.i = 0x43f0000000000000ull; l64 = f64.d; l128 = l64 * l64; lx = (double) C256.w[3] * l64 * l128; l2 = (double) C256.w[2] * l128; lx = FENCE (lx + l2); l1 = (double) C256.w[1] * l64; lx = FENCE (lx + l1); l0 = (double) C256.w[0]; lx = FENCE (lx + l0); // sqrt(C256) ly.d = 1.0 / sqrt (lx); MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull; ey = 0x3ff - (ly.i >> 52); // A10*RS^2, scaled by 2^(2*ey+104) __mul_64x256_to_320 (ARS0, MY, C256); __mul_64x320_to_384 (ARS, MY, ARS0); // shr by k=(2*ey+104)-128 // expect k is in the range (192, 256) if result in [10^33, 10^34) // apply an additional signed shift by 1 at the same time (to get eps=eps0/2) k = (ey << 1) + 104 - 128 - 192; k2 = 64 - k; ES.w[0] = (ARS.w[3] >> (k + 1)) | (ARS.w[4] << (k2 - 1)); ES.w[1] = (ARS.w[4] >> k) | (ARS.w[5] << k2); ES.w[1] = ((SINT64) ES.w[1]) >> 1; // A*RS >> 192 (for error term computation) ARS1.w[0] = ARS0.w[3]; ARS1.w[1] = ARS0.w[4]; // A*RS>>64 ARS00.w[0] = ARS0.w[1]; ARS00.w[1] = ARS0.w[2]; ARS00.w[2] = ARS0.w[3]; ARS00.w[3] = ARS0.w[4]; if (((SINT64) ES.w[1]) < 0) { ES.w[0] = -ES.w[0]; ES.w[1] = -ES.w[1]; if (ES.w[0]) ES.w[1]--; // A*RS*eps __mul_128x128_to_256 (AE, ES, ARS1); __add_carry_out (S.w[0], CY, ARS00.w[0], AE.w[0]); __add_carry_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY); __add_carry_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY); S.w[3] = ARS00.w[3] + AE.w[3] + CY; } else { // A*RS*eps __mul_128x128_to_256 (AE, ES, ARS1); __sub_borrow_out (S.w[0], CY, ARS00.w[0], AE.w[0]); __sub_borrow_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY); __sub_borrow_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY); S.w[3] = ARS00.w[3] - AE.w[3] - CY; } // 3/2*eps^2, scaled by 2^128 ES32 = ES.w[1] + (ES.w[1] >> 1); __mul_64x64_to_128 (ES2, ES32, ES.w[1]); // A*RS*3/2*eps^2 __mul_128x128_to_256 (AE2, ES2, ARS1); // result, scaled by 2^(ey+52-64) __add_carry_out (S.w[0], CY, S.w[0], AE2.w[0]); __add_carry_in_out (S.w[1], CY, S.w[1], AE2.w[1], CY); __add_carry_in_out (S.w[2], CY, S.w[2], AE2.w[2], CY); S.w[3] = S.w[3] + AE2.w[3] + CY; // k in (0, 64) k = ey + 51 - 128; k2 = 64 - k; S.w[0] = (S.w[1] >> k) | (S.w[2] << k2); S.w[1] = (S.w[2] >> k) | (S.w[3] << k2); // round to nearest S.w[0]++; if (!S.w[0]) S.w[1]++; pCS->w[0] = (S.w[1] << 63) | (S.w[0] >> 1); pCS->w[1] = S.w[1] >> 1; } #endif #endif
Go to most recent revision | Compare with Previous | Blame | View Log