libgcc/config/libbid/bid_sqrt_macros.h

   1 /* Copyright (C) 2007-2016 Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 3, or (at your option) any later
   8 version.
   9
  10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13 for more details.
  14
  15 Under Section 7 of GPL version 3, you are granted additional
  16 permissions described in the GCC Runtime Library Exception, version
  17 3.1, as published by the Free Software Foundation.
  18
  19 You should have received a copy of the GNU General Public License and
  20 a copy of the GCC Runtime Library Exception along with this program;
  21 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  22 <http://www.gnu.org/licenses/>.  */
  23
  24 #ifndef _SQRT_MACROS_H_
  25 #define _SQRT_MACROS_H_
  26
  27 #define FENCE __fence
  28
  29 #if DOUBLE_EXTENDED_ON
  30
  31 extern BINARY80 SQRT80 (BINARY80);
  32
  33
  34 __BID_INLINE__ UINT64
  35 short_sqrt128 (UINT128 A10) {
  36   BINARY80 lx, ly, l64;
  37   int_float f64;
  38
  39   // 2^64
  40   f64.i = 0x5f800000;
  41   l64 = (BINARY80) f64.d;
  42   lx = (BINARY80) A10.w[1] * l64 + (BINARY80) A10.w[0];
  43   ly = SQRT80 (lx);
  44   return (UINT64) ly;
  45 }
  46
  47
  48 __BID_INLINE__ void
  49 long_sqrt128 (UINT128 * pCS, UINT256 C256) {
  50   UINT256 C4;
  51   UINT128 CS;
  52   UINT64 X;
  53   SINT64 SE;
  54   BINARY80 l64, lm64, l128, lxL, lx, ly, lS, lSH, lSL, lE, l3, l2,
  55     l1, l0, lp, lCl;
  56   int_float fx, f64, fm64;
  57   int *ple = (int *) &lx;
  58
  59   // 2^64
  60   f64.i = 0x5f800000;
  61   l64 = (BINARY80) f64.d;
  62
  63   l128 = l64 * l64;
  64   lx = l3 = (BINARY80) C256.w[3] * l64 * l128;
  65   l2 = (BINARY80) C256.w[2] * l128;
  66   lx = FENCE (lx + l2);
  67   l1 = (BINARY80) C256.w[1] * l64;
  68   lx = FENCE (lx + l1);
  69   l0 = (BINARY80) C256.w[0];
  70   lx = FENCE (lx + l0);
  71   // sqrt(C256)
  72   lS = SQRT80 (lx);
  73
  74   // get coefficient
  75   // 2^(-64)
  76   fm64.i = 0x1f800000;
  77   lm64 = (BINARY80) fm64.d;
  78   CS.w[1] = (UINT64) (lS * lm64);
  79   CS.w[0] = (UINT64) (lS - (BINARY80) CS.w[1] * l64);
  80
  81   ///////////////////////////////////////
  82   //  CAUTION!
  83   //  little endian code only
  84   //  add solution for big endian
  85   //////////////////////////////////////
  86   lSH = lS;
  87   *((UINT64 *) & lSH) &= 0xffffffff00000000ull;
  88
  89   // correction for C256 rounding
  90   lCl = FENCE (l3 - lx);
  91   lCl = FENCE (lCl + l2);
  92   lCl = FENCE (lCl + l1);
  93   lCl = FENCE (lCl + l0);
  94
  95   lSL = lS - lSH;
  96
  97   //////////////////////////////////////////
  98   //   Watch for compiler re-ordering
  99   //
 100   /////////////////////////////////////////
 101   // C256-S^2
 102   lxL = FENCE (lx - lSH * lSH);
 103   lp = lSH * lSL;
 104   lp += lp;
 105   lxL = FENCE (lxL - lp);
 106   lSL *= lSL;
 107   lxL = FENCE (lxL - lSL);
 108   lCl += lxL;
 109
 110   // correction term
 111   lE = lCl / (lS + lS);
 112
 113   // get low part of coefficient
 114   X = CS.w[0];
 115   if (lCl >= 0) {
 116     SE = (SINT64) (lE);
 117     CS.w[0] += SE;
 118     if (CS.w[0] < X)
 119       CS.w[1]++;
 120   } else {
 121     SE = (SINT64) (-lE);
 122     CS.w[0] -= SE;
 123     if (CS.w[0] > X)
 124       CS.w[1]--;
 125   }
 126
 127   pCS->w[0] = CS.w[0];
 128   pCS->w[1] = CS.w[1];
 129 }
 130
 131 #else
 132
 133 extern double sqrt (double);
 134
 135 __BID_INLINE__ UINT64
 136 short_sqrt128 (UINT128 A10) {
 137   UINT256 ARS, ARS0, AE0, AE, S;
 138
 139   UINT64 MY, ES, CY;
 140   double lx, l64;
 141   int_double f64, ly;
 142   int ey, k;
 143
 144   // 2^64
 145   f64.i = 0x43f0000000000000ull;
 146   l64 = f64.d;
 147   lx = (double) A10.w[1] * l64 + (double) A10.w[0];
 148   ly.d = 1.0 / sqrt (lx);
 149
 150   MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
 151   ey = 0x3ff - (ly.i >> 52);
 152
 153   // A10*RS^2
 154   __mul_64x128_to_192 (ARS0, MY, A10);
 155   __mul_64x192_to_256 (ARS, MY, ARS0);
 156
 157   // shr by 2*ey+40, to get a 64-bit value
 158   k = (ey << 1) + 104 - 64;
 159   if (k >= 128) {
 160     if (k > 128)
 161       ES = (ARS.w[2] >> (k - 128)) | (ARS.w[3] << (192 - k));
 162     else
 163       ES = ARS.w[2];
 164   } else {
 165     if (k >= 64) {
 166       ARS.w[0] = ARS.w[1];
 167       ARS.w[1] = ARS.w[2];
 168       k -= 64;
 169     }
 170     if (k) {
 171       __shr_128 (ARS, ARS, k);
 172     }
 173     ES = ARS.w[0];
 174   }
 175
 176   ES = ((SINT64) ES) >> 1;
 177
 178   if (((SINT64) ES) < 0) {
 179     ES = -ES;
 180
 181     // A*RS*eps (scaled by 2^64)
 182     __mul_64x192_to_256 (AE0, ES, ARS0);
 183
 184     AE.w[0] = AE0.w[1];
 185     AE.w[1] = AE0.w[2];
 186     AE.w[2] = AE0.w[3];
 187
 188     __add_carry_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
 189     __add_carry_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
 190     S.w[2] = ARS0.w[2] + AE.w[2] + CY;
 191   } else {
 192     // A*RS*eps (scaled by 2^64)
 193     __mul_64x192_to_256 (AE0, ES, ARS0);
 194
 195     AE.w[0] = AE0.w[1];
 196     AE.w[1] = AE0.w[2];
 197     AE.w[2] = AE0.w[3];
 198
 199     __sub_borrow_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
 200     __sub_borrow_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
 201     S.w[2] = ARS0.w[2] - AE.w[2] - CY;
 202   }
 203
 204   k = ey + 51;
 205
 206   if (k >= 64) {
 207     if (k >= 128) {
 208       S.w[0] = S.w[2];
 209       S.w[1] = 0;
 210       k -= 128;
 211     } else {
 212       S.w[0] = S.w[1];
 213       S.w[1] = S.w[2];
 214     }
 215     k -= 64;
 216   }
 217   if (k) {
 218     __shr_128 (S, S, k);
 219   }
 220
 221
 222   return (UINT64) ((S.w[0] + 1) >> 1);
 223
 224 }
 225
 226
 227
 228 __BID_INLINE__ void
 229 long_sqrt128 (UINT128 * pCS, UINT256 C256) {
 230   UINT512 ARS0, ARS;
 231   UINT256 ARS00, AE, AE2, S;
 232   UINT128 ES, ES2, ARS1;
 233   UINT64 ES32, CY, MY;
 234   double l64, l128, lx, l2, l1, l0;
 235   int_double f64, ly;
 236   int ey, k, k2;
 237
 238   // 2^64
 239   f64.i = 0x43f0000000000000ull;
 240   l64 = f64.d;
 241
 242   l128 = l64 * l64;
 243   lx = (double) C256.w[3] * l64 * l128;
 244   l2 = (double) C256.w[2] * l128;
 245   lx = FENCE (lx + l2);
 246   l1 = (double) C256.w[1] * l64;
 247   lx = FENCE (lx + l1);
 248   l0 = (double) C256.w[0];
 249   lx = FENCE (lx + l0);
 250   // sqrt(C256)
 251   ly.d = 1.0 / sqrt (lx);
 252
 253   MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
 254   ey = 0x3ff - (ly.i >> 52);
 255
 256   // A10*RS^2, scaled by 2^(2*ey+104)
 257   __mul_64x256_to_320 (ARS0, MY, C256);
 258   __mul_64x320_to_384 (ARS, MY, ARS0);
 259
 260   // shr by k=(2*ey+104)-128
 261   // expect k is in the range (192, 256) if result in [10^33, 10^34)
 262   // apply an additional signed shift by 1 at the same time (to get eps=eps0/2)
 263   k = (ey << 1) + 104 - 128 - 192;
 264   k2 = 64 - k;
 265   ES.w[0] = (ARS.w[3] >> (k + 1)) | (ARS.w[4] << (k2 - 1));
 266   ES.w[1] = (ARS.w[4] >> k) | (ARS.w[5] << k2);
 267   ES.w[1] = ((SINT64) ES.w[1]) >> 1;
 268
 269   // A*RS >> 192 (for error term computation)
 270   ARS1.w[0] = ARS0.w[3];
 271   ARS1.w[1] = ARS0.w[4];
 272
 273   // A*RS>>64
 274   ARS00.w[0] = ARS0.w[1];
 275   ARS00.w[1] = ARS0.w[2];
 276   ARS00.w[2] = ARS0.w[3];
 277   ARS00.w[3] = ARS0.w[4];
 278
 279   if (((SINT64) ES.w[1]) < 0) {
 280     ES.w[0] = -ES.w[0];
 281     ES.w[1] = -ES.w[1];
 282     if (ES.w[0])
 283       ES.w[1]--;
 284
 285     // A*RS*eps
 286     __mul_128x128_to_256 (AE, ES, ARS1);
 287
 288     __add_carry_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
 289     __add_carry_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
 290     __add_carry_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
 291     S.w[3] = ARS00.w[3] + AE.w[3] + CY;
 292   } else {
 293     // A*RS*eps
 294     __mul_128x128_to_256 (AE, ES, ARS1);
 295
 296     __sub_borrow_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
 297     __sub_borrow_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
 298     __sub_borrow_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
 299     S.w[3] = ARS00.w[3] - AE.w[3] - CY;
 300   }
 301
 302   // 3/2*eps^2, scaled by 2^128
 303   ES32 = ES.w[1] + (ES.w[1] >> 1);
 304   __mul_64x64_to_128 (ES2, ES32, ES.w[1]);
 305   // A*RS*3/2*eps^2
 306   __mul_128x128_to_256 (AE2, ES2, ARS1);
 307
 308   // result, scaled by 2^(ey+52-64)
 309   __add_carry_out (S.w[0], CY, S.w[0], AE2.w[0]);
 310   __add_carry_in_out (S.w[1], CY, S.w[1], AE2.w[1], CY);
 311   __add_carry_in_out (S.w[2], CY, S.w[2], AE2.w[2], CY);
 312   S.w[3] = S.w[3] + AE2.w[3] + CY;
 313
 314   // k in (0, 64)
 315   k = ey + 51 - 128;
 316   k2 = 64 - k;
 317   S.w[0] = (S.w[1] >> k) | (S.w[2] << k2);
 318   S.w[1] = (S.w[2] >> k) | (S.w[3] << k2);
 319
 320   // round to nearest
 321   S.w[0]++;
 322   if (!S.w[0])
 323     S.w[1]++;
 324
 325   pCS->w[0] = (S.w[1] << 63) | (S.w[0] >> 1);
 326   pCS->w[1] = S.w[1] >> 1;
 327
 328 }
 329
 330 #endif
 331 #endif