Added x86_64 math functions to libroot.

These have been pulled from glibc 2.10, the last version before a
bunch of changes were made that would have made porting more difficult.
None of this has been tested yet, it is not currently possible to do
so: I'm just trying to get libroot compiling so that I can work on the
runtime loader. I will test them when I am able to.
This commit is contained in:
Alex Smith
2012-07-28 12:23:47 +01:00
parent 04c1e6d017
commit a3210bbc76
112 changed files with 9139 additions and 9 deletions
@@ -41,9 +41,6 @@ const fenv_t __fe_dfl_env = {
__INITIAL_MXCSR__
};
extern inline int feclearexcept(int __excepts);
extern inline int fegetexceptflag(fexcept_t *__flagp, int __excepts);
int
fesetexceptflag(const fexcept_t *flagp, int excepts)
{
@@ -72,10 +69,6 @@ feraiseexcept(int excepts)
return (0);
}
extern inline int fetestexcept(int __excepts);
extern inline int fegetround(void);
extern inline int fesetround(int __round);
int
fegetenv(fenv_t *envp)
{
@@ -105,8 +98,6 @@ feholdexcept(fenv_t *envp)
return (0);
}
extern inline int fesetenv(const fenv_t *__envp);
int
feupdateenv(const fenv_t *envp)
{
@@ -0,0 +1,79 @@
/* s_atanhl.c -- long double version of s_atan.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* __ieee754_atanhl(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
*
* Special cases:
* atanhl(x) is NaN if |x| > 1 with signal;
* atanhl(NaN) is that NaN with no signal;
* atanhl(+-1) is +-INF with signal.
*
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double one = 1.0, huge = 1e4900L;
#else
static long double one = 1.0, huge = 1e4900L;
#endif
#ifdef __STDC__
static const long double zero = 0.0;
#else
static double long zero = 0.0;
#endif
#ifdef __STDC__
long double __ieee754_atanhl(long double x)
#else
long double __ieee754_atanhl(x)
long double x;
#endif
{
long double t;
int32_t ix;
u_int32_t se,i0,i1;
GET_LDOUBLE_WORDS(se,i0,i1,x);
ix = se&0x7fff;
if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff)
/* |x|>1 */
return (x-x)/(x-x);
if(ix==0x3fff)
return x/zero;
if(ix<0x3fe3&&(huge+x)>zero) return x; /* x<2**-28 */
SET_LDOUBLE_EXP(x,ix);
if(ix<0x3ffe) { /* x < 0.5 */
t = x+x;
t = 0.5*__log1pl(t+t*x/(one-x));
} else
t = 0.5*__log1pl((x+x)/(one-x));
if(se<=0x7fff) return t; else return -t;
}
@@ -0,0 +1,644 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Long double expansions are
Copyright (C) 2001 Stephen L. Moshier <[email protected]>
and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under
the following terms:
This 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.
This 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 this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j0(x):
* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
* 2. Reduce x to |x| since j0(x)=j0(-x), and
* for x in (0,2)
* j0(x) = 1 - z/4 + z^2*R0/S0, where z = x*x;
* for x in (2,inf)
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* as follow:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (cos(x) + sin(x))
* sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j0(nan)= nan
* j0(0) = 1
* j0(inf) = 0
*
* Method -- y0(x):
* 1. For x<2.
* Since
* y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
* therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
* We use the following function to approximate y0,
* y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
*
* Note: For tiny x, U/V = u0 and j0(x)~1, hence
* y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
* 2. For x>=2.
* y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* by the method mentioned above.
* 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static long double pzero (long double), qzero (long double);
#else
static long double pzero (), qzero ();
#endif
#ifdef __STDC__
static const long double
#else
static long double
#endif
huge = 1e4930L,
one = 1.0L,
invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
tpi = 6.3661977236758134307553505349005744813784e-1L,
/* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
0 <= x <= 2
peak relative error 1.41e-22 */
R[5] = {
4.287176872744686992880841716723478740566E7L,
-6.652058897474241627570911531740907185772E5L,
7.011848381719789863458364584613651091175E3L,
-3.168040850193372408702135490809516253693E1L,
6.030778552661102450545394348845599300939E-2L,
},
S[4] = {
2.743793198556599677955266341699130654342E9L,
3.364330079384816249840086842058954076201E7L,
1.924119649412510777584684927494642526573E5L,
6.239282256012734914211715620088714856494E2L,
/* 1.000000000000000000000000000000000000000E0L,*/
};
#ifdef __STDC__
static const long double zero = 0.0;
#else
static long double zero = 0.0;
#endif
#ifdef __STDC__
long double
__ieee754_j0l (long double x)
#else
long double
__ieee754_j0l (x)
long double x;
#endif
{
long double z, s, c, ss, cc, r, u, v;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff)
return one / (x * x);
x = fabsl (x);
if (ix >= 0x4000) /* |x| >= 2.0 */
{
__sincosl (x, &s, &c);
ss = s - c;
cc = s + c;
if (ix < 0x7ffe)
{ /* make sure x+x not overflow */
z = -__cosl (x + x);
if ((s * c) < zero)
cc = z / ss;
else
ss = z / cc;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix > 0x4080) /* 2^129 */
z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
else
{
u = pzero (x);
v = qzero (x);
z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x);
}
return z;
}
if (ix < 0x3fef) /* |x| < 2**-16 */
{
if (huge + x > one)
{ /* raise inexact if x != 0 */
if (ix < 0x3fde) /* |x| < 2^-33 */
return one;
else
return one - 0.25 * x * x;
}
}
z = x * x;
r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4]))));
s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
if (ix < 0x3fff)
{ /* |x| < 1.00 */
return (one - 0.25 * z + z * (r / s));
}
else
{
u = 0.5 * x;
return ((one + u) * (one - u) + z * (r / s));
}
}
/* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
0 < x <= 2
peak relative error 1.7e-21 */
#ifdef __STDC__
static const long double
#else
static long double
#endif
U[6] = {
-1.054912306975785573710813351985351350861E10L,
2.520192609749295139432773849576523636127E10L,
-1.856426071075602001239955451329519093395E9L,
4.079209129698891442683267466276785956784E7L,
-3.440684087134286610316661166492641011539E5L,
1.005524356159130626192144663414848383774E3L,
};
#ifdef __STDC__
static const long double
#else
static long double
#endif
V[5] = {
1.429337283720789610137291929228082613676E11L,
2.492593075325119157558811370165695013002E9L,
2.186077620785925464237324417623665138376E7L,
1.238407896366385175196515057064384929222E5L,
4.693924035211032457494368947123233101664E2L,
/* 1.000000000000000000000000000000000000000E0L */
};
#ifdef __STDC__
long double
__ieee754_y0l (long double x)
#else
long double
__ieee754_y0l (x)
long double x;
#endif
{
long double z, s, c, ss, cc, u, v;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (se & 0x8000)
return zero / (zero * x);
if (ix >= 0x7fff)
return one / (x + x * x);
if ((i0 | i1) == 0)
return -HUGE_VALL + x; /* -inf and overflow exception. */
if (ix >= 0x4000)
{ /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) + cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
__sincosl (x, &s, &c);
ss = s - c;
cc = s + c;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix < 0x7ffe)
{ /* make sure x+x not overflow */
z = -__cosl (x + x);
if ((s * c) < zero)
cc = z / ss;
else
ss = z / cc;
}
if (ix > 0x4080) /* 1e39 */
z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
else
{
u = pzero (x);
v = qzero (x);
z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
}
return z;
}
if (ix <= 0x3fde) /* x < 2^-33 */
{
z = -7.380429510868722527629822444004602747322E-2L
+ tpi * __ieee754_logl (x);
return z;
}
z = x * x;
u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
}
/* The asymptotic expansions of pzero is
* 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
* For x >= 2, We approximate pzero by
* pzero(x) = 1 + s^2 R(s^2) / S(s^2)
*/
#ifdef __STDC__
static const long double pR8[7] = {
#else
static long double pR8[7] = {
#endif
/* 8 <= x <= inf
Peak relative error 4.62 */
-4.094398895124198016684337960227780260127E-9L,
-8.929643669432412640061946338524096893089E-7L,
-6.281267456906136703868258380673108109256E-5L,
-1.736902783620362966354814353559382399665E-3L,
-1.831506216290984960532230842266070146847E-2L,
-5.827178869301452892963280214772398135283E-2L,
-2.087563267939546435460286895807046616992E-2L,
};
#ifdef __STDC__
static const long double pS8[6] = {
#else
static long double pS8[6] = {
#endif
5.823145095287749230197031108839653988393E-8L,
1.279281986035060320477759999428992730280E-5L,
9.132668954726626677174825517150228961304E-4L,
2.606019379433060585351880541545146252534E-2L,
2.956262215119520464228467583516287175244E-1L,
1.149498145388256448535563278632697465675E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double pR5[7] = {
#else
static long double pR5[7] = {
#endif
/* 4.54541015625 <= x <= 8
Peak relative error 6.51E-22 */
-2.041226787870240954326915847282179737987E-7L,
-2.255373879859413325570636768224534428156E-5L,
-7.957485746440825353553537274569102059990E-4L,
-1.093205102486816696940149222095559439425E-2L,
-5.657957849316537477657603125260701114646E-2L,
-8.641175552716402616180994954177818461588E-2L,
-1.354654710097134007437166939230619726157E-2L,
};
#ifdef __STDC__
static const long double pS5[6] = {
#else
static long double pS5[6] = {
#endif
2.903078099681108697057258628212823545290E-6L,
3.253948449946735405975737677123673867321E-4L,
1.181269751723085006534147920481582279979E-2L,
1.719212057790143888884745200257619469363E-1L,
1.006306498779212467670654535430694221924E0L,
2.069568808688074324555596301126375951502E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double pR3[7] = {
#else
static long double pR3[7] = {
#endif
/* 2.85711669921875 <= x <= 4.54541015625
peak relative error 5.25e-21 */
-5.755732156848468345557663552240816066802E-6L,
-3.703675625855715998827966962258113034767E-4L,
-7.390893350679637611641350096842846433236E-3L,
-5.571922144490038765024591058478043873253E-2L,
-1.531290690378157869291151002472627396088E-1L,
-1.193350853469302941921647487062620011042E-1L,
-8.567802507331578894302991505331963782905E-3L,
};
#ifdef __STDC__
static const long double pS3[6] = {
#else
static long double pS3[6] = {
#endif
8.185931139070086158103309281525036712419E-5L,
5.398016943778891093520574483111255476787E-3L,
1.130589193590489566669164765853409621081E-1L,
9.358652328786413274673192987670237145071E-1L,
3.091711512598349056276917907005098085273E0L,
3.594602474737921977972586821673124231111E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double pR2[7] = {
#else
static long double pR2[7] = {
#endif
/* 2 <= x <= 2.85711669921875
peak relative error 2.64e-21 */
-1.219525235804532014243621104365384992623E-4L,
-4.838597135805578919601088680065298763049E-3L,
-5.732223181683569266223306197751407418301E-2L,
-2.472947430526425064982909699406646503758E-1L,
-3.753373645974077960207588073975976327695E-1L,
-1.556241316844728872406672349347137975495E-1L,
-5.355423239526452209595316733635519506958E-3L,
};
#ifdef __STDC__
static const long double pS2[6] = {
#else
static long double pS2[6] = {
#endif
1.734442793664291412489066256138894953823E-3L,
7.158111826468626405416300895617986926008E-2L,
9.153839713992138340197264669867993552641E-1L,
4.539209519433011393525841956702487797582E0L,
8.868932430625331650266067101752626253644E0L,
6.067161890196324146320763844772857713502E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static long double
pzero (long double x)
#else
static long double
pzero (x)
long double x;
#endif
{
#ifdef __STDC__
const long double *p, *q;
#else
long double *p, *q;
#endif
long double z, r, s;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x4002)
{
p = pR8;
q = pS8;
} /* x >= 8 */
else
{
i1 = (ix << 16) | (i0 >> 16);
if (i1 >= 0x40019174) /* x >= 4.54541015625 */
{
p = pR5;
q = pS5;
}
else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
{
p = pR3;
q = pS3;
}
else if (ix >= 0x4000) /* x better be >= 2 */
{
p = pR2;
q = pS2;
}
}
z = one / (x * x);
r =
p[0] + z * (p[1] +
z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
s =
q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
return (one + z * r / s);
}
/* For x >= 8, the asymptotic expansions of qzero is
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
* We approximate qzero by
* qzero(x) = s*(-.125 + R(s^2) / S(s^2))
*/
#ifdef __STDC__
static const long double qR8[7] = {
#else
static long double qR8[7] = {
#endif
/* 8 <= x <= inf
peak relative error 2.23e-21 */
3.001267180483191397885272640777189348008E-10L,
8.693186311430836495238494289942413810121E-8L,
8.496875536711266039522937037850596580686E-6L,
3.482702869915288984296602449543513958409E-4L,
6.036378380706107692863811938221290851352E-3L,
3.881970028476167836382607922840452192636E-2L,
6.132191514516237371140841765561219149638E-2L,
};
#ifdef __STDC__
static const long double qS8[7] = {
#else
static long double qS8[7] = {
#endif
4.097730123753051126914971174076227600212E-9L,
1.199615869122646109596153392152131139306E-6L,
1.196337580514532207793107149088168946451E-4L,
5.099074440112045094341500497767181211104E-3L,
9.577420799632372483249761659674764460583E-2L,
7.385243015344292267061953461563695918646E-1L,
1.917266424391428937962682301561699055943E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double qR5[7] = {
#else
static long double qR5[7] = {
#endif
/* 4.54541015625 <= x <= 8
peak relative error 1.03e-21 */
3.406256556438974327309660241748106352137E-8L,
4.855492710552705436943630087976121021980E-6L,
2.301011739663737780613356017352912281980E-4L,
4.500470249273129953870234803596619899226E-3L,
3.651376459725695502726921248173637054828E-2L,
1.071578819056574524416060138514508609805E-1L,
7.458950172851611673015774675225656063757E-2L,
};
#ifdef __STDC__
static const long double qS5[7] = {
#else
static long double qS5[7] = {
#endif
4.650675622764245276538207123618745150785E-7L,
6.773573292521412265840260065635377164455E-5L,
3.340711249876192721980146877577806687714E-3L,
7.036218046856839214741678375536970613501E-2L,
6.569599559163872573895171876511377891143E-1L,
2.557525022583599204591036677199171155186E0L,
3.457237396120935674982927714210361269133E0L,
/* 1.000000000000000000000000000000000000000E0L,*/
};
#ifdef __STDC__
static const long double qR3[7] = {
#else
static long double qR3[7] = {
#endif
/* 2.85711669921875 <= x <= 4.54541015625
peak relative error 5.24e-21 */
1.749459596550816915639829017724249805242E-6L,
1.446252487543383683621692672078376929437E-4L,
3.842084087362410664036704812125005761859E-3L,
4.066369994699462547896426554180954233581E-2L,
1.721093619117980251295234795188992722447E-1L,
2.538595333972857367655146949093055405072E-1L,
8.560591367256769038905328596020118877936E-2L,
};
#ifdef __STDC__
static const long double qS3[7] = {
#else
static long double qS3[7] = {
#endif
2.388596091707517488372313710647510488042E-5L,
2.048679968058758616370095132104333998147E-3L,
5.824663198201417760864458765259945181513E-2L,
6.953906394693328750931617748038994763958E-1L,
3.638186936390881159685868764832961092476E0L,
7.900169524705757837298990558459547842607E0L,
5.992718532451026507552820701127504582907E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double qR2[7] = {
#else
static long double qR2[7] = {
#endif
/* 2 <= x <= 2.85711669921875
peak relative error 1.58e-21 */
6.306524405520048545426928892276696949540E-5L,
3.209606155709930950935893996591576624054E-3L,
5.027828775702022732912321378866797059604E-2L,
3.012705561838718956481911477587757845163E-1L,
6.960544893905752937420734884995688523815E-1L,
5.431871999743531634887107835372232030655E-1L,
9.447736151202905471899259026430157211949E-2L,
};
#ifdef __STDC__
static const long double qS2[7] = {
#else
static long double qS2[7] = {
#endif
8.610579901936193494609755345106129102676E-4L,
4.649054352710496997203474853066665869047E-2L,
8.104282924459837407218042945106320388339E-1L,
5.807730930825886427048038146088828206852E0L,
1.795310145936848873627710102199881642939E1L,
2.281313316875375733663657188888110605044E1L,
1.011242067883822301487154844458322200143E1L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static long double
qzero (long double x)
#else
static long double
qzero (x)
long double x;
#endif
{
#ifdef __STDC__
const long double *p, *q;
#else
long double *p, *q;
#endif
long double s, r, z;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x4002) /* x >= 8 */
{
p = qR8;
q = qS8;
}
else
{
i1 = (ix << 16) | (i0 >> 16);
if (i1 >= 0x40019174) /* x >= 4.54541015625 */
{
p = qR5;
q = qS5;
}
else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
{
p = qR3;
q = qS3;
}
else if (ix >= 0x4000) /* x better be >= 2 */
{
p = qR2;
q = qS2;
}
}
z = one / (x * x);
r =
p[0] + z * (p[1] +
z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
s =
q[0] + z * (q[1] +
z * (q[2] +
z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
return (-.125 + z * r / s) / x;
}
@@ -0,0 +1,651 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Long double expansions are
Copyright (C) 2001 Stephen L. Moshier <[email protected]>
and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under
the following terms:
This 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.
This 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 this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j1(x):
* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
* 2. Reduce x to |x| since j1(x)=-j1(-x), and
* for x in (0,2)
* j1(x) = x/2 + x*z*R0/S0, where z = x*x;
* for x in (2,inf)
* j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* as follow:
* cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (sin(x) + cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j1(nan)= nan
* j1(0) = 0
* j1(inf) = 0
*
* Method -- y1(x):
* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
* 2. For x<2.
* Since
* y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
* therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
* We use the following function to approximate y1,
* y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
* Note: For tiny x, 1/x dominate y1 and hence
* y1(tiny) = -2/pi/tiny
* 3. For x>=2.
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* by method mentioned above.
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static long double pone (long double), qone (long double);
#else
static long double pone (), qone ();
#endif
#ifdef __STDC__
static const long double
#else
static long double
#endif
huge = 1e4930L,
one = 1.0L,
invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
tpi = 6.3661977236758134307553505349005744813784e-1L,
/* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
0 <= x <= 2
Peak relative error 4.5e-21 */
R[5] = {
-9.647406112428107954753770469290757756814E7L,
2.686288565865230690166454005558203955564E6L,
-3.689682683905671185891885948692283776081E4L,
2.195031194229176602851429567792676658146E2L,
-5.124499848728030297902028238597308971319E-1L,
},
S[4] =
{
1.543584977988497274437410333029029035089E9L,
2.133542369567701244002565983150952549520E7L,
1.394077011298227346483732156167414670520E5L,
5.252401789085732428842871556112108446506E2L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static const long double zero = 0.0;
#else
static long double zero = 0.0;
#endif
#ifdef __STDC__
long double
__ieee754_j1l (long double x)
#else
long double
__ieee754_j1l (x)
long double x;
#endif
{
long double z, c, r, s, ss, cc, u, v, y;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff)
return one / x;
y = fabsl (x);
if (ix >= 0x4000)
{ /* |x| >= 2.0 */
__sincosl (y, &s, &c);
ss = -s - c;
cc = s - c;
if (ix < 0x7ffe)
{ /* make sure y+y not overflow */
z = __cosl (y + y);
if ((s * c) > zero)
cc = z / ss;
else
ss = z / cc;
}
/*
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
*/
if (ix > 0x4080)
z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
else
{
u = pone (y);
v = qone (y);
z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
}
if (se & 0x8000)
return -z;
else
return z;
}
if (ix < 0x3fde) /* |x| < 2^-33 */
{
if (huge + x > one)
return 0.5 * x; /* inexact if x!=0 necessary */
}
z = x * x;
r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
r *= x;
return (x * 0.5 + r / s);
}
/* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
0 <= x <= 2
Peak relative error 2.3e-23 */
#ifdef __STDC__
static const long double U0[6] = {
#else
static long double U0[6] = {
#endif
-5.908077186259914699178903164682444848615E10L,
1.546219327181478013495975514375773435962E10L,
-6.438303331169223128870035584107053228235E8L,
9.708540045657182600665968063824819371216E6L,
-6.138043997084355564619377183564196265471E4L,
1.418503228220927321096904291501161800215E2L,
};
#ifdef __STDC__
static const long double V0[5] = {
#else
static long double V0[5] = {
#endif
3.013447341682896694781964795373783679861E11L,
4.669546565705981649470005402243136124523E9L,
3.595056091631351184676890179233695857260E7L,
1.761554028569108722903944659933744317994E5L,
5.668480419646516568875555062047234534863E2L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
long double
__ieee754_y1l (long double x)
#else
long double
__ieee754_y1l (x)
long double x;
#endif
{
long double z, s, c, ss, cc, u, v;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (se & 0x8000)
return zero / (zero * x);
if (ix >= 0x7fff)
return one / (x + x * x);
if ((i0 | i1) == 0)
return -HUGE_VALL + x; /* -inf and overflow exception. */
if (ix >= 0x4000)
{ /* |x| >= 2.0 */
__sincosl (x, &s, &c);
ss = -s - c;
cc = s - c;
if (ix < 0x7fe00000)
{ /* make sure x+x not overflow */
z = __cosl (x + x);
if ((s * c) > zero)
cc = z / ss;
else
ss = z / cc;
}
/* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
* where x0 = x-3pi/4
* Better formula:
* cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (cos(x) + sin(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if (ix > 0x4080)
z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
else
{
u = pone (x);
v = qone (x);
z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
}
return z;
}
if (ix <= 0x3fbe)
{ /* x < 2**-65 */
return (-tpi / x);
}
z = x * x;
u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
return (x * (u / v) +
tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
}
/* For x >= 8, the asymptotic expansions of pone is
* 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
* We approximate pone by
* pone(x) = 1 + (R/S)
*/
/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
P1(x) = 1 + z^2 R(z^2), z=1/x
8 <= x <= inf (0 <= z <= 0.125)
Peak relative error 5.2e-22 */
#ifdef __STDC__
static const long double pr8[7] = {
#else
static long double pr8[7] = {
#endif
8.402048819032978959298664869941375143163E-9L,
1.813743245316438056192649247507255996036E-6L,
1.260704554112906152344932388588243836276E-4L,
3.439294839869103014614229832700986965110E-3L,
3.576910849712074184504430254290179501209E-2L,
1.131111483254318243139953003461511308672E-1L,
4.480715825681029711521286449131671880953E-2L,
};
#ifdef __STDC__
static const long double ps8[6] = {
#else
static long double ps8[6] = {
#endif
7.169748325574809484893888315707824924354E-8L,
1.556549720596672576431813934184403614817E-5L,
1.094540125521337139209062035774174565882E-3L,
3.060978962596642798560894375281428805840E-2L,
3.374146536087205506032643098619414507024E-1L,
1.253830208588979001991901126393231302559E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
P1(x) = 1 + z^2 R(z^2), z=1/x
4.54541015625 <= x <= 8
Peak relative error 7.7e-22 */
#ifdef __STDC__
static const long double pr5[7] = {
#else
static long double pr5[7] = {
#endif
4.318486887948814529950980396300969247900E-7L,
4.715341880798817230333360497524173929315E-5L,
1.642719430496086618401091544113220340094E-3L,
2.228688005300803935928733750456396149104E-2L,
1.142773760804150921573259605730018327162E-1L,
1.755576530055079253910829652698703791957E-1L,
3.218803858282095929559165965353784980613E-2L,
};
#ifdef __STDC__
static const long double ps5[6] = {
#else
static long double ps5[6] = {
#endif
3.685108812227721334719884358034713967557E-6L,
4.069102509511177498808856515005792027639E-4L,
1.449728676496155025507893322405597039816E-2L,
2.058869213229520086582695850441194363103E-1L,
1.164890985918737148968424972072751066553E0L,
2.274776933457009446573027260373361586841E0L,
/* 1.000000000000000000000000000000000000000E0L,*/
};
/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
P1(x) = 1 + z^2 R(z^2), z=1/x
2.85711669921875 <= x <= 4.54541015625
Peak relative error 6.5e-21 */
#ifdef __STDC__
static const long double pr3[7] = {
#else
static long double pr3[7] = {
#endif
1.265251153957366716825382654273326407972E-5L,
8.031057269201324914127680782288352574567E-4L,
1.581648121115028333661412169396282881035E-2L,
1.179534658087796321928362981518645033967E-1L,
3.227936912780465219246440724502790727866E-1L,
2.559223765418386621748404398017602935764E-1L,
2.277136933287817911091370397134882441046E-2L,
};
#ifdef __STDC__
static const long double ps3[6] = {
#else
static long double ps3[6] = {
#endif
1.079681071833391818661952793568345057548E-4L,
6.986017817100477138417481463810841529026E-3L,
1.429403701146942509913198539100230540503E-1L,
1.148392024337075609460312658938700765074E0L,
3.643663015091248720208251490291968840882E0L,
3.990702269032018282145100741746633960737E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
P1(x) = 1 + z^2 R(z^2), z=1/x
2 <= x <= 2.85711669921875
Peak relative error 3.5e-21 */
#ifdef __STDC__
static const long double pr2[7] = {
#else
static long double pr2[7] = {
#endif
2.795623248568412225239401141338714516445E-4L,
1.092578168441856711925254839815430061135E-2L,
1.278024620468953761154963591853679640560E-1L,
5.469680473691500673112904286228351988583E-1L,
8.313769490922351300461498619045639016059E-1L,
3.544176317308370086415403567097130611468E-1L,
1.604142674802373041247957048801599740644E-2L,
};
#ifdef __STDC__
static const long double ps2[6] = {
#else
static long double ps2[6] = {
#endif
2.385605161555183386205027000675875235980E-3L,
9.616778294482695283928617708206967248579E-2L,
1.195215570959693572089824415393951258510E0L,
5.718412857897054829999458736064922974662E0L,
1.065626298505499086386584642761602177568E1L,
6.809140730053382188468983548092322151791E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static long double
pone (long double x)
#else
static long double
pone (x)
long double x;
#endif
{
#ifdef __STDC__
const long double *p, *q;
#else
long double *p, *q;
#endif
long double z, r, s;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x4002) /* x >= 8 */
{
p = pr8;
q = ps8;
}
else
{
i1 = (ix << 16) | (i0 >> 16);
if (i1 >= 0x40019174) /* x >= 4.54541015625 */
{
p = pr5;
q = ps5;
}
else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
{
p = pr3;
q = ps3;
}
else if (ix >= 0x4000) /* x better be >= 2 */
{
p = pr2;
q = ps2;
}
}
z = one / (x * x);
r = p[0] + z * (p[1] +
z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
return one + z * r / s;
}
/* For x >= 8, the asymptotic expansions of qone is
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
* We approximate pone by
* qone(x) = s*(0.375 + (R/S))
*/
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
8 <= x <= inf
Peak relative error 8.3e-22 */
#ifdef __STDC__
static const long double qr8[7] = {
#else
static long double qr8[7] = {
#endif
-5.691925079044209246015366919809404457380E-10L,
-1.632587664706999307871963065396218379137E-7L,
-1.577424682764651970003637263552027114600E-5L,
-6.377627959241053914770158336842725291713E-4L,
-1.087408516779972735197277149494929568768E-2L,
-6.854943629378084419631926076882330494217E-2L,
-1.055448290469180032312893377152490183203E-1L,
};
#ifdef __STDC__
static const long double qs8[7] = {
#else
static long double qs8[7] = {
#endif
5.550982172325019811119223916998393907513E-9L,
1.607188366646736068460131091130644192244E-6L,
1.580792530091386496626494138334505893599E-4L,
6.617859900815747303032860443855006056595E-3L,
1.212840547336984859952597488863037659161E-1L,
9.017885953937234900458186716154005541075E-1L,
2.201114489712243262000939120146436167178E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
4.54541015625 <= x <= 8
Peak relative error 4.1e-22 */
#ifdef __STDC__
static const long double qr5[7] = {
#else
static long double qr5[7] = {
#endif
-6.719134139179190546324213696633564965983E-8L,
-9.467871458774950479909851595678622044140E-6L,
-4.429341875348286176950914275723051452838E-4L,
-8.539898021757342531563866270278505014487E-3L,
-6.818691805848737010422337101409276287170E-2L,
-1.964432669771684034858848142418228214855E-1L,
-1.333896496989238600119596538299938520726E-1L,
};
#ifdef __STDC__
static const long double qs5[7] = {
#else
static long double qs5[7] = {
#endif
6.552755584474634766937589285426911075101E-7L,
9.410814032118155978663509073200494000589E-5L,
4.561677087286518359461609153655021253238E-3L,
9.397742096177905170800336715661091535805E-2L,
8.518538116671013902180962914473967738771E-1L,
3.177729183645800174212539541058292579009E0L,
4.006745668510308096259753538973038902990E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
2.85711669921875 <= x <= 4.54541015625
Peak relative error 2.2e-21 */
#ifdef __STDC__
static const long double qr3[7] = {
#else
static long double qr3[7] = {
#endif
-3.618746299358445926506719188614570588404E-6L,
-2.951146018465419674063882650970344502798E-4L,
-7.728518171262562194043409753656506795258E-3L,
-8.058010968753999435006488158237984014883E-2L,
-3.356232856677966691703904770937143483472E-1L,
-4.858192581793118040782557808823460276452E-1L,
-1.592399251246473643510898335746432479373E-1L,
};
#ifdef __STDC__
static const long double qs3[7] = {
#else
static long double qs3[7] = {
#endif
3.529139957987837084554591421329876744262E-5L,
2.973602667215766676998703687065066180115E-3L,
8.273534546240864308494062287908662592100E-2L,
9.613359842126507198241321110649974032726E-1L,
4.853923697093974370118387947065402707519E0L,
1.002671608961669247462020977417828796933E1L,
7.028927383922483728931327850683151410267E0L,
/* 1.000000000000000000000000000000000000000E0L, */
};
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
2 <= x <= 2.85711669921875
Peak relative error 6.9e-22 */
#ifdef __STDC__
static const long double qr2[7] = {
#else
static long double qr2[7] = {
#endif
-1.372751603025230017220666013816502528318E-4L,
-6.879190253347766576229143006767218972834E-3L,
-1.061253572090925414598304855316280077828E-1L,
-6.262164224345471241219408329354943337214E-1L,
-1.423149636514768476376254324731437473915E0L,
-1.087955310491078933531734062917489870754E0L,
-1.826821119773182847861406108689273719137E-1L,
};
#ifdef __STDC__
static const long double qs2[7] = {
#else
static long double qs2[7] = {
#endif
1.338768933634451601814048220627185324007E-3L,
7.071099998918497559736318523932241901810E-2L,
1.200511429784048632105295629933382142221E0L,
8.327301713640367079030141077172031825276E0L,
2.468479301872299311658145549931764426840E1L,
2.961179686096262083509383820557051621644E1L,
1.201402313144305153005639494661767354977E1L,
/* 1.000000000000000000000000000000000000000E0L, */
};
#ifdef __STDC__
static long double
qone (long double x)
#else
static long double
qone (x)
long double x;
#endif
{
#ifdef __STDC__
const long double *p, *q;
#else
long double *p, *q;
#endif
static long double s, r, z;
int32_t ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x4002) /* x >= 8 */
{
p = qr8;
q = qs8;
}
else
{
i1 = (ix << 16) | (i0 >> 16);
if (i1 >= 0x40019174) /* x >= 4.54541015625 */
{
p = qr5;
q = qs5;
}
else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
{
p = qr3;
q = qs3;
}
else if (ix >= 0x4000) /* x better be >= 2 */
{
p = qr2;
q = qs2;
}
}
z = one / (x * x);
r =
p[0] + z * (p[1] +
z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
s =
q[0] + z * (q[1] +
z * (q[2] +
z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
return (.375 + z * r / s) / x;
}
@@ -0,0 +1,386 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Modifications for long double are
Copyright (C) 2001 Stephen L. Moshier <[email protected]>
and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under
the following terms:
This 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.
This 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 this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
*
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with overflow signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
* Note 2. About jn(n,x), yn(n,x)
* For n=0, j0(x) is called,
* for n=1, j1(x) is called,
* for n<x, forward recursion us used starting
* from values of j0(x) and j1(x).
* for n>x, a continued fraction approximation to
* j(n,x)/j(n-1,x) is evaluated and then backward
* recursion is used starting from a supposed value
* for j(n,x). The resulting value of j(0,x) is
* compared with the actual value to correct the
* supposed value of j(n,x).
*
* yn(n,x) is similar in all respects, except
* that forward recursion is used for all
* values of n>1.
*
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double
#else
static long double
#endif
invsqrtpi = 5.64189583547756286948079e-1L, two = 2.0e0L, one = 1.0e0L;
#ifdef __STDC__
static const long double zero = 0.0L;
#else
static long double zero = 0.0L;
#endif
#ifdef __STDC__
long double
__ieee754_jnl (int n, long double x)
#else
long double
__ieee754_jnl (n, x)
int n;
long double x;
#endif
{
u_int32_t se, i0, i1;
int32_t i, ix, sgn;
long double a, b, temp, di;
long double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
/* if J(n,NaN) is NaN */
if ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0))
return x + x;
if (n < 0)
{
n = -n;
x = -x;
se ^= 0x8000;
}
if (n == 0)
return (__ieee754_j0l (x));
if (n == 1)
return (__ieee754_j1l (x));
sgn = (n & 1) & (se >> 15); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x);
if ((ix | i0 | i1) == 0 || ix >= 0x7fff) /* if x is 0 or inf */
b = zero;
else if ((long double) n <= x)
{
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x412D)
{ /* x > 2**302 */
/* ??? This might be a futile gesture.
If x exceeds X_TLOSS anyway, the wrapper function
will set the result to zero. */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
long double s;
long double c;
__sincosl (x, &s, &c);
switch (n & 3)
{
case 0:
temp = c + s;
break;
case 1:
temp = -c + s;
break;
case 2:
temp = -c - s;
break;
case 3:
temp = c - s;
break;
}
b = invsqrtpi * temp / __ieee754_sqrtl (x);
}
else
{
a = __ieee754_j0l (x);
b = __ieee754_j1l (x);
for (i = 1; i < n; i++)
{
temp = b;
b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
a = temp;
}
}
}
else
{
if (ix < 0x3fde)
{ /* x < 2**-33 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n >= 400) /* underflow, result < 10^-4952 */
b = zero;
else
{
temp = x * 0.5;
b = temp;
for (a = one, i = 2; i <= n; i++)
{
a *= (long double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b = b / a;
}
}
else
{
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
/* determine k */
long double t, v;
long double q0, q1, h, tmp;
int32_t k, m;
w = (n + n) / (long double) x;
h = 2.0L / (long double) x;
q0 = w;
z = w + h;
q1 = w * z - 1.0L;
k = 1;
while (q1 < 1.0e11L)
{
k += 1;
z += h;
tmp = z * q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t);
a = t;
b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = n;
v = two / x;
tmp = tmp * __ieee754_logl (fabsl (v * tmp));
if (tmp < 1.1356523406294143949491931077970765006170e+04L)
{
for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{
temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= two;
}
}
else
{
for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{
temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= two;
/* scale b to avoid spurious overflow */
if (b > 1e100L)
{
a /= b;
t /= b;
b = one;
}
}
}
b = (t * __ieee754_j0l (x) / b);
}
}
if (sgn == 1)
return -b;
else
return b;
}
#ifdef __STDC__
long double
__ieee754_ynl (int n, long double x)
#else
long double
__ieee754_ynl (n, x)
int n;
long double x;
#endif
{
u_int32_t se, i0, i1;
int32_t i, ix;
int32_t sign;
long double a, b, temp;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
/* if Y(n,NaN) is NaN */
if ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0))
return x + x;
if ((ix | i0 | i1) == 0)
return -HUGE_VALL + x; /* -inf and overflow exception. */
if (se & 0x8000)
return zero / (zero * x);
sign = 1;
if (n < 0)
{
n = -n;
sign = 1 - ((n & 1) << 1);
}
if (n == 0)
return (__ieee754_y0l (x));
if (n == 1)
return (sign * __ieee754_y1l (x));
if (ix == 0x7fff)
return zero;
if (ix >= 0x412D)
{ /* x > 2**302 */
/* ??? See comment above on the possible futility of this. */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
long double s;
long double c;
__sincosl (x, &s, &c);
switch (n & 3)
{
case 0:
temp = s - c;
break;
case 1:
temp = -s - c;
break;
case 2:
temp = -s + c;
break;
case 3:
temp = s + c;
break;
}
b = invsqrtpi * temp / __ieee754_sqrtl (x);
}
else
{
a = __ieee754_y0l (x);
b = __ieee754_y1l (x);
/* quit if b is -inf */
GET_LDOUBLE_WORDS (se, i0, i1, b);
for (i = 1; i < n && se != 0xffff; i++)
{
temp = b;
b = ((long double) (i + i) / x) * b - a;
GET_LDOUBLE_WORDS (se, i0, i1, b);
a = temp;
}
}
if (sign > 0)
return b;
else
return -b;
}
@@ -0,0 +1,183 @@
/* @(#)e_rem_pio2.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
#endif
/* __ieee754_rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2()
*/
#include "math.h"
#include "math_private.h"
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
#ifdef __STDC__
static const int32_t two_over_pi[] = {
#else
static int32_t two_over_pi[] = {
#endif
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
#ifdef __STDC__
static const int32_t npio2_hw[] = {
#else
static int32_t npio2_hw[] = {
#endif
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x404858EB, 0x404921FB,
};
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
#ifdef __STDC__
static const double
#else
static double
#endif
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
#ifdef __STDC__
int32_t __ieee754_rem_pio2(double x, double *y)
#else
int32_t __ieee754_rem_pio2(x,y)
double x,y[];
#endif
{
double z,w,t,r,fn;
double tx[3];
int32_t e0,i,j,nx,n,ix,hx;
u_int32_t low;
GET_HIGH_WORD(hx,x); /* high word of x */
ix = hx&0x7fffffff;
if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{y[0] = x; y[1] = 0; return 0;}
if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
if(hx>0) {
z = x - pio2_1;
if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z - pio2_1t;
y[1] = (z-y[0])-pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z-y[0])-pio2_2t;
}
return 1;
} else { /* negative x */
z = x + pio2_1;
if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z + pio2_1t;
y[1] = (z-y[0])+pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z-y[0])+pio2_2t;
}
return -1;
}
}
if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
t = fabs(x);
n = (int32_t) (t*invpio2+half);
fn = (double)n;
r = t-fn*pio2_1;
w = fn*pio2_1t; /* 1st round good to 85 bit */
if(n<32&&ix!=npio2_hw[n-1]) {
y[0] = r-w; /* quick check no cancellation */
} else {
u_int32_t high;
j = ix>>20;
y[0] = r-w;
GET_HIGH_WORD(high,y[0]);
i = j-((high>>20)&0x7ff);
if(i>16) { /* 2nd iteration needed, good to 118 */
t = r;
w = fn*pio2_2;
r = t-w;
w = fn*pio2_2t-((t-r)-w);
y[0] = r-w;
GET_HIGH_WORD(high,y[0]);
i = j-((high>>20)&0x7ff);
if(i>49) { /* 3rd iteration need, 151 bits acc */
t = r; /* will cover all possible cases */
w = fn*pio2_3;
r = t-w;
w = fn*pio2_3t-((t-r)-w);
y[0] = r-w;
}
}
}
y[1] = (r-y[0])-w;
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
else return n;
}
/*
* all other (large) arguments
*/
if(ix>=0x7ff00000) { /* x is inf or NaN */
y[0]=y[1]=x-x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD(low,x);
SET_LOW_WORD(z,low);
e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
for(i=0;i<2;i++) {
tx[i] = (double)((int32_t)(z));
z = (z-tx[i])*two24;
}
tx[2] = z;
nx = 3;
while(tx[nx-1]==zero) nx--; /* skip zero term */
n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
return n;
}
@@ -0,0 +1,70 @@
/* s_asinhl.c -- long double version of s_asinh.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* asinhl(x)
* Method :
* Based on
* asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ]
* we have
* asinhl(x) := x if 1+x*x=1,
* := signl(x)*(logl(x)+ln2)) for large |x|, else
* := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else
* := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double
#else
static long double
#endif
one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge= 1.000000000000000000e+4900L;
#ifdef __STDC__
long double __asinhl(long double x)
#else
long double __asinhl(x)
long double x;
#endif
{
long double t,w;
int32_t hx,ix;
GET_LDOUBLE_EXP(hx,x);
ix = hx&0x7fff;
if(ix==0x7fff) return x+x; /* x is inf or NaN */
if(ix< 0x3fde) { /* |x|<2**-34 */
if(huge+x>one) return x; /* return x inexact except 0 */
}
if(ix>0x4020) { /* |x| > 2**34 */
w = __ieee754_logl(fabsl(x))+ln2;
} else if (ix>0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
w = __ieee754_logl(2.0*t+one/(__ieee754_sqrtl(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =__log1pl(fabsl(x)+t/(one+__ieee754_sqrtl(one+t)));
}
if(hx&0x8000) return -w; else return w;
}
weak_alias (__asinhl, asinhl)
@@ -0,0 +1,71 @@
/* Compute cubic root of double value.
Copyright (C) 1997 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Dirk Alboth <[email protected]> and
Ulrich Drepper <[email protected]>, 1997.
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"
#define CBRT2 1.2599210498948731648 /* 2^(1/3) */
#define SQR_CBRT2 1.5874010519681994748 /* 2^(2/3) */
/* We don't use long double values here since U need not be computed
with full precision. */
static const double factor[5] =
{
1.0 / SQR_CBRT2,
1.0 / CBRT2,
1.0,
CBRT2,
SQR_CBRT2
};
static const long double third = 0.3333333333333333333333333L;
long double
__cbrtl (long double x)
{
long double xm, u;
int xe;
/* Reduce X. XM now is an range 1.0 to 0.5. */
xm = __frexpl (fabs (x), &xe);
/* If X is not finite or is null return it (with raising exceptions
if necessary.
Note: *Our* version of `frexp' sets XE to zero if the argument is
Inf or NaN. This is not portable but faster. */
if (xe == 0 && fpclassify (x) <= FP_ZERO)
return x + x;
u = (((-1.34661104733595206551E-1 * xm
+ 5.46646013663955245034E-1) * xm
- 9.54382247715094465250E-1) * xm
+ 1.13999833547172932737E0) * xm
+ 4.02389795645447521269E-1;
u *= factor[2 + xe % 3];
u = __ldexpl (x > 0.0 ? u : -u, xe / 3);
u -= (u - (x / (u * u))) * third;
u -= (u - (x / (u * u))) * third;
return u;
}
weak_alias (__cbrtl, cbrtl)
@@ -0,0 +1,40 @@
/* Return positive difference between arguments.
Copyright (C) 1997, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <[email protected]>, 1997.
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>
double
__fdim (double x, double y)
{
int clsx = fpclassify (x);
int clsy = fpclassify (y);
if (clsx == FP_NAN || clsy == FP_NAN
|| (y < 0 && clsx == FP_INFINITE && clsy == FP_INFINITE))
/* Raise invalid flag. */
return x - y;
return x <= y ? 0 : x - y;
}
weak_alias (__fdim, fdim)
#ifdef NO_LONG_DOUBLE
strong_alias (__fdim, __fdiml)
weak_alias (__fdim, fdiml)
#endif
@@ -0,0 +1,36 @@
/* Return positive difference between arguments.
Copyright (C) 1997, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <[email protected]>, 1997.
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>
float
__fdimf (float x, float y)
{
int clsx = fpclassify (x);
int clsy = fpclassify (y);
if (clsx == FP_NAN || clsy == FP_NAN
|| (y < 0 && clsx == FP_INFINITE && clsy == FP_INFINITE))
/* Raise invalid flag. */
return x - y;
return x <= y ? 0 : x - y;
}
weak_alias (__fdimf, fdimf)
@@ -0,0 +1,70 @@
/* s_frexpl.c -- long double version of s_frexp.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* for non-zero x
* x = frexpl(arg,&exp);
* return a long double fp quantity x such that 0.5 <= |x| <1.0
* and the corresponding binary exponent "exp". That is
* arg = x*2^exp.
* If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg
* with *exp=0.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double
#else
static long double
#endif
#if LDBL_MANT_DIG == 64
two65 = 3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */
#else
# error "Cannot handle this MANT_DIG"
#endif
#ifdef __STDC__
long double __frexpl(long double x, int *eptr)
#else
long double __frexpl(x, eptr)
long double x; int *eptr;
#endif
{
u_int32_t se, hx, ix, lx;
GET_LDOUBLE_WORDS(se,hx,lx,x);
ix = 0x7fff&se;
*eptr = 0;
if(ix==0x7fff||((ix|hx|lx)==0)) return x; /* 0,inf,nan */
if (ix==0x0000) { /* subnormal */
x *= two65;
GET_LDOUBLE_EXP(se,x);
ix = se&0x7fff;
*eptr = -65;
}
*eptr += ix-16382;
se = (se & 0x8000) | 0x3ffe;
SET_LDOUBLE_EXP(x,se);
return x;
}
weak_alias (__frexpl, frexpl)
@@ -0,0 +1,37 @@
/* s_ldexpl.c -- long double version of s_ldexp.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
#include <math.h>
#include "math_private.h"
#include <errno.h>
#ifdef __STDC__
long double __ldexpl(long double value, int exp)
#else
long double __ldexpl(value, exp)
long double value; int exp;
#endif
{
if(!__finitel(value)||value==0.0) return value;
value = __scalbnl(value,exp);
if(!__finitel(value)||value==0.0) __set_errno (ERANGE);
return value;
}
weak_alias (__ldexpl, ldexpl)
@@ -0,0 +1,105 @@
/* Adapted for use as nearbyint by Ulrich Drepper <[email protected]>. */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
#endif
/*
* rint(x)
* Return x rounded to integral value according to the prevailing
* rounding mode.
* Method:
* Using floating addition.
* Exception:
* Inexact flag raised if x not equal to rint(x).
*/
#include <fenv.h>
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const double
#else
static double
#endif
TWO52[2]={
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
};
#ifdef __STDC__
double __nearbyint(double x)
#else
double __nearbyint(x)
double x;
#endif
{
fenv_t env;
int32_t i0,j0,sx;
u_int32_t i,i1;
double w,t;
EXTRACT_WORDS(i0,i1,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) {
if(((i0&0x7fffffff)|i1)==0) return x;
i1 |= (i0&0x0fffff);
i0 &= 0xfffe0000;
i0 |= ((i1|-i1)>>12)&0x80000;
SET_HIGH_WORD(x,i0);
feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
fesetenv (&env);
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) return x; /* x is integral */
i>>=1;
if(((i0&i)|i1)!=0) {
if (j0==19)
i1 = 0x40000000;
else if (j0<18)
i0 = (i0&(~i))|((0x20000)>>j0);
else
{
i0 &= ~i;
i1 = 0x80000000;
}
}
}
} else if (j0>51) {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
i>>=1;
if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
}
INSERT_WORDS(x,i0,i1);
feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
fesetenv (&env);
return t;
}
weak_alias (__nearbyint, nearbyint)
#ifdef NO_LONG_DOUBLE
strong_alias (__nearbyint, __nearbyintl)
weak_alias (__nearbyint, nearbyintl)
#endif
@@ -0,0 +1,77 @@
/* s_rintf.c -- float version of s_rint.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
*/
/* Adapted for use as nearbyint by Ulrich Drepper <[email protected]>. */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <fenv.h>
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const float
#else
static float
#endif
TWO23[2]={
8.3886080000e+06, /* 0x4b000000 */
-8.3886080000e+06, /* 0xcb000000 */
};
#ifdef __STDC__
float __nearbyintf(float x)
#else
float __nearbyintf(x)
float x;
#endif
{
fenv_t env;
int32_t i0,j0,sx;
u_int32_t i,i1;
float w,t;
GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>23)&0xff)-0x7f;
if(j0<23) {
if(j0<0) {
if((i0&0x7fffffff)==0) return x;
i1 = (i0&0x07fffff);
i0 &= 0xfff00000;
i0 |= ((i1|-i1)>>9)&0x400000;
SET_FLOAT_WORD(x,i0);
feholdexcept (&env);
w = TWO23[sx]+x;
t = w-TWO23[sx];
fesetenv (&env);
GET_FLOAT_WORD(i0,t);
SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
i = (0x007fffff)>>j0;
if((i0&i)==0) return x; /* x is integral */
i>>=1;
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
}
} else {
if(j0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
feholdexcept (&env);
w = TWO23[sx]+x;
t = w-TWO23[sx];
fesetenv (&env);
return t;
}
weak_alias (__nearbyintf, nearbyintf)
@@ -0,0 +1,98 @@
/* @(#)s_rint.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
#endif
/*
* rint(x)
* Return x rounded to integral value according to the prevailing
* rounding mode.
* Method:
* Using floating addition.
* Exception:
* Inexact flag raised if x not equal to rint(x).
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const double
#else
static double
#endif
TWO52[2]={
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
};
#ifdef __STDC__
double __rint(double x)
#else
double __rint(x)
double x;
#endif
{
int32_t i0,j0,sx;
u_int32_t i,i1;
double w,t;
EXTRACT_WORDS(i0,i1,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) {
if(((i0&0x7fffffff)|i1)==0) return x;
i1 |= (i0&0x0fffff);
i0 &= 0xfffe0000;
i0 |= ((i1|-i1)>>12)&0x80000;
SET_HIGH_WORD(x,i0);
w = TWO52[sx]+x;
t = w-TWO52[sx];
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) return x; /* x is integral */
i>>=1;
if(((i0&i)|i1)!=0) {
if (j0==19)
i1 = 0x40000000;
else if (j0<18)
i0 = (i0&(~i))|((0x20000)>>j0);
else
{
i0 &= ~i;
i1 = 0x80000000;
}
}
}
} else if (j0>51) {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
i>>=1;
if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
}
INSERT_WORDS(x,i0,i1);
w = TWO52[sx]+x;
return w-TWO52[sx];
}
weak_alias (__rint, rint)
#ifdef NO_LONG_DOUBLE
strong_alias (__rint, __rintl)
weak_alias (__rint, rintl)
#endif
@@ -0,0 +1,72 @@
/* s_rintf.c -- float version of s_rint.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_rintf.c,v 1.4 1995/05/10 20:48:06 jtc Exp $";
#endif
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const float
#else
static float
#endif
TWO23[2]={
8.3886080000e+06, /* 0x4b000000 */
-8.3886080000e+06, /* 0xcb000000 */
};
#ifdef __STDC__
float __rintf(float x)
#else
float __rintf(x)
float x;
#endif
{
int32_t i0,j0,sx;
u_int32_t i,i1;
float w,t;
GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>23)&0xff)-0x7f;
if(j0<23) {
if(j0<0) {
if((i0&0x7fffffff)==0) return x;
i1 = (i0&0x07fffff);
i0 &= 0xfff00000;
i0 |= ((i1|-i1)>>9)&0x400000;
SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
t = w-TWO23[sx];
GET_FLOAT_WORD(i0,t);
SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
i = (0x007fffff)>>j0;
if((i0&i)==0) return x; /* x is integral */
i>>=1;
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
}
} else {
if(j0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
return w-TWO23[sx];
}
weak_alias (__rintf, rintf)
@@ -0,0 +1,70 @@
/* @(#)s_scalbn.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
#endif
/*
* scalbn (double x, int n)
* scalbn(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const double
#else
static double
#endif
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
#ifdef __STDC__
double __scalbln (double x, long int n)
#else
double __scalbln (x,n)
double x; long int n;
#endif
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
}
if (k==0x7ff) return x+x; /* NaN or Inf */
k = k+n;
if (n> 50000 || k > 0x7fe)
return huge*__copysign(huge,x); /* overflow */
if (n< -50000) return tiny*__copysign(tiny,x); /*underflow*/
if (k > 0) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
if (k <= -54)
return tiny*__copysign(tiny,x); /*underflow*/
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
}
weak_alias (__scalbln, scalbln)
#ifdef NO_LONG_DOUBLE
strong_alias (__scalbln, __scalblnl)
weak_alias (__scalbln, scalblnl)
#endif
@@ -0,0 +1,63 @@
/* s_scalbnf.c -- float version of s_scalbn.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_scalbnf.c,v 1.4 1995/05/10 20:48:10 jtc Exp $";
#endif
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const float
#else
static float
#endif
two25 = 3.355443200e+07, /* 0x4c000000 */
twom25 = 2.9802322388e-08, /* 0x33000000 */
huge = 1.0e+30,
tiny = 1.0e-30;
#ifdef __STDC__
float __scalblnf (float x, long int n)
#else
float __scalblnf (x,n)
float x; long int n;
#endif
{
int32_t k,ix;
GET_FLOAT_WORD(ix,x);
k = (ix&0x7f800000)>>23; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((ix&0x7fffffff)==0) return x; /* +-0 */
x *= two25;
GET_FLOAT_WORD(ix,x);
k = ((ix&0x7f800000)>>23) - 25;
}
if (k==0xff) return x+x; /* NaN or Inf */
k = k+n;
if (n> 50000 || k > 0xfe)
return huge*copysignf(huge,x); /* overflow */
if (n< -50000)
return tiny*copysignf(tiny,x); /*underflow*/
if (k > 0) /* normal result */
{SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
if (k <= -25)
return tiny*copysignf(tiny,x); /*underflow*/
k += 25; /* subnormal result */
SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
return x*twom25;
}
weak_alias (__scalblnf, scalblnf)
@@ -0,0 +1,71 @@
/* s_scalbnl.c -- long double version of s_scalbn.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* scalbnl (long double x, int n)
* scalbnl(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double
#else
static long double
#endif
two63 = 4.50359962737049600000e+15,
twom63 = 1.08420217248550443400e-19,
huge = 1.0e+4900L,
tiny = 1.0e-4900L;
#ifdef __STDC__
long double __scalblnl (long double x, long int n)
#else
long double __scalblnl (x,n)
long double x; long int n;
#endif
{
int32_t k,es,hx,lx;
GET_LDOUBLE_WORDS(es,hx,lx,x);
k = es&0x7fff; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two63;
GET_LDOUBLE_EXP(es,x);
k = (hx&0x7fff) - 63;
}
if (k==0x7fff) return x+x; /* NaN or Inf */
k = k+n;
if (n> 50000 || k > 0x7ffe)
return huge*__copysignl(huge,x); /* overflow */
if (n< -50000)
return tiny*__copysignl(tiny,x);
if (k > 0) /* normal result */
{SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;}
if (k <= -63)
return tiny*__copysignl(tiny,x); /*underflow*/
k += 63; /* subnormal result */
SET_LDOUBLE_EXP(x,(es&0x8000)|k);
return x*twom63;
}
weak_alias (__scalblnl, scalblnl)
@@ -0,0 +1,95 @@
/* s_tanhl.c -- long double version of s_tanh.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, [email protected].
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* tanhl(x)
* Return the Hyperbolic Tangent of x
*
* Method :
* x -x
* e - e
* 0. tanhl(x) is defined to be -----------
* x -x
* e + e
* 1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
* 2. 0 <= x <= 2**-55 : tanhl(x) := x*(one+x)
* -t
* 2**-55 < x <= 1 : tanhl(x) := -----; t = expm1l(-2x)
* t + 2
* 2
* 1 <= x <= 23.0 : tanhl(x) := 1- ----- ; t=expm1l(2x)
* t + 2
* 23.0 < x <= INF : tanhl(x) := 1.
*
* Special cases:
* tanhl(NaN) is NaN;
* only tanhl(0)=0 is exact for finite argument.
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
#else
static long double one=1.0, two=2.0, tiny = 1.0e-4900L;
#endif
#ifdef __STDC__
long double __tanhl(long double x)
#else
long double __tanhl(x)
long double x;
#endif
{
long double t,z;
int32_t se;
u_int32_t j0,j1,ix;
/* High word of |x|. */
GET_LDOUBLE_WORDS(se,j0,j1,x);
ix = se&0x7fff;
/* x is INF or NaN */
if(ix==0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se&0x8000) return one/x-one; /* tanhl(-inf)= -1; */
else return one/x+one; /* tanhl(+inf)=+1 */
}
/* |x| < 23 */
if (ix < 0x4003 || (ix == 0x4003 && j0 < 0xb8000000u)) {/* |x|<23 */
if ((ix|j0|j1) == 0)
return x; /* x == +- 0 */
if (ix<0x3fc8) /* |x|<2**-55 */
return x*(one+tiny); /* tanh(small) = small */
if (ix>=0x3fff) { /* |x|>=1 */
t = __expm1l(two*fabsl(x));
z = one - two/(t+two);
} else {
t = __expm1l(-two*fabsl(x));
z= -t/(t+two);
}
/* |x| > 23, return +-1 */
} else {
z = one - tiny; /* raised inexact flag */
}
return (se&0x8000)? -z: z;
}
weak_alias (__tanhl, tanhl)
@@ -0,0 +1,585 @@
/* These values are accurate to 52+12 bits when represented as
a double. */
static const double exp2_accuratetable[512] = {
0.707106781187802013759 /* 0x0.b504f333fb3f80007 */,
0.708064712808760599040 /* 0x0.b543baa0f71b38000 */,
0.709023942160304065938 /* 0x0.b58297d3a8d518002 */,
0.709984470998547667624 /* 0x0.b5c18ad39b4ba0001 */,
0.710946301084324217006 /* 0x0.b60093a85e8d30001 */,
0.711909434180505784637 /* 0x0.b63fb25984e628005 */,
0.712873872052760648733 /* 0x0.b67ee6eea3b5f8003 */,
0.713839616467838999908 /* 0x0.b6be316f518c98001 */,
0.714806669195984345523 /* 0x0.b6fd91e328d148007 */,
0.715775032009894562898 /* 0x0.b73d0851c69e20002 */,
0.716744706683768884058 /* 0x0.b77c94c2c9b3d0003 */,
0.717715694995770148178 /* 0x0.b7bc373dd52eb0003 */,
0.718687998724665488852 /* 0x0.b7fbefca8cd530004 */,
0.719661619652575468291 /* 0x0.b83bbe70981da8001 */,
0.720636559564428180758 /* 0x0.b87ba337a194b0006 */,
0.721612820246623098989 /* 0x0.b8bb9e27556508004 */,
0.722590403488338473025 /* 0x0.b8fbaf4762c798006 */,
0.723569311081411870036 /* 0x0.b93bd69f7be1d0000 */,
0.724549544820974333906 /* 0x0.b97c1437567828007 */,
0.725531106502312561633 /* 0x0.b9bc6816a87ae8002 */,
0.726513997924421062181 /* 0x0.b9fcd2452bee00000 */,
0.727498220889519875430 /* 0x0.ba3d52ca9e6148002 */,
0.728483777200401694265 /* 0x0.ba7de9aebe05c8003 */,
0.729470668664712662563 /* 0x0.babe96f94e62a8002 */,
0.730458897090379144517 /* 0x0.baff5ab2134df0004 */,
0.731448464287988597833 /* 0x0.bb4034e0d38ab0000 */,
0.732439372072965166897 /* 0x0.bb81258d5b2d60001 */,
0.733431622260458326859 /* 0x0.bbc22cbf75fd28001 */,
0.734425216668725511232 /* 0x0.bc034a7ef32c00001 */,
0.735420157118880535324 /* 0x0.bc447ed3a50fe0005 */,
0.736416445434497690674 /* 0x0.bc85c9c560b350001 */,
0.737414083433310718618 /* 0x0.bcc72b5bf4b4e0000 */,
0.738413072966152328496 /* 0x0.bd08a39f5417a8007 */,
0.739413415848264365956 /* 0x0.bd4a32974abcd0002 */,
0.740415113911250699637 /* 0x0.bd8bd84bb68300002 */,
0.741418168994518067562 /* 0x0.bdcd94c47ddd30003 */,
0.742422582936659858376 /* 0x0.be0f6809865968006 */,
0.743428357577745613238 /* 0x0.be515222b72530003 */,
0.744435494762383687126 /* 0x0.be935317fc6ba0002 */,
0.745443996335090397492 /* 0x0.bed56af1423de8001 */,
0.746453864145572798553 /* 0x0.bf1799b67a6248007 */,
0.747465100043933849969 /* 0x0.bf59df6f970e70002 */,
0.748477705883256683178 /* 0x0.bf9c3c248dbee8001 */,
0.749491683518965001732 /* 0x0.bfdeafdd568308000 */,
0.750507034813367890373 /* 0x0.c0213aa1f0fc38004 */,
0.751523761622240105153 /* 0x0.c063dc7a559ca0003 */,
0.752541865811731880422 /* 0x0.c0a6956e883ed8000 */,
0.753561349247157341600 /* 0x0.c0e965868bd220006 */,
0.754582213796583967110 /* 0x0.c12c4cca664cb8002 */,
0.755604461332336940791 /* 0x0.c16f4b42225350006 */,
0.756628093726406381068 /* 0x0.c1b260f5ca2c48002 */,
0.757653112855631305506 /* 0x0.c1f58ded6d72d8001 */,
0.758679520599333412360 /* 0x0.c238d2311e7d08001 */,
0.759707318837184453227 /* 0x0.c27c2dc8f00368005 */,
0.760736509456435783249 /* 0x0.c2bfa0bcfd1400000 */,
0.761767094336480043995 /* 0x0.c3032b155818d0000 */,
0.762799075372231349951 /* 0x0.c346ccda248cc0001 */,
0.763832454453522768941 /* 0x0.c38a8613805488005 */,
0.764867233473625618441 /* 0x0.c3ce56c98d1ca8005 */,
0.765903414329434539816 /* 0x0.c4123f04708d80002 */,
0.766940998920452976510 /* 0x0.c4563ecc532dc0001 */,
0.767979989148100838946 /* 0x0.c49a56295f9f88006 */,
0.769020386915772125040 /* 0x0.c4de8523c2b0a0001 */,
0.770062194131770905170 /* 0x0.c522cbc3ae94e0003 */,
0.771105412703856241146 /* 0x0.c5672a1154e6b8004 */,
0.772150044545352520777 /* 0x0.c5aba014ed5f18003 */,
0.773196091570364285606 /* 0x0.c5f02dd6b09288003 */,
0.774243555696622731700 /* 0x0.c634d35edb1260003 */,
0.775292438842697939641 /* 0x0.c67990b5aa5c18004 */,
0.776342742931542928455 /* 0x0.c6be65e360bed8000 */,
0.777394469888802008854 /* 0x0.c70352f0437f50004 */,
0.778447621641124243320 /* 0x0.c74857e498fd00006 */,
0.779502200118583399303 /* 0x0.c78d74c8ab5b60000 */,
0.780558207255445668515 /* 0x0.c7d2a9a4c959f8000 */,
0.781615644985491186966 /* 0x0.c817f681412f80002 */,
0.782674515247667956808 /* 0x0.c85d5b6666c150006 */,
0.783734819983036512536 /* 0x0.c8a2d85c904760003 */,
0.784796561133562109454 /* 0x0.c8e86d6c14f850002 */,
0.785859740645942328471 /* 0x0.c92e1a9d513ec8002 */,
0.786924360469767103536 /* 0x0.c973dff8a4b390007 */,
0.787990422552312885808 /* 0x0.c9b9bd866c6440007 */,
0.789057928854407064640 /* 0x0.c9ffb34f1444b0001 */,
0.790126881326406182996 /* 0x0.ca45c15afcc570001 */,
0.791197281930050233534 /* 0x0.ca8be7b292db38000 */,
0.792269132620954885659 /* 0x0.cad2265e3cbee8000 */,
0.793342435380726906957 /* 0x0.cb187d667d3d38006 */,
0.794417192158282659010 /* 0x0.cb5eecd3b33158006 */,
0.795493404931386649540 /* 0x0.cba574ae5d2e80001 */,
0.796571075671306805268 /* 0x0.cbec14fef2a348004 */,
0.797650206352955137846 /* 0x0.cc32cdcdef0000000 */,
0.798730798954342069432 /* 0x0.cc799f23d11d18000 */,
0.799812855456121796232 /* 0x0.ccc089091abb28004 */,
0.800896377841454287795 /* 0x0.cd078b86505c18003 */,
0.801981368096190028208 /* 0x0.cd4ea6a3f97720007 */,
0.803067828208752554378 /* 0x0.cd95da6aa057b8007 */,
0.804155760170129796375 /* 0x0.cddd26e2d21b28001 */,
0.805245165974338261710 /* 0x0.ce248c151f3330001 */,
0.806336047619038653883 /* 0x0.ce6c0a0a1c1350001 */,
0.807428407102107836855 /* 0x0.ceb3a0ca5d6be0006 */,
0.808522246427078927792 /* 0x0.cefb505e7e2550007 */,
0.809617567597010201484 /* 0x0.cf4318cf18a268002 */,
0.810714372621179513182 /* 0x0.cf8afa24ce1c98004 */,
0.811812663508675536069 /* 0x0.cfd2f4683f9810005 */,
0.812912442272482604912 /* 0x0.d01b07a2126188003 */,
0.814013710929394895825 /* 0x0.d06333daeff618001 */,
0.815116471495287542325 /* 0x0.d0ab791b80d028006 */,
0.816220725993571205593 /* 0x0.d0f3d76c75b330000 */,
0.817326476447408967199 /* 0x0.d13c4ed67f1cf8000 */,
0.818433724883006474832 /* 0x0.d184df6250e3b0001 */,
0.819542473330909460055 /* 0x0.d1cd8918a3a328004 */,
0.820652723822034690935 /* 0x0.d2164c02305fa0002 */,
0.821764478391968422618 /* 0x0.d25f2827b53fb0005 */,
0.822877739077315761840 /* 0x0.d2a81d91f188b8000 */,
0.823992507918612782109 /* 0x0.d2f12c49a8d290005 */,
0.825108786960634610365 /* 0x0.d33a5457a35e40003 */,
0.826226578247117093869 /* 0x0.d38395c4a84848007 */,
0.827345883828319528258 /* 0x0.d3ccf09985d958004 */,
0.828466705754248966560 /* 0x0.d41664df0a1320005 */,
0.829589046080638992111 /* 0x0.d45ff29e094330000 */,
0.830712906863802391671 /* 0x0.d4a999df585a20005 */,
0.831838290163696481037 /* 0x0.d4f35aabd04a60006 */,
0.832965198041969556729 /* 0x0.d53d350c4be258002 */,
0.834093632565442222342 /* 0x0.d5872909aba050007 */,
0.835223595802037643865 /* 0x0.d5d136acd138e8006 */,
0.836355089820669306292 /* 0x0.d61b5dfe9f7780004 */,
0.837488116698010487424 /* 0x0.d6659f0801afa8005 */,
0.838622678508982644113 /* 0x0.d6aff9d1e147d8004 */,
0.839758777333464490056 /* 0x0.d6fa6e652d19e0000 */,
0.840896415254110962690 /* 0x0.d744fccad70d00003 */,
0.842035594355151628676 /* 0x0.d78fa50bd2c3b0000 */,
0.843176316724478125433 /* 0x0.d7da673117e730007 */,
0.844318584453106590905 /* 0x0.d8254343a19038003 */,
0.845462399634695271912 /* 0x0.d870394c6dbf30003 */,
0.846607764365415071965 /* 0x0.d8bb49547d37c0004 */,
0.847754680744707056494 /* 0x0.d9067364d45608003 */,
0.848903150873708822763 /* 0x0.d951b7867953b0006 */,
0.850053176859071113491 /* 0x0.d99d15c2787a30006 */,
0.851204760807439786431 /* 0x0.d9e88e21de11a0003 */,
0.852357904828824897169 /* 0x0.da3420adba1508003 */,
0.853512611037803181642 /* 0x0.da7fcd6f2184d8005 */,
0.854668881550406100980 /* 0x0.dacb946f2afaf8000 */,
0.855826718478671755185 /* 0x0.db1775b6e8ad48000 */,
0.856986123964844970247 /* 0x0.db63714f8e0818006 */,
0.858147100114499461478 /* 0x0.dbaf87422625b8000 */,
0.859309649060962410524 /* 0x0.dbfbb797daa460002 */,
0.860473772936213743282 /* 0x0.dc480259d3a710001 */,
0.861639473872910177676 /* 0x0.dc9467913a0f48006 */,
0.862806754008130227807 /* 0x0.dce0e7473b9b28003 */,
0.863975615481124226159 /* 0x0.dd2d8185086c20006 */,
0.865146060433749419813 /* 0x0.dd7a3653d38168005 */,
0.866318091005120138881 /* 0x0.ddc705bcccd628000 */,
0.867491709362415264210 /* 0x0.de13efc9434100004 */,
0.868666917636779056818 /* 0x0.de60f4825df9b8005 */,
0.869843717989716047624 /* 0x0.deae13f16599c0003 */,
0.871022112578215268471 /* 0x0.defb4e1f9dc388002 */,
0.872202103559697183859 /* 0x0.df48a3164a92f0001 */,
0.873383693097737778847 /* 0x0.df9612deb6e878007 */,
0.874566883362160263365 /* 0x0.dfe39d82348310001 */,
0.875751676517234511901 /* 0x0.e031430a0f0688000 */,
0.876938074732511840819 /* 0x0.e07f037f97e548001 */,
0.878126080186539592654 /* 0x0.e0ccdeec2a75e0006 */,
0.879315695055312818168 /* 0x0.e11ad5591f4078001 */,
0.880506921518618312932 /* 0x0.e168e6cfd2f880004 */,
0.881699761760385225541 /* 0x0.e1b71359a6df60003 */,
0.882894217964411143207 /* 0x0.e2055afffc1178000 */,
0.884090292325693805080 /* 0x0.e253bdcc3ffbb8001 */,
0.885287987031581180559 /* 0x0.e2a23bc7d7a1d8002 */,
0.886487304278189114386 /* 0x0.e2f0d4fc31ab80004 */,
0.887688246263368285778 /* 0x0.e33f8972bea8a8005 */,
0.888890815189881999840 /* 0x0.e38e5934f49010007 */,
0.890095013257492739835 /* 0x0.e3dd444c460bd0007 */,
0.891300842677948068626 /* 0x0.e42c4ac232f380000 */,
0.892508305659222567226 /* 0x0.e47b6ca036f8b8005 */,
0.893717404414979710310 /* 0x0.e4caa9efd40e58002 */,
0.894928141160697743242 /* 0x0.e51a02ba8e2610007 */,
0.896140518115016826430 /* 0x0.e5697709ecab90000 */,
0.897354537501434679237 /* 0x0.e5b906e77c61d0006 */,
0.898570201543732793877 /* 0x0.e608b25cca5ba8005 */,
0.899787512470129891014 /* 0x0.e6587973688ce8002 */,
0.901006472512270728537 /* 0x0.e6a85c34ecadb8000 */,
0.902227083902570559127 /* 0x0.e6f85aaaed4f20006 */,
0.903449348881299796343 /* 0x0.e74874df09a530003 */,
0.904673269686823378091 /* 0x0.e798aadadecba0007 */,
0.905898848559668845585 /* 0x0.e7e8fca80c3ee0001 */,
0.907126087750156795426 /* 0x0.e8396a503c3fe0005 */,
0.908354989505901100354 /* 0x0.e889f3dd1615b0002 */,
0.909585556079328783087 /* 0x0.e8da9958465228007 */,
0.910817789726044213523 /* 0x0.e92b5acb7d0578001 */,
0.912051692703457872481 /* 0x0.e97c38406c3c30003 */,
0.913287267274154990210 /* 0x0.e9cd31c0cbb370001 */,
0.914524515702244578108 /* 0x0.ea1e475654d540000 */,
0.915763440256158633982 /* 0x0.ea6f790ac5cc78001 */,
0.917004043205012497909 /* 0x0.eac0c6e7dd8448007 */,
0.918246326823137892807 /* 0x0.eb1230f760a428007 */,
0.919490293387826285200 /* 0x0.eb63b7431714a8007 */,
0.920735945178816406225 /* 0x0.ebb559d4cb6f30007 */,
0.921983284479243714322 /* 0x0.ec0718b64c0940002 */,
0.923232313574974705626 /* 0x0.ec58f3f16a3910002 */,
0.924483034755387955725 /* 0x0.ecaaeb8ffb3168005 */,
0.925735450311948926408 /* 0x0.ecfcff9bd67078000 */,
0.926989562542820610982 /* 0x0.ed4f301edad1a0007 */,
0.928245373740515189457 /* 0x0.eda17d22e0f9b0001 */,
0.929502886213858126045 /* 0x0.edf3e6b1d37d40001 */,
0.930762102264245716494 /* 0x0.ee466cd594c5c8005 */,
0.932023024199046146183 /* 0x0.ee990f980dcdb0005 */,
0.933285654329454095216 /* 0x0.eeebcf032bc470007 */,
0.934549994971191289044 /* 0x0.ef3eab20e0d3c0001 */,
0.935816048439005676599 /* 0x0.ef91a3fb1e1340004 */,
0.937083817055075818404 /* 0x0.efe4b99bdcc618006 */,
0.938353303143720007819 /* 0x0.f037ec0d1889b8000 */,
0.939624509028518128972 /* 0x0.f08b3b58cc2bb8006 */,
0.940897437041863904384 /* 0x0.f0dea788fc2a90000 */,
0.942172089516254085427 /* 0x0.f13230a7ad21b8003 */,
0.943448468787511540534 /* 0x0.f185d6bee754e0006 */,
0.944726577195256100890 /* 0x0.f1d999d8b73478005 */,
0.946006417082291717338 /* 0x0.f22d79ff2cb130000 */,
0.947287990793413858827 /* 0x0.f281773c59ec48007 */,
0.948571300678290207925 /* 0x0.f2d5919a566268001 */,
0.949856349088629370320 /* 0x0.f329c9233bceb0001 */,
0.951143138379053731954 /* 0x0.f37e1de1272068002 */,
0.952431670908847949364 /* 0x0.f3d28fde3a6728006 */,
0.953721949039916472305 /* 0x0.f4271f249a93f0001 */,
0.955013975135367898520 /* 0x0.f47bcbbe6deab0001 */,
0.956307751564417496418 /* 0x0.f4d095b5e16638004 */,
0.957603280698967163097 /* 0x0.f5257d1524f590006 */,
0.958900564911197350604 /* 0x0.f57a81e668d628000 */,
0.960199606581278120057 /* 0x0.f5cfa433e60e50007 */,
0.961500408088936442422 /* 0x0.f624e407d527a0007 */,
0.962802971817578789903 /* 0x0.f67a416c72b760006 */,
0.964107300155846558292 /* 0x0.f6cfbc6c011458004 */,
0.965413395493874504368 /* 0x0.f7255510c439a8002 */,
0.966721260225105960572 /* 0x0.f77b0b6503c5b8006 */,
0.968030896745834645873 /* 0x0.f7d0df730a7940005 */,
0.969342307458006424716 /* 0x0.f826d145294be8003 */,
0.970655494764855020231 /* 0x0.f87ce0e5b29fd8000 */,
0.971970461071268720958 /* 0x0.f8d30e5efaa8f0004 */,
0.973287208789983648852 /* 0x0.f92959bb5e3c08001 */,
0.974605740331924708124 /* 0x0.f97fc305383028004 */,
0.975926058115625383329 /* 0x0.f9d64a46ebb9f8004 */,
0.977248164559556209435 /* 0x0.fa2cef8adbfc68004 */,
0.978572062087848637573 /* 0x0.fa83b2db7253d0007 */,
0.979897753126343307191 /* 0x0.fada944319fda0005 */,
0.981225240104636631254 /* 0x0.fb3193cc425870002 */,
0.982554525455618277276 /* 0x0.fb88b1815e61d0003 */,
0.983885611617111077747 /* 0x0.fbdfed6ce683e0007 */,
0.985218501026348891812 /* 0x0.fc3747995282f8006 */,
0.986553196127724962867 /* 0x0.fc8ec0112202a0005 */,
0.987889699367056062238 /* 0x0.fce656ded63710002 */,
0.989228013193998778636 /* 0x0.fd3e0c0cf48d50005 */,
0.990568140061241164686 /* 0x0.fd95dfa605c7b0003 */,
0.991910082424819927754 /* 0x0.fdedd1b4965710004 */,
0.993253842749249660216 /* 0x0.fe45e2433bfea0000 */,
0.994599423484053835071 /* 0x0.fe9e115c7c05f0005 */,
0.995946827107488830167 /* 0x0.fef65f0afb4c28006 */,
0.997296056085008264529 /* 0x0.ff4ecb59509cc8001 */,
0.998647112892057764479 /* 0x0.ffa756521dbfd0007 */,
1.000000000000000000000 /* 0x1.00000000000000000 */,
1.001354719891689004659 /* 0x1.0058c86da14aa0005 */,
1.002711275050312211844 /* 0x1.00b1afa5abead0003 */,
1.004069667960743483835 /* 0x1.010ab5b2cc0660009 */,
1.005429901112333324093 /* 0x1.0163da9fb2af30008 */,
1.006791976999887428009 /* 0x1.01bd1e7716f6a0008 */,
1.008155898118476168101 /* 0x1.02168143b03890006 */,
1.009521666967782227439 /* 0x1.027003103ae320002 */,
1.010889286051850133326 /* 0x1.02c9a3e7783030002 */,
1.012258757875921233497 /* 0x1.032363d42aaa8000e */,
1.013630084952214405194 /* 0x1.037d42e11c88d0000 */,
1.015003269791313389451 /* 0x1.03d741191635a0001 */,
1.016378314911229763267 /* 0x1.04315e86e84630008 */,
1.017755222831652872635 /* 0x1.048b9b35652800002 */,
1.019133996077934645224 /* 0x1.04e5f72f65827000b */,
1.020514637175266248212 /* 0x1.0540727fc1cfa0006 */,
1.021897148653734488385 /* 0x1.059b0d3157ebb0002 */,
1.023281533050062419584 /* 0x1.05f5c74f0cfeb0002 */,
1.024667792897328677539 /* 0x1.0650a0e3c22ee0003 */,
1.026055930738840826806 /* 0x1.06ab99fa63e1b0008 */,
1.027445949118511947550 /* 0x1.0706b29ddf2700009 */,
1.028837850584049418178 /* 0x1.0761ead9253ab0009 */,
1.030231637685799839262 /* 0x1.07bd42b72a3f80008 */,
1.031627312979383592802 /* 0x1.0818ba42e824a000c */,
1.033024879021186448496 /* 0x1.0874518759b0b0008 */,
1.034424338374263729911 /* 0x1.08d0088f80ffa0006 */,
1.035825693601787333992 /* 0x1.092bdf66604e30005 */,
1.037228947273990842283 /* 0x1.0987d617019cd000a */,
1.038634101961269928846 /* 0x1.09e3ecac6f199000f */,
1.040041160239590700707 /* 0x1.0a402331b91270002 */,
1.041450124688240164200 /* 0x1.0a9c79b1f37c3000b */,
1.042860997889083929381 /* 0x1.0af8f038352160000 */,
1.044273782427270314011 /* 0x1.0b5586cf986890006 */,
1.045688480893644856116 /* 0x1.0bb23d833dfbf0006 */,
1.047105095879385272564 /* 0x1.0c0f145e46e330007 */,
1.048523629981608529302 /* 0x1.0c6c0b6bdaadc000f */,
1.049944085800634585634 /* 0x1.0cc922b72470a000f */,
1.051366465939483019223 /* 0x1.0d265a4b5238b0007 */,
1.052790773004648849929 /* 0x1.0d83b23395e510002 */,
1.054217009607077093512 /* 0x1.0de12a7b263970006 */,
1.055645178360430591625 /* 0x1.0e3ec32d3cf680000 */,
1.057075281882416506511 /* 0x1.0e9c7c55184f5000e */,
1.058507322794714378170 /* 0x1.0efa55fdfad51000a */,
1.059941303721639416236 /* 0x1.0f58503329fed0003 */,
1.061377227289284297385 /* 0x1.0fb66affed37f0000 */,
1.062815096132297298980 /* 0x1.1014a66f95540000c */,
1.064254912884593951029 /* 0x1.1073028d725850007 */,
1.065696680185205469411 /* 0x1.10d17f64d9ea2000b */,
1.067140400676658718053 /* 0x1.11301d012586a0007 */,
1.068586077004890055886 /* 0x1.118edb6db26ab0003 */,
1.070033711820396415998 /* 0x1.11edbab5e2d6e000b */,
1.071483307775789262099 /* 0x1.124cbae51b5ef0001 */,
1.072934867526001312439 /* 0x1.12abdc06c3240000c */,
1.074388393734249103080 /* 0x1.130b1e264a62e0005 */,
1.075843889063253344684 /* 0x1.136a814f20ccd0003 */,
1.077301356179926061823 /* 0x1.13ca058cbaaed000b */,
1.078760797756675327056 /* 0x1.1429aaea9260e000e */,
1.080222216468626150775 /* 0x1.148971742537c0009 */,
1.081685614993597610617 /* 0x1.14e95934f37e8000b */,
1.083150996013011013776 /* 0x1.1549623881762000d */,
1.084618362213087383633 /* 0x1.15a98c8a58a6a000b */,
1.086087716284427351384 /* 0x1.1609d8360768c0008 */,
1.087559060917626885283 /* 0x1.166a45471c13f0008 */,
1.089032398810997337465 /* 0x1.16cad3c92d7b50009 */,
1.090507732647478578212 /* 0x1.172b83c7c18b5000f */,
1.091985065182095926460 /* 0x1.178c554ead72a000c */,
1.093464399073070136880 /* 0x1.17ed48695befe000c */,
1.094945737045367906172 /* 0x1.184e5d23812500007 */,
1.096429081816546080591 /* 0x1.18af9388c90e40005 */,
1.097914436104650892651 /* 0x1.1910eba4e031a0001 */,
1.099401802629782043408 /* 0x1.19726583755720003 */,
1.100891184121537858001 /* 0x1.19d4013041b860007 */,
1.102382583308144647940 /* 0x1.1a35beb6fd0cd0007 */,
1.103876002922312915544 /* 0x1.1a979e2363fa10000 */,
1.105371445702084232160 /* 0x1.1af99f8139025000e */,
1.106868914387219016199 /* 0x1.1b5bc2dc408b9000e */,
1.108368411723785085252 /* 0x1.1bbe084045eb30002 */,
1.109869940458469095340 /* 0x1.1c206fb91524c000e */,
1.111373503344554869449 /* 0x1.1c82f952817cc0001 */,
1.112879103137133007859 /* 0x1.1ce5a51860344000f */,
1.114386742595953938610 /* 0x1.1d4873168babf000e */,
1.115896424484008608911 /* 0x1.1dab6358e1d4a000f */,
1.117408151567338414664 /* 0x1.1e0e75eb43f9c000c */,
1.118921926613465345265 /* 0x1.1e71aad995078000f */,
1.120437752409564780022 /* 0x1.1ed5022fcd8600003 */,
1.121955631720569668277 /* 0x1.1f387bf9cd88b0000 */,
1.123475567332998359439 /* 0x1.1f9c18438cdec000a */,
1.124997562033035469759 /* 0x1.1fffd71902f970002 */,
1.126521618608448571713 /* 0x1.2063b88629079000e */,
1.128047739853580200284 /* 0x1.20c7bc96ff72a0002 */,
1.129575928566289189112 /* 0x1.212be3578a81e0006 */,
1.131106187546149888259 /* 0x1.21902cd3d05f70007 */,
1.132638519598779369743 /* 0x1.21f49917ddda5000c */,
1.134172927531616359481 /* 0x1.2259282fc1c24000e */,
1.135709414157753949251 /* 0x1.22bdda27911e90007 */,
1.137247982292643566662 /* 0x1.2322af0b638e60007 */,
1.138788634756517259562 /* 0x1.2387a6e755f270000 */,
1.140331374372893558110 /* 0x1.23ecc1c788c890006 */,
1.141876203969685699176 /* 0x1.2451ffb821639000c */,
1.143423126377846266197 /* 0x1.24b760c5486dc0009 */,
1.144972144431494420774 /* 0x1.251ce4fb2a0cc0005 */,
1.146523260971646252006 /* 0x1.25828c65f9fb8000d */,
1.148076478839068270690 /* 0x1.25e85711ebaeb0000 */,
1.149631800883562204903 /* 0x1.264e450b3c8a30008 */,
1.151189229953253789786 /* 0x1.26b4565e281a20003 */,
1.152748768902654319399 /* 0x1.271a8b16f0f000002 */,
1.154310420590433317050 /* 0x1.2780e341de2fc0001 */,
1.155874187878668246681 /* 0x1.27e75eeb3abc90007 */,
1.157440073633736243899 /* 0x1.284dfe1f5633e000a */,
1.159008080725518974322 /* 0x1.28b4c0ea840d90001 */,
1.160578212048386514965 /* 0x1.291ba75932ae60000 */,
1.162150470417516290340 /* 0x1.2982b177796850008 */,
1.163724858777502646494 /* 0x1.29e9df51fdd900001 */,
1.165301379991388053320 /* 0x1.2a5130f50bf34000e */,
1.166880036952526289469 /* 0x1.2ab8a66d10fdc0008 */,
1.168460832550151540268 /* 0x1.2b203fc675b7a000a */,
1.170043769683112966389 /* 0x1.2b87fd0dad7260008 */,
1.171628851252754177681 /* 0x1.2befde4f2e3da000d */,
1.173216080163546060084 /* 0x1.2c57e397719940002 */,
1.174805459325657830448 /* 0x1.2cc00cf2f7491000c */,
1.176396991650083379037 /* 0x1.2d285a6e3ff90000b */,
1.177990680055698513602 /* 0x1.2d90cc15d4ff90005 */,
1.179586527463262646306 /* 0x1.2df961f641c57000c */,
1.181184536796979545103 /* 0x1.2e621c1c157cd000d */,
1.182784710984701836994 /* 0x1.2ecafa93e35af0004 */,
1.184387052960675701386 /* 0x1.2f33fd6a459cb0000 */,
1.185991565661414393112 /* 0x1.2f9d24abd8fd1000e */,
1.187598252026902612178 /* 0x1.300670653e083000a */,
1.189207115003001469262 /* 0x1.306fe0a31bc040008 */,
1.190818157535919796833 /* 0x1.30d9757219895000e */,
1.192431382587621380206 /* 0x1.31432edef01a1000f */,
1.194046793097208292195 /* 0x1.31ad0cf63f0630008 */,
1.195664392040319823392 /* 0x1.32170fc4ce0db000c */,
1.197284182375793593084 /* 0x1.32813757527750005 */,
1.198906167074650808198 /* 0x1.32eb83ba8eef3000f */,
1.200530349107333139048 /* 0x1.3355f4fb457e5000d */,
1.202156731453099647353 /* 0x1.33c08b2641df9000c */,
1.203785317090505513368 /* 0x1.342b46484f07b0005 */,
1.205416109005122526928 /* 0x1.3496266e3fa270005 */,
1.207049110184904572310 /* 0x1.35012ba4e8fa10000 */,
1.208684323627194912036 /* 0x1.356c55f92aabb0004 */,
1.210321752322854882437 /* 0x1.35d7a577dd33f0004 */,
1.211961399276747286580 /* 0x1.36431a2de8748000d */,
1.213603267492579629347 /* 0x1.36aeb4283309e000c */,
1.215247359985374142610 /* 0x1.371a7373b00160000 */,
1.216893679753690671322 /* 0x1.3786581d404e90000 */,
1.218542229828181611183 /* 0x1.37f26231e82e4000c */,
1.220193013225231215567 /* 0x1.385e91be9c2d20002 */,
1.221846032973555429280 /* 0x1.38cae6d05e66f0000 */,
1.223501292099485437962 /* 0x1.393761742e5830001 */,
1.225158793636904830441 /* 0x1.39a401b713cb3000e */,
1.226818540625497444577 /* 0x1.3a10c7a61ceae0007 */,
1.228480536107136034131 /* 0x1.3a7db34e5a4a50003 */,
1.230144783126481566885 /* 0x1.3aeac4bcdf8d60001 */,
1.231811284734168454619 /* 0x1.3b57fbfec6e950008 */,
1.233480043984379381835 /* 0x1.3bc559212e7a2000f */,
1.235151063936380300149 /* 0x1.3c32dc3139f2a0004 */,
1.236824347652524913647 /* 0x1.3ca0853c106ac000e */,
1.238499898199571624970 /* 0x1.3d0e544eddd240003 */,
1.240177718649636107175 /* 0x1.3d7c4976d3fcd0000 */,
1.241857812073360767273 /* 0x1.3dea64c1231f70004 */,
1.243540181554270152039 /* 0x1.3e58a63b099920005 */,
1.245224830175077013244 /* 0x1.3ec70df1c4e46000e */,
1.246911761022835740725 /* 0x1.3f359bf29741c000e */,
1.248600977188942806639 /* 0x1.3fa4504ac7b800009 */,
1.250292481770148400634 /* 0x1.40132b07a330d000a */,
1.251986277866492969263 /* 0x1.40822c367a340000b */,
1.253682368581898742876 /* 0x1.40f153e4a18e0000d */,
1.255380757024939564249 /* 0x1.4160a21f73289000d */,
1.257081446308726757662 /* 0x1.41d016f44deaa000c */,
1.258784439550028944083 /* 0x1.423fb27094c090008 */,
1.260489739869405489991 /* 0x1.42af74a1aec1c0006 */,
1.262197350394008266193 /* 0x1.431f5d950a453000c */,
1.263907274252603851764 /* 0x1.438f6d58176860004 */,
1.265619514578811388761 /* 0x1.43ffa3f84b9eb000d */,
1.267334074511444086425 /* 0x1.44700183221180008 */,
1.269050957191869555296 /* 0x1.44e0860618b930006 */,
1.270770165768063009230 /* 0x1.4551318eb4d20000e */,
1.272491703389059036805 /* 0x1.45c2042a7cc26000b */,
1.274215573211836316547 /* 0x1.4632fde6ffacd000d */,
1.275941778396075143580 /* 0x1.46a41ed1cfac40001 */,
1.277670322103555911043 /* 0x1.471566f8812ac0000 */,
1.279401207505722393185 /* 0x1.4786d668b33260005 */,
1.281134437771823675369 /* 0x1.47f86d3002637000a */,
1.282870016078732078362 /* 0x1.486a2b5c13c00000e */,
1.284607945607987078432 /* 0x1.48dc10fa916bd0004 */,
1.286348229545787758022 /* 0x1.494e1e192aaa30007 */,
1.288090871080605159846 /* 0x1.49c052c5913df000c */,
1.289835873406902644341 /* 0x1.4a32af0d7d8090002 */,
1.291583239722392528754 /* 0x1.4aa532feab5e10002 */,
1.293332973229098792374 /* 0x1.4b17dea6db8010008 */,
1.295085077135345708087 /* 0x1.4b8ab213d57d9000d */,
1.296839554650994097442 /* 0x1.4bfdad53629e10003 */,
1.298596408992440220988 /* 0x1.4c70d0735358a000d */,
1.300355643380135983739 /* 0x1.4ce41b817c99e0001 */,
1.302117261036232376282 /* 0x1.4d578e8bb52cb0003 */,
1.303881265192249561154 /* 0x1.4dcb299fde2920008 */,
1.305647659079073541490 /* 0x1.4e3eeccbd7f4c0003 */,
1.307416445934474813521 /* 0x1.4eb2d81d8a86f000b */,
1.309187629001237640529 /* 0x1.4f26eba2e35a5000e */,
1.310961211525240921493 /* 0x1.4f9b2769d35090009 */,
1.312737196755087820678 /* 0x1.500f8b804e4a30000 */,
1.314515587949291131086 /* 0x1.508417f4530d00009 */,
1.316296388365203462468 /* 0x1.50f8ccd3df1840003 */,
1.318079601265708777911 /* 0x1.516daa2cf60020002 */,
1.319865229921343141607 /* 0x1.51e2b00da3c2b0007 */,
1.321653277603506371251 /* 0x1.5257de83f5512000d */,
1.323443747588034513690 /* 0x1.52cd359dfc7d5000e */,
1.325236643161341820781 /* 0x1.5342b569d6baa000f */,
1.327031967602244177939 /* 0x1.53b85df59921b0000 */,
1.328829724206201046165 /* 0x1.542e2f4f6b17e0006 */,
1.330629916266568235675 /* 0x1.54a4298571b27000e */,
1.332432547083447937938 /* 0x1.551a4ca5d97190009 */,
1.334237619959296017340 /* 0x1.559098bed16bf0008 */,
1.336045138203900251029 /* 0x1.56070dde90c800000 */,
1.337855105129210686631 /* 0x1.567dac13510cd0009 */,
1.339667524053662184301 /* 0x1.56f4736b52e2c000c */,
1.341482398296830025383 /* 0x1.576b63f4d8333000f */,
1.343299731186792467254 /* 0x1.57e27dbe2c40e0003 */,
1.345119526053918823702 /* 0x1.5859c0d59cd37000f */,
1.346941786233264881662 /* 0x1.58d12d497cd9a0005 */,
1.348766515064854010261 /* 0x1.5948c32824b87000c */,
1.350593715891792223641 /* 0x1.59c0827ff03890007 */,
1.352423392064920459908 /* 0x1.5a386b5f43a3e0006 */,
1.354255546937278120764 /* 0x1.5ab07dd485af1000c */,
1.356090183865519494030 /* 0x1.5b28b9ee21085000f */,
1.357927306213322804534 /* 0x1.5ba11fba8816e000b */,
1.359766917346459269620 /* 0x1.5c19af482f8f2000f */,
1.361609020638567812980 /* 0x1.5c9268a594cc00004 */,
1.363453619463660171403 /* 0x1.5d0b4be135916000c */,
1.365300717204201985683 /* 0x1.5d84590998eeb0005 */,
1.367150317245710233754 /* 0x1.5dfd902d494e40001 */,
1.369002422974674892971 /* 0x1.5e76f15ad22c40008 */,
1.370857037789471544224 /* 0x1.5ef07ca0cc166000b */,
1.372714165088220639199 /* 0x1.5f6a320dcf5280006 */,
1.374573808273481745378 /* 0x1.5fe411b0790800009 */,
1.376435970755022220096 /* 0x1.605e1b976e4b1000e */,
1.378300655944092456600 /* 0x1.60d84fd155d15000e */,
1.380167867259843417228 /* 0x1.6152ae6cdf0030003 */,
1.382037608124419003675 /* 0x1.61cd3778bc879000d */,
1.383909881963391264069 /* 0x1.6247eb03a4dc40009 */,
1.385784692209972801544 /* 0x1.62c2c91c56d9b0002 */,
1.387662042298923203992 /* 0x1.633dd1d1930ec0001 */,
1.389541935670444372533 /* 0x1.63b90532200630004 */,
1.391424375772021271329 /* 0x1.6434634ccc4cc0007 */,
1.393309366052102982208 /* 0x1.64afec30677e90008 */,
1.395196909966106124701 /* 0x1.652b9febc8e0f000d */,
1.397087010973788290271 /* 0x1.65a77e8dcc7f10004 */,
1.398979672539331309267 /* 0x1.66238825534170000 */,
1.400874898129892187656 /* 0x1.669fbcc1415600008 */,
1.402772691220124823310 /* 0x1.671c1c708328e000a */,
1.404673055288671035301 /* 0x1.6798a7420988b000d */,
1.406575993818903302975 /* 0x1.68155d44ca77a000f */,
1.408481510297352468121 /* 0x1.68923e87bf70e000a */,
1.410389608216942924956 /* 0x1.690f4b19e8f74000c */,
1.412300291075172076232 /* 0x1.698c830a4c94c0008 */
};
#define S (1.0/4503599627370496.0) /* 2^-52 */
static const float exp2_deltatable[512] = {
11527*S, -963*S, 884*S, -781*S, -2363*S, -3441*S, 123*S, 526*S,
-6*S, 1254*S, -1138*S, 1519*S, 1576*S, -65*S, 1040*S, 793*S,
-1662*S, -5063*S, -387*S, 968*S, -941*S, 984*S, -2856*S, -545*S,
495*S, -5246*S, -2109*S, 1281*S, 2075*S, 909*S, -1642*S,-78233*S,
-31653*S, -265*S, 130*S, 430*S, 2482*S, -742*S, 1616*S, -2213*S,
-519*S, 20*S, -3134*S,-13981*S, 1343*S, -1740*S, 247*S, 1679*S,
-1097*S, 3131*S, 871*S, -1480*S, 1936*S, -1827*S, 17325*S, 528*S,
-322*S, 1404*S, -152*S, -1845*S, -212*S, 2639*S, -476*S, 2960*S,
-962*S, -1012*S, -1231*S, 3030*S, 1659*S, -486*S, 2154*S, 1728*S,
-2793*S, 699*S, -1560*S, -2125*S, 2156*S, 142*S, -1888*S, 4426*S,
-13443*S, 1970*S, -50*S, 1771*S,-43399*S, 4979*S, -2448*S, -370*S,
1414*S, 1075*S, 232*S, 206*S, 873*S, 2141*S, 2970*S, 1279*S,
-2331*S, 336*S, -2595*S, 753*S, -3384*S, -616*S, 89*S, -818*S,
5755*S, -241*S, -528*S, -661*S, -3777*S, -354*S, 250*S, 3881*S,
2632*S, -2131*S, 2565*S, -316*S, 1746*S, -2541*S, -1324*S, -50*S,
2564*S, -782*S, 1176*S, 6452*S, -1002*S, 1288*S, 336*S, -185*S,
3063*S, 3784*S, 2169*S, 686*S, 328*S, -400*S, 312*S, -4517*S,
-1457*S, 1046*S, -1530*S, -685*S, 1328*S,-49815*S, -895*S, 1063*S,
-2091*S, -672*S, -1710*S, -665*S, 1545*S, 1819*S,-45265*S, 3548*S,
-554*S, -568*S, 4752*S, -1907*S,-13738*S, 675*S, 9611*S, -1115*S,
-815*S, 408*S, -1281*S, -937*S,-16376*S, -4772*S, -1440*S, 992*S,
788*S, 10364*S, -1602*S, -661*S, -1783*S, -265*S, -20*S, -3781*S,
-861*S, -345*S, -994*S, 1364*S, -5339*S, 1620*S, 9390*S, -1066*S,
-305*S, -170*S, 175*S, 2461*S, -490*S, -769*S, -1450*S, 3315*S,
2418*S, -45*S, -852*S, -1295*S, -488*S, -96*S, 1142*S, -2639*S,
7905*S, -9306*S, -3859*S, 760*S, 1057*S, -1570*S, 3977*S, 209*S,
-514*S, 7151*S, 1646*S, 627*S, 599*S, -774*S, -1468*S, 633*S,
-473*S, 851*S, 2406*S, 143*S, 74*S, 4260*S, 1177*S, -913*S,
2670*S, -3298*S, -1662*S, -120*S, -3264*S, -2148*S, 410*S, 2078*S,
-2098*S, -926*S, 3580*S, -1289*S, 2450*S, -1158*S, 907*S, -590*S,
986*S, 1801*S, 1145*S, -1677*S, 3455*S, 956*S, 710*S, 144*S,
153*S, -255*S, -1898*S, 28102*S, 2748*S, 1194*S, -3009*S, 7076*S,
0*S, -2720*S, 711*S, 1225*S, -3034*S, -473*S, 378*S, -1046*S,
962*S, -2006*S, 4647*S, 3206*S, 1769*S, -2665*S, 1254*S, 2025*S,
-2430*S, 6193*S, 1224*S, -856*S, -1592*S, -325*S, -1521*S, 1827*S,
-264*S, 2403*S, -1065*S, 967*S, -681*S, -2106*S, -474*S, 1333*S,
-893*S, 2296*S, 592*S, -1220*S, -326*S, 990*S, 139*S, 206*S,
-779*S, -1683*S, 1238*S, 6098*S, 136*S, 1197*S, 790*S, -107*S,
-1004*S, -2449*S, 939*S, 5568*S, 156*S, 1812*S, 2792*S, -1094*S,
-2677*S, -251*S, 2297*S, 943*S, -1329*S, 2883*S, -853*S, -2626*S,
-105929*S, -6552*S, 1095*S, -1508*S, 1003*S, 5039*S, -2600*S, -749*S,
1790*S, 890*S, 2016*S, -1073*S, 624*S, -2084*S, -1536*S, -1330*S,
358*S, 2444*S, -179*S,-25759*S, -243*S, -552*S, -124*S, 3766*S,
1192*S, -1614*S, 6*S, -1227*S, 345*S, -981*S, -295*S, -1006*S,
-995*S, -1195*S, 706*S, 2512*S, -1758*S, -734*S, -6286*S, -922*S,
1530*S, 1542*S, 1223*S, 61*S, -83*S, 522*S,116937*S, -914*S,
-418*S, -7339*S, 249*S, -520*S, -762*S, 426*S, -505*S, 2664*S,
-1093*S, -1035*S, 2130*S, 4878*S, 1982*S, 1551*S, 2304*S, 193*S,
1532*S, -7268*S, 24357*S, 531*S, 2676*S, -1170*S, 1465*S, -1917*S,
2143*S, 1466*S, -7*S, -7300*S, 3297*S, -1197*S, -289*S, -1548*S,
26226*S, 4401*S, 4123*S, -1588*S, 4243*S, 4069*S, -1276*S, -2010*S,
1407*S, 1478*S, 488*S, -2366*S, -2909*S, -2534*S, -1285*S, 7095*S,
-645*S, -2089*S, -944*S, -40*S, -1363*S, -833*S, 917*S, 1609*S,
1286*S, 1677*S, 1613*S, -2295*S, -1248*S, 40*S, 26*S, 2038*S,
698*S, 2675*S, -1755*S, -3522*S, -1614*S, -6111*S, 270*S, 1822*S,
-234*S, -2844*S, -1201*S, -830*S, 1193*S, 2354*S, 47*S, 1522*S,
-78*S, -640*S, 2425*S, -1596*S, 1563*S, 1169*S, -1006*S, -83*S,
2362*S, -3521*S, -314*S, 1814*S, -1751*S, 305*S, 1715*S, -3741*S,
7847*S, 1291*S, 1206*S, 36*S, 1397*S, -1419*S, -1194*S, -2014*S,
1742*S, -578*S, -207*S, 875*S, 1539*S, 2826*S, -1165*S, -909*S,
1849*S, 927*S, 2018*S, -981*S, 1637*S, -463*S, 905*S, 6618*S,
400*S, 630*S, 2614*S, 900*S, 2323*S, -1094*S, -1858*S, -212*S,
-2069*S, 747*S, 1845*S, -1450*S, 444*S, -213*S, -438*S, 1158*S,
4738*S, 2497*S, -370*S, -2016*S, -518*S, -1160*S, -1510*S, 123*S
};
/* Maximum magnitude in above table: 116937 */
#undef S
@@ -0,0 +1,352 @@
/* Accurate tables for exp2f().
Copyright (C) 1998 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Geoffrey Keating <[email protected]>
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. */
/* This table has the property that, for all integers -128 <= i <= 127,
exp(i/256.0 + __exp2f_deltatable[i-128]) == __exp2f_atable[i+128] + r
for some -2^-35 < r < 2^-35 (abs(r) < 2^-36 if i <= 0); and that
__exp2f_deltatable[i+128] == t * 2^-30
for integer t so that abs(t) <= 43447 * 2^0. */
#define W30 (9.31322575e-10)
static const float __exp2f_deltatable[256] = {
-810*W30, 283*W30, -1514*W30, 1304*W30,
-1148*W30, -98*W30, -744*W30, -156*W30,
-419*W30, -155*W30, 474*W30, 167*W30,
-1984*W30, -826*W30, 692*W30, 781*W30,
-578*W30, -411*W30, -129*W30, -1500*W30,
654*W30, -141*W30, -816*W30, -53*W30,
148*W30, 493*W30, -2214*W30, 760*W30,
260*W30, 750*W30, -1300*W30, 1424*W30,
-1445*W30, -339*W30, -680*W30, -349*W30,
-922*W30, 531*W30, 193*W30, -2892*W30,
290*W30, -2145*W30, -276*W30, 485*W30,
-695*W30, 215*W30, -7093*W30, 412*W30,
-4596*W30, 367*W30, 592*W30, -615*W30,
-97*W30, -1066*W30, 972*W30, -226*W30,
-625*W30, -374*W30, -5647*W30, -180*W30,
20349*W30, -447*W30, 111*W30, -4164*W30,
-87*W30, -21*W30, -251*W30, 66*W30,
-517*W30, 2093*W30, -263*W30, 182*W30,
-601*W30, 475*W30, -483*W30, -1251*W30,
-373*W30, 1471*W30, -92*W30, -215*W30,
-97*W30, -190*W30, 0*W30, -290*W30,
-2647*W30, 1940*W30, -582*W30, 28*W30,
833*W30, 1493*W30, 34*W30, 321*W30,
3327*W30, -35*W30, 177*W30, -135*W30,
-796*W30, -428*W30, 129*W30, 9332*W30,
-12*W30, -69*W30, -1743*W30, 6508*W30,
-60*W30, 359*W30, 43447*W30, 15*W30,
-23*W30, -305*W30, -375*W30, -652*W30,
667*W30, 269*W30, -1575*W30, 185*W30,
-329*W30, 200*W30, 6002*W30, 163*W30,
-647*W30, 19*W30, -603*W30, -755*W30,
742*W30, -438*W30, 3587*W30, 2560*W30,
0*W30, -520*W30, -241*W30, -299*W30,
-1270*W30, -991*W30, -1138*W30, 255*W30,
-1192*W30, 1722*W30, 1023*W30, 3700*W30,
-1388*W30, -1551*W30, -2549*W30, 27*W30,
282*W30, 673*W30, 113*W30, 1561*W30,
72*W30, 873*W30, 87*W30, -395*W30,
-433*W30, 629*W30, 3440*W30, -284*W30,
-592*W30, -103*W30, -46*W30, -3844*W30,
1712*W30, 303*W30, 1555*W30, -631*W30,
-1400*W30, -961*W30, -854*W30, -276*W30,
407*W30, 833*W30, -345*W30, -1501*W30,
121*W30, -1581*W30, 400*W30, 150*W30,
1224*W30, -139*W30, -563*W30, 879*W30,
933*W30, 2939*W30, 788*W30, 211*W30,
530*W30, -192*W30, 706*W30, -13347*W30,
1065*W30, 3*W30, 111*W30, -208*W30,
-360*W30, -532*W30, -291*W30, 483*W30,
987*W30, -33*W30, -1373*W30, -166*W30,
-1174*W30, -3955*W30, 1601*W30, -280*W30,
1405*W30, 600*W30, -1659*W30, -23*W30,
390*W30, 449*W30, 570*W30, -13143*W30,
-9*W30, -1646*W30, 1201*W30, 294*W30,
2181*W30, -1173*W30, 1388*W30, -4504*W30,
190*W30, -2304*W30, 211*W30, 239*W30,
48*W30, -817*W30, 1018*W30, 1828*W30,
-663*W30, 1408*W30, 408*W30, -36*W30,
1295*W30, -230*W30, 1341*W30, 9*W30,
40*W30, 705*W30, 186*W30, 376*W30,
557*W30, 5866*W30, 363*W30, -1558*W30,
718*W30, 669*W30, 1369*W30, -2972*W30,
-468*W30, -121*W30, -219*W30, 667*W30,
29954*W30, 366*W30, 48*W30, -203*W30
};
static const float __exp2f_atable[256] /* __attribute__((mode(SF))) */ = {
0.707106411447, /* 0x0.b504ecfff */
0.709024071690, /* 0x0.b58299fff */
0.710945606239, /* 0x0.b60088000 */
0.712874472142, /* 0x0.b67ef1000 */
0.714806139464, /* 0x0.b6fd88fff */
0.716744661340, /* 0x0.b77c94000 */
0.718687653549, /* 0x0.b7fbea000 */
0.720636486992, /* 0x0.b87ba1fff */
0.722590208040, /* 0x0.b8fbabfff */
0.724549472323, /* 0x0.b97c12fff */
0.726514220228, /* 0x0.b9fcd5fff */
0.728483855735, /* 0x0.ba7deb000 */
0.730457961549, /* 0x0.baff4afff */
0.732438981522, /* 0x0.bb811efff */
0.734425544748, /* 0x0.bc0350000 */
0.736416816713, /* 0x0.bc85d0000 */
0.738412797450, /* 0x0.bd089efff */
0.740414917465, /* 0x0.bd8bd4fff */
0.742422521111, /* 0x0.be0f66fff */
0.744434773914, /* 0x0.be9346fff */
0.746454179287, /* 0x0.bf179f000 */
0.748477637755, /* 0x0.bf9c3afff */
0.750506639473, /* 0x0.c02133fff */
0.752541840064, /* 0x0.c0a694fff */
0.754582285889, /* 0x0.c12c4e000 */
0.756628334525, /* 0x0.c1b265000 */
0.758678436269, /* 0x0.c238bffff */
0.760736882681, /* 0x0.c2bfa6fff */
0.762799203401, /* 0x0.c346cf000 */
0.764867603790, /* 0x0.c3ce5d000 */
0.766940355298, /* 0x0.c45633fff */
0.769021093841, /* 0x0.c4de90fff */
0.771104693409, /* 0x0.c5671dfff */
0.773195922364, /* 0x0.c5f02afff */
0.775292098512, /* 0x0.c6798afff */
0.777394294745, /* 0x0.c70350000 */
0.779501736166, /* 0x0.c78d6d000 */
0.781615912910, /* 0x0.c817fafff */
0.783734917628, /* 0x0.c8a2d9fff */
0.785858273516, /* 0x0.c92e02000 */
0.787990570071, /* 0x0.c9b9c0000 */
0.790125787245, /* 0x0.ca45aefff */
0.792268991467, /* 0x0.cad223fff */
0.794417440881, /* 0x0.cb5ef0fff */
0.796570718287, /* 0x0.cbec0efff */
0.798730909811, /* 0x0.cc79a0fff */
0.800892710672, /* 0x0.cd074dfff */
0.803068041795, /* 0x0.cd95ddfff */
0.805242776881, /* 0x0.ce2464000 */
0.807428598393, /* 0x0.ceb3a3fff */
0.809617877002, /* 0x0.cf431dfff */
0.811812341211, /* 0x0.cfd2eefff */
0.814013659956, /* 0x0.d06333000 */
0.816220164311, /* 0x0.d0f3ce000 */
0.818434238424, /* 0x0.d184e7fff */
0.820652604094, /* 0x0.d21649fff */
0.822877407074, /* 0x0.d2a818000 */
0.825108587751, /* 0x0.d33a51000 */
0.827342867839, /* 0x0.d3ccbdfff */
0.829588949684, /* 0x0.d45ff1000 */
0.831849217401, /* 0x0.d4f411fff */
0.834093391880, /* 0x0.d58724fff */
0.836355149750, /* 0x0.d61b5f000 */
0.838620424257, /* 0x0.d6afd3fff */
0.840896368027, /* 0x0.d744fc000 */
0.843176305293, /* 0x0.d7da66fff */
0.845462262643, /* 0x0.d87037000 */
0.847754716864, /* 0x0.d90673fff */
0.850052893157, /* 0x0.d99d10fff */
0.852359056469, /* 0x0.da3433fff */
0.854668736446, /* 0x0.dacb91fff */
0.856986224651, /* 0x0.db6373000 */
0.859309315673, /* 0x0.dbfbb1fff */
0.861639738080, /* 0x0.dc946bfff */
0.863975346095, /* 0x0.dd2d7d000 */
0.866317391394, /* 0x0.ddc6f9fff */
0.868666708472, /* 0x0.de60f1000 */
0.871022939695, /* 0x0.defb5c000 */
0.873383641229, /* 0x0.df9611fff */
0.875751554968, /* 0x0.e03141000 */
0.878126025200, /* 0x0.e0ccde000 */
0.880506813521, /* 0x0.e168e4fff */
0.882894217966, /* 0x0.e2055afff */
0.885287821299, /* 0x0.e2a239000 */
0.887686729423, /* 0x0.e33f6ffff */
0.890096127973, /* 0x0.e3dd56fff */
0.892507970338, /* 0x0.e47b67000 */
0.894928157336, /* 0x0.e51a03000 */
0.897355020043, /* 0x0.e5b90efff */
0.899788379682, /* 0x0.e65888000 */
0.902227103705, /* 0x0.e6f85afff */
0.904673457151, /* 0x0.e798ae000 */
0.907128036008, /* 0x0.e8398afff */
0.909585535528, /* 0x0.e8da99000 */
0.912051796915, /* 0x0.e97c3a000 */
0.914524436003, /* 0x0.ea1e46000 */
0.917003571999, /* 0x0.eac0bf000 */
0.919490039339, /* 0x0.eb63b2fff */
0.921983361257, /* 0x0.ec071a000 */
0.924488604054, /* 0x0.ecab48fff */
0.926989555360, /* 0x0.ed4f30000 */
0.929502844812, /* 0x0.edf3e6000 */
0.932021975503, /* 0x0.ee98fdfff */
0.934553921208, /* 0x0.ef3eecfff */
0.937083780759, /* 0x0.efe4b8fff */
0.939624726786, /* 0x0.f08b3f000 */
0.942198514924, /* 0x0.f133ebfff */
0.944726586343, /* 0x0.f1d99a000 */
0.947287976728, /* 0x0.f28176fff */
0.949856162070, /* 0x0.f329c5fff */
0.952431440345, /* 0x0.f3d28bfff */
0.955013573175, /* 0x0.f47bc5000 */
0.957603693021, /* 0x0.f52584000 */
0.960199773321, /* 0x0.f5cfa7000 */
0.962801992906, /* 0x0.f67a31000 */
0.965413510788, /* 0x0.f72556fff */
0.968030691152, /* 0x0.f7d0dc000 */
0.970655620084, /* 0x0.f87ce2fff */
0.973290979849, /* 0x0.f92998fff */
0.975926160805, /* 0x0.f9d64bfff */
0.978571653370, /* 0x0.fa83ac000 */
0.981225252139, /* 0x0.fb3193fff */
0.983885228626, /* 0x0.fbdfe6fff */
0.986552715296, /* 0x0.fc8eb7fff */
0.989228487027, /* 0x0.fd3e14000 */
0.991909801964, /* 0x0.fdedcd000 */
0.994601726545, /* 0x0.fe9e38000 */
0.997297704209, /* 0x0.ff4ee6fff */
1.000000000000, /* 0x1.000000000 */
1.002710938457, /* 0x1.00b1aa000 */
1.005429744692, /* 0x1.0163d7ffe */
1.008155703526, /* 0x1.02167dffe */
1.010888457284, /* 0x1.02c995fff */
1.013629436498, /* 0x1.037d38000 */
1.016377568250, /* 0x1.043152000 */
1.019134163841, /* 0x1.04e5f9ffe */
1.021896362316, /* 0x1.059b00000 */
1.024668931945, /* 0x1.0650b3ffe */
1.027446627635, /* 0x1.0706be001 */
1.030234098408, /* 0x1.07bd6bffe */
1.033023953416, /* 0x1.087441ffe */
1.035824656494, /* 0x1.092bce000 */
1.038632392900, /* 0x1.09e3d0001 */
1.041450142840, /* 0x1.0a9c79ffe */
1.044273972530, /* 0x1.0b558a001 */
1.047105550795, /* 0x1.0c0f1c001 */
1.049944162390, /* 0x1.0cc924001 */
1.052791833895, /* 0x1.0d83c4001 */
1.055645227426, /* 0x1.0e3ec3fff */
1.058507919326, /* 0x1.0efa60001 */
1.061377286898, /* 0x1.0fb66bfff */
1.064254641510, /* 0x1.1072fdffe */
1.067140102389, /* 0x1.113018000 */
1.070034146304, /* 0x1.11edc1fff */
1.072937250162, /* 0x1.12ac04001 */
1.075843691823, /* 0x1.136a7dfff */
1.078760385496, /* 0x1.1429a3ffe */
1.081685543070, /* 0x1.14e958000 */
1.084618330005, /* 0x1.15a98c000 */
1.087556362176, /* 0x1.166a18001 */
1.090508937863, /* 0x1.172b98001 */
1.093464612954, /* 0x1.17ed4bfff */
1.096430182434, /* 0x1.18afa5ffe */
1.099401354802, /* 0x1.19725e000 */
1.102381587017, /* 0x1.1a35adfff */
1.105370759965, /* 0x1.1af994000 */
1.108367800686, /* 0x1.1bbdfdffe */
1.111373305331, /* 0x1.1c82f6000 */
1.114387035385, /* 0x1.1d4878001 */
1.117408752440, /* 0x1.1e0e7ffff */
1.120437502874, /* 0x1.1ed4fe000 */
1.123474478729, /* 0x1.1f9c06000 */
1.126521706601, /* 0x1.2063ba001 */
1.129574775716, /* 0x1.212bd0001 */
1.132638812065, /* 0x1.21f49e000 */
1.135709524130, /* 0x1.22bddbffe */
1.138789534565, /* 0x1.2387b5fff */
1.141876101508, /* 0x1.2451fe000 */
1.144971728301, /* 0x1.251cddffe */
1.148077130296, /* 0x1.25e861ffe */
1.151189923305, /* 0x1.26b462001 */
1.154312610610, /* 0x1.278107ffe */
1.157440662410, /* 0x1.284e08001 */
1.160578370109, /* 0x1.291baa001 */
1.163725256932, /* 0x1.29e9e6000 */
1.166879892324, /* 0x1.2ab8a3ffe */
1.170044302935, /* 0x1.2b8805fff */
1.173205971694, /* 0x1.2c5739ffe */
1.176397800428, /* 0x1.2d2867ffe */
1.179586529747, /* 0x1.2df962001 */
1.182784795737, /* 0x1.2ecafbffe */
1.185991406414, /* 0x1.2f9d21ffe */
1.189206838636, /* 0x1.306fdc001 */
1.192430973067, /* 0x1.314328000 */
1.195664167430, /* 0x1.32170c001 */
1.198906540890, /* 0x1.32eb8a001 */
1.202157497408, /* 0x1.33c098000 */
1.205416083326, /* 0x1.349625fff */
1.208683252332, /* 0x1.356c43fff */
1.211961269402, /* 0x1.364318001 */
1.215246438983, /* 0x1.371a64000 */
1.218539118740, /* 0x1.37f22dffe */
1.221847295770, /* 0x1.38cafc000 */
1.225158572187, /* 0x1.39a3fdfff */
1.228481650325, /* 0x1.3a7dc5ffe */
1.231811761846, /* 0x1.3b5803fff */
1.235149741144, /* 0x1.3c32c5ffe */
1.238499879811, /* 0x1.3d0e53ffe */
1.241858124726, /* 0x1.3dea69fff */
1.245225191102, /* 0x1.3ec713fff */
1.248601436624, /* 0x1.3fa458000 */
1.251975655584, /* 0x1.40817a001 */
1.255380749731, /* 0x1.4160a2001 */
1.258783102010, /* 0x1.423f9bffe */
1.262198328973, /* 0x1.431f6e000 */
1.265619754780, /* 0x1.43ffa7fff */
1.269052743928, /* 0x1.44e0a4001 */
1.272490739830, /* 0x1.45c1f4000 */
1.275942921659, /* 0x1.46a432001 */
1.279397487615, /* 0x1.478697ffe */
1.282870173427, /* 0x1.486a2dffe */
1.286346316319, /* 0x1.494dfdffe */
1.289836049094, /* 0x1.4a32b2001 */
1.293333172770, /* 0x1.4b17e1ffe */
1.296839594835, /* 0x1.4bfdadfff */
1.300354957560, /* 0x1.4ce40fffe */
1.303882122055, /* 0x1.4dcb38001 */
1.307417988757, /* 0x1.4eb2f1ffe */
1.310960650439, /* 0x1.4f9b1dfff */
1.314516782746, /* 0x1.50842bfff */
1.318079948424, /* 0x1.516daffff */
1.321653246888, /* 0x1.5257de000 */
1.325237751030, /* 0x1.5342c8001 */
1.328829526907, /* 0x1.542e2c000 */
1.332433700535, /* 0x1.551a5fffe */
1.336045145966, /* 0x1.56070dffe */
1.339667558645, /* 0x1.56f473ffe */
1.343300342533, /* 0x1.57e287ffe */
1.346941947961, /* 0x1.58d130001 */
1.350594043714, /* 0x1.59c087ffe */
1.354256033883, /* 0x1.5ab085fff */
1.357932448365, /* 0x1.5ba175ffe */
1.361609339707, /* 0x1.5c926dfff */
1.365299344044, /* 0x1.5d8441ffe */
1.369003057507, /* 0x1.5e76fc001 */
1.372714757920, /* 0x1.5f6a3c000 */
1.376437187179, /* 0x1.605e2fffe */
1.380165219333, /* 0x1.615282001 */
1.383909463864, /* 0x1.6247e3ffe */
1.387661933907, /* 0x1.633dd0000 */
1.391424179060, /* 0x1.64345fffe */
1.395197510706, /* 0x1.652ba9fff */
1.399006724329, /* 0x1.66254dffe */
1.402773022651, /* 0x1.671c22000 */
1.406576037403, /* 0x1.68155dfff */
1.410389423392, /* 0x1.690f48001 */
};
@@ -0,0 +1,136 @@
SubDir HAIKU_TOP src system libroot posix glibc arch x86_64 ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc include arch $(TARGET_ARCH) ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc include arch generic ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc include ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc libio ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc stdlib ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc math ;
SubDirSysHdrs $(HAIKU_TOP) src system libroot posix glibc ;
SubDirHdrs $(HAIKU_TOP) src system libroot posix glibc arch generic ;
UsePrivateHeaders libroot ;
if $(OPTIM) = -O0 {
OPTIM = -O ;
}
# don't compile with debugging
DEBUG = 0 ;
SubDirCcFlags -D_GNU_SOURCE -D_IEEE_LIBM -DPIC ;
SubDirAsFlags -DPIC ;
local genericSources =
cmp.c dbl2mpn.c divrem.c
memrchr.c
mpn2dbl.c mpn2flt.c mpn2ldbl.c
mul.c mul_n.c
branred.c doasin.c dosincos.c halfulp.c mpa.c mpatan.c mpatan2.c mpexp.c
mplog.c mpsqrt.c mptan.c sincos32.c slowexp.c slowpow.c
e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_acoshl.c e_asin.c e_asinf.c
e_asinl.c e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_atanhl.c e_cosh.c
e_coshf.c e_coshl.c e_exp.c e_exp10.c e_exp10f.c e_exp10l.c e_exp2.c
e_exp2f.c e_expf.c e_fmod.c e_fmodf.c e_gamma_r.c e_gammaf_r.c e_gammal_r.c
e_hypot.c e_hypotf.c e_hypotl.c e_j0.c e_j0f.c e_j0l.c e_j1.c e_j1f.c
e_j1l.c e_jn.c e_jnf.c e_jnl.c e_lgamma_r.c e_lgammaf_r.c e_lgammal_r.c
e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c e_pow.c e_powf.c
e_rem_pio2.c e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c
e_scalbf.c e_sinh.c e_sinhf.c e_sinhl.c
k_cos.c k_cosf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c k_sinf.c k_tan.c
k_tanf.c
s_asinh.c s_asinhf.c s_asinhl.c s_atan.c s_atanf.c s_cacos.c s_cacosf.c
s_cacosh.c s_cacoshf.c s_cacoshl.c s_cacosl.c s_casin.c s_casinf.c
s_casinh.c s_casinhf.c s_casinhl.c s_casinl.c s_catan.c s_catanf.c
s_catanh.c s_catanhf.c s_catanhl.c s_catanl.c s_cbrt.c s_cbrtf.c
s_cbrtl.c s_ccos.c s_ccosf.c s_ccosh.c s_ccoshf.c s_ccoshl.c s_ccosl.c
s_ceil.c s_ceilf.c s_cexp.c s_cexpf.c s_cexpl.c s_clog.c s_clog10.c
s_clog10f.c s_clog10l.c s_clogf.c s_clogl.c s_cos.c s_cosf.c s_cpow.c
s_cpowf.c s_cpowl.c s_cproj.c s_cprojf.c s_cprojl.c s_csin.c s_csinf.c
s_csinh.c s_csinhf.c s_csinhl.c s_csinl.c s_csqrt.c s_csqrtf.c s_csqrtl.c
s_ctan.c s_ctanf.c s_ctanh.c s_ctanhf.c s_ctanhl.c s_ctanl.c s_erf.c
s_erff.c s_erfl.c s_expm1.c s_expm1f.c s_fdim.c s_fdimf.c s_finite.c
s_finitef.c s_floor.c s_floorf.c s_fma.c s_fmaf.c s_fmal.c s_fpclassify.c
s_fpclassifyf.c s_frexp.c s_frexpf.c s_frexpl.c s_ilogb.c s_ilogbf.c
s_isinf.c s_isinff.c s_isnan.c s_isnanf.c s_ldexp.c s_ldexpf.c s_ldexpl.c
s_llround.c s_llroundf.c s_llroundl.c s_log1p.c s_log1pf.c s_logb.c
s_logbf.c s_lround.c s_lroundf.c s_lroundl.c s_matherr.c s_modf.c s_modff.c
s_modfl.c s_nan.c s_nanf.c s_nanl.c s_nearbyint.c s_nearbyintf.c
s_nextafter.c s_nextafterf.c s_nexttowardl.c s_remquo.c s_remquof.c
s_remquol.c s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c s_scalbln.c
s_scalblnf.c s_scalblnl.c s_scalbn.c s_scalbnf.c s_signbit.c s_signbitf.c
s_signbitl.c s_signgam.c s_significand.c s_significandf.c s_sin.c
s_sincosf.c s_sinf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tanhl.c
s_trunc.c s_truncf.c
t_exp.c
w_acos.c w_acosf.c w_acosh.c w_acoshf.c w_acoshl.c w_acosl.c w_asin.c
w_asinf.c w_asinl.c w_atan2.c w_atan2f.c w_atan2l.c w_atanh.c w_atanhf.c
w_atanhl.c w_cosh.c w_coshf.c w_coshl.c w_drem.c w_dremf.c w_dreml.c
w_exp.c w_exp10.c w_exp10f.c w_exp10l.c w_exp2.c w_exp2f.c w_exp2l.c
w_expf.c w_expl.c w_fmod.c w_fmodf.c w_fmodl.c w_hypot.c w_hypotf.c
w_hypotl.c w_j0.c w_j0f.c w_j0l.c w_j1.c w_j1f.c w_j1l.c w_jn.c w_jnf.c
w_jnl.c w_lgamma.c w_lgamma_r.c w_lgammaf.c w_lgammaf_r.c w_lgammal.c
w_lgammal_r.c w_log.c w_log10.c w_log10f.c w_log10l.c w_log2.c w_log2f.c
w_log2l.c w_logf.c w_logl.c w_pow.c w_powf.c w_powl.c w_remainder.c
w_remainderf.c w_remainderl.c w_scalb.c w_scalbf.c w_scalbl.c w_sinh.c
w_sinhf.c w_sinhl.c w_sqrt.c w_sqrtf.c w_sqrtl.c w_tgamma.c w_tgammaf.c
w_tgammal.c
;
MergeObject posix_gnu_arch_$(TARGET_ARCH)_e.o :
e_acosl.c e_atan2l.c e_exp2l.S e_expl.c e_fmodl.S e_log10l.S e_log2l.S
e_logl.S e_powl.S e_remainderl.S e_rem_pio2l.c e_scalbl.S e_sqrt.c
e_sqrtf.c e_sqrtl.c
;
MergeObject posix_gnu_arch_$(TARGET_ARCH)_k.o :
k_cosl.c k_rem_pio2l.c k_sinl.c k_tanl.c
;
MergeObject posix_gnu_arch_$(TARGET_ARCH)_s.o :
s_atanl.c s_ceill.S s_copysignf.S s_copysignl.S s_copysign.S s_cosl.S
s_expm1l.S s_fabs.c s_fabsf.c s_fabsl.S s_fdiml.S s_finitel.S s_floorl.S
s_fmaxf.S s_fmaxl.S s_fmax.S s_fminf.S s_fminl.S s_fmin.S s_fpclassifyl.c
s_ilogbl.S s_isinfl.c s_isnanl.c s_llrintf.S s_llrintl.S s_llrint.S
s_log1pl.S s_logbl.c s_lrintf.S s_lrintl.S s_lrint.S s_nearbyintl.S
s_nextafterl.c s_nexttoward.c s_nexttowardf.c s_rintl.c s_scalbnl.S
s_significandl.c s_sincosl.S s_sincos.S s_sinl.S s_tanl.S s_truncl.S
;
MergeObject posix_gnu_arch_$(TARGET_ARCH)_other.o :
add_n.S
addmul_1.S
fegetround.c
fesetenv.c
fesetround.c
ldbl2mpn.c
mul_1.S
lshift.S rshift.S
sub_n.S
submul_1.S
;
MergeObject posix_gnu_arch_$(TARGET_ARCH)_generic.o :
$(genericSources)
;
MergeObjectFromObjects posix_gnu_arch_$(TARGET_ARCH).o : :
posix_gnu_arch_$(TARGET_ARCH)_e.o
posix_gnu_arch_$(TARGET_ARCH)_k.o
posix_gnu_arch_$(TARGET_ARCH)_s.o
posix_gnu_arch_$(TARGET_ARCH)_other.o
posix_gnu_arch_$(TARGET_ARCH)_generic.o
;
SEARCH on [ FGristFiles $(genericSources) ]
= [ FDirName $(HAIKU_TOP) src system libroot posix glibc arch generic ] ;
@@ -0,0 +1,42 @@
/* Add two limb vectors of the same length > 0 and store sum in a third
limb vector.
Copyright (C) 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_add_n)
leaq (%rsi,%rcx,8), %rsi
leaq (%rdi,%rcx,8), %rdi
leaq (%rdx,%rcx,8), %rdx
negq %rcx
xorl %eax, %eax # clear cy
.p2align 2
L(loop):
movq (%rsi,%rcx,8), %rax
movq (%rdx,%rcx,8), %r10
adcq %r10, %rax
movq %rax, (%rdi,%rcx,8)
incq %rcx
jne L(loop)
movq %rcx, %rax # zero %rax
adcq %rax, %rax
ret
END (__mpn_add_n)
@@ -0,0 +1,46 @@
/* AMD64 __mpn_addmul_1 -- Multiply a limb vector with a limb and add
the result to a second limb vector.
Copyright (C) 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_addmul_1)
movq %rdx, %r11
leaq (%rsi,%rdx,8), %rsi
leaq (%rdi,%rdx,8), %rdi
negq %r11
xorl %r8d, %r8d
xorl %r10d, %r10d
.p2align 2
L(loop):
movq (%rsi,%r11,8), %rax
mulq %rcx
addq (%rdi,%r11,8), %rax
adcq %r10, %rdx
addq %r8, %rax
movq %r10, %r8
movq %rax, (%rdi,%r11,8)
adcq %rdx, %r8
incq %r11
jne L(loop)
movq %r8, %rax
ret
END (__mpn_addmul_1)
@@ -0,0 +1,25 @@
/*
* Written by J.T. Conklin <[email protected]>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <[email protected]>.
*/
#include "math_private.h"
long double
__ieee754_acosl (long double x)
{
long double res;
/* acosl = atanl (sqrtl(1 - x^2) / x) */
asm ( "fld %%st\n"
"fmul %%st(0)\n" /* x^2 */
"fld1\n"
"fsubp\n" /* 1 - x^2 */
"fsqrt\n" /* sqrtl (1 - x^2) */
"fxch %%st(1)\n"
"fpatan"
: "=t" (res) : "0" (x) : "st(1)");
return res;
}
@@ -0,0 +1,18 @@
/*
* Written by J.T. Conklin <[email protected]>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <[email protected]>.
*/
#include "math_private.h"
long double
__ieee754_atan2l (long double y, long double x)
{
long double res;
asm ("fpatan" : "=t" (res) : "u" (y), "0" (x) : "st(1)");
return res;
}
@@ -0,0 +1,38 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Public domain.
*/
#include <machine/asm.h>
ENTRY(__ieee754_exp2l)
fldt 8(%rsp)
/* I added the following ugly construct because exp(+-Inf) resulted
in NaN. The ugliness results from the bright minds at Intel.
For the i686 the code can be written better.
-- drepper@cygnus.com. */
fxam /* Is NaN or +-Inf? */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fld %st
frndint /* int(x) */
fsubr %st,%st(1) /* fract(x) */
fxch
f2xm1 /* 2^(fract(x)) - 1 */
fld1
faddp /* 2^(fract(x)) */
fscale /* e^x */
fstp %st(1)
ret
1: testl $0x200, %eax /* Test sign. */
jz 2f /* If positive, jump. */
fstp %st
fldz /* Set result to 0. */
2: ret
END (__ieee754_exp2l)
@@ -0,0 +1,77 @@
/*
* Written by J.T. Conklin <[email protected]>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <[email protected]>.
*/
/*
* The 8087 method for the exponential function is to calculate
* exp(x) = 2^(x log2(e))
* after separating integer and fractional parts
* x log2(e) = i + f, |f| <= .5
* 2^i is immediate but f needs to be precise for long double accuracy.
* Suppress range reduction error in computing f by the following.
* Separate x into integer and fractional parts
* x = xi + xf, |xf| <= .5
* Separate log2(e) into the sum of an exact number c0 and small part c1.
* c0 + c1 = log2(e) to extra precision
* Then
* f = (c0 xi - i) + c0 xf + c1 x
* where c0 xi is exact and so also is (c0 xi - i).
* -- [email protected]
*/
#include "math_private.h"
static const long double c0 = 1.44268798828125L;
static const long double c1 = 7.05260771340735992468e-6L;
long double
__ieee754_expl (long double x)
{
long double res;
/* I added the following ugly construct because expl(+-Inf) resulted
in NaN. The ugliness results from the bright minds at Intel.
For the i686 the code can be written better.
-- [email protected]. */
asm ("fxam\n\t" /* Is NaN or +-Inf? */
"fstsw %%ax\n\t"
"movb $0x45, %%dh\n\t"
"andb %%ah, %%dh\n\t"
"cmpb $0x05, %%dh\n\t"
"je 1f\n\t" /* Is +-Inf, jump. */
"fldl2e\n\t" /* 1 log2(e) */
"fmul %%st(1),%%st\n\t" /* 1 x log2(e) */
"frndint\n\t" /* 1 i */
"fld %%st(1)\n\t" /* 2 x */
"frndint\n\t" /* 2 xi */
"fld %%st(1)\n\t" /* 3 i */
"fldt %2\n\t" /* 4 c0 */
"fld %%st(2)\n\t" /* 5 xi */
"fmul %%st(1),%%st\n\t" /* 5 c0 xi */
"fsubp %%st,%%st(2)\n\t" /* 4 f = c0 xi - i */
"fld %%st(4)\n\t" /* 5 x */
"fsub %%st(3),%%st\n\t" /* 5 xf = x - xi */
"fmulp %%st,%%st(1)\n\t" /* 4 c0 xf */
"faddp %%st,%%st(1)\n\t" /* 3 f = f + c0 xf */
"fldt %3\n\t" /* 4 */
"fmul %%st(4),%%st\n\t" /* 4 c1 * x */
"faddp %%st,%%st(1)\n\t" /* 3 f = f + c1 * x */
"f2xm1\n\t" /* 3 2^(fract(x * log2(e))) - 1 */
"fld1\n\t" /* 4 1.0 */
"faddp\n\t" /* 3 2^(fract(x * log2(e))) */
"fstp %%st(1)\n\t" /* 2 */
"fscale\n\t" /* 2 scale factor is st(1); e^x */
"fstp %%st(1)\n\t" /* 1 */
"fstp %%st(1)\n\t" /* 0 */
"jmp 2f\n\t"
"1:\ttestl $0x200, %%eax\n\t" /* Test sign. */
"jz 2f\n\t" /* If positive, jump. */
"fstp %%st\n\t"
"fldz\n\t" /* Set result to 0. */
"2:\t\n"
: "=t" (res) : "0" (x), "m" (c0), "m" (c1) : "ax", "dx");
return res;
}
@@ -0,0 +1,22 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
ENTRY(__ieee754_fmodl)
fldt 24(%rsp)
fldt 8(%rsp)
1: fprem
fstsw %ax
and $04,%ah
jnz 1b
fstp %st(1)
ret
END (__ieee754_fmodl)
@@ -0,0 +1,67 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
*
* Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
/* It is not important that this constant is precise. It is only
a value which is known to be on the safe side for using the
fyl2xp1 instruction. */
ASM_TYPE_DIRECTIVE(limit,@object)
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_log10l)
fldlg2 // log10(2)
fldt 8(%rsp) // x : log10(2)
fxam
fnstsw
fld %st // x : x : log10(2)
testb $1, %ah
jnz 3f // in case x is NaN or ±Inf
4: fsubl MO(one) // x-1 : x : log10(2)
fld %st // x-1 : x-1 : x : log10(2)
fabs // |x-1| : x-1 : x : log10(2)
fcompl MO(limit) // x-1 : x : log10(2)
fnstsw // x-1 : x : log10(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log10(2)
fyl2xp1 // log10(x)
ret
2: fstp %st(0) // x : log10(2)
fyl2x // log10(x)
ret
3: testb $4, %ah
jnz 4b // in case x is ±Inf
fstp %st(1)
fstp %st(1)
ret
END(__ieee754_log10l)
@@ -0,0 +1,64 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
* Public domain.
*
* Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
/* It is not important that this constant is precise. It is only
a value which is known to be on the safe side for using the
fyl2xp1 instruction. */
ASM_TYPE_DIRECTIVE(limit,@object)
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_log2l)
fldl MO(one)
fldt 8(%rsp) // x : 1
fxam
fnstsw
fld %st // x : x : 1
testb $1, %ah
jnz 3f // in case x is NaN or ±Inf
4: fsub %st(2), %st // x-1 : x : 1
fld %st // x-1 : x-1 : x : 1
fabs // |x-1| : x-1 : x : 1
fcompl MO(limit) // x-1 : x : 1
fnstsw // x-1 : x : 1
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : 1
fyl2xp1 // log(x)
ret
2: fstp %st(0) // x : 1
fyl2x // log(x)
ret
3: testb $4, %ah
jnz 4b // in case x is ±Inf
fstp %st(1)
fstp %st(1)
ret
END (__ieee754_log2l)
@@ -0,0 +1,66 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
/* It is not important that this constant is precise. It is only
a value which is known to be on the safe side for using the
fyl2xp1 instruction. */
ASM_TYPE_DIRECTIVE(limit,@object)
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_logl)
fldln2 // log(2)
fldt 8(%rsp) // x : log(2)
fxam
fnstsw
fld %st // x : x : log(2)
testb $1, %ah
jnz 3f // in case x is NaN or +-Inf
4: fsubl MO(one) // x-1 : x : log(2)
fld %st // x-1 : x-1 : x : log(2)
fabs // |x-1| : x-1 : x : log(2)
fcompl MO(limit) // x-1 : x : log(2)
fnstsw // x-1 : x : log(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret
2: fstp %st(0) // x : log(2)
fyl2x // log(x)
ret
3: testb $4, %ah
jnz 4b // in case x is +-Inf
fstp %st(1)
fstp %st(1)
ret
END (__ieee754_logl)
@@ -0,0 +1,341 @@
/* ix87 specific implementation of pow function.
Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2007
Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
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 <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(infinity,@object)
inf_zero:
infinity:
.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
ASM_SIZE_DIRECTIVE(infinity)
ASM_TYPE_DIRECTIVE(zero,@object)
zero: .double 0.0
ASM_SIZE_DIRECTIVE(zero)
ASM_TYPE_DIRECTIVE(minf_mzero,@object)
minf_mzero:
minfinity:
.byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
ASM_TYPE_DIRECTIVE(limit,@object)
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
ASM_TYPE_DIRECTIVE(p63,@object)
p63:
.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_SIZE_DIRECTIVE(p63)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_powl)
fldt 24(%rsp) // y
fxam
fnstsw
movb %ah, %dl
andb $0x45, %ah
cmpb $0x40, %ah // is y == 0 ?
je 11f
cmpb $0x05, %ah // is y == ±inf ?
je 12f
cmpb $0x01, %ah // is y == NaN ?
je 30f
fldt 8(%rsp) // x : y
fxam
fnstsw
movb %ah, %dh
andb $0x45, %ah
cmpb $0x40, %ah
je 20f // x is ±0
cmpb $0x05, %ah
je 15f // x is ±inf
fxch // y : x
/* fistpll raises invalid exception for |y| >= 1L<<63. */
fldl MO(p63) // 1L<<63 : y : x
fld %st(1) // y : 1L<<63 : y : x
fabs // |y| : 1L<<63 : y : x
fcomip %st(1), %st // 1L<<63 : y : x
fstp %st(0) // y : x
jnc 2f
/* First see whether `y' is a natural number. In this case we
can use a more precise algorithm. */
fld %st // y : y : x
fistpll -8(%rsp) // y : x
fildll -8(%rsp) // int(y) : y : x
fucomip %st(1),%st // y : x
jne 2f
/* OK, we have an integer value for y. */
mov -8(%rsp),%eax
mov -4(%rsp),%edx
orl $0, %edx
fstp %st(0) // x
jns 4f // y >= 0, jump
fdivrl MO(one) // 1/x (now referred to as x)
negl %eax
adcl $0, %edx
negl %edx
4: fldl MO(one) // 1 : x
fxch
6: shrdl $1, %edx, %eax
jnc 5f
fxch
fmul %st(1) // x : ST*x
fxch
5: fmul %st(0), %st // x*x : ST*x
shrl $1, %edx
movl %eax, %ecx
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
ret
/* y is ±NAN */
30: fldt 8(%rsp) // x : y
fldl MO(one) // 1.0 : x : y
fucomip %st(1),%st // x : y
je 31f
fxch // y : x
31: fstp %st(1)
ret
.align ALIGNARG(4)
2: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fldl MO(limit) // 0.29 : 1.0 : x : y
fld %st(2) // x : 0.29 : 1.0 : x : y
fsub %st(2) // x-1 : 0.29 : 1.0 : x : y
fabs // |x-1| : 0.29 : 1.0 : x : y
fucompp // 1.0 : x : y
fnstsw
fxch // x : 1.0 : y
test $4500,%eax
jz 7f
fsub %st(1) // x-1 : 1.0 : y
fyl2xp1 // log2(x) : y
jmp 8f
7: fyl2x // log2(x) : y
8: fmul %st(1) // y*log2(x) : y
fxam
fnstsw
andb $0x45, %ah
cmpb $0x05, %ah // is y*log2(x) == ±inf ?
je 28f
fst %st(1) // y*log2(x) : y*log2(x)
frndint // int(y*log2(x)) : y*log2(x)
fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
fxch // fract(y*log2(x)) : int(y*log2(x))
f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
ret
28: fstp %st(1) // y*log2(x)
fldl MO(one) // 1 : y*log2(x)
fscale // 2^(y*log2(x)) : y*log2(x)
fstp %st(1) // 2^(y*log2(x))
ret
// pow(x,±0) = 1
.align ALIGNARG(4)
11: fstp %st(0) // pop y
fldl MO(one)
ret
// y == ±inf
.align ALIGNARG(4)
12: fstp %st(0) // pop y
fldl MO(one) // 1
fldt 8(%rsp) // x : 1
fabs // abs(x) : 1
fucompp // < 1, == 1, or > 1
fnstsw
andb $0x45, %ah
cmpb $0x45, %ah
je 13f // jump if x is NaN
cmpb $0x40, %ah
je 14f // jump if |x| == 1
shlb $1, %ah
xorb %ah, %dl
andl $2, %edx
#ifdef PIC
lea inf_zero(%rip),%rcx
fldl (%rcx, %rdx, 4)
#else
fldl inf_zero(,%rdx, 4)
#endif
ret
.align ALIGNARG(4)
14: fldl MO(one)
ret
.align ALIGNARG(4)
13: fldt 8(%rsp) // load x == NaN
ret
.align ALIGNARG(4)
// x is ±inf
15: fstp %st(0) // y
testb $2, %dh
jz 16f // jump if x == +inf
// We must find out whether y is an odd integer.
fld %st // y : y
fistpll -8(%rsp) // y
fildll -8(%rsp) // int(y) : y
fucomip %st(1),%st
ffreep %st // <empty>
jne 17f
// OK, the value is an integer, but is it odd?
mov -8(%rsp), %eax
mov -4(%rsp), %edx
andb $1, %al
jz 18f // jump if not odd
// It's an odd integer.
shrl $31, %edx
#ifdef PIC
lea minf_mzero(%rip),%rcx
fldl (%rcx, %rdx, 8)
#else
fldl minf_mzero(,%rdx, 8)
#endif
ret
.align ALIGNARG(4)
16: fcompl MO(zero)
fnstsw
shrl $5, %eax
andl $8, %eax
#ifdef PIC
lea inf_zero(%rip),%rcx
fldl (%rcx, %rax, 1)
#else
fldl inf_zero(,%rax, 1)
#endif
ret
.align ALIGNARG(4)
17: shll $30, %edx // sign bit for y in right position
18: shrl $31, %edx
#ifdef PIC
lea inf_zero(%rip),%rcx
fldl (%rcx, %rdx, 8)
#else
fldl inf_zero(,%rdx, 8)
#endif
ret
.align ALIGNARG(4)
// x is ±0
20: fstp %st(0) // y
testb $2, %dl
jz 21f // y > 0
// x is ±0 and y is < 0. We must find out whether y is an odd integer.
testb $2, %dh
jz 25f
fld %st // y : y
fistpll -8(%rsp) // y
fildll -8(%rsp) // int(y) : y
fucomip %st(1),%st
ffreep %st // <empty>
jne 26f
// OK, the value is an integer, but is it odd?
mov -8(%rsp),%eax
mov -4(%rsp),%edx
andb $1, %al
jz 27f // jump if not odd
// It's an odd integer.
// Raise divide-by-zero exception and get minus infinity value.
fldl MO(one)
fdivl MO(zero)
fchs
ret
25: fstp %st(0)
26:
27: // Raise divide-by-zero exception and get infinity value.
fldl MO(one)
fdivl MO(zero)
ret
.align ALIGNARG(4)
// x is ±0 and y is > 0. We must find out whether y is an odd integer.
21: testb $2, %dh
jz 22f
fld %st // y : y
fistpll -8(%rsp) // y
fildll -8(%rsp) // int(y) : y
fucomip %st(1),%st
ffreep %st // <empty>
jne 23f
// OK, the value is an integer, but is it odd?
mov -8(%rsp),%eax
mov -4(%rsp),%edx
andb $1, %al
jz 24f // jump if not odd
// It's an odd integer.
fldl MO(mzero)
ret
22: fstp %st(0)
23:
24: fldl MO(zero)
ret
END(__ieee754_powl)
@@ -0,0 +1,3 @@
/* Empty. This file is only meant to avoid compiling the file with the
same name in the libm-ieee754 directory. The code is not used since
there is an assembler version for all users of this file. */
@@ -0,0 +1,20 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
ENTRY(__ieee754_remainderl)
fldt 24(%rsp)
fldt 8(%rsp)
1: fprem1
fstsw %ax
testl $0x400,%eax
jnz 1b
fstp %st(1)
ret
END (__ieee754_remainderl)
@@ -0,0 +1,100 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>
*
* Correct handling of y==-inf <drepper@gnu>
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(zero_nan,@object)
zero_nan:
.double 0.0
nan: .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
minus_zero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
ASM_SIZE_DIRECTIVE(zero_nan)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_scalbl)
fldt 24(%rsp)
fxam
fnstsw
fldt 8(%rsp)
andl $0x4700, %eax
cmpl $0x0700, %eax
je 1f
andl $0x4500, %eax
cmpl $0x0100, %eax
je 2f
fxam
fnstsw
andl $0x4500, %eax
cmpl $0x0100, %eax
je 3f
fld %st(1)
frndint
fcomip %st(2), %st
jne 4f
fscale
fstp %st(1)
ret
/* y is -inf */
1: fxam
fnstsw
movl 16(%rsp), %edx
shrl $5, %eax
fstp %st
fstp %st
andl $0x8000, %edx
andl $8, %eax
jnz 4f
shrl $11, %edx
addl %edx, %eax
#ifdef PIC
lea zero_nan(%rip),%rdx
fldl (%rdx,%rax,1)
#else
fldl zero_nan(%rax, 1)
#endif
ret
/* The result is NaN, but we must not raise an exception.
So use a variable. */
2: fstp %st
fstp %st
fldl MO(nan)
ret
/* The first parameter is a NaN. Return it. */
3: fstp %st(1)
ret
/* Return NaN and raise the invalid exception. */
4: fstp %st
fstp %st
fldz
fdiv %st
ret
END(__ieee754_scalbl)
@@ -0,0 +1,30 @@
/* Square root of floating point number.
Copyright (C) 2002 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_private.h"
double
__ieee754_sqrt (double x)
{
double res;
asm ("sqrtsd %0, %1" : "=x" (res) : "x" (x));
return res;
}
@@ -0,0 +1,30 @@
/* Square root of floating point number.
Copyright (C) 2002 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_private.h"
float
__ieee754_sqrtf (float x)
{
float res;
asm ("sqrtss %0, %1" : "=x" (res) : "x" (x));
return res;
}
@@ -0,0 +1,18 @@
/*
* Written by J.T. Conklin <[email protected]>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <[email protected]>.
*/
#include "math_private.h"
long double
__ieee754_sqrtl (long double x)
{
long double res;
asm ("fsqrt" : "=t" (res) : "0" (x));
return res;
}
@@ -0,0 +1,33 @@
/* Return current rounding direction.
Copyright (C) 1997, 2001 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <[email protected]>, 1997.
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 <fenv.h>
int
fegetround (void)
{
int cw;
/* We only check the x87 FPU unit. The SSE unit should be the same
- and if it's not the same there's no way to signal it. */
__asm__ ("fnstcw %0" : "=m" (*&cw));
return cw & 0xc00;
}
@@ -0,0 +1,89 @@
/* Install given floating-point environment.
Copyright (C) 2001, 2002, 2003 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 <fenv.h>
#include <assert.h>
int
fesetenv (const fenv_t *envp)
{
fenv_t temp;
/* Install the environment specified by ENVP. But there are a few
values which we do not want to come from the saved environment.
Therefore, we get the current environment and replace the values
we want to use from the environment specified by the parameter. */
__asm__ ("fnstenv %0\n"
"stmxcsr %1" : "=m" (*&temp), "=m" (*&temp.__mxcsr));
if (envp == FE_DFL_ENV)
{
temp.__control_word |= FE_ALL_EXCEPT;
temp.__control_word &= ~FE_TOWARDZERO;
temp.__status_word &= ~FE_ALL_EXCEPT;
temp.__eip = 0;
temp.__cs_selector = 0;
temp.__opcode = 0;
temp.__data_offset = 0;
temp.__data_selector = 0;
/* Set mask for SSE MXCSR. */
temp.__mxcsr |= (FE_ALL_EXCEPT << 7);
/* Set rounding to FE_TONEAREST. */
temp.__mxcsr &= ~ 0x6000;
temp.__mxcsr |= (FE_TONEAREST << 3);
}
else if (envp == FE_NOMASK_ENV)
{
temp.__control_word &= ~(FE_ALL_EXCEPT | FE_TOWARDZERO);
temp.__status_word &= ~FE_ALL_EXCEPT;
temp.__eip = 0;
temp.__cs_selector = 0;
temp.__opcode = 0;
temp.__data_offset = 0;
temp.__data_selector = 0;
/* Set mask for SSE MXCSR. */
/* Set rounding to FE_TONEAREST. */
temp.__mxcsr &= ~ 0x6000;
temp.__mxcsr |= (FE_TONEAREST << 3);
/* Do not mask exceptions. */
temp.__mxcsr &= ~(FE_ALL_EXCEPT << 7);
}
else
{
temp.__control_word &= ~(FE_ALL_EXCEPT | FE_TOWARDZERO);
temp.__control_word |= (envp->__control_word
& (FE_ALL_EXCEPT | FE_TOWARDZERO));
temp.__status_word &= ~FE_ALL_EXCEPT;
temp.__status_word |= envp->__status_word & FE_ALL_EXCEPT;
temp.__eip = envp->__eip;
temp.__cs_selector = envp->__cs_selector;
temp.__opcode = envp->__opcode;
temp.__data_offset = envp->__data_offset;
temp.__data_selector = envp->__data_selector;
temp.__mxcsr = envp->__mxcsr;
}
__asm__ ("fldenv %0\n"
"ldmxcsr %1" : : "m" (temp), "m" (temp.__mxcsr));
/* Success. */
return 0;
}
libm_hidden_def (fesetenv)
@@ -0,0 +1,47 @@
/* Set current rounding direction.
Copyright (C) 2001, 2005 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 <fenv.h>
int
fesetround (int round)
{
unsigned short int cw;
int mxcsr;
if ((round & ~0xc00) != 0)
/* ROUND is no valid rounding mode. */
return 1;
/* First set the x87 FPU. */
asm ("fnstcw %0" : "=m" (*&cw));
cw &= ~0xc00;
cw |= round;
asm ("fldcw %0" : : "m" (*&cw));
/* And now the MSCSR register for SSE, the precision is at different bit
positions in the different units, we need to shift it 3 bits. */
asm ("stmxcsr %0" : "=m" (*&mxcsr));
mxcsr &= ~ 0x6000;
mxcsr |= round << 3;
asm ("ldmxcsr %0" : : "m" (*&mxcsr));
return 0;
}
libm_hidden_def (fesetround)
@@ -0,0 +1 @@
/* Not needed. */
@@ -0,0 +1 @@
/* Not needed. */
@@ -0,0 +1 @@
/* Not needed. */
@@ -0,0 +1 @@
/* Not needed. */
@@ -0,0 +1 @@
#include "../x86/ldbl2mpn.c"
@@ -0,0 +1,60 @@
/* AMD64 __mpn_lshift --
Copyright 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_lshift)
movq -8(%rsi,%rdx,8), %mm7
movd %ecx, %mm1
movl $64, %eax
subl %ecx, %eax
movd %eax, %mm0
movq %mm7, %mm3
psrlq %mm0, %mm7
movd %mm7, %rax
subq $2, %rdx
jl L(endo)
.p2align 2
L(loop):
movq (%rsi,%rdx,8), %mm6
movq %mm6, %mm2
psrlq %mm0, %mm6
psllq %mm1, %mm3
por %mm6, %mm3
movq %mm3, 8(%rdi,%rdx,8)
je L(ende)
movq -8(%rsi,%rdx,8), %mm7
movq %mm7, %mm3
psrlq %mm0, %mm7
psllq %mm1, %mm2
por %mm7, %mm2
movq %mm2, (%rdi,%rdx,8)
subq $2, %rdx
jge L(loop)
L(endo):
movq %mm3, %mm2
L(ende):
psllq %mm1, %mm2
movq %mm2, (%rdi)
emms
ret
END (__mpn_lshift)
@@ -0,0 +1,21 @@
#ifndef _MATH_PRIVATE_H
#define math_opt_barrier(x) \
({ __typeof(x) __x; \
if (sizeof (x) <= sizeof (double)) \
__asm ("" : "=x" (__x) : "0" (x)); \
else \
__asm ("" : "=t" (__x) : "0" (x)); \
__x; })
#define math_force_eval(x) \
do \
{ \
if (sizeof (x) <= sizeof (double)) \
__asm __volatile ("" : : "x" (x)); \
else \
__asm __volatile ("" : : "f" (x)); \
} \
while (0)
#include <math/math_private.h>
#endif
@@ -0,0 +1,42 @@
/* AMD64 __mpn_mul_1 -- Multiply a limb vector with a limb and store
the result in a second limb vector.
Copyright (C) 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include <sysdep.h>
#include "asm-syntax.h"
.text
ENTRY (__mpn_mul_1)
movq %rdx, %r11
leaq (%rsi,%rdx,8), %rsi
leaq (%rdi,%rdx,8), %rdi
negq %r11
xorl %r8d, %r8d
L(loop):
movq (%rsi,%r11,8), %rax
mulq %rcx
addq %r8, %rax
movl $0, %r8d
adcq %rdx, %r8
movq %rax, (%rdi,%r11,8)
incq %r11
jne L(loop)
movq %r8, %rax
ret
END (__mpn_mul_1)
@@ -0,0 +1,62 @@
/* AMD64 __mpn_rshift --
Copyright (C) 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_rshift)
movq (%rsi), %mm7
movd %ecx, %mm1
movl $64, %eax
subl %ecx, %eax
movd %eax, %mm0
movq %mm7, %mm3
psllq %mm0, %mm7
movd %mm7, %rax
leaq (%rsi,%rdx,8), %rsi
leaq (%rdi,%rdx,8), %rdi
negq %rdx
addq $2, %rdx
jg L(endo)
.p2align 2
L(loop):
movq -8(%rsi,%rdx,8), %mm6
movq %mm6, %mm2
psllq %mm0, %mm6
psrlq %mm1, %mm3
por %mm6, %mm3
movq %mm3, -16(%rdi,%rdx,8)
je L(ende)
movq (%rsi,%rdx,8), %mm7
movq %mm7, %mm3
psllq %mm0, %mm7
psrlq %mm1, %mm2
por %mm7, %mm2
movq %mm2, -8(%rdi,%rdx,8)
addq $2, %rdx
jle L(loop)
L(endo):
movq %mm3, %mm2
L(ende):
psrlq %mm1, %mm2
movq %mm2, -8(%rdi)
emms
ret
END (__mpn_rshift)
@@ -0,0 +1,22 @@
/*
* Written by J.T. Conklin <[email protected]>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <[email protected]>.
*/
#include "math_private.h"
long double
__atanl (long double x)
{
long double res;
asm ("fld1\n"
"fpatan"
: "=t" (res) : "0" (x));
return res;
}
weak_alias (__atanl, atanl)
@@ -0,0 +1,31 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Changes for x86-64 by Andreas Jaeger <aj@suse.de>
* Public domain.
*/
#include <machine/asm.h>
ENTRY(__ceill)
fldt 8(%rsp)
fstcw -4(%rsp) /* store fpu control word */
/* We use here %edx although only the low 1 bits are defined.
But none of the operations should care and they are faster
than the 16 bit operations. */
movl $0x0800,%edx /* round towards +oo */
orl -4(%rsp),%edx
andl $0xfbff,%edx
movl %edx,-8(%rsp)
fldcw -8(%rsp) /* load modified control word */
frndint /* round */
fldcw -4(%rsp) /* restore original control word */
ret
END (__ceill)
weak_alias (__ceill, ceill)
@@ -0,0 +1,55 @@
/* copy sign, double version.
Copyright (C) 2002, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(signmask,@object)
signmask:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
.byte 0, 0, 0, 0, 0, 0, 0, 0
ASM_SIZE_DIRECTIVE(signmask)
ASM_TYPE_DIRECTIVE(othermask,@object)
othermask:
.byte 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f
.byte 0, 0, 0, 0, 0, 0, 0, 0
ASM_SIZE_DIRECTIVE(othermask)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__copysign)
andpd MO(othermask),%xmm0
andpd MO(signmask),%xmm1
orpd %xmm1,%xmm0
ret
END (__copysign)
weak_alias (__copysign, copysign)
@@ -0,0 +1,50 @@
/* copy sign, double version.
Copyright (C) 2002, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(mask,@object)
mask:
.byte 0xff, 0xff, 0xff, 0x7f
ASM_SIZE_DIRECTIVE(mask)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__copysignf)
movss MO(mask),%xmm3
andps %xmm3,%xmm0
andnps %xmm1,%xmm3
orps %xmm3,%xmm0
retq
END (__copysignf)
weak_alias (__copysignf, copysignf)
@@ -0,0 +1,22 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>.
* Adopted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Public domain.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
ENTRY(__copysignl)
movl 32(%rsp),%edx
movl 16(%rsp),%eax
andl $0x8000,%edx
andl $0x7fff,%eax
orl %edx,%eax
movl %eax,16(%rsp)
fldt 8(%rsp)
ret
END (__copysignl)
weak_alias (__copysignl, copysignl)
@@ -0,0 +1,40 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Fixed errno handling by Ulrich Drepper <drepper@redhat.com>.
*/
#include <machine/asm.h>
ENTRY(__cosl)
fldt 8(%rsp)
fxam
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 3f
4: fcos
fnstsw %ax
testl $0x400,%eax
jnz 1f
ret
.align ALIGNARG(4)
1: fldpi
fadd %st(0)
fxch %st(1)
2: fprem1
fnstsw %ax
testl $0x400,%eax
jnz 2b
fstp %st(1)
fcos
ret
3: //call __errno_location@PLT
//movl $EDOM, (%rax)
jmp 4b
END (__cosl)
weak_alias (__cosl, cosl)
@@ -0,0 +1,89 @@
/* ix87 specific implementation of exp(x)-1.
Copyright (C) 1996,1997,2001,2002,2008,2009 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
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. */
/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
#include <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(minus1,@object)
minus1: .double -1.0
ASM_SIZE_DIRECTIVE(minus1)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
ASM_TYPE_DIRECTIVE(l2e,@object)
l2e: .tfloat 1.442695040888963407359924681002
ASM_SIZE_DIRECTIVE(l2e)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__expm1l)
movzwl 8+8(%rsp), %eax // load sign bit and 15-bit exponent
xorb $0x80, %ah // invert sign bit (now 1 is "positive")
cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)?
jae __expl@PLT // (if num is denormal, it is at least >= 64.0)
fldt 8(%rsp) // x
fxam // Is NaN or +-Inf?
fstsw %ax
movb $0x45, %ch
andb %ah, %ch
cmpb $0x40, %ch
je 3f // If +-0, jump.
cmpb $0x05, %ch
je 2f // If +-Inf, jump.
fldt MO(l2e) // log2(e) : x
fmulp // log2(e)*x
fld %st // log2(e)*x : log2(e)*x
frndint // int(log2(e)*x) : log2(e)*x
fsubr %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x)
fxch // fract(log2(e)*x) : int(log2(e)*x)
f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x)
fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fscale // 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrp %st, %st(1) // 2^(log2(e)*x)-1
ret
2: testl $0x200, %eax // Test sign.
jz 3f // If positive, jump.
fstp %st
fldl MO(minus1) // Set result to -1.0.
3: ret
END(__expm1l)
libm_hidden_def (__expm1l)
weak_alias (__expm1l, expm1l)
@@ -0,0 +1,27 @@
/* Absolute value of floating point number.
Copyright (C) 2002 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>
double
__fabs (double x)
{
return __builtin_fabs (x);
}
weak_alias (__fabs, fabs)
@@ -0,0 +1,27 @@
/* Absolute value of floating point number.
Copyright (C) 2002 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>
float
__fabsf (float x)
{
return __builtin_fabsf (x);
}
weak_alias (__fabsf, fabsf)
@@ -0,0 +1,28 @@
/* Absolute value of floating point number.
Copyright (C) 2002 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 <sysdep.h>
.text
ENTRY(__fabsl)
fldt 8(%rsp)
fabs
ret
END(__fabsl)
weak_alias (__fabsl, fabsl)
@@ -0,0 +1,44 @@
/* Compute positive difference.
Copyright (C) 1997, 1998, 2002, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <sysdep.h>
.text
ENTRY(__fdiml)
fldt 8(%rsp) // x
fldt 24(%rsp) // x : y
fucomi %st(1), %st
jp 1f
jc 3f
fstp %st(1)
fldz
jmp 2f
3: fsubrp %st, %st(1)
ret
1: fucomi %st(0), %st
fcmovnu %st(1), %st
2: fstp %st(1)
ret
END(__fdiml)
weak_alias (__fdiml, fdiml)
@@ -0,0 +1,16 @@
/*
* Written by Joe Keane <jgk@jgk.org>.
* Adopted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
ENTRY(__finitel)
movl 16(%rsp),%eax
orl $0xffff8000, %eax
incl %eax
shrl $31, %eax
ret
END (__finitel)
weak_alias (__finitel, finitel)
hidden_def (__finitel)
@@ -0,0 +1,30 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Changes for x86-64 by Andreas Jaeger <aj@suse.de>=09
* Public domain.
*/
#include <machine/asm.h>
ENTRY(__floorl)
fldt 8(%rsp)
fstcw -4(%rsp) /* store fpu control word */
/* We use here %edx although only the low 1 bits are defined.
But none of the operations should care and they are faster
than the 16 bit operations. */
movl $0x400,%edx /* round towards -oo */
orl -4(%rsp),%edx
andl $0xf7ff,%edx
movl %edx,-8(%rsp)
fldcw -8(%rsp) /* load modified control word */
frndint /* round */
fldcw -4(%rsp) /* restore original control word */
ret
END (__floorl)
weak_alias (__floorl, floorl)
@@ -0,0 +1,36 @@
/* Compute maximum of two numbers, regarding NaN as missing argument.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <sysdep.h>
.text
ENTRY(__fmax)
ucomisd %xmm0, %xmm1
jp 1f // jump if unordered
maxsd %xmm1, %xmm0
jmp 2f
1: ucomisd %xmm1, %xmm1 // Is xmm1 a NaN?
jp 2f // then return xmm0
movsd %xmm1, %xmm0 // otherwise return xmm1
2: ret
END(__fmax)
weak_alias (__fmax, fmax)
@@ -0,0 +1,36 @@
/* Compute maximum of two numbers, regarding NaN as missing argument.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <sysdep.h>
.text
ENTRY(__fmaxf)
ucomiss %xmm0, %xmm1
jp 1f // jump if unordered
maxss %xmm1, %xmm0
jmp 2f
1: ucomiss %xmm1, %xmm1 // Is xmm1 a NaN?
jp 2f // then return xmm0
movss %xmm1, %xmm0 // otherwise return xmm1
2: ret
END(__fmaxf)
weak_alias (__fmaxf, fmaxf)
@@ -0,0 +1,40 @@
/* Compute maximum of two numbers, regarding NaN as missing argument.
Copyright (C) 1997, 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <sysdep.h>
.text
ENTRY(__fmaxl)
fldt 8(%rsp) // x
fldt 24(%rsp) // x : y
fucomi %st(0), %st
fcmovu %st(1), %st // now %st contains y if not NaN, x otherwise
fxch
fucomi %st(1), %st
fcmovb %st(1), %st
fstp %st(1)
ret
END(__fmaxl)
weak_alias (__fmaxl, fmaxl)
@@ -0,0 +1,36 @@
/* Compute minimum of two numbers, regarding NaN as missing argument.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <sysdep.h>
.text
ENTRY(__fmin)
ucomisd %xmm0, %xmm1
jp 1f // jump if unordered
minsd %xmm1, %xmm0
jmp 2f
1: ucomisd %xmm1, %xmm1 // Is xmm1 a NaN?
jp 2f // then return xmm0
movsd %xmm1, %xmm0 // otherwise return xmm1
2: ret
END(__fmin)
weak_alias (__fmin, fmin)
@@ -0,0 +1,36 @@
/* Compute minimum of two numbers, regarding NaN as missing argument.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.de>, 2002.
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 <sysdep.h>
.text
ENTRY(__fminf)
ucomiss %xmm0, %xmm1
jp 1f // jump if unordered
minss %xmm1, %xmm0
jmp 2f
1: ucomiss %xmm1, %xmm1 // Is xmm1 a NaN?
jp 2f // then return xmm0
movss %xmm1, %xmm0 // otherwise return xmm1
2: ret
END(__fminf)
weak_alias (__fminf, fminf)
@@ -0,0 +1,38 @@
/* Compute minimum of two numbers, regarding NaN as missing argument.
Copyright (C) 1997, 1998, 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <sysdep.h>
.text
ENTRY(__fminl)
fldt 8(%rsp) // x
fldt 24(%rsp) // x : y
fucomi %st(0), %st
fcmovu %st(1), %st // now %st contains y if not NaN, x otherwise
fucomi %st(1), %st
fcmovnb %st(1), %st
fstp %st(1)
ret
END(__fminl)
weak_alias (__fminl, fminl)
@@ -0,0 +1,43 @@
/* Return classification value corresponding to argument.
Copyright (C) 1997, 2000, 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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"
int
__fpclassifyl (long double x)
{
u_int32_t ex, hx, lx;
int retval = FP_NORMAL;
GET_LDOUBLE_WORDS (ex, hx, lx, x);
ex &= 0x7fff;
if ((ex | lx | hx) == 0)
retval = FP_ZERO;
else if (ex == 0 && (hx & 0x80000000) == 0)
retval = FP_SUBNORMAL;
else if (ex == 0x7fff)
retval = ((hx & 0x7fffffff) | lx) != 0 ? FP_NAN : FP_INFINITE;
return retval;
}
libm_hidden_def (__fpclassifyl)
@@ -0,0 +1,35 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Public domain.
*/
#include <machine/asm.h>
ENTRY(__ilogbl)
fldt 8(%rsp)
/* I added the following ugly construct because ilogb(+-Inf) is
required to return INT_MAX in ISO C99.
-- jakub@redhat.com. */
fxam /* Is NaN or +-Inf? */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fxtract
fstp %st
fistpl -4(%rsp)
fwait
movl -4(%rsp),%eax
ret
1: fstp %st
movl $0x7fffffff, %eax
ret
END (__ilogbl)
weak_alias (__ilogbl, ilogbl)
@@ -0,0 +1,37 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Change for long double by Ulrich Drepper <drepper@cygnus.com>.
* Intel i387 specific version.
* Public domain.
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0;
* no branching!
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
int __isinfl(long double x)
#else
int __isinfl(x)
long double x;
#endif
{
int32_t se,hx,lx;
GET_LDOUBLE_WORDS(se,hx,lx,x);
/* This additional ^ 0x80000000 is necessary because in Intel's
internal representation of the implicit one is explicit. */
lx |= (hx ^ 0x80000000) | ((se & 0x7fff) ^ 0x7fff);
lx |= -lx;
se &= 0x8000;
return ~(lx >> 31) & (1 - (se >> 14));
}
hidden_def (__isinfl)
weak_alias (__isinfl, isinfl)
@@ -0,0 +1,48 @@
/* s_isnanl.c -- long double version for i387 of s_isnan.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* isnanl(x) returns 1 is x is nan, else 0;
* no branching!
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
int __isnanl(long double x)
#else
int __isnanl(x)
long double x;
#endif
{
int32_t se,hx,lx;
GET_LDOUBLE_WORDS(se,hx,lx,x);
se = (se & 0x7fff) << 1;
/* The additional & 0x7fffffff is required because Intel's
extended format has the normally implicit 1 explicit
present. Sigh! */
lx |= hx & 0x7fffffff;
se |= (u_int32_t)(lx|(-lx))>>31;
se = 0xfffe - se;
return (int)((u_int32_t)(se))>>16;
}
hidden_def (__isnanl)
weak_alias (__isnanl, isnanl)
@@ -0,0 +1,31 @@
/* Round argument to nearest integral value according to current rounding
direction.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.d>, 2002.
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 <sysdep.h>
.text
ENTRY(__llrint)
cvtsd2si %xmm0,%rax
ret
END(__llrint)
weak_alias (__llrint, llrint)
strong_alias (__llrint, __lrint)
weak_alias (__llrint, lrint)
@@ -0,0 +1,31 @@
/* Round argument to nearest integral value according to current rounding
direction.
Copyright (C) 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Andreas Jaeger <aj@suse.d>, 2002.
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 <sysdep.h>
.text
ENTRY(__llrintf)
cvtss2si %xmm0,%rax
ret
END(__llrintf)
weak_alias (__llrintf, llrintf)
strong_alias (__llrintf, __lrintf)
weak_alias (__llrintf, lrintf)
@@ -0,0 +1,34 @@
/* Round argument to nearest integral value according to current rounding
direction.
Copyright (C) 1997, 2002 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 <sysdep.h>
.text
ENTRY(__llrintl)
fldt 8(%rsp)
fistpll -8(%rsp)
fwait
movq -8(%rsp),%rax
ret
END(__llrintl)
weak_alias (__llrintl, llrintl)
strong_alias (__llrintl, __lrintl)
weak_alias (__llrintl, lrintl)
@@ -0,0 +1,72 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $")
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
/* The fyl2xp1 can only be used for values in
-1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2
0.29 is a safe value.
*/
limit: .tfloat 0.29
/* Please note: we use a double value here. Since 1.0 has
an exact representation this does not effect the accuracy
but it helps to optimize the code. */
one: .double 1.0
/*
* Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29,
* otherwise fyl2x with the needed extra computation.
*/
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__log1pl)
fldln2
fldt 8(%rsp)
fxam
fnstsw
fld %st
testb $1, %ah
jnz 3f // in case x is NaN or ±Inf
4:
fabs
fldt MO(limit)
fcompp
fnstsw
andb $1,%ah
jz 2f
faddl MO(one)
fyl2x
ret
2: fyl2xp1
ret
3: testb $4, %ah
jnz 4b // in case x is ±Inf
fstp %st(1)
fstp %st(1)
ret
END (__log1pl)
weak_alias (__log1pl, log1pl)
@@ -0,0 +1,19 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Public domain.
*/
#include "math_private.h"
long double
__logbl (long double x)
{
long double res;
asm ("fxtract\n"
"fstp %%st" : "=t" (res) : "0" (x));
return res;
}
weak_alias (__logbl, logbl)
@@ -0,0 +1 @@
/* Not needed, see s_llrint.S. */
@@ -0,0 +1 @@
/* Not needed, see s_llrintf.S. */
@@ -0,0 +1 @@
/* Not needed, see s_llrintl.S. */
@@ -0,0 +1,21 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*/
/* Adapted for use as nearbyint by Ulrich Drepper <[email protected]>. */
#include <machine/asm.h>
ENTRY(__nearbyintl)
fldt 8(%rsp)
fnstcw -4(%rsp)
movl -4(%rsp), %eax
orl $0x20, %eax
movl %eax, -8(%rsp)
fldcw -8(%rsp)
frndint
fclex
fldcw -4(%rsp)
ret
END (__nearbyintl)
weak_alias (__nearbyintl, nearbyintl)
@@ -0,0 +1,124 @@
/* s_nextafterl.c -- long double version of s_nextafter.c.
* Special version for i387.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* IEEE functions
* nextafterl(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
long double __nextafterl(long double x, long double y)
#else
long double __nextafterl(x,y)
long double x,y;
#endif
{
u_int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
int32_t esx,esy;
GET_LDOUBLE_WORDS(esx,hx,lx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = esx&0x7fff; /* |x| */
iy = esy&0x7fff; /* |y| */
/* Intel's extended format has the normally implicit 1 explicit
present. Sigh! */
if(((ix==0x7fff)&&(((hx&0x7fffffff)|lx)!=0)) || /* x is nan */
((iy==0x7fff)&&(((hy&0x7fffffff)|ly)!=0))) /* y is nan */
return x+y;
if(x==y) return y; /* x=y, return y */
if((ix|hx|lx)==0) { /* x == 0 */
long double u;
SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */
u = math_opt_barrier (x);
u = u * u;
math_force_eval (u); /* raise underflow flag */
return x;
}
if(esx>=0) { /* x > 0 */
if(esx>esy||((esx==esy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
/* x > y, x -= ulp */
if(lx==0) {
if (hx <= 0x80000000) {
if (esx == 0) {
--hx;
} else {
esx -= 1;
hx = hx - 1;
if (esx > 0)
hx |= 0x80000000;
}
} else
hx -= 1;
}
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) {
hx += 1;
if (hx==0 || (esx == 0 && hx == 0x80000000)) {
esx += 1;
hx |= 0x80000000;
}
}
}
} else { /* x < 0 */
if(esy>=0||(esx>esy||((esx==esy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){
/* x < y, x -= ulp */
if(lx==0) {
if (hx <= 0x80000000) {
esx -= 1;
hx = hx - 1;
if ((esx&0x7fff) > 0)
hx |= 0x80000000;
} else
hx -= 1;
}
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) {
hx += 1;
if (hx==0 || (esx == 0xffff8000 && hx == 0x80000000)) {
esx += 1;
hx |= 0x80000000;
}
}
}
}
esy = esx&0x7fff;
if(esy==0x7fff) return x+x; /* overflow */
if(esy==0) {
long double u = x*x; /* underflow */
math_force_eval (u); /* raise underflow flag */
}
SET_LDOUBLE_WORDS(x,esx,hx,lx);
return x;
}
weak_alias (__nextafterl, nextafterl)
strong_alias (__nextafterl, __nexttowardl)
weak_alias (__nextafterl, nexttowardl)
@@ -0,0 +1,107 @@
/* s_nexttoward.c
* Special i387 version
* Conversion from s_nextafter.c by Ulrich Drepper, Cygnus Support,
* drepper@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* IEEE functions
* nexttoward(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include "math.h"
#include "math_private.h"
#include <float.h>
#ifdef __STDC__
double __nexttoward(double x, long double y)
#else
double __nexttoward(x,y)
double x;
long double y;
#endif
{
int32_t hx,ix,iy;
u_int32_t lx,hy,ly,esy;
EXTRACT_WORDS(hx,lx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = esy&0x7fff; /* |y| */
/* Intel's extended format has the normally implicit 1 explicit
present. Sigh! */
if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
((iy>=0x7fff)&&((hy&0x7fffffff)|ly)!=0)) /* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
double u;
INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
u = math_opt_barrier (x);
u = u * u;
math_force_eval (u); /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00
|| (((ix>>20)&0x7ff)==iy-0x3c00
&& (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
&& (lx<<11)>ly)))) { /* x > y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
} else { /* x < 0 */
if (esy<0x8000||((ix>>20)&0x7ff)>iy-0x3c00
|| (((ix>>20)&0x7ff)==iy-0x3c00
&& (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
&& (lx<<11)>ly)))) {/* x < y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
}
hy = hx&0x7ff00000;
if(hy>=0x7ff00000) {
x = x+x; /* overflow */
if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
/* Force conversion to double. */
asm ("" : "+m"(x));
return x;
}
if(hy<0x00100000) {
double u = x*x; /* underflow */
math_force_eval (u); /* raise underflow flag */
}
INSERT_WORDS(x,hx,lx);
return x;
}
weak_alias (__nexttoward, nexttoward)
#ifdef NO_LONG_DOUBLE
strong_alias (__nexttoward, __nexttowardl)
weak_alias (__nexttoward, nexttowardl)
#endif
@@ -0,0 +1,87 @@
/* s_nexttowardf.c -- float version of s_nextafter.c.
* Special i387 version.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
#include "math.h"
#include "math_private.h"
#include <float.h>
#ifdef __STDC__
float __nexttowardf(float x, long double y)
#else
float __nexttowardf(x,y)
float x;
long double y;
#endif
{
int32_t hx,ix,iy;
u_int32_t hy,ly,esy;
GET_FLOAT_WORD(hx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = esy&0x7fff; /* |y| */
/* Intel's extended format has the normally implicit 1 explicit
present. Sigh! */
if((ix>0x7f800000) || /* x is nan */
(iy>=0x7fff&&(((hy&0x7fffffff)|ly)!=0))) /* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
float u;
SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
u = math_opt_barrier (x);
u = u * u;
math_force_eval (u); /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80
|| (((ix>>23)&0xff)==iy-0x3f80
&& ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x > y, x -= ulp */
hx -= 1;
} else { /* x < y, x += ulp */
hx += 1;
}
} else { /* x < 0 */
if(esy<0x8000||((ix>>23)&0xff)>iy-0x3f80
|| (((ix>>23)&0xff)==iy-0x3f80
&& ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x < y, x -= ulp */
hx -= 1;
} else { /* x > y, x += ulp */
hx += 1;
}
}
hy = hx&0x7f800000;
if(hy>=0x7f800000) {
x = x+x; /* overflow */
if (FLT_EVAL_METHOD != 0)
/* Force conversion to float. */
asm ("" : "+m"(x));
return x;
}
if(hy<0x00800000) {
float u = x*x; /* underflow */
math_force_eval (u); /* raise underflow flag */
}
SET_FLOAT_WORD(x,hx);
return x;
}
weak_alias (__nexttowardf, nexttowardf)
@@ -0,0 +1,18 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Public domain.
*/
#include "math_private.h"
long double
__rintl (long double x)
{
long double res;
asm ("frndint" : "=t" (res) : "0" (x));
return res;
}
weak_alias (__rintl, rintl)
@@ -0,0 +1,18 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Changes for x86-64 by Andreas Jaeger <aj@suse.de>=09
* Public domain.
*/
#include <machine/asm.h>
ENTRY(__scalbnl)
movl %edi,-4(%rsp)
fildl -4(%rsp)
fldt 8(%rsp)
fscale
fstp %st(1)
ret
END (__scalbnl)
weak_alias (__scalbnl, scalbnl)
@@ -0,0 +1,19 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Changes for long double by Ulrich Drepper <drepper@cygnus.com>
* Public domain.
*/
#include "math_private.h"
long double
__significandl (long double x)
{
long double res;
asm ("fxtract\n"
"fstp %%st(1)" : "=t" (res) : "0" (x));
return res;
}
weak_alias (__significandl, significandl)
@@ -0,0 +1,61 @@
/* Compute sine and cosine of argument.
Copyright (C) 1997, 2000, 2001, 2005 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <machine/asm.h>
#include "bp-sym.h"
#include "bp-asm.h"
#define PARMS LINKAGE /* no space for saved regs */
#define ANGLE PARMS
#define SINP ANGLE+12
#define COSP SINP+PTR_SIZE
.text
ENTRY (BP_SYM (__sincos))
ENTER
movsd %xmm0, -8(%rsp)
fldl -8(%rsp)
fsincos
fnstsw %ax
testl $0x400,%eax
jnz 1f
fstpl (%rsi)
fstpl (%rdi)
LEAVE
retq
1: fldpi
fadd %st(0)
fxch %st(1)
2: fprem1
fnstsw %ax
testl $0x400,%eax
jnz 2b
fstp %st(1)
fsincos
fstpl (%rsi)
fstpl (%rdi)
LEAVE
retq
END (BP_SYM (__sincos))
weak_alias (BP_SYM (__sincos), BP_SYM (sincos))
@@ -0,0 +1,60 @@
/* Compute sine and cosine of argument.
Copyright (C) 1997, 2000, 2001, 2005 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <machine/asm.h>
#include "bp-sym.h"
#include "bp-asm.h"
#define PARMS LINKAGE /* no space for saved regs */
#define ANGLE PARMS
#define SINP ANGLE+12
#define COSP SINP+PTR_SIZE
.text
ENTRY (BP_SYM (__sincosl))
ENTER
fldt 8(%rsp)
fsincos
fnstsw %ax
testl $0x400,%eax
jnz 1f
fstpt (%rsi)
fstpt (%rdi)
LEAVE
retq
1: fldpi
fadd %st(0)
fxch %st(1)
2: fprem1
fnstsw %ax
testl $0x400,%eax
jnz 2b
fstp %st(1)
fsincos
fstpt (%rsi)
fstpt (%rdi)
LEAVE
retq
END (BP_SYM (__sincosl))
weak_alias (BP_SYM (__sincosl), BP_SYM (sincosl))
@@ -0,0 +1,40 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Fixed errno handling by Ulrich Drepper <drepper@redhat.com>.
*/
#include <machine/asm.h>
ENTRY(__sinl)
fldt 8(%rsp)
fxam
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 3f
4: fsin
fnstsw %ax
testl $0x400,%eax
jnz 1f
ret
.align ALIGNARG(4)
1: fldpi
fadd %st(0)
fxch %st(1)
2: fprem1
fnstsw %ax
testl $0x400,%eax
jnz 2b
fstp %st(1)
fsin
ret
3: //call __errno_location@PLT
//movl $EDOM, (%rax)
jmp 4b
END (__sinl)
weak_alias (__sinl, sinl)
@@ -0,0 +1,43 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
* Fixed errno handling by Ulrich Drepper <drepper@redhat.com>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
ENTRY(__tanl)
fldt 8(%rsp)
fxam
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 3f
4: fptan
fnstsw %ax
testl $0x400,%eax
jnz 1f
fstp %st(0)
ret
1: fldpi
fadd %st(0)
fxch %st(1)
2: fprem1
fstsw %ax
testl $0x400,%eax
jnz 2b
fstp %st(1)
fptan
fstp %st(0)
ret
3: //call __errno_location@PLT
//movl $EDOM, (%rax)
jmp 4b
END (__tanl)
weak_alias (__tanl, tanl)
@@ -0,0 +1,34 @@
/* Truncate long double value.
Copyright (C) 1997, 2003 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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 <machine/asm.h>
ENTRY(__truncl)
fldt 8(%rsp)
fstcw -4(%rsp)
movl $0xc00, %edx
orl -4(%rsp), %edx
movl %edx, -8(%rsp)
fldcw -8(%rsp)
frndint
fldcw -4(%rsp)
ret
END(__truncl)
weak_alias (__truncl, truncl)
@@ -0,0 +1,42 @@
/* AMD64 __mpn_sub_n -- Add two limb vectors of the same length > 0 and store
sum in a third limb vector.
Copyright (C) 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_sub_n)
leaq (%rsi,%rcx,8), %rsi
leaq (%rdi,%rcx,8), %rdi
leaq (%rdx,%rcx,8), %rdx
negq %rcx
xorl %eax, %eax # clear cy
.p2align 2
L(loop):
movq (%rsi,%rcx,8), %rax
movq (%rdx,%rcx,8), %r10
sbbq %r10, %rax
movq %rax, (%rdi,%rcx,8)
incq %rcx
jne L(loop)
movq %rcx, %rax # zero %rax
adcq %rax, %rax
ret
END (__mpn_sub_n)
@@ -0,0 +1,46 @@
/* AMD64 __mpn_submul_1 -- Multiply a limb vector with a limb and subtract
the result from a second limb vector.
Copyright (C) 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP 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 MP 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 MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "sysdep.h"
#include "asm-syntax.h"
.text
ENTRY (__mpn_submul_1)
movq %rdx, %r11
leaq (%rsi,%r11,8), %rsi
leaq (%rdi,%r11,8), %rdi
negq %r11
xorl %r8d, %r8d
.p2align 3
L(loop):
movq (%rsi,%r11,8), %rax
movq (%rdi,%r11,8), %r10
mulq %rcx
subq %r8, %r10
movl $0, %r8d
adcl %r8d, %r8d
subq %rax, %r10
adcq %rdx, %r8
movq %r10, (%rdi,%r11,8)
incq %r11
jne L(loop)
movq %r8, %rax
ret
END (__mpn_submul_1)
@@ -0,0 +1,24 @@
/* Definitions for x86 syntax variations.
Copyright (C) 1992,1994,1995,1997,2000,2012 Free Software Foundation, Inc.
This file is part of the GNU C Library. Its master source is NOT part of
the C library, however. The master source lives in the GNU MP 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, see
<http://www.gnu.org/licenses/>. */
#undef ALIGN
#define ALIGN(log) .align 1<<log
#undef L
#define L(body) .L##body
@@ -0,0 +1,133 @@
/* Macros to swap the order of bytes in integer values.
Copyright (C) 1997, 1998, 2000, 2002, 2003 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. */
#if !defined _BYTESWAP_H && !defined _NETINET_IN_H
# error "Never use <bits/byteswap.h> directly; include <byteswap.h> instead."
#endif
#ifndef _BITS_BYTESWAP_H
#define _BITS_BYTESWAP_H 1
#include <bits/wordsize.h>
/* Swap bytes in 16 bit value. */
#define __bswap_constant_16(x) \
((((x) >> 8) & 0xff) | (((x) & 0xff) << 8))
#if defined __GNUC__ && __GNUC__ >= 2
# define __bswap_16(x) \
(__extension__ \
({ register unsigned short int __v, __x = (x); \
if (__builtin_constant_p (__x)) \
__v = __bswap_constant_16 (__x); \
else \
__asm__ ("rorw $8, %w0" \
: "=r" (__v) \
: "0" (__x) \
: "cc"); \
__v; }))
#else
/* This is better than nothing. */
# define __bswap_16(x) \
(__extension__ \
({ register unsigned short int __x = (x); __bswap_constant_16 (__x); }))
#endif
/* Swap bytes in 32 bit value. */
#define __bswap_constant_32(x) \
((((x) & 0xff000000) >> 24) | (((x) & 0x00ff0000) >> 8) | \
(((x) & 0x0000ff00) << 8) | (((x) & 0x000000ff) << 24))
#if defined __GNUC__ && __GNUC__ >= 2
# if __WORDSIZE == 64 || (defined __i486__ || defined __pentium__ \
|| defined __pentiumpro__ || defined __pentium4__ \
|| defined __k8__ || defined __athlon__ \
|| defined __k6__)
/* To swap the bytes in a word the i486 processors and up provide the
`bswap' opcode. On i386 we have to use three instructions. */
# define __bswap_32(x) \
(__extension__ \
({ register unsigned int __v, __x = (x); \
if (__builtin_constant_p (__x)) \
__v = __bswap_constant_32 (__x); \
else \
__asm__ ("bswap %0" : "=r" (__v) : "0" (__x)); \
__v; }))
# else
# define __bswap_32(x) \
(__extension__ \
({ register unsigned int __v, __x = (x); \
if (__builtin_constant_p (__x)) \
__v = __bswap_constant_32 (__x); \
else \
__asm__ ("rorw $8, %w0;" \
"rorl $16, %0;" \
"rorw $8, %w0" \
: "=r" (__v) \
: "0" (__x) \
: "cc"); \
__v; }))
# endif
#else
# define __bswap_32(x) \
(__extension__ \
({ register unsigned int __x = (x); __bswap_constant_32 (__x); }))
#endif
#if defined __GNUC__ && __GNUC__ >= 2
/* Swap bytes in 64 bit value. */
# define __bswap_constant_64(x) \
((((x) & 0xff00000000000000ull) >> 56) \
| (((x) & 0x00ff000000000000ull) >> 40) \
| (((x) & 0x0000ff0000000000ull) >> 24) \
| (((x) & 0x000000ff00000000ull) >> 8) \
| (((x) & 0x00000000ff000000ull) << 8) \
| (((x) & 0x0000000000ff0000ull) << 24) \
| (((x) & 0x000000000000ff00ull) << 40) \
| (((x) & 0x00000000000000ffull) << 56))
# if __WORDSIZE == 64
# define __bswap_64(x) \
(__extension__ \
({ register unsigned long __v, __x = (x); \
if (__builtin_constant_p (__x)) \
__v = __bswap_constant_64 (__x); \
else \
__asm__ ("bswap %q0" : "=r" (__v) : "0" (__x)); \
__v; }))
# else
# define __bswap_64(x) \
(__extension__ \
({ union { __extension__ unsigned long long int __ll; \
unsigned int __l[2]; } __w, __r; \
if (__builtin_constant_p (x)) \
__r.__ll = __bswap_constant_64 (x); \
else \
{ \
__w.__ll = (x); \
__r.__l[0] = __bswap_32 (__w.__l[1]); \
__r.__l[1] = __bswap_32 (__w.__l[0]); \
} \
__r.__ll; }))
# endif
#endif
#endif /* _BITS_BYTESWAP_H */
@@ -0,0 +1,7 @@
/* x86_64 is little-endian. */
#ifndef _ENDIAN_H
# error "Never use <bits/endian.h> directly; include <endian.h> instead."
#endif
#define __BYTE_ORDER __LITTLE_ENDIAN

Some files were not shown because too many files have changed in this diff Show More