0009-fsl-e500-e5500-e6500-603e-fsqrt-implementation.patch 55 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581
  1. From 74923ca4b1ae0ed5a2478e7d265b37534f6815d7 Mon Sep 17 00:00:00 2001
  2. From: Khem Raj <raj.khem@gmail.com>
  3. Date: Wed, 18 Mar 2015 00:01:50 +0000
  4. Subject: [PATCH] fsl e500/e5500/e6500/603e fsqrt implementation
  5. Upstream-Status: Pending
  6. Signed-off-by: Edmar Wienskoski <edmar@freescale.com>
  7. Signed-off-by: Khem Raj <raj.khem@gmail.com>
  8. ---
  9. sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c | 134 ++++++++++++++++++
  10. sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c | 101 +++++++++++++
  11. sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c | 134 ++++++++++++++++++
  12. .../powerpc/powerpc32/e500mc/fpu/e_sqrtf.c | 101 +++++++++++++
  13. sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c | 134 ++++++++++++++++++
  14. sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c | 101 +++++++++++++
  15. sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c | 134 ++++++++++++++++++
  16. sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c | 101 +++++++++++++
  17. sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c | 134 ++++++++++++++++++
  18. sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c | 101 +++++++++++++
  19. sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c | 134 ++++++++++++++++++
  20. sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c | 101 +++++++++++++
  21. .../linux/powerpc/powerpc32/603e/fpu/Implies | 1 +
  22. .../powerpc/powerpc32/e300c3/fpu/Implies | 2 +
  23. .../powerpc/powerpc32/e500mc/fpu/Implies | 1 +
  24. .../linux/powerpc/powerpc32/e5500/fpu/Implies | 1 +
  25. .../linux/powerpc/powerpc32/e6500/fpu/Implies | 1 +
  26. .../linux/powerpc/powerpc64/e5500/fpu/Implies | 1 +
  27. .../linux/powerpc/powerpc64/e6500/fpu/Implies | 1 +
  28. 19 files changed, 1418 insertions(+)
  29. create mode 100644 sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
  30. create mode 100644 sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
  31. create mode 100644 sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
  32. create mode 100644 sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
  33. create mode 100644 sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
  34. create mode 100644 sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
  35. create mode 100644 sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
  36. create mode 100644 sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
  37. create mode 100644 sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
  38. create mode 100644 sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
  39. create mode 100644 sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
  40. create mode 100644 sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
  41. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
  42. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
  43. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
  44. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
  45. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
  46. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
  47. create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
  48. diff --git a/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
  49. new file mode 100644
  50. index 0000000000..71e516d1c8
  51. --- /dev/null
  52. +++ b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
  53. @@ -0,0 +1,134 @@
  54. +/* Double-precision floating point square root.
  55. + Copyright (C) 2010 Free Software Foundation, Inc.
  56. + This file is part of the GNU C Library.
  57. +
  58. + The GNU C Library is free software; you can redistribute it and/or
  59. + modify it under the terms of the GNU Lesser General Public
  60. + License as published by the Free Software Foundation; either
  61. + version 2.1 of the License, or (at your option) any later version.
  62. +
  63. + The GNU C Library is distributed in the hope that it will be useful,
  64. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  65. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  66. + Lesser General Public License for more details.
  67. +
  68. + You should have received a copy of the GNU Lesser General Public
  69. + License along with the GNU C Library; if not, write to the Free
  70. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  71. + 02111-1307 USA. */
  72. +
  73. +#include <math.h>
  74. +#include <math_private.h>
  75. +#include <fenv_libc.h>
  76. +#include <inttypes.h>
  77. +
  78. +#include <sysdep.h>
  79. +#include <ldsodefs.h>
  80. +
  81. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  82. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  83. +static const float two108 = 3.245185536584267269e+32;
  84. +static const float twom54 = 5.551115123125782702e-17;
  85. +static const float half = 0.5;
  86. +
  87. +/* The method is based on the descriptions in:
  88. +
  89. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  90. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  91. +
  92. + We find the actual square root and half of its reciprocal
  93. + simultaneously. */
  94. +
  95. +#ifdef __STDC__
  96. +double
  97. +__ieee754_sqrt (double b)
  98. +#else
  99. +double
  100. +__ieee754_sqrt (b)
  101. + double b;
  102. +#endif
  103. +{
  104. + if (__builtin_expect (b > 0, 1))
  105. + {
  106. + double y, g, h, d, r;
  107. + ieee_double_shape_type u;
  108. +
  109. + if (__builtin_expect (b != a_inf.value, 1))
  110. + {
  111. + fenv_t fe;
  112. +
  113. + fe = fegetenv_register ();
  114. +
  115. + u.value = b;
  116. +
  117. + relax_fenv_state ();
  118. +
  119. + __asm__ ("frsqrte %[estimate], %[x]\n"
  120. + : [estimate] "=f" (y) : [x] "f" (b));
  121. +
  122. + /* Following Muller et al, page 168, equation 5.20.
  123. +
  124. + h goes to 1/(2*sqrt(b))
  125. + g goes to sqrt(b).
  126. +
  127. + We need three iterations to get within 1ulp. */
  128. +
  129. + /* Indicate that these can be performed prior to the branch. GCC
  130. + insists on sinking them below the branch, however; it seems like
  131. + they'd be better before the branch so that we can cover any latency
  132. + from storing the argument and loading its high word. Oh well. */
  133. +
  134. + g = b * y;
  135. + h = 0.5 * y;
  136. +
  137. + /* Handle small numbers by scaling. */
  138. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  139. + return __ieee754_sqrt (b * two108) * twom54;
  140. +
  141. +#define FMADD(a_, c_, b_) \
  142. + ({ double __r; \
  143. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  144. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  145. + __r;})
  146. +#define FNMSUB(a_, c_, b_) \
  147. + ({ double __r; \
  148. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  149. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  150. + __r;})
  151. +
  152. + r = FNMSUB (g, h, half);
  153. + g = FMADD (g, r, g);
  154. + h = FMADD (h, r, h);
  155. +
  156. + r = FNMSUB (g, h, half);
  157. + g = FMADD (g, r, g);
  158. + h = FMADD (h, r, h);
  159. +
  160. + r = FNMSUB (g, h, half);
  161. + g = FMADD (g, r, g);
  162. + h = FMADD (h, r, h);
  163. +
  164. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  165. +
  166. + /* Final refinement. */
  167. + d = FNMSUB (g, g, b);
  168. +
  169. + fesetenv_register (fe);
  170. + return FMADD (d, h, g);
  171. + }
  172. + }
  173. + else if (b < 0)
  174. + {
  175. + /* For some reason, some PowerPC32 processors don't implement
  176. + FE_INVALID_SQRT. */
  177. +#ifdef FE_INVALID_SQRT
  178. + feraiseexcept (FE_INVALID_SQRT);
  179. +
  180. + fenv_union_t u = { .fenv = fegetenv_register () };
  181. + if ((u.l & FE_INVALID) == 0)
  182. +#endif
  183. + feraiseexcept (FE_INVALID);
  184. + b = a_nan.value;
  185. + }
  186. + return f_wash (b);
  187. +}
  188. diff --git a/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
  189. new file mode 100644
  190. index 0000000000..26fa067abf
  191. --- /dev/null
  192. +++ b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
  193. @@ -0,0 +1,101 @@
  194. +/* Single-precision floating point square root.
  195. + Copyright (C) 2010 Free Software Foundation, Inc.
  196. + This file is part of the GNU C Library.
  197. +
  198. + The GNU C Library is free software; you can redistribute it and/or
  199. + modify it under the terms of the GNU Lesser General Public
  200. + License as published by the Free Software Foundation; either
  201. + version 2.1 of the License, or (at your option) any later version.
  202. +
  203. + The GNU C Library is distributed in the hope that it will be useful,
  204. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  205. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  206. + Lesser General Public License for more details.
  207. +
  208. + You should have received a copy of the GNU Lesser General Public
  209. + License along with the GNU C Library; if not, write to the Free
  210. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  211. + 02111-1307 USA. */
  212. +
  213. +#include <math.h>
  214. +#include <math_private.h>
  215. +#include <fenv_libc.h>
  216. +#include <inttypes.h>
  217. +
  218. +#include <sysdep.h>
  219. +#include <ldsodefs.h>
  220. +
  221. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  222. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  223. +static const float threehalf = 1.5;
  224. +
  225. +/* The method is based on the descriptions in:
  226. +
  227. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  228. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  229. +
  230. + We find the reciprocal square root and use that to compute the actual
  231. + square root. */
  232. +
  233. +#ifdef __STDC__
  234. +float
  235. +__ieee754_sqrtf (float b)
  236. +#else
  237. +float
  238. +__ieee754_sqrtf (b)
  239. + float b;
  240. +#endif
  241. +{
  242. + if (__builtin_expect (b > 0, 1))
  243. + {
  244. +#define FMSUB(a_, c_, b_) \
  245. + ({ double __r; \
  246. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  247. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  248. + __r;})
  249. +#define FNMSUB(a_, c_, b_) \
  250. + ({ double __r; \
  251. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  252. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  253. + __r;})
  254. +
  255. + if (__builtin_expect (b != a_inf.value, 1))
  256. + {
  257. + double y, x;
  258. + fenv_t fe;
  259. +
  260. + fe = fegetenv_register ();
  261. +
  262. + relax_fenv_state ();
  263. +
  264. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  265. + y = FMSUB (threehalf, b, b);
  266. +
  267. + /* Initial estimate. */
  268. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  269. +
  270. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  271. + x = x * FNMSUB (y, x * x, threehalf);
  272. + x = x * FNMSUB (y, x * x, threehalf);
  273. + x = x * FNMSUB (y, x * x, threehalf);
  274. +
  275. + /* All done. */
  276. + fesetenv_register (fe);
  277. + return x * b;
  278. + }
  279. + }
  280. + else if (b < 0)
  281. + {
  282. + /* For some reason, some PowerPC32 processors don't implement
  283. + FE_INVALID_SQRT. */
  284. +#ifdef FE_INVALID_SQRT
  285. + feraiseexcept (FE_INVALID_SQRT);
  286. +
  287. + fenv_union_t u = { .fenv = fegetenv_register () };
  288. + if ((u.l & FE_INVALID) == 0)
  289. +#endif
  290. + feraiseexcept (FE_INVALID);
  291. + b = a_nan.value;
  292. + }
  293. + return f_washf (b);
  294. +}
  295. diff --git a/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
  296. new file mode 100644
  297. index 0000000000..71e516d1c8
  298. --- /dev/null
  299. +++ b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
  300. @@ -0,0 +1,134 @@
  301. +/* Double-precision floating point square root.
  302. + Copyright (C) 2010 Free Software Foundation, Inc.
  303. + This file is part of the GNU C Library.
  304. +
  305. + The GNU C Library is free software; you can redistribute it and/or
  306. + modify it under the terms of the GNU Lesser General Public
  307. + License as published by the Free Software Foundation; either
  308. + version 2.1 of the License, or (at your option) any later version.
  309. +
  310. + The GNU C Library is distributed in the hope that it will be useful,
  311. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  312. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  313. + Lesser General Public License for more details.
  314. +
  315. + You should have received a copy of the GNU Lesser General Public
  316. + License along with the GNU C Library; if not, write to the Free
  317. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  318. + 02111-1307 USA. */
  319. +
  320. +#include <math.h>
  321. +#include <math_private.h>
  322. +#include <fenv_libc.h>
  323. +#include <inttypes.h>
  324. +
  325. +#include <sysdep.h>
  326. +#include <ldsodefs.h>
  327. +
  328. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  329. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  330. +static const float two108 = 3.245185536584267269e+32;
  331. +static const float twom54 = 5.551115123125782702e-17;
  332. +static const float half = 0.5;
  333. +
  334. +/* The method is based on the descriptions in:
  335. +
  336. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  337. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  338. +
  339. + We find the actual square root and half of its reciprocal
  340. + simultaneously. */
  341. +
  342. +#ifdef __STDC__
  343. +double
  344. +__ieee754_sqrt (double b)
  345. +#else
  346. +double
  347. +__ieee754_sqrt (b)
  348. + double b;
  349. +#endif
  350. +{
  351. + if (__builtin_expect (b > 0, 1))
  352. + {
  353. + double y, g, h, d, r;
  354. + ieee_double_shape_type u;
  355. +
  356. + if (__builtin_expect (b != a_inf.value, 1))
  357. + {
  358. + fenv_t fe;
  359. +
  360. + fe = fegetenv_register ();
  361. +
  362. + u.value = b;
  363. +
  364. + relax_fenv_state ();
  365. +
  366. + __asm__ ("frsqrte %[estimate], %[x]\n"
  367. + : [estimate] "=f" (y) : [x] "f" (b));
  368. +
  369. + /* Following Muller et al, page 168, equation 5.20.
  370. +
  371. + h goes to 1/(2*sqrt(b))
  372. + g goes to sqrt(b).
  373. +
  374. + We need three iterations to get within 1ulp. */
  375. +
  376. + /* Indicate that these can be performed prior to the branch. GCC
  377. + insists on sinking them below the branch, however; it seems like
  378. + they'd be better before the branch so that we can cover any latency
  379. + from storing the argument and loading its high word. Oh well. */
  380. +
  381. + g = b * y;
  382. + h = 0.5 * y;
  383. +
  384. + /* Handle small numbers by scaling. */
  385. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  386. + return __ieee754_sqrt (b * two108) * twom54;
  387. +
  388. +#define FMADD(a_, c_, b_) \
  389. + ({ double __r; \
  390. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  391. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  392. + __r;})
  393. +#define FNMSUB(a_, c_, b_) \
  394. + ({ double __r; \
  395. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  396. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  397. + __r;})
  398. +
  399. + r = FNMSUB (g, h, half);
  400. + g = FMADD (g, r, g);
  401. + h = FMADD (h, r, h);
  402. +
  403. + r = FNMSUB (g, h, half);
  404. + g = FMADD (g, r, g);
  405. + h = FMADD (h, r, h);
  406. +
  407. + r = FNMSUB (g, h, half);
  408. + g = FMADD (g, r, g);
  409. + h = FMADD (h, r, h);
  410. +
  411. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  412. +
  413. + /* Final refinement. */
  414. + d = FNMSUB (g, g, b);
  415. +
  416. + fesetenv_register (fe);
  417. + return FMADD (d, h, g);
  418. + }
  419. + }
  420. + else if (b < 0)
  421. + {
  422. + /* For some reason, some PowerPC32 processors don't implement
  423. + FE_INVALID_SQRT. */
  424. +#ifdef FE_INVALID_SQRT
  425. + feraiseexcept (FE_INVALID_SQRT);
  426. +
  427. + fenv_union_t u = { .fenv = fegetenv_register () };
  428. + if ((u.l & FE_INVALID) == 0)
  429. +#endif
  430. + feraiseexcept (FE_INVALID);
  431. + b = a_nan.value;
  432. + }
  433. + return f_wash (b);
  434. +}
  435. diff --git a/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
  436. new file mode 100644
  437. index 0000000000..26fa067abf
  438. --- /dev/null
  439. +++ b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
  440. @@ -0,0 +1,101 @@
  441. +/* Single-precision floating point square root.
  442. + Copyright (C) 2010 Free Software Foundation, Inc.
  443. + This file is part of the GNU C Library.
  444. +
  445. + The GNU C Library is free software; you can redistribute it and/or
  446. + modify it under the terms of the GNU Lesser General Public
  447. + License as published by the Free Software Foundation; either
  448. + version 2.1 of the License, or (at your option) any later version.
  449. +
  450. + The GNU C Library is distributed in the hope that it will be useful,
  451. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  452. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  453. + Lesser General Public License for more details.
  454. +
  455. + You should have received a copy of the GNU Lesser General Public
  456. + License along with the GNU C Library; if not, write to the Free
  457. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  458. + 02111-1307 USA. */
  459. +
  460. +#include <math.h>
  461. +#include <math_private.h>
  462. +#include <fenv_libc.h>
  463. +#include <inttypes.h>
  464. +
  465. +#include <sysdep.h>
  466. +#include <ldsodefs.h>
  467. +
  468. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  469. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  470. +static const float threehalf = 1.5;
  471. +
  472. +/* The method is based on the descriptions in:
  473. +
  474. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  475. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  476. +
  477. + We find the reciprocal square root and use that to compute the actual
  478. + square root. */
  479. +
  480. +#ifdef __STDC__
  481. +float
  482. +__ieee754_sqrtf (float b)
  483. +#else
  484. +float
  485. +__ieee754_sqrtf (b)
  486. + float b;
  487. +#endif
  488. +{
  489. + if (__builtin_expect (b > 0, 1))
  490. + {
  491. +#define FMSUB(a_, c_, b_) \
  492. + ({ double __r; \
  493. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  494. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  495. + __r;})
  496. +#define FNMSUB(a_, c_, b_) \
  497. + ({ double __r; \
  498. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  499. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  500. + __r;})
  501. +
  502. + if (__builtin_expect (b != a_inf.value, 1))
  503. + {
  504. + double y, x;
  505. + fenv_t fe;
  506. +
  507. + fe = fegetenv_register ();
  508. +
  509. + relax_fenv_state ();
  510. +
  511. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  512. + y = FMSUB (threehalf, b, b);
  513. +
  514. + /* Initial estimate. */
  515. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  516. +
  517. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  518. + x = x * FNMSUB (y, x * x, threehalf);
  519. + x = x * FNMSUB (y, x * x, threehalf);
  520. + x = x * FNMSUB (y, x * x, threehalf);
  521. +
  522. + /* All done. */
  523. + fesetenv_register (fe);
  524. + return x * b;
  525. + }
  526. + }
  527. + else if (b < 0)
  528. + {
  529. + /* For some reason, some PowerPC32 processors don't implement
  530. + FE_INVALID_SQRT. */
  531. +#ifdef FE_INVALID_SQRT
  532. + feraiseexcept (FE_INVALID_SQRT);
  533. +
  534. + fenv_union_t u = { .fenv = fegetenv_register () };
  535. + if ((u.l & FE_INVALID) == 0)
  536. +#endif
  537. + feraiseexcept (FE_INVALID);
  538. + b = a_nan.value;
  539. + }
  540. + return f_washf (b);
  541. +}
  542. diff --git a/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
  543. new file mode 100644
  544. index 0000000000..71e516d1c8
  545. --- /dev/null
  546. +++ b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
  547. @@ -0,0 +1,134 @@
  548. +/* Double-precision floating point square root.
  549. + Copyright (C) 2010 Free Software Foundation, Inc.
  550. + This file is part of the GNU C Library.
  551. +
  552. + The GNU C Library is free software; you can redistribute it and/or
  553. + modify it under the terms of the GNU Lesser General Public
  554. + License as published by the Free Software Foundation; either
  555. + version 2.1 of the License, or (at your option) any later version.
  556. +
  557. + The GNU C Library is distributed in the hope that it will be useful,
  558. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  559. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  560. + Lesser General Public License for more details.
  561. +
  562. + You should have received a copy of the GNU Lesser General Public
  563. + License along with the GNU C Library; if not, write to the Free
  564. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  565. + 02111-1307 USA. */
  566. +
  567. +#include <math.h>
  568. +#include <math_private.h>
  569. +#include <fenv_libc.h>
  570. +#include <inttypes.h>
  571. +
  572. +#include <sysdep.h>
  573. +#include <ldsodefs.h>
  574. +
  575. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  576. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  577. +static const float two108 = 3.245185536584267269e+32;
  578. +static const float twom54 = 5.551115123125782702e-17;
  579. +static const float half = 0.5;
  580. +
  581. +/* The method is based on the descriptions in:
  582. +
  583. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  584. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  585. +
  586. + We find the actual square root and half of its reciprocal
  587. + simultaneously. */
  588. +
  589. +#ifdef __STDC__
  590. +double
  591. +__ieee754_sqrt (double b)
  592. +#else
  593. +double
  594. +__ieee754_sqrt (b)
  595. + double b;
  596. +#endif
  597. +{
  598. + if (__builtin_expect (b > 0, 1))
  599. + {
  600. + double y, g, h, d, r;
  601. + ieee_double_shape_type u;
  602. +
  603. + if (__builtin_expect (b != a_inf.value, 1))
  604. + {
  605. + fenv_t fe;
  606. +
  607. + fe = fegetenv_register ();
  608. +
  609. + u.value = b;
  610. +
  611. + relax_fenv_state ();
  612. +
  613. + __asm__ ("frsqrte %[estimate], %[x]\n"
  614. + : [estimate] "=f" (y) : [x] "f" (b));
  615. +
  616. + /* Following Muller et al, page 168, equation 5.20.
  617. +
  618. + h goes to 1/(2*sqrt(b))
  619. + g goes to sqrt(b).
  620. +
  621. + We need three iterations to get within 1ulp. */
  622. +
  623. + /* Indicate that these can be performed prior to the branch. GCC
  624. + insists on sinking them below the branch, however; it seems like
  625. + they'd be better before the branch so that we can cover any latency
  626. + from storing the argument and loading its high word. Oh well. */
  627. +
  628. + g = b * y;
  629. + h = 0.5 * y;
  630. +
  631. + /* Handle small numbers by scaling. */
  632. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  633. + return __ieee754_sqrt (b * two108) * twom54;
  634. +
  635. +#define FMADD(a_, c_, b_) \
  636. + ({ double __r; \
  637. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  638. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  639. + __r;})
  640. +#define FNMSUB(a_, c_, b_) \
  641. + ({ double __r; \
  642. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  643. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  644. + __r;})
  645. +
  646. + r = FNMSUB (g, h, half);
  647. + g = FMADD (g, r, g);
  648. + h = FMADD (h, r, h);
  649. +
  650. + r = FNMSUB (g, h, half);
  651. + g = FMADD (g, r, g);
  652. + h = FMADD (h, r, h);
  653. +
  654. + r = FNMSUB (g, h, half);
  655. + g = FMADD (g, r, g);
  656. + h = FMADD (h, r, h);
  657. +
  658. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  659. +
  660. + /* Final refinement. */
  661. + d = FNMSUB (g, g, b);
  662. +
  663. + fesetenv_register (fe);
  664. + return FMADD (d, h, g);
  665. + }
  666. + }
  667. + else if (b < 0)
  668. + {
  669. + /* For some reason, some PowerPC32 processors don't implement
  670. + FE_INVALID_SQRT. */
  671. +#ifdef FE_INVALID_SQRT
  672. + feraiseexcept (FE_INVALID_SQRT);
  673. +
  674. + fenv_union_t u = { .fenv = fegetenv_register () };
  675. + if ((u.l & FE_INVALID) == 0)
  676. +#endif
  677. + feraiseexcept (FE_INVALID);
  678. + b = a_nan.value;
  679. + }
  680. + return f_wash (b);
  681. +}
  682. diff --git a/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
  683. new file mode 100644
  684. index 0000000000..26fa067abf
  685. --- /dev/null
  686. +++ b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
  687. @@ -0,0 +1,101 @@
  688. +/* Single-precision floating point square root.
  689. + Copyright (C) 2010 Free Software Foundation, Inc.
  690. + This file is part of the GNU C Library.
  691. +
  692. + The GNU C Library is free software; you can redistribute it and/or
  693. + modify it under the terms of the GNU Lesser General Public
  694. + License as published by the Free Software Foundation; either
  695. + version 2.1 of the License, or (at your option) any later version.
  696. +
  697. + The GNU C Library is distributed in the hope that it will be useful,
  698. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  699. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  700. + Lesser General Public License for more details.
  701. +
  702. + You should have received a copy of the GNU Lesser General Public
  703. + License along with the GNU C Library; if not, write to the Free
  704. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  705. + 02111-1307 USA. */
  706. +
  707. +#include <math.h>
  708. +#include <math_private.h>
  709. +#include <fenv_libc.h>
  710. +#include <inttypes.h>
  711. +
  712. +#include <sysdep.h>
  713. +#include <ldsodefs.h>
  714. +
  715. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  716. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  717. +static const float threehalf = 1.5;
  718. +
  719. +/* The method is based on the descriptions in:
  720. +
  721. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  722. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  723. +
  724. + We find the reciprocal square root and use that to compute the actual
  725. + square root. */
  726. +
  727. +#ifdef __STDC__
  728. +float
  729. +__ieee754_sqrtf (float b)
  730. +#else
  731. +float
  732. +__ieee754_sqrtf (b)
  733. + float b;
  734. +#endif
  735. +{
  736. + if (__builtin_expect (b > 0, 1))
  737. + {
  738. +#define FMSUB(a_, c_, b_) \
  739. + ({ double __r; \
  740. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  741. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  742. + __r;})
  743. +#define FNMSUB(a_, c_, b_) \
  744. + ({ double __r; \
  745. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  746. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  747. + __r;})
  748. +
  749. + if (__builtin_expect (b != a_inf.value, 1))
  750. + {
  751. + double y, x;
  752. + fenv_t fe;
  753. +
  754. + fe = fegetenv_register ();
  755. +
  756. + relax_fenv_state ();
  757. +
  758. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  759. + y = FMSUB (threehalf, b, b);
  760. +
  761. + /* Initial estimate. */
  762. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  763. +
  764. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  765. + x = x * FNMSUB (y, x * x, threehalf);
  766. + x = x * FNMSUB (y, x * x, threehalf);
  767. + x = x * FNMSUB (y, x * x, threehalf);
  768. +
  769. + /* All done. */
  770. + fesetenv_register (fe);
  771. + return x * b;
  772. + }
  773. + }
  774. + else if (b < 0)
  775. + {
  776. + /* For some reason, some PowerPC32 processors don't implement
  777. + FE_INVALID_SQRT. */
  778. +#ifdef FE_INVALID_SQRT
  779. + feraiseexcept (FE_INVALID_SQRT);
  780. +
  781. + fenv_union_t u = { .fenv = fegetenv_register () };
  782. + if ((u.l & FE_INVALID) == 0)
  783. +#endif
  784. + feraiseexcept (FE_INVALID);
  785. + b = a_nan.value;
  786. + }
  787. + return f_washf (b);
  788. +}
  789. diff --git a/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
  790. new file mode 100644
  791. index 0000000000..71e516d1c8
  792. --- /dev/null
  793. +++ b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
  794. @@ -0,0 +1,134 @@
  795. +/* Double-precision floating point square root.
  796. + Copyright (C) 2010 Free Software Foundation, Inc.
  797. + This file is part of the GNU C Library.
  798. +
  799. + The GNU C Library is free software; you can redistribute it and/or
  800. + modify it under the terms of the GNU Lesser General Public
  801. + License as published by the Free Software Foundation; either
  802. + version 2.1 of the License, or (at your option) any later version.
  803. +
  804. + The GNU C Library is distributed in the hope that it will be useful,
  805. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  806. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  807. + Lesser General Public License for more details.
  808. +
  809. + You should have received a copy of the GNU Lesser General Public
  810. + License along with the GNU C Library; if not, write to the Free
  811. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  812. + 02111-1307 USA. */
  813. +
  814. +#include <math.h>
  815. +#include <math_private.h>
  816. +#include <fenv_libc.h>
  817. +#include <inttypes.h>
  818. +
  819. +#include <sysdep.h>
  820. +#include <ldsodefs.h>
  821. +
  822. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  823. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  824. +static const float two108 = 3.245185536584267269e+32;
  825. +static const float twom54 = 5.551115123125782702e-17;
  826. +static const float half = 0.5;
  827. +
  828. +/* The method is based on the descriptions in:
  829. +
  830. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  831. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  832. +
  833. + We find the actual square root and half of its reciprocal
  834. + simultaneously. */
  835. +
  836. +#ifdef __STDC__
  837. +double
  838. +__ieee754_sqrt (double b)
  839. +#else
  840. +double
  841. +__ieee754_sqrt (b)
  842. + double b;
  843. +#endif
  844. +{
  845. + if (__builtin_expect (b > 0, 1))
  846. + {
  847. + double y, g, h, d, r;
  848. + ieee_double_shape_type u;
  849. +
  850. + if (__builtin_expect (b != a_inf.value, 1))
  851. + {
  852. + fenv_t fe;
  853. +
  854. + fe = fegetenv_register ();
  855. +
  856. + u.value = b;
  857. +
  858. + relax_fenv_state ();
  859. +
  860. + __asm__ ("frsqrte %[estimate], %[x]\n"
  861. + : [estimate] "=f" (y) : [x] "f" (b));
  862. +
  863. + /* Following Muller et al, page 168, equation 5.20.
  864. +
  865. + h goes to 1/(2*sqrt(b))
  866. + g goes to sqrt(b).
  867. +
  868. + We need three iterations to get within 1ulp. */
  869. +
  870. + /* Indicate that these can be performed prior to the branch. GCC
  871. + insists on sinking them below the branch, however; it seems like
  872. + they'd be better before the branch so that we can cover any latency
  873. + from storing the argument and loading its high word. Oh well. */
  874. +
  875. + g = b * y;
  876. + h = 0.5 * y;
  877. +
  878. + /* Handle small numbers by scaling. */
  879. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  880. + return __ieee754_sqrt (b * two108) * twom54;
  881. +
  882. +#define FMADD(a_, c_, b_) \
  883. + ({ double __r; \
  884. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  885. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  886. + __r;})
  887. +#define FNMSUB(a_, c_, b_) \
  888. + ({ double __r; \
  889. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  890. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  891. + __r;})
  892. +
  893. + r = FNMSUB (g, h, half);
  894. + g = FMADD (g, r, g);
  895. + h = FMADD (h, r, h);
  896. +
  897. + r = FNMSUB (g, h, half);
  898. + g = FMADD (g, r, g);
  899. + h = FMADD (h, r, h);
  900. +
  901. + r = FNMSUB (g, h, half);
  902. + g = FMADD (g, r, g);
  903. + h = FMADD (h, r, h);
  904. +
  905. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  906. +
  907. + /* Final refinement. */
  908. + d = FNMSUB (g, g, b);
  909. +
  910. + fesetenv_register (fe);
  911. + return FMADD (d, h, g);
  912. + }
  913. + }
  914. + else if (b < 0)
  915. + {
  916. + /* For some reason, some PowerPC32 processors don't implement
  917. + FE_INVALID_SQRT. */
  918. +#ifdef FE_INVALID_SQRT
  919. + feraiseexcept (FE_INVALID_SQRT);
  920. +
  921. + fenv_union_t u = { .fenv = fegetenv_register () };
  922. + if ((u.l & FE_INVALID) == 0)
  923. +#endif
  924. + feraiseexcept (FE_INVALID);
  925. + b = a_nan.value;
  926. + }
  927. + return f_wash (b);
  928. +}
  929. diff --git a/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
  930. new file mode 100644
  931. index 0000000000..26fa067abf
  932. --- /dev/null
  933. +++ b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
  934. @@ -0,0 +1,101 @@
  935. +/* Single-precision floating point square root.
  936. + Copyright (C) 2010 Free Software Foundation, Inc.
  937. + This file is part of the GNU C Library.
  938. +
  939. + The GNU C Library is free software; you can redistribute it and/or
  940. + modify it under the terms of the GNU Lesser General Public
  941. + License as published by the Free Software Foundation; either
  942. + version 2.1 of the License, or (at your option) any later version.
  943. +
  944. + The GNU C Library is distributed in the hope that it will be useful,
  945. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  946. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  947. + Lesser General Public License for more details.
  948. +
  949. + You should have received a copy of the GNU Lesser General Public
  950. + License along with the GNU C Library; if not, write to the Free
  951. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  952. + 02111-1307 USA. */
  953. +
  954. +#include <math.h>
  955. +#include <math_private.h>
  956. +#include <fenv_libc.h>
  957. +#include <inttypes.h>
  958. +
  959. +#include <sysdep.h>
  960. +#include <ldsodefs.h>
  961. +
  962. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  963. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  964. +static const float threehalf = 1.5;
  965. +
  966. +/* The method is based on the descriptions in:
  967. +
  968. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  969. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  970. +
  971. + We find the reciprocal square root and use that to compute the actual
  972. + square root. */
  973. +
  974. +#ifdef __STDC__
  975. +float
  976. +__ieee754_sqrtf (float b)
  977. +#else
  978. +float
  979. +__ieee754_sqrtf (b)
  980. + float b;
  981. +#endif
  982. +{
  983. + if (__builtin_expect (b > 0, 1))
  984. + {
  985. +#define FMSUB(a_, c_, b_) \
  986. + ({ double __r; \
  987. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  988. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  989. + __r;})
  990. +#define FNMSUB(a_, c_, b_) \
  991. + ({ double __r; \
  992. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  993. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  994. + __r;})
  995. +
  996. + if (__builtin_expect (b != a_inf.value, 1))
  997. + {
  998. + double y, x;
  999. + fenv_t fe;
  1000. +
  1001. + fe = fegetenv_register ();
  1002. +
  1003. + relax_fenv_state ();
  1004. +
  1005. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  1006. + y = FMSUB (threehalf, b, b);
  1007. +
  1008. + /* Initial estimate. */
  1009. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  1010. +
  1011. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  1012. + x = x * FNMSUB (y, x * x, threehalf);
  1013. + x = x * FNMSUB (y, x * x, threehalf);
  1014. + x = x * FNMSUB (y, x * x, threehalf);
  1015. +
  1016. + /* All done. */
  1017. + fesetenv_register (fe);
  1018. + return x * b;
  1019. + }
  1020. + }
  1021. + else if (b < 0)
  1022. + {
  1023. + /* For some reason, some PowerPC32 processors don't implement
  1024. + FE_INVALID_SQRT. */
  1025. +#ifdef FE_INVALID_SQRT
  1026. + feraiseexcept (FE_INVALID_SQRT);
  1027. +
  1028. + fenv_union_t u = { .fenv = fegetenv_register () };
  1029. + if ((u.l & FE_INVALID) == 0)
  1030. +#endif
  1031. + feraiseexcept (FE_INVALID);
  1032. + b = a_nan.value;
  1033. + }
  1034. + return f_washf (b);
  1035. +}
  1036. diff --git a/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
  1037. new file mode 100644
  1038. index 0000000000..71e516d1c8
  1039. --- /dev/null
  1040. +++ b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
  1041. @@ -0,0 +1,134 @@
  1042. +/* Double-precision floating point square root.
  1043. + Copyright (C) 2010 Free Software Foundation, Inc.
  1044. + This file is part of the GNU C Library.
  1045. +
  1046. + The GNU C Library is free software; you can redistribute it and/or
  1047. + modify it under the terms of the GNU Lesser General Public
  1048. + License as published by the Free Software Foundation; either
  1049. + version 2.1 of the License, or (at your option) any later version.
  1050. +
  1051. + The GNU C Library is distributed in the hope that it will be useful,
  1052. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  1053. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  1054. + Lesser General Public License for more details.
  1055. +
  1056. + You should have received a copy of the GNU Lesser General Public
  1057. + License along with the GNU C Library; if not, write to the Free
  1058. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  1059. + 02111-1307 USA. */
  1060. +
  1061. +#include <math.h>
  1062. +#include <math_private.h>
  1063. +#include <fenv_libc.h>
  1064. +#include <inttypes.h>
  1065. +
  1066. +#include <sysdep.h>
  1067. +#include <ldsodefs.h>
  1068. +
  1069. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  1070. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  1071. +static const float two108 = 3.245185536584267269e+32;
  1072. +static const float twom54 = 5.551115123125782702e-17;
  1073. +static const float half = 0.5;
  1074. +
  1075. +/* The method is based on the descriptions in:
  1076. +
  1077. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  1078. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  1079. +
  1080. + We find the actual square root and half of its reciprocal
  1081. + simultaneously. */
  1082. +
  1083. +#ifdef __STDC__
  1084. +double
  1085. +__ieee754_sqrt (double b)
  1086. +#else
  1087. +double
  1088. +__ieee754_sqrt (b)
  1089. + double b;
  1090. +#endif
  1091. +{
  1092. + if (__builtin_expect (b > 0, 1))
  1093. + {
  1094. + double y, g, h, d, r;
  1095. + ieee_double_shape_type u;
  1096. +
  1097. + if (__builtin_expect (b != a_inf.value, 1))
  1098. + {
  1099. + fenv_t fe;
  1100. +
  1101. + fe = fegetenv_register ();
  1102. +
  1103. + u.value = b;
  1104. +
  1105. + relax_fenv_state ();
  1106. +
  1107. + __asm__ ("frsqrte %[estimate], %[x]\n"
  1108. + : [estimate] "=f" (y) : [x] "f" (b));
  1109. +
  1110. + /* Following Muller et al, page 168, equation 5.20.
  1111. +
  1112. + h goes to 1/(2*sqrt(b))
  1113. + g goes to sqrt(b).
  1114. +
  1115. + We need three iterations to get within 1ulp. */
  1116. +
  1117. + /* Indicate that these can be performed prior to the branch. GCC
  1118. + insists on sinking them below the branch, however; it seems like
  1119. + they'd be better before the branch so that we can cover any latency
  1120. + from storing the argument and loading its high word. Oh well. */
  1121. +
  1122. + g = b * y;
  1123. + h = 0.5 * y;
  1124. +
  1125. + /* Handle small numbers by scaling. */
  1126. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  1127. + return __ieee754_sqrt (b * two108) * twom54;
  1128. +
  1129. +#define FMADD(a_, c_, b_) \
  1130. + ({ double __r; \
  1131. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  1132. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1133. + __r;})
  1134. +#define FNMSUB(a_, c_, b_) \
  1135. + ({ double __r; \
  1136. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  1137. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1138. + __r;})
  1139. +
  1140. + r = FNMSUB (g, h, half);
  1141. + g = FMADD (g, r, g);
  1142. + h = FMADD (h, r, h);
  1143. +
  1144. + r = FNMSUB (g, h, half);
  1145. + g = FMADD (g, r, g);
  1146. + h = FMADD (h, r, h);
  1147. +
  1148. + r = FNMSUB (g, h, half);
  1149. + g = FMADD (g, r, g);
  1150. + h = FMADD (h, r, h);
  1151. +
  1152. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  1153. +
  1154. + /* Final refinement. */
  1155. + d = FNMSUB (g, g, b);
  1156. +
  1157. + fesetenv_register (fe);
  1158. + return FMADD (d, h, g);
  1159. + }
  1160. + }
  1161. + else if (b < 0)
  1162. + {
  1163. + /* For some reason, some PowerPC32 processors don't implement
  1164. + FE_INVALID_SQRT. */
  1165. +#ifdef FE_INVALID_SQRT
  1166. + feraiseexcept (FE_INVALID_SQRT);
  1167. +
  1168. + fenv_union_t u = { .fenv = fegetenv_register () };
  1169. + if ((u.l & FE_INVALID) == 0)
  1170. +#endif
  1171. + feraiseexcept (FE_INVALID);
  1172. + b = a_nan.value;
  1173. + }
  1174. + return f_wash (b);
  1175. +}
  1176. diff --git a/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
  1177. new file mode 100644
  1178. index 0000000000..26fa067abf
  1179. --- /dev/null
  1180. +++ b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
  1181. @@ -0,0 +1,101 @@
  1182. +/* Single-precision floating point square root.
  1183. + Copyright (C) 2010 Free Software Foundation, Inc.
  1184. + This file is part of the GNU C Library.
  1185. +
  1186. + The GNU C Library is free software; you can redistribute it and/or
  1187. + modify it under the terms of the GNU Lesser General Public
  1188. + License as published by the Free Software Foundation; either
  1189. + version 2.1 of the License, or (at your option) any later version.
  1190. +
  1191. + The GNU C Library is distributed in the hope that it will be useful,
  1192. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  1193. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  1194. + Lesser General Public License for more details.
  1195. +
  1196. + You should have received a copy of the GNU Lesser General Public
  1197. + License along with the GNU C Library; if not, write to the Free
  1198. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  1199. + 02111-1307 USA. */
  1200. +
  1201. +#include <math.h>
  1202. +#include <math_private.h>
  1203. +#include <fenv_libc.h>
  1204. +#include <inttypes.h>
  1205. +
  1206. +#include <sysdep.h>
  1207. +#include <ldsodefs.h>
  1208. +
  1209. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  1210. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  1211. +static const float threehalf = 1.5;
  1212. +
  1213. +/* The method is based on the descriptions in:
  1214. +
  1215. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  1216. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  1217. +
  1218. + We find the reciprocal square root and use that to compute the actual
  1219. + square root. */
  1220. +
  1221. +#ifdef __STDC__
  1222. +float
  1223. +__ieee754_sqrtf (float b)
  1224. +#else
  1225. +float
  1226. +__ieee754_sqrtf (b)
  1227. + float b;
  1228. +#endif
  1229. +{
  1230. + if (__builtin_expect (b > 0, 1))
  1231. + {
  1232. +#define FMSUB(a_, c_, b_) \
  1233. + ({ double __r; \
  1234. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  1235. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1236. + __r;})
  1237. +#define FNMSUB(a_, c_, b_) \
  1238. + ({ double __r; \
  1239. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  1240. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1241. + __r;})
  1242. +
  1243. + if (__builtin_expect (b != a_inf.value, 1))
  1244. + {
  1245. + double y, x;
  1246. + fenv_t fe;
  1247. +
  1248. + fe = fegetenv_register ();
  1249. +
  1250. + relax_fenv_state ();
  1251. +
  1252. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  1253. + y = FMSUB (threehalf, b, b);
  1254. +
  1255. + /* Initial estimate. */
  1256. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  1257. +
  1258. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  1259. + x = x * FNMSUB (y, x * x, threehalf);
  1260. + x = x * FNMSUB (y, x * x, threehalf);
  1261. + x = x * FNMSUB (y, x * x, threehalf);
  1262. +
  1263. + /* All done. */
  1264. + fesetenv_register (fe);
  1265. + return x * b;
  1266. + }
  1267. + }
  1268. + else if (b < 0)
  1269. + {
  1270. + /* For some reason, some PowerPC32 processors don't implement
  1271. + FE_INVALID_SQRT. */
  1272. +#ifdef FE_INVALID_SQRT
  1273. + feraiseexcept (FE_INVALID_SQRT);
  1274. +
  1275. + fenv_union_t u = { .fenv = fegetenv_register () };
  1276. + if ((u.l & FE_INVALID) == 0)
  1277. +#endif
  1278. + feraiseexcept (FE_INVALID);
  1279. + b = a_nan.value;
  1280. + }
  1281. + return f_washf (b);
  1282. +}
  1283. diff --git a/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
  1284. new file mode 100644
  1285. index 0000000000..71e516d1c8
  1286. --- /dev/null
  1287. +++ b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
  1288. @@ -0,0 +1,134 @@
  1289. +/* Double-precision floating point square root.
  1290. + Copyright (C) 2010 Free Software Foundation, Inc.
  1291. + This file is part of the GNU C Library.
  1292. +
  1293. + The GNU C Library is free software; you can redistribute it and/or
  1294. + modify it under the terms of the GNU Lesser General Public
  1295. + License as published by the Free Software Foundation; either
  1296. + version 2.1 of the License, or (at your option) any later version.
  1297. +
  1298. + The GNU C Library is distributed in the hope that it will be useful,
  1299. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  1300. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  1301. + Lesser General Public License for more details.
  1302. +
  1303. + You should have received a copy of the GNU Lesser General Public
  1304. + License along with the GNU C Library; if not, write to the Free
  1305. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  1306. + 02111-1307 USA. */
  1307. +
  1308. +#include <math.h>
  1309. +#include <math_private.h>
  1310. +#include <fenv_libc.h>
  1311. +#include <inttypes.h>
  1312. +
  1313. +#include <sysdep.h>
  1314. +#include <ldsodefs.h>
  1315. +
  1316. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  1317. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  1318. +static const float two108 = 3.245185536584267269e+32;
  1319. +static const float twom54 = 5.551115123125782702e-17;
  1320. +static const float half = 0.5;
  1321. +
  1322. +/* The method is based on the descriptions in:
  1323. +
  1324. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  1325. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  1326. +
  1327. + We find the actual square root and half of its reciprocal
  1328. + simultaneously. */
  1329. +
  1330. +#ifdef __STDC__
  1331. +double
  1332. +__ieee754_sqrt (double b)
  1333. +#else
  1334. +double
  1335. +__ieee754_sqrt (b)
  1336. + double b;
  1337. +#endif
  1338. +{
  1339. + if (__builtin_expect (b > 0, 1))
  1340. + {
  1341. + double y, g, h, d, r;
  1342. + ieee_double_shape_type u;
  1343. +
  1344. + if (__builtin_expect (b != a_inf.value, 1))
  1345. + {
  1346. + fenv_t fe;
  1347. +
  1348. + fe = fegetenv_register ();
  1349. +
  1350. + u.value = b;
  1351. +
  1352. + relax_fenv_state ();
  1353. +
  1354. + __asm__ ("frsqrte %[estimate], %[x]\n"
  1355. + : [estimate] "=f" (y) : [x] "f" (b));
  1356. +
  1357. + /* Following Muller et al, page 168, equation 5.20.
  1358. +
  1359. + h goes to 1/(2*sqrt(b))
  1360. + g goes to sqrt(b).
  1361. +
  1362. + We need three iterations to get within 1ulp. */
  1363. +
  1364. + /* Indicate that these can be performed prior to the branch. GCC
  1365. + insists on sinking them below the branch, however; it seems like
  1366. + they'd be better before the branch so that we can cover any latency
  1367. + from storing the argument and loading its high word. Oh well. */
  1368. +
  1369. + g = b * y;
  1370. + h = 0.5 * y;
  1371. +
  1372. + /* Handle small numbers by scaling. */
  1373. + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
  1374. + return __ieee754_sqrt (b * two108) * twom54;
  1375. +
  1376. +#define FMADD(a_, c_, b_) \
  1377. + ({ double __r; \
  1378. + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
  1379. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1380. + __r;})
  1381. +#define FNMSUB(a_, c_, b_) \
  1382. + ({ double __r; \
  1383. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  1384. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1385. + __r;})
  1386. +
  1387. + r = FNMSUB (g, h, half);
  1388. + g = FMADD (g, r, g);
  1389. + h = FMADD (h, r, h);
  1390. +
  1391. + r = FNMSUB (g, h, half);
  1392. + g = FMADD (g, r, g);
  1393. + h = FMADD (h, r, h);
  1394. +
  1395. + r = FNMSUB (g, h, half);
  1396. + g = FMADD (g, r, g);
  1397. + h = FMADD (h, r, h);
  1398. +
  1399. + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
  1400. +
  1401. + /* Final refinement. */
  1402. + d = FNMSUB (g, g, b);
  1403. +
  1404. + fesetenv_register (fe);
  1405. + return FMADD (d, h, g);
  1406. + }
  1407. + }
  1408. + else if (b < 0)
  1409. + {
  1410. + /* For some reason, some PowerPC32 processors don't implement
  1411. + FE_INVALID_SQRT. */
  1412. +#ifdef FE_INVALID_SQRT
  1413. + feraiseexcept (FE_INVALID_SQRT);
  1414. +
  1415. + fenv_union_t u = { .fenv = fegetenv_register () };
  1416. + if ((u.l & FE_INVALID) == 0)
  1417. +#endif
  1418. + feraiseexcept (FE_INVALID);
  1419. + b = a_nan.value;
  1420. + }
  1421. + return f_wash (b);
  1422. +}
  1423. diff --git a/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
  1424. new file mode 100644
  1425. index 0000000000..26fa067abf
  1426. --- /dev/null
  1427. +++ b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
  1428. @@ -0,0 +1,101 @@
  1429. +/* Single-precision floating point square root.
  1430. + Copyright (C) 2010 Free Software Foundation, Inc.
  1431. + This file is part of the GNU C Library.
  1432. +
  1433. + The GNU C Library is free software; you can redistribute it and/or
  1434. + modify it under the terms of the GNU Lesser General Public
  1435. + License as published by the Free Software Foundation; either
  1436. + version 2.1 of the License, or (at your option) any later version.
  1437. +
  1438. + The GNU C Library is distributed in the hope that it will be useful,
  1439. + but WITHOUT ANY WARRANTY; without even the implied warranty of
  1440. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  1441. + Lesser General Public License for more details.
  1442. +
  1443. + You should have received a copy of the GNU Lesser General Public
  1444. + License along with the GNU C Library; if not, write to the Free
  1445. + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  1446. + 02111-1307 USA. */
  1447. +
  1448. +#include <math.h>
  1449. +#include <math_private.h>
  1450. +#include <fenv_libc.h>
  1451. +#include <inttypes.h>
  1452. +
  1453. +#include <sysdep.h>
  1454. +#include <ldsodefs.h>
  1455. +
  1456. +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
  1457. +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
  1458. +static const float threehalf = 1.5;
  1459. +
  1460. +/* The method is based on the descriptions in:
  1461. +
  1462. + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
  1463. + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
  1464. +
  1465. + We find the reciprocal square root and use that to compute the actual
  1466. + square root. */
  1467. +
  1468. +#ifdef __STDC__
  1469. +float
  1470. +__ieee754_sqrtf (float b)
  1471. +#else
  1472. +float
  1473. +__ieee754_sqrtf (b)
  1474. + float b;
  1475. +#endif
  1476. +{
  1477. + if (__builtin_expect (b > 0, 1))
  1478. + {
  1479. +#define FMSUB(a_, c_, b_) \
  1480. + ({ double __r; \
  1481. + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
  1482. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1483. + __r;})
  1484. +#define FNMSUB(a_, c_, b_) \
  1485. + ({ double __r; \
  1486. + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
  1487. + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
  1488. + __r;})
  1489. +
  1490. + if (__builtin_expect (b != a_inf.value, 1))
  1491. + {
  1492. + double y, x;
  1493. + fenv_t fe;
  1494. +
  1495. + fe = fegetenv_register ();
  1496. +
  1497. + relax_fenv_state ();
  1498. +
  1499. + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
  1500. + y = FMSUB (threehalf, b, b);
  1501. +
  1502. + /* Initial estimate. */
  1503. + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
  1504. +
  1505. + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
  1506. + x = x * FNMSUB (y, x * x, threehalf);
  1507. + x = x * FNMSUB (y, x * x, threehalf);
  1508. + x = x * FNMSUB (y, x * x, threehalf);
  1509. +
  1510. + /* All done. */
  1511. + fesetenv_register (fe);
  1512. + return x * b;
  1513. + }
  1514. + }
  1515. + else if (b < 0)
  1516. + {
  1517. + /* For some reason, some PowerPC32 processors don't implement
  1518. + FE_INVALID_SQRT. */
  1519. +#ifdef FE_INVALID_SQRT
  1520. + feraiseexcept (FE_INVALID_SQRT);
  1521. +
  1522. + fenv_union_t u = { .fenv = fegetenv_register () };
  1523. + if ((u.l & FE_INVALID) == 0)
  1524. +#endif
  1525. + feraiseexcept (FE_INVALID);
  1526. + b = a_nan.value;
  1527. + }
  1528. + return f_washf (b);
  1529. +}
  1530. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
  1531. new file mode 100644
  1532. index 0000000000..b103b4dea5
  1533. --- /dev/null
  1534. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
  1535. @@ -0,0 +1 @@
  1536. +powerpc/powerpc32/603e/fpu
  1537. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
  1538. new file mode 100644
  1539. index 0000000000..64db17fada
  1540. --- /dev/null
  1541. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
  1542. @@ -0,0 +1,2 @@
  1543. +# e300c3 is a variant of 603e so use the same optimizations for sqrt
  1544. +powerpc/powerpc32/603e/fpu
  1545. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
  1546. new file mode 100644
  1547. index 0000000000..7eac5fcf02
  1548. --- /dev/null
  1549. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
  1550. @@ -0,0 +1 @@
  1551. +powerpc/powerpc32/e500mc/fpu
  1552. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
  1553. new file mode 100644
  1554. index 0000000000..264b2a7700
  1555. --- /dev/null
  1556. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
  1557. @@ -0,0 +1 @@
  1558. +powerpc/powerpc32/e5500/fpu
  1559. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
  1560. new file mode 100644
  1561. index 0000000000..a25934467b
  1562. --- /dev/null
  1563. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
  1564. @@ -0,0 +1 @@
  1565. +powerpc/powerpc32/e6500/fpu
  1566. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
  1567. new file mode 100644
  1568. index 0000000000..a7bc854be8
  1569. --- /dev/null
  1570. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
  1571. @@ -0,0 +1 @@
  1572. +powerpc/powerpc64/e5500/fpu
  1573. diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
  1574. new file mode 100644
  1575. index 0000000000..04ff8cc181
  1576. --- /dev/null
  1577. +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
  1578. @@ -0,0 +1 @@
  1579. +powerpc/powerpc64/e6500/fpu