diff options
Diffstat (limited to 'recipes-devtools/eglibc/files/glibc.fix_sqrt.patch')
| -rw-r--r-- | recipes-devtools/eglibc/files/glibc.fix_sqrt.patch | 500 |
1 files changed, 0 insertions, 500 deletions
diff --git a/recipes-devtools/eglibc/files/glibc.fix_sqrt.patch b/recipes-devtools/eglibc/files/glibc.fix_sqrt.patch deleted file mode 100644 index a648b4f..0000000 --- a/recipes-devtools/eglibc/files/glibc.fix_sqrt.patch +++ /dev/null | |||
| @@ -1,500 +0,0 @@ | |||
| 1 | Implements sqrt with fsqrte. | ||
| 2 | Adds functions for 603e (also used for e300c3, e600 and e500mc) | ||
| 3 | Adds functions for e5500 in 64 bit mode (also used for e6500) | ||
| 4 | |||
| 5 | diff -ruN libc-orig/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c libc-e5500/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c | ||
| 6 | --- libc-orig/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c 1969-12-31 18:00:00.000000000 -0600 | ||
| 7 | +++ libc-e5500/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c 2011-05-18 13:00:40.423986787 -0500 | ||
| 8 | @@ -0,0 +1,134 @@ | ||
| 9 | +/* Double-precision floating point square root. | ||
| 10 | + Copyright (C) 2010 Free Software Foundation, Inc. | ||
| 11 | + This file is part of the GNU C Library. | ||
| 12 | + | ||
| 13 | + The GNU C Library is free software; you can redistribute it and/or | ||
| 14 | + modify it under the terms of the GNU Lesser General Public | ||
| 15 | + License as published by the Free Software Foundation; either | ||
| 16 | + version 2.1 of the License, or (at your option) any later version. | ||
| 17 | + | ||
| 18 | + The GNU C Library is distributed in the hope that it will be useful, | ||
| 19 | + but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 20 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
| 21 | + Lesser General Public License for more details. | ||
| 22 | + | ||
| 23 | + You should have received a copy of the GNU Lesser General Public | ||
| 24 | + License along with the GNU C Library; if not, write to the Free | ||
| 25 | + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA | ||
| 26 | + 02111-1307 USA. */ | ||
| 27 | + | ||
| 28 | +#include <math.h> | ||
| 29 | +#include <math_private.h> | ||
| 30 | +#include <fenv_libc.h> | ||
| 31 | +#include <inttypes.h> | ||
| 32 | + | ||
| 33 | +#include <sysdep.h> | ||
| 34 | +#include <ldsodefs.h> | ||
| 35 | + | ||
| 36 | +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 }; | ||
| 37 | +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 }; | ||
| 38 | +static const float two108 = 3.245185536584267269e+32; | ||
| 39 | +static const float twom54 = 5.551115123125782702e-17; | ||
| 40 | +static const float half = 0.5; | ||
| 41 | + | ||
| 42 | +/* The method is based on the descriptions in: | ||
| 43 | + | ||
| 44 | + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5; | ||
| 45 | + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9 | ||
| 46 | + | ||
| 47 | + We find the actual square root and half of its reciprocal | ||
| 48 | + simultaneously. */ | ||
| 49 | + | ||
| 50 | +#ifdef __STDC__ | ||
| 51 | +double | ||
| 52 | +__ieee754_sqrt (double b) | ||
| 53 | +#else | ||
| 54 | +double | ||
| 55 | +__ieee754_sqrt (b) | ||
| 56 | + double b; | ||
| 57 | +#endif | ||
| 58 | +{ | ||
| 59 | + if (__builtin_expect (b > 0, 1)) | ||
| 60 | + { | ||
| 61 | + double y, g, h, d, r; | ||
| 62 | + ieee_double_shape_type u; | ||
| 63 | + | ||
| 64 | + if (__builtin_expect (b != a_inf.value, 1)) | ||
| 65 | + { | ||
| 66 | + fenv_t fe; | ||
| 67 | + | ||
| 68 | + fe = fegetenv_register (); | ||
| 69 | + | ||
| 70 | + u.value = b; | ||
| 71 | + | ||
| 72 | + relax_fenv_state (); | ||
| 73 | + | ||
| 74 | + __asm__ ("frsqrte %[estimate], %[x]\n" | ||
| 75 | + : [estimate] "=f" (y) : [x] "f" (b)); | ||
| 76 | + | ||
| 77 | + /* Following Muller et al, page 168, equation 5.20. | ||
| 78 | + | ||
| 79 | + h goes to 1/(2*sqrt(b)) | ||
| 80 | + g goes to sqrt(b). | ||
| 81 | + | ||
| 82 | + We need three iterations to get within 1ulp. */ | ||
| 83 | + | ||
| 84 | + /* Indicate that these can be performed prior to the branch. GCC | ||
| 85 | + insists on sinking them below the branch, however; it seems like | ||
| 86 | + they'd be better before the branch so that we can cover any latency | ||
| 87 | + from storing the argument and loading its high word. Oh well. */ | ||
| 88 | + | ||
| 89 | + g = b * y; | ||
| 90 | + h = 0.5 * y; | ||
| 91 | + | ||
| 92 | + /* Handle small numbers by scaling. */ | ||
| 93 | + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0)) | ||
| 94 | + return __ieee754_sqrt (b * two108) * twom54; | ||
| 95 | + | ||
| 96 | +#define FMADD(a_, c_, b_) \ | ||
| 97 | + ({ double __r; \ | ||
| 98 | + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \ | ||
| 99 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 100 | + __r;}) | ||
| 101 | +#define FNMSUB(a_, c_, b_) \ | ||
| 102 | + ({ double __r; \ | ||
| 103 | + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 104 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 105 | + __r;}) | ||
| 106 | + | ||
| 107 | + r = FNMSUB (g, h, half); | ||
| 108 | + g = FMADD (g, r, g); | ||
| 109 | + h = FMADD (h, r, h); | ||
| 110 | + | ||
| 111 | + r = FNMSUB (g, h, half); | ||
| 112 | + g = FMADD (g, r, g); | ||
| 113 | + h = FMADD (h, r, h); | ||
| 114 | + | ||
| 115 | + r = FNMSUB (g, h, half); | ||
| 116 | + g = FMADD (g, r, g); | ||
| 117 | + h = FMADD (h, r, h); | ||
| 118 | + | ||
| 119 | + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */ | ||
| 120 | + | ||
| 121 | + /* Final refinement. */ | ||
| 122 | + d = FNMSUB (g, g, b); | ||
| 123 | + | ||
| 124 | + fesetenv_register (fe); | ||
| 125 | + return FMADD (d, h, g); | ||
| 126 | + } | ||
| 127 | + } | ||
| 128 | + else if (b < 0) | ||
| 129 | + { | ||
| 130 | + /* For some reason, some PowerPC32 processors don't implement | ||
| 131 | + FE_INVALID_SQRT. */ | ||
| 132 | +#ifdef FE_INVALID_SQRT | ||
| 133 | + feraiseexcept (FE_INVALID_SQRT); | ||
| 134 | + | ||
| 135 | + fenv_union_t u = { .fenv = fegetenv_register () }; | ||
| 136 | + if ((u.l[1] & FE_INVALID) == 0) | ||
| 137 | +#endif | ||
| 138 | + feraiseexcept (FE_INVALID); | ||
| 139 | + b = a_nan.value; | ||
| 140 | + } | ||
| 141 | + return f_wash (b); | ||
| 142 | +} | ||
| 143 | diff -ruN libc-orig/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c libc-e5500/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c | ||
| 144 | --- libc-orig/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c 1969-12-31 18:00:00.000000000 -0600 | ||
| 145 | +++ libc-e5500/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c 2011-05-18 13:00:40.423986787 -0500 | ||
| 146 | @@ -0,0 +1,101 @@ | ||
| 147 | +/* Single-precision floating point square root. | ||
| 148 | + Copyright (C) 2010 Free Software Foundation, Inc. | ||
| 149 | + This file is part of the GNU C Library. | ||
| 150 | + | ||
| 151 | + The GNU C Library is free software; you can redistribute it and/or | ||
| 152 | + modify it under the terms of the GNU Lesser General Public | ||
| 153 | + License as published by the Free Software Foundation; either | ||
| 154 | + version 2.1 of the License, or (at your option) any later version. | ||
| 155 | + | ||
| 156 | + The GNU C Library is distributed in the hope that it will be useful, | ||
| 157 | + but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 158 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
| 159 | + Lesser General Public License for more details. | ||
| 160 | + | ||
| 161 | + You should have received a copy of the GNU Lesser General Public | ||
| 162 | + License along with the GNU C Library; if not, write to the Free | ||
| 163 | + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA | ||
| 164 | + 02111-1307 USA. */ | ||
| 165 | + | ||
| 166 | +#include <math.h> | ||
| 167 | +#include <math_private.h> | ||
| 168 | +#include <fenv_libc.h> | ||
| 169 | +#include <inttypes.h> | ||
| 170 | + | ||
| 171 | +#include <sysdep.h> | ||
| 172 | +#include <ldsodefs.h> | ||
| 173 | + | ||
| 174 | +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 }; | ||
| 175 | +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 }; | ||
| 176 | +static const float threehalf = 1.5; | ||
| 177 | + | ||
| 178 | +/* The method is based on the descriptions in: | ||
| 179 | + | ||
| 180 | + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5; | ||
| 181 | + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9 | ||
| 182 | + | ||
| 183 | + We find the reciprocal square root and use that to compute the actual | ||
| 184 | + square root. */ | ||
| 185 | + | ||
| 186 | +#ifdef __STDC__ | ||
| 187 | +float | ||
| 188 | +__ieee754_sqrtf (float b) | ||
| 189 | +#else | ||
| 190 | +float | ||
| 191 | +__ieee754_sqrtf (b) | ||
| 192 | + float b; | ||
| 193 | +#endif | ||
| 194 | +{ | ||
| 195 | + if (__builtin_expect (b > 0, 1)) | ||
| 196 | + { | ||
| 197 | +#define FMSUB(a_, c_, b_) \ | ||
| 198 | + ({ double __r; \ | ||
| 199 | + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 200 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 201 | + __r;}) | ||
| 202 | +#define FNMSUB(a_, c_, b_) \ | ||
| 203 | + ({ double __r; \ | ||
| 204 | + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 205 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 206 | + __r;}) | ||
| 207 | + | ||
| 208 | + if (__builtin_expect (b != a_inf.value, 1)) | ||
| 209 | + { | ||
| 210 | + double y, x; | ||
| 211 | + fenv_t fe; | ||
| 212 | + | ||
| 213 | + fe = fegetenv_register (); | ||
| 214 | + | ||
| 215 | + relax_fenv_state (); | ||
| 216 | + | ||
| 217 | + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */ | ||
| 218 | + y = FMSUB (threehalf, b, b); | ||
| 219 | + | ||
| 220 | + /* Initial estimate. */ | ||
| 221 | + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b)); | ||
| 222 | + | ||
| 223 | + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */ | ||
| 224 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 225 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 226 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 227 | + | ||
| 228 | + /* All done. */ | ||
| 229 | + fesetenv_register (fe); | ||
| 230 | + return x * b; | ||
| 231 | + } | ||
| 232 | + } | ||
| 233 | + else if (b < 0) | ||
| 234 | + { | ||
| 235 | + /* For some reason, some PowerPC32 processors don't implement | ||
| 236 | + FE_INVALID_SQRT. */ | ||
| 237 | +#ifdef FE_INVALID_SQRT | ||
| 238 | + feraiseexcept (FE_INVALID_SQRT); | ||
| 239 | + | ||
| 240 | + fenv_union_t u = { .fenv = fegetenv_register () }; | ||
| 241 | + if ((u.l[1] & FE_INVALID) == 0) | ||
| 242 | +#endif | ||
| 243 | + feraiseexcept (FE_INVALID); | ||
| 244 | + b = a_nan.value; | ||
| 245 | + } | ||
| 246 | + return f_washf (b); | ||
| 247 | +} | ||
| 248 | diff -ruN libc-orig/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c libc-e5500/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c | ||
| 249 | --- libc-orig/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c 1969-12-31 18:00:00.000000000 -0600 | ||
| 250 | +++ libc-e5500/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c 2011-05-18 11:15:22.467987000 -0500 | ||
| 251 | @@ -0,0 +1,134 @@ | ||
| 252 | +/* Double-precision floating point square root. | ||
| 253 | + Copyright (C) 2010 Free Software Foundation, Inc. | ||
| 254 | + This file is part of the GNU C Library. | ||
| 255 | + | ||
| 256 | + The GNU C Library is free software; you can redistribute it and/or | ||
| 257 | + modify it under the terms of the GNU Lesser General Public | ||
| 258 | + License as published by the Free Software Foundation; either | ||
| 259 | + version 2.1 of the License, or (at your option) any later version. | ||
| 260 | + | ||
| 261 | + The GNU C Library is distributed in the hope that it will be useful, | ||
| 262 | + but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 263 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
| 264 | + Lesser General Public License for more details. | ||
| 265 | + | ||
| 266 | + You should have received a copy of the GNU Lesser General Public | ||
| 267 | + License along with the GNU C Library; if not, write to the Free | ||
| 268 | + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA | ||
| 269 | + 02111-1307 USA. */ | ||
| 270 | + | ||
| 271 | +#include <math.h> | ||
| 272 | +#include <math_private.h> | ||
| 273 | +#include <fenv_libc.h> | ||
| 274 | +#include <inttypes.h> | ||
| 275 | + | ||
| 276 | +#include <sysdep.h> | ||
| 277 | +#include <ldsodefs.h> | ||
| 278 | + | ||
| 279 | +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 }; | ||
| 280 | +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 }; | ||
| 281 | +static const float two108 = 3.245185536584267269e+32; | ||
| 282 | +static const float twom54 = 5.551115123125782702e-17; | ||
| 283 | +static const float half = 0.5; | ||
| 284 | + | ||
| 285 | +/* The method is based on the descriptions in: | ||
| 286 | + | ||
| 287 | + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5; | ||
| 288 | + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9 | ||
| 289 | + | ||
| 290 | + We find the actual square root and half of its reciprocal | ||
| 291 | + simultaneously. */ | ||
| 292 | + | ||
| 293 | +#ifdef __STDC__ | ||
| 294 | +double | ||
| 295 | +__ieee754_sqrt (double b) | ||
| 296 | +#else | ||
| 297 | +double | ||
| 298 | +__ieee754_sqrt (b) | ||
| 299 | + double b; | ||
| 300 | +#endif | ||
| 301 | +{ | ||
| 302 | + if (__builtin_expect (b > 0, 1)) | ||
| 303 | + { | ||
| 304 | + double y, g, h, d, r; | ||
| 305 | + ieee_double_shape_type u; | ||
| 306 | + | ||
| 307 | + if (__builtin_expect (b != a_inf.value, 1)) | ||
| 308 | + { | ||
| 309 | + fenv_t fe; | ||
| 310 | + | ||
| 311 | + fe = fegetenv_register (); | ||
| 312 | + | ||
| 313 | + u.value = b; | ||
| 314 | + | ||
| 315 | + relax_fenv_state (); | ||
| 316 | + | ||
| 317 | + __asm__ ("frsqrte %[estimate], %[x]\n" | ||
| 318 | + : [estimate] "=f" (y) : [x] "f" (b)); | ||
| 319 | + | ||
| 320 | + /* Following Muller et al, page 168, equation 5.20. | ||
| 321 | + | ||
| 322 | + h goes to 1/(2*sqrt(b)) | ||
| 323 | + g goes to sqrt(b). | ||
| 324 | + | ||
| 325 | + We need three iterations to get within 1ulp. */ | ||
| 326 | + | ||
| 327 | + /* Indicate that these can be performed prior to the branch. GCC | ||
| 328 | + insists on sinking them below the branch, however; it seems like | ||
| 329 | + they'd be better before the branch so that we can cover any latency | ||
| 330 | + from storing the argument and loading its high word. Oh well. */ | ||
| 331 | + | ||
| 332 | + g = b * y; | ||
| 333 | + h = 0.5 * y; | ||
| 334 | + | ||
| 335 | + /* Handle small numbers by scaling. */ | ||
| 336 | + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0)) | ||
| 337 | + return __ieee754_sqrt (b * two108) * twom54; | ||
| 338 | + | ||
| 339 | +#define FMADD(a_, c_, b_) \ | ||
| 340 | + ({ double __r; \ | ||
| 341 | + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \ | ||
| 342 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 343 | + __r;}) | ||
| 344 | +#define FNMSUB(a_, c_, b_) \ | ||
| 345 | + ({ double __r; \ | ||
| 346 | + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 347 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 348 | + __r;}) | ||
| 349 | + | ||
| 350 | + r = FNMSUB (g, h, half); | ||
| 351 | + g = FMADD (g, r, g); | ||
| 352 | + h = FMADD (h, r, h); | ||
| 353 | + | ||
| 354 | + r = FNMSUB (g, h, half); | ||
| 355 | + g = FMADD (g, r, g); | ||
| 356 | + h = FMADD (h, r, h); | ||
| 357 | + | ||
| 358 | + r = FNMSUB (g, h, half); | ||
| 359 | + g = FMADD (g, r, g); | ||
| 360 | + h = FMADD (h, r, h); | ||
| 361 | + | ||
| 362 | + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */ | ||
| 363 | + | ||
| 364 | + /* Final refinement. */ | ||
| 365 | + d = FNMSUB (g, g, b); | ||
| 366 | + | ||
| 367 | + fesetenv_register (fe); | ||
| 368 | + return FMADD (d, h, g); | ||
| 369 | + } | ||
| 370 | + } | ||
| 371 | + else if (b < 0) | ||
| 372 | + { | ||
| 373 | + /* For some reason, some PowerPC32 processors don't implement | ||
| 374 | + FE_INVALID_SQRT. */ | ||
| 375 | +#ifdef FE_INVALID_SQRT | ||
| 376 | + feraiseexcept (FE_INVALID_SQRT); | ||
| 377 | + | ||
| 378 | + fenv_union_t u = { .fenv = fegetenv_register () }; | ||
| 379 | + if ((u.l[1] & FE_INVALID) == 0) | ||
| 380 | +#endif | ||
| 381 | + feraiseexcept (FE_INVALID); | ||
| 382 | + b = a_nan.value; | ||
| 383 | + } | ||
| 384 | + return f_wash (b); | ||
| 385 | +} | ||
| 386 | diff -ruN libc-orig/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c libc-e5500/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c | ||
| 387 | --- libc-orig/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c 1969-12-31 18:00:00.000000000 -0600 | ||
| 388 | +++ libc-e5500/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c 2011-05-18 11:15:22.467987000 -0500 | ||
| 389 | @@ -0,0 +1,101 @@ | ||
| 390 | +/* Single-precision floating point square root. | ||
| 391 | + Copyright (C) 2010 Free Software Foundation, Inc. | ||
| 392 | + This file is part of the GNU C Library. | ||
| 393 | + | ||
| 394 | + The GNU C Library is free software; you can redistribute it and/or | ||
| 395 | + modify it under the terms of the GNU Lesser General Public | ||
| 396 | + License as published by the Free Software Foundation; either | ||
| 397 | + version 2.1 of the License, or (at your option) any later version. | ||
| 398 | + | ||
| 399 | + The GNU C Library is distributed in the hope that it will be useful, | ||
| 400 | + but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 401 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
| 402 | + Lesser General Public License for more details. | ||
| 403 | + | ||
| 404 | + You should have received a copy of the GNU Lesser General Public | ||
| 405 | + License along with the GNU C Library; if not, write to the Free | ||
| 406 | + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA | ||
| 407 | + 02111-1307 USA. */ | ||
| 408 | + | ||
| 409 | +#include <math.h> | ||
| 410 | +#include <math_private.h> | ||
| 411 | +#include <fenv_libc.h> | ||
| 412 | +#include <inttypes.h> | ||
| 413 | + | ||
| 414 | +#include <sysdep.h> | ||
| 415 | +#include <ldsodefs.h> | ||
| 416 | + | ||
| 417 | +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 }; | ||
| 418 | +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 }; | ||
| 419 | +static const float threehalf = 1.5; | ||
| 420 | + | ||
| 421 | +/* The method is based on the descriptions in: | ||
| 422 | + | ||
| 423 | + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5; | ||
| 424 | + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9 | ||
| 425 | + | ||
| 426 | + We find the reciprocal square root and use that to compute the actual | ||
| 427 | + square root. */ | ||
| 428 | + | ||
| 429 | +#ifdef __STDC__ | ||
| 430 | +float | ||
| 431 | +__ieee754_sqrtf (float b) | ||
| 432 | +#else | ||
| 433 | +float | ||
| 434 | +__ieee754_sqrtf (b) | ||
| 435 | + float b; | ||
| 436 | +#endif | ||
| 437 | +{ | ||
| 438 | + if (__builtin_expect (b > 0, 1)) | ||
| 439 | + { | ||
| 440 | +#define FMSUB(a_, c_, b_) \ | ||
| 441 | + ({ double __r; \ | ||
| 442 | + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 443 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 444 | + __r;}) | ||
| 445 | +#define FNMSUB(a_, c_, b_) \ | ||
| 446 | + ({ double __r; \ | ||
| 447 | + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \ | ||
| 448 | + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \ | ||
| 449 | + __r;}) | ||
| 450 | + | ||
| 451 | + if (__builtin_expect (b != a_inf.value, 1)) | ||
| 452 | + { | ||
| 453 | + double y, x; | ||
| 454 | + fenv_t fe; | ||
| 455 | + | ||
| 456 | + fe = fegetenv_register (); | ||
| 457 | + | ||
| 458 | + relax_fenv_state (); | ||
| 459 | + | ||
| 460 | + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */ | ||
| 461 | + y = FMSUB (threehalf, b, b); | ||
| 462 | + | ||
| 463 | + /* Initial estimate. */ | ||
| 464 | + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b)); | ||
| 465 | + | ||
| 466 | + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */ | ||
| 467 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 468 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 469 | + x = x * FNMSUB (y, x * x, threehalf); | ||
| 470 | + | ||
| 471 | + /* All done. */ | ||
| 472 | + fesetenv_register (fe); | ||
| 473 | + return x * b; | ||
| 474 | + } | ||
| 475 | + } | ||
| 476 | + else if (b < 0) | ||
| 477 | + { | ||
| 478 | + /* For some reason, some PowerPC32 processors don't implement | ||
| 479 | + FE_INVALID_SQRT. */ | ||
| 480 | +#ifdef FE_INVALID_SQRT | ||
| 481 | + feraiseexcept (FE_INVALID_SQRT); | ||
| 482 | + | ||
| 483 | + fenv_union_t u = { .fenv = fegetenv_register () }; | ||
| 484 | + if ((u.l[1] & FE_INVALID) == 0) | ||
| 485 | +#endif | ||
| 486 | + feraiseexcept (FE_INVALID); | ||
| 487 | + b = a_nan.value; | ||
| 488 | + } | ||
| 489 | + return f_washf (b); | ||
| 490 | +} | ||
| 491 | diff -ruN libc-orig/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies libc-e5500/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies | ||
| 492 | --- libc-orig/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies 1969-12-31 18:00:00.000000000 -0600 | ||
| 493 | +++ libc-e5500/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies 2011-05-18 13:01:33.987988922 -0500 | ||
| 494 | @@ -0,0 +1 @@ | ||
| 495 | +powerpc/powerpc32/603e/fpu | ||
| 496 | diff -ruN libc-orig/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies libc-e5500/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies | ||
| 497 | --- libc-orig/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies 1969-12-31 18:00:00.000000000 -0600 | ||
| 498 | +++ libc-e5500/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies 2011-05-18 11:15:22.412987000 -0500 | ||
| 499 | @@ -0,0 +1 @@ | ||
| 500 | +powerpc/powerpc64/e5500/fpu | ||
