123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581 |
- From 74923ca4b1ae0ed5a2478e7d265b37534f6815d7 Mon Sep 17 00:00:00 2001
- From: Khem Raj <raj.khem@gmail.com>
- Date: Wed, 18 Mar 2015 00:01:50 +0000
- Subject: [PATCH] fsl e500/e5500/e6500/603e fsqrt implementation
- Upstream-Status: Pending
- Signed-off-by: Edmar Wienskoski <edmar@freescale.com>
- Signed-off-by: Khem Raj <raj.khem@gmail.com>
- ---
- sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c | 134 ++++++++++++++++++
- sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c | 101 +++++++++++++
- sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c | 134 ++++++++++++++++++
- .../powerpc/powerpc32/e500mc/fpu/e_sqrtf.c | 101 +++++++++++++
- sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c | 134 ++++++++++++++++++
- sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c | 101 +++++++++++++
- sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c | 134 ++++++++++++++++++
- sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c | 101 +++++++++++++
- sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c | 134 ++++++++++++++++++
- sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c | 101 +++++++++++++
- sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c | 134 ++++++++++++++++++
- sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c | 101 +++++++++++++
- .../linux/powerpc/powerpc32/603e/fpu/Implies | 1 +
- .../powerpc/powerpc32/e300c3/fpu/Implies | 2 +
- .../powerpc/powerpc32/e500mc/fpu/Implies | 1 +
- .../linux/powerpc/powerpc32/e5500/fpu/Implies | 1 +
- .../linux/powerpc/powerpc32/e6500/fpu/Implies | 1 +
- .../linux/powerpc/powerpc64/e5500/fpu/Implies | 1 +
- .../linux/powerpc/powerpc64/e6500/fpu/Implies | 1 +
- 19 files changed, 1418 insertions(+)
- create mode 100644 sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
- create mode 100644 sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
- create mode 100644 sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
- create mode 100644 sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
- create mode 100644 sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
- create mode 100644 sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
- create mode 100644 sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
- create mode 100644 sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
- diff --git a/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/603e/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e500mc/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e5500/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc32/e6500/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc64/e5500/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
- new file mode 100644
- index 0000000000..71e516d1c8
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrt.c
- @@ -0,0 +1,134 @@
- +/* Double-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float two108 = 3.245185536584267269e+32;
- +static const float twom54 = 5.551115123125782702e-17;
- +static const float half = 0.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the actual square root and half of its reciprocal
- + simultaneously. */
- +
- +#ifdef __STDC__
- +double
- +__ieee754_sqrt (double b)
- +#else
- +double
- +__ieee754_sqrt (b)
- + double b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- + double y, g, h, d, r;
- + ieee_double_shape_type u;
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + u.value = b;
- +
- + relax_fenv_state ();
- +
- + __asm__ ("frsqrte %[estimate], %[x]\n"
- + : [estimate] "=f" (y) : [x] "f" (b));
- +
- + /* Following Muller et al, page 168, equation 5.20.
- +
- + h goes to 1/(2*sqrt(b))
- + g goes to sqrt(b).
- +
- + We need three iterations to get within 1ulp. */
- +
- + /* Indicate that these can be performed prior to the branch. GCC
- + insists on sinking them below the branch, however; it seems like
- + they'd be better before the branch so that we can cover any latency
- + from storing the argument and loading its high word. Oh well. */
- +
- + g = b * y;
- + h = 0.5 * y;
- +
- + /* Handle small numbers by scaling. */
- + if (__builtin_expect ((u.parts.msw & 0x7ff00000) <= 0x02000000, 0))
- + return __ieee754_sqrt (b * two108) * twom54;
- +
- +#define FMADD(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmadd %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + r = FNMSUB (g, h, half);
- + g = FMADD (g, r, g);
- + h = FMADD (h, r, h);
- +
- + /* g is now +/- 1ulp, or exactly equal to, the square root of b. */
- +
- + /* Final refinement. */
- + d = FNMSUB (g, g, b);
- +
- + fesetenv_register (fe);
- + return FMADD (d, h, g);
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_wash (b);
- +}
- diff --git a/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
- new file mode 100644
- index 0000000000..26fa067abf
- --- /dev/null
- +++ b/sysdeps/powerpc/powerpc64/e6500/fpu/e_sqrtf.c
- @@ -0,0 +1,101 @@
- +/* Single-precision floating point square root.
- + Copyright (C) 2010 Free Software Foundation, Inc.
- + This file is part of the GNU C Library.
- +
- + The GNU C Library is free software; you can redistribute it and/or
- + modify it under the terms of the GNU Lesser General Public
- + License as published by the Free Software Foundation; either
- + version 2.1 of the License, or (at your option) any later version.
- +
- + The GNU C Library is distributed in the hope that it will be useful,
- + but WITHOUT ANY WARRANTY; without even the implied warranty of
- + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- + Lesser General Public License for more details.
- +
- + You should have received a copy of the GNU Lesser General Public
- + License along with the GNU C Library; if not, write to the Free
- + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- + 02111-1307 USA. */
- +
- +#include <math.h>
- +#include <math_private.h>
- +#include <fenv_libc.h>
- +#include <inttypes.h>
- +
- +#include <sysdep.h>
- +#include <ldsodefs.h>
- +
- +static const ieee_float_shape_type a_nan = {.word = 0x7fc00000 };
- +static const ieee_float_shape_type a_inf = {.word = 0x7f800000 };
- +static const float threehalf = 1.5;
- +
- +/* The method is based on the descriptions in:
- +
- + _The Handbook of Floating-Pointer Arithmetic_ by Muller et al., chapter 5;
- + _IA-64 and Elementary Functions: Speed and Precision_ by Markstein, chapter 9
- +
- + We find the reciprocal square root and use that to compute the actual
- + square root. */
- +
- +#ifdef __STDC__
- +float
- +__ieee754_sqrtf (float b)
- +#else
- +float
- +__ieee754_sqrtf (b)
- + float b;
- +#endif
- +{
- + if (__builtin_expect (b > 0, 1))
- + {
- +#define FMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +#define FNMSUB(a_, c_, b_) \
- + ({ double __r; \
- + __asm__ ("fnmsub %[r], %[a], %[c], %[b]\n" \
- + : [r] "=f" (__r) : [a] "f" (a_), [c] "f" (c_), [b] "f" (b_)); \
- + __r;})
- +
- + if (__builtin_expect (b != a_inf.value, 1))
- + {
- + double y, x;
- + fenv_t fe;
- +
- + fe = fegetenv_register ();
- +
- + relax_fenv_state ();
- +
- + /* Compute y = 1.5 * b - b. Uses fewer constants than y = 0.5 * b. */
- + y = FMSUB (threehalf, b, b);
- +
- + /* Initial estimate. */
- + __asm__ ("frsqrte %[x], %[b]\n" : [x] "=f" (x) : [b] "f" (b));
- +
- + /* Iterate. x_{n+1} = x_n * (1.5 - y * (x_n * x_n)). */
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- + x = x * FNMSUB (y, x * x, threehalf);
- +
- + /* All done. */
- + fesetenv_register (fe);
- + return x * b;
- + }
- + }
- + else if (b < 0)
- + {
- + /* For some reason, some PowerPC32 processors don't implement
- + FE_INVALID_SQRT. */
- +#ifdef FE_INVALID_SQRT
- + feraiseexcept (FE_INVALID_SQRT);
- +
- + fenv_union_t u = { .fenv = fegetenv_register () };
- + if ((u.l & FE_INVALID) == 0)
- +#endif
- + feraiseexcept (FE_INVALID);
- + b = a_nan.value;
- + }
- + return f_washf (b);
- +}
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
- new file mode 100644
- index 0000000000..b103b4dea5
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/603e/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc32/603e/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
- new file mode 100644
- index 0000000000..64db17fada
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e300c3/fpu/Implies
- @@ -0,0 +1,2 @@
- +# e300c3 is a variant of 603e so use the same optimizations for sqrt
- +powerpc/powerpc32/603e/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
- new file mode 100644
- index 0000000000..7eac5fcf02
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e500mc/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc32/e500mc/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
- new file mode 100644
- index 0000000000..264b2a7700
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e5500/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc32/e5500/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
- new file mode 100644
- index 0000000000..a25934467b
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/e6500/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc32/e6500/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
- new file mode 100644
- index 0000000000..a7bc854be8
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e5500/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc64/e5500/fpu
- diff --git a/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
- new file mode 100644
- index 0000000000..04ff8cc181
- --- /dev/null
- +++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/e6500/fpu/Implies
- @@ -0,0 +1 @@
- +powerpc/powerpc64/e6500/fpu
|